
<<TableOfContents()>>

= 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:

[[attachment:cdse.clt]]

{{{
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.

{{{#!wiki  warning
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.
}}}

{{{#!wiki  important
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 [[ajm/camcasp/multipoles/mulfit|Mulfit: Reducing the rank of your multipole model]]

Here's the command file you need:

[[attachment:cdse.mulfit]]

{{{
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 [[https://app.ph.qmul.ac.uk/wiki/ajm:orient:display#differences_against_references_calculated_with_orient|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.

[[attachment:cdse_difference_map.ornt]]

{{{
! 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.

{{{#!wiki  important
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?
}}}