Using CamCASP with ECPs

For Lei

Interfacing with NWChem: Translate the movecs file into ASCII format

Create the Cluster input file

The Cluster input file is needed to generate the CamCASP command file. Of course you can create the CamCASP command file directly, but this way is more fail-safe, and also generally much faster.

Here's an example Cluster input file:


Title CdSe cluster

  Units Bohr Degrees

Molecule cdse
   Atom-charges Auto
   Charge  0
   Units Angstrom
   Cd     -0.29494191    -3.09007130     0.00000000 
   Cd     -0.94198592    -0.26700523     0.00000000 
   Cd      0.98794852     1.97545329     0.00000000 
   Se      1.71643919    -0.49259764     0.00000000 
   Se     -1.58061571     2.24485365     0.00000000 


  Molecule   cdse

  Basis      User-Def
  Aux-basis  User-Def
  SCFcode    NWChem
  MO-File    cdse-asc.movecs  Format ASCII

  #method dma-only



Some points:

To run this use:

  $ cluster < cdse.clt  

It will create a file called cdse.cks.

Editing the CamCASP command file

The CamCASP command file will need to be edited to give CamCASP the correct number of electrons on the atoms. Recall that you are using ECPs so the number of electrons are reduced. For the LANL2 ECP you have to make the changes:

Now run CamCASP:

  $ camcasp < cdse.cks > cdse.out &  

This will take only a few minutes, but could take longer for the larger systems. The outputs will include the and files.


Check the CamCASP output file for the Total multipoles: the Q00 term should be equal to the total charge on the system. If it is not (nearly) equal to that, then something has probably gone wrong.


To run CamCASP on comanche you will need to place the following commands in your .bashrc file:

  • export CAMCASP=/home/alston/CamCASP/current

    export PATH=$CAMCASP/bin:$PATH export CORES=12 export SCRATCH=/scratch/SSD_1TB/${USER}

Save the file and then type

  • $ source ~/.bashrc

You won't have to do this again as the shell will automatically execute this file when you login.

Mulfit: reduce ranks of multipoles

See the more complete tutorial at Mulfit: Reducing the rank of your multipole model

Here's the command file you need:


  Cd-Se clusters
Units Bohr


! DMA moments in Orient format : these are the moments read by Mulfit.

Total rank 5

! Radii of atoms in Bohr
! Mulfit does not have the radii of atoms larger than Argon so we need 
! to provide them. I have used 8 Bohr here. This is clearly nonsense, but it 
! was just a test. You should try more realistic radii.
Cd  8.0
Se  8.0

! This will be the name of the output multipole moments


! list of sites and ranks
! WARNING : End ranks with one blank line
Cd   0
Se   0

Keep track of the errors in the reduced rank models reported by Mulfit. These should not be taken seriously, but can be informative. For example, here's part of the output for one example:

Fitted multipoles

Cd             -0.557359   -5.839388    0.000000  Type Cd    Rank 0
                    Q00  =  0.19014

Cd             -1.780095   -0.504567    0.000000  Type Cd    Rank 0
                    Q00  =  0.42917

Cd              1.866952    3.733065    0.000000  Type Cd    Rank 0
                    Q00  =  0.65323

Se              3.243600   -0.930875    0.000000  Type Se    Rank 0
                    Q00  = -0.64112

Se             -2.986931    4.242158    0.000000  Type Se    Rank 0
                    Q00  = -0.63141


Goodness of fit parameter for the total fit:    6.3084 kJ/mol

Total multipoles
Coordinates:      0.000     0.000     0.000 bohr
                            reference   fitted
                    00  =      0.000     0.000
                   |0|  =      0.000     0.000

                    10  =      0.000     0.000
                    11c =      0.334     0.156
                    11s =     -1.009    -0.970
                   |1|  =      1.063     0.982

                    20  =     -2.716     2.453
                    21c =      0.000     0.000

You can see that Mulfit reports a Goodness of the fit and, more importantly, it also says how the fitted model reproduces the total multipoles of various ranks. Monitor these.

Potential Maps with Orient

For details see the tutorial on using Orient to compare multipole models.

Here we will compute difference maps of the potential from the DMA multipoles (rank 5) and reduced-rank multipoles obtained from Mulfit (rank 0). For a good model the potential difference should be small compared with the original potential. Say, the errors should be 10% or so. To do this step you need to have Orient installed on your computer and need to have it compiled with OPEN_GL.

Here's a sample input file. You will find a description on the tutorial site listed above.


! ORIENT display commands


      Sites     50 polarizable     11
      S-functions 50000
      Alphas 50000
      Parameter-sets 50000
      Pairs 100000

      Cd         Z     48    colour Red    Radius   4.5
      Se         Z     34    colour Blue

Molecule  cdse-ref at  0.0 0.0 0.0
Edit cdse-ref
   Bonds Auto

Molecule  cdse at  0.0 0.0 0.0
Edit cdse
   Bonds Auto

Units Bohr kJ/mol

  ! Generate the grid around one of the cdse molecules
  ! Here we generate  a grid on a 1.5 x vdW surface. 
    Name 1.5vdW
    Title "Elst : L5"
    Molecule cdse-ref
    Smoothed  0.5
    Step  0.75 B
    Radius scale 1.5  Add 0.0
  ! And the potential on this grid for the reference cdse 
  Map L5
    Molecule cdse-ref
    Title "Elst map of potential"
  ! And the potential on the same grid for the test cdse
  Map L0
    Molecule cdse
    Title "Elst map of potential"

  ! Now take the difference of the two:
  Map Diff
    Molecule cdse
    Title "Diff Elst map of potential"
    Difference L5 minus L0
  ! Write out some useful information about the difference map:
  Describe Map Diff

    0   210  0.25  1
    6   240  0.75  1
   12   300  1.0   0
   18   360  0.75  1
   24   390  0.25  1

  ! And finally, display the difference map:
  Show Diff
    Viewport 14
    Colour-scale min -30 max +30 top +30

  ! To view the original potential map comment out the Show Diff...End block and 
  ! un-comment the following lines:
  ! Show L5
  !   Viewport 14
  !   Colour-scale min -100 max +100 top +100
  !   Ball-and-stick
  ! End



To run this use:

  $ orient < cdse_difference_map.ornt   

If all goes well you will see some output on your screen and also an image in a separate window.

Look for the following (at the end of the output):

Minimum scalar value  -27.01155 kJ/mol maximum   26.80488 kJ/mol
R.m.s. difference =  12.87270 kJ/mol

This gives you a quantitative measure of the error made by the L0 model compared with the L5 model.


Things to do:

  • Here we have used a 1.5vdW surface. This can be changed. You may want to check the potential on a 2.0vdW or a 1.0vdW surface.
  • Orient does not have radii for cadmium so I have set a radius of 4.5 Bohr in the Types...End block. I think this is OK, but we should make sure it is.

  • View the potential of the L5 model. What does it tell you?

AJMPublic/camcasp/ecp (last edited 2021-04-14 11:47:27 by apw109)