# Using CamCASP with ECPs

## For Lei

### Interfacing with NWChem: Translate the movecs file into ASCII format

• The first step is to translate the JOB.movecs file (which is binary) into ASCII format.
• For NWChem this is done using the command: $readNWCHEMmocs JOB.movecs --ascii JOB-asc.movecs • Here JOB.movecs is the file created by NWChem and JOB-acs.movecs will be the ASCII form of this file. • Check to see that this file has been correctly read. • First look for the number of basis functions in BFNS. This should be the same as given by NWChem. • Next make sure that the number of molecular orbitals (NMOS) is the same. • Finally, check some energies to make sure that they are the same as in NWChem. ### 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 Global Units Bohr Degrees End 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 End Run-Type CamCASP-Only Molecule cdse Basis User-Def Aux-basis User-Def SCFcode NWChem MO-File cdse-asc.movecs Format ASCII #method dma-only End Finish Some points: • Get the units correct! • Basis: is set to 'user-def' (case does not matter). This means that CamCASP will look for the basis sets in$CAMCASP/basis/nwchem/user-def/.

• Aux-basis: is also 'user-def' so CamCASP will look in $CAMCASP/auxiliary/user-def/ • If there are no appropriate basis sets in these directories the CamCASP calculation will fail. • METHOD: is dma-only for just the calculation of the DMA multipoles. • The MO-file should be the name of the file in which the molecular orbitals are in ASCII format. 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:

• Cd: 48 to 12 electrons
• Se: 34 to 6 electrons

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 DMA_Z0_L5.mom and DMA_Z4_L5.mom files. Warning 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. Important 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: Title Cd-Se clusters Units Bohr Verbose ! DMA moments in Orient format : these are the moments read by Mulfit. DMA DMA_Z4_L5.mom 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. Radii Cd 8.0 Se 8.0 ! This will be the name of the output multipole moments Orient cdse1205-DMA_Z4_L5-to-L0.mom ! list of sites and ranks ! WARNING : End ranks with one blank line Ranks 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 End 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 UNITS BOHR Parameters Sites 50 polarizable 11 S-functions 50000 Alphas 50000 Parameter-sets 50000 Pairs 100000 End Types Cd Z 48 colour Red Radius 4.5 Se Z 34 colour Blue End Molecule cdse-ref at 0.0 0.0 0.0 #include DMA_Z0_L5.mom End Edit cdse-ref Bonds Auto End Molecule cdse at 0.0 0.0 0.0 #include cdse1205-DMA_Z4_L5-to-L0.mom End Edit cdse Bonds Auto End Units Bohr kJ/mol Display ! Generate the grid around one of the cdse molecules ! Here we generate a grid on a 1.5 x vdW surface. Grid Name 1.5vdW Title "Elst : L5" Molecule cdse-ref Smoothed 0.5 Step 0.75 B Radius scale 1.5 Add 0.0 End ! And the potential on this grid for the reference cdse Map L5 Molecule cdse-ref Title "Elst map of potential" Potential End ! And the potential on the same grid for the test cdse Map L0 Molecule cdse Title "Elst map of potential" Potential End ! Now take the difference of the two: Map Diff Molecule cdse Title "Diff Elst map of potential" Difference L5 minus L0 End ! Write out some useful information about the difference map: Describe Map Diff Colour-map 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 End ! And finally, display the difference map: Show Diff Viewport 14 Colour-scale min -30 max +30 top +30 Ball-and-stick End ! 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 End Finish 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.

Important

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)