Navigation:

Distributed molecular properties

Introduction

CamCASP makes it very easy to calculate molecular properties like:

The multipoles are usually calculated using Anthony Stone's GDMA2 method. The GDMA code is strongly interfaced with CamCASP and is able to use wavefunctions read in via the CamCASP interfaces to various SCF/DFT codes like GAMESS(US), NWChem, DALTON2 and Gaussian (this is still a partial interface only).

Localized distributed polarizabilities are calculated using the Williams-Stone-Misquitta algorithm, and from these, the two-centerd dispersion models can be calculated.

As part of the polarizability calculation CamCASP calculates non-local (distributed) polarizabilities, and these can be use directly to calculate dispersion energies using the Dispersion programme. While this approach is very accurate and includes charge-flow effects that are important in low-dimensional delocalized systems (1-D wires), it is also quite computationally expensive to use these non-local models (the dispersion energy involves a quadruple sum over sites!). We will describe this sort of calculation elsewhere. For now we focus on the creation of local models.

Warning

CamCASP is in a state of change (goes without saying). This tutorial was written for version 5.9.x of the code, but from 6.0.x the format of the command files has changed - for the better.

Consequently, I have included both sets of files in the examples below. Please ensure that you are using files consistent with your version of the code!

Papers and Theory

I'll put in the relevant theoretical details later; for now, here are the papers that decsribe the methods:

  1. Distributed Multipole Analysis:  Stability for Large Basis Sets, Stone, J. Chem. Theory Comput., 1, 1128–1132 (2005). This is the GDMA2 method for distributed multipoles.

  2. Distributed polarizabilities obtained using a constrained density-fitting algorithm, Misquitta and Stone, J. Chem. Phys. 124, 024111 (2006). These non-local polarizabilities form the basis for the WSM polarizabilities.

  3. Accurate Induction Energies for Small Organic Molecules: 1. Theory, Misquitta and Stone, J. Chem. Theory Comput., 4, 7-18 (2008). Describes the Williams-Stone-Misquitta (WSM) polarizability models

  4. Accurate Induction Energies for Small Organic Molecules. 2. Development and Testing of Distributed Polarizability Models against SAPT(DFT) Energies//, Misquitta, Stone and Price, J. Chem. Theory Comput., 4, 19-32 (2008). Examples of the WSM polarizabilities.

  5. Dispersion energies for small organic molecules: first row atoms, Misquitta and Stone, Mol. Phys., 106, 1631-1643 (2008). Here we use the WSM models to create dispersion models.

Goals

Multipole models

These will be of the form: $Q^{a}_{lm}$, where $a$ is a site and $lm$ are the angular momentum labels. These are normally defined in the laboratory frame. They can be rotated into a local frame using the ORIENT program, but this is not normally done.

WSM polarizabilities

These are of the form : $\alpha^{a}_{lm,l'm'}$, where $a$ is a site and $lm$/$l'm'$ are the angular momentum labels. Since these are local polarizabilities, they depend on single sites only. These are normally defined in a local axis frame that is defined before the localization step in the WSM procedure. The axis frame should reflect the local symmetry, typically with the $z$-axis defined to be along terminal bonds.

To calculate polarization energies we need the distributed multipole moments in addition to the distributed polarizabilities. Additionally, we need to determine the damping. For this see the sections on Polarization and Charge-Transfer.

The WSM polarizabilities are always an approximation. It is the non-local polarizabilities that provide the full polarizability description. In localizing these we necessarily make an approximation. For many systems (those with large HOMO-LUMO gaps) the approximation may be a very good one, but for others (low dimensional systems with small gaps) the true physical polarizability description may indeed be non-local. In the latter case, the WSM description will fail.

However, when the WSM localization works, the WSM method gives us, in a well-defined sense, the best polarizability description given the limitations of the model. This is important. So, if we wish to have only rank 1 terms (dipole-dipole polarizabilities) on all atoms, then the WSM rank 1 model will be the best rank 1 model possible in the sense that it would be best able to reproduce the point-to-point polarizabilities. Likewise, if an isotropic polarizability description was needed, we'd best the most accurate model possible.

Important

The WSM models do not favour any particular type of interaction or orientation: The models will seek to best reproduce all the point-to-point polarizabilities which are uniformly distributed on a pseudo-random set of points around the molecule. The only way to bias the model is to bias the choice of point-to-point polarizabilities. So, for very specific application, when a model is desired to work for very specific interactions - say to favour just the hydrogen-bonded, or stacking orientations - the WSM models may not be the most appropriate.

WSM dispersion models

The dispersion models we need are in the form (using the SAPT notation; here the distributed multipole expansion is denoted by //DM//):

$$E^{(2)}_{\rm disp}({\rm DM}) = - \sum_{a,b>a} \sum_{n=6}^{N_{\rm max}^{ab}} f_{n}(\beta^{ab} r_{ab}) \frac{C_{n}^{ab}(\Omega_{ab})}{r_{ab}^{n}}$$

where $a$ and $b$ label sites, $f_{n}$ is a damping function of order $n$, $r_{ab}$ is the inter-site distance, $\Omega_{ab}$ are angular variables describing the orientation of the local axes on sites $a$ and $b$ relative the distance vector joining them, $N_{\rm max}^{ab}$ is the maximum order of the expansion for the particular site pair, and $C_{n}^{ab}(\Omega_{ab})$ are the angular-dependent dispersion coefficients.

We will typically use the Tang-Toennies incomplete Gamma functions for the damping:

$\begin{equation} f_{n}(\beta^{ab} r_{ab}) = 1 - \exp(-\beta^{ab} r_{ab}) \sum_{k=0}^{n}\frac{(\beta^{ab} r_{ab})^{k}}{k!} \end{equation}$

where $\beta^{ab}$ is the damping constant which may be angle-dependent, but is often assumed to be a simple constant dependent on the pair of sites $(a,b)$.

Often, even the site dependence is dropped and $\beta^{ab}$ is assumed to depend on the molecular types alone. Based on a simple idea, we have demonstrated that a good choice for the damping parameter for the //dispersion expansion only// is to use

$$ \beta^{AB}_{\rm disp} = \sqrt{2 I_{\rm A}} + \sqrt{2 I_{\rm B}}$$

where $I_{\rm A}$ and $I_{\rm B}$ are the vertical (first) ionization energies for the interacting molecules.

Our goal here is to use the WSM frequency-denpendent polarizabilities to calculate the $C_{n}^{ab}(\Omega_{ab})$ dispersion coefficients.

Important

We have recently started using damping models based on the ISA atomic ionization energies (those estimated from the ISA atomic domains). More on this later.

Distributed molecular properties with CamCASP 5.9.x and earlier

Warning

This section is for CamCASP 5.9.x and earlier versions only!

Files and Quantities

You will need:

  • Molecular geometry in either XYZ format or the Z-matrix format.
  • This should be a closed-shell system!!!
  • For calculations with an asymptotically corrected XC functional you will need the first vertical IP. Either use an experimental value or calculate it using the PBE0 functional in a $\Delta\mbox{DFT}$ procedure.

Files needed:

  • Cluster file for the properties calculation.
  • Axis file for local axis definitions (if not used, all properties will be calculated in the global (lab) axis system).

We will illustrate this with the water molecule. Here is the cluster file for water (file name: H2O.clt):

H2O.clt

Title H2O properties calculation

Global
  Units Bohr Degrees
  Overwrite Yes
End

Molecule H2O
  I.P. 0.4638 a.u. ( NIST Chem Webbook )
  ! Vibrationally averaged geometry from
  ! Mas and Szalewicz [J Chem Phys 104, 7606 (1996)].
  Units Bohr
  O   8.0    0.0000000000    0.0000000000    0.0000000000 Type O
  H1  1.0   -1.4536519600    0.0000000000   -1.1216873200 Type H
  H2  1.0    1.4536519600    0.0000000000   -1.1216873200 Type H
End

Run-Type
  Molecule    H2O
  File-prefix H2O_aTZ
  !
  Basis     aVTZ
  Aux-Basis aVTZ
  !
  Functional PBE0
  AC         MULTPOLE
  Kernel     ALDA+CHF
  SCF-Code   DALTON2006  ( As I still use the old DALTON )
  !
  Orient files for localization and display
  Process file
  Sites file
  Interface file
End

Finish

Some of the options are defaults (PBE0, DALTON) and some would be automatically set given others (the aug-cc-pVTZ auxiliary basis would automatically be set with the aug-cc-pVTZ main basis). See the CamCASP User's Guide for details of the Cluster input.

The axis file takes the form (file name H2O.axes):

H2O.axes

Axes
  O  z between H1 and H2  y from H1 to H2
  H1 z from H1 to O       y from H1 to H2
  H2 z from H2 to O       y from H2 to H1
End

See the ORIENT manual for detail of the syntax.

This job finishes rather quickly. The results will be in a directory //H2O_aTZ/OUT/ //. Have a look at the CamCASP output in file //H2O_aTZ/OUT/H2O_aTZ_camcasp.out//. The total polarizability tensor should be something like:

 Order:    1 by    1
     9.391882     0.000000     0.000000
     0.000000    10.034214     0.000000
     0.000000     0.000000     8.750848
    Isotropic polarizability:     9.39231483
  Anisotropic polarizability:     1.11142748

The formatting will not be the same as above.

The //OUT/ // directory will contain a few more files:

$ ls H2O_aTZ/OUT/
H2O_aTZ-data-summary.data       
H2O_aTZ.out                     
H2O_aTZ_0.0005_1000_f11_NL4.pol 
H2O_aTZ_DMA2_L4.mom
H2O_aTZ_camcasp.out
H2O_aTZ_dalton.out
H2O_daTZ_lim2.0_4.0_p2000_f11.p2p

The non-local polarizabilities are in file //H2O_aTZ_0.0005_1000_f11_NL4.pol//, the distributed multipole moments in //H2O_aTZ_DMA2_L4.mom// and the point-to-point polarizabilities in file //H2O_daTZ_lim2.0_4.0_p2000_f11.p2p//.

Localizing the non-local polarizabilities

In this step we will localize the non-local polarizabilities and create a dispersion model for the water dimer. We will create a rank 3 model. By default the scripts will limit polarizabilities on hydrogen atoms to rank 1. If you wish to increase this use the //--hlimit// option to the localize script.

$ mkdir L3
$ cd L3
$ cp ../H2O.clt .
$ cluster < H2O.clt  
  ( we need the various interface files Cluster creates)
$ ln -s ../H2O_aTZ/OUT .
  ( we need to make this link as the localize script need to find
    the point-to-point and the non-local polarizabilities )
$ localize.py --limit 3 --axes ../H2O.axes H2O_aTZ

Using polarizability file OUT/H2O_aTZ_0.0005_1000_f11_NL4.pol
Using p2p file OUT/H2O_aTZ_lim2.0_4.0_p2000_f11.p2p
Splitting OUT/H2O_aTZ_0.0005_1000_f11_NL4.pol into individual frequencies
.../usr/bin/csplit --prefix=H2O_aTZ_NL4_ --suffix-format=%03d.pol  
--elide-empty-files --quiet OUT/H2O_aTZ_0.0005_1000_f11_NL4.pol /# INDEX/ '{*}'
 done
Localizing polarizabilities using LW procedure ...
000 001 002 003 004 005 006 007 008 009 010  ... done
Splitting the point-to-point polarizability file ...
 Splitting point-to-point pols for H2O_aTZ
 Frequencies :            0  to           10
 Number of points in grid :         2000
 Total p2p pols per frequency :      2001000
 Split p2p pols in files with prefix : H2O_aTZ_010
done
Refining the local polarizabilities
000 001 002 003 004 005 006 007 008 009 010  ... done
Calculating the dispersion coefficients ...
...done.
Dispersion coefficients are in H2O_aTZ_ref_wt4_L3_casimir.out
The dispersion potential, in Orient form, is in H2O_aTZ_wt4_L3_C12.pot
$

The localize.py script does quite a number of things. In brief:

  1. Split the point-to-point polarizability file into files for each of the 11 frequencies (Static + 10 imaginary frequencies. The latter for the integration to get the dispersion coefficients using the Casimir-Polder integrals).
  2. Run Cluster to generate the ORIENT input files for localization.
  3. Run ORIENT to localize the non-local polarizabilities (using, in this example and by default, the Lillestolen-Wheatley localization scheme).
  4. Run the Process code to write out the PFIT commands for the refinement stage.
  5. Run PFIT to refine the localized polarizabilities.
  6. Run Process to write out the Casimir commands.
  7. Run Casimir to calculate the dispersion coefficients.

And the quantities we want are:

  • The localized and refined WSM rank 3 polarizabilities in file //H2O_daTZ_ref_wt4_L3_static.pol//
  • and the dispersion model in file //H2O_daTZ_wt4_L3_C12.pot//

Important

The localize.py script attempts to do a lot of steps and things can go wrong. If something does, use the //--keep// option to get it to keep all intermediate files so you get an idea of when and where it failed. Start by looking into the results of the localization and refinement steps. The latter will be in files of the sort //H2O_daTZ_ref_wt4_L3_000.out//. The total polarizabilities are printed out at the end. Check to see if they make sense.

Warning

Make sure your local axis definition reflects the symmetry of the system. If it doesn't, the total molecular polarizability will come out wrong. As yet there isn't a fool-proof way of doing this. Best to draw a diagram and make sure your axis system satisfies the point group (full or sub-group) of the system.

For a rank 1 isotropic model we would use (in this case we use a d-aug-cc-pVTZ basis):

$ localize.py --axes H2O.axes --limit 1 H2O_daTZ

This gives the following WSM static polarizability model (done with a d-aug-cc-pVTZ basis):

$ more H2O_daTZ_ref_wt4_L1_static.pol
# Static polarizabilities

O O
  0.000  0.0000000      0.0000000      0.0000000    
  0.000  6.6033732      0.0000000      0.0000000    
  0.000  0.0000000      6.2713511      0.0000000    
  0.000  0.0000000      0.0000000      6.5586589    

H1 H1
  0.000  0.0000000      0.0000000      0.0000000    
  0.000  2.2630552     0.42434356E-01  0.0000000    
  0.000 0.42434356E-01  1.0983324      0.0000000    
  0.000  0.0000000      0.0000000      1.3045396    

H2 H2
  0.000  0.0000000      0.0000000      0.0000000    
  0.000  2.2630552     0.42434356E-01  0.0000000    
  0.000 0.42434356E-01  1.0983324      0.0000000    
  0.000  0.0000000      0.0000000      1.3045396

For an //isotropic// rank 3 model use

$ localize.py --limit 3 --isotropic  H2O_daTZ

In the latter case, the dispersion model is

$ more H2O_daTZ_wt4_L3iso_C12iso.pot
  O  O     C6      C7       C8       C9      C10      C11      C12
    00   00   0   24.34135  0.0  537.4367  0.0  15579.54   0.0  332073.2
  End
  H  O     C6      C7       C8       C9      C10      C11      C12
    00   00   0   4.381638   0.0  48.00146   0.0    900.4451
  End
  H  H     C6      C7       C8       C9      C10      C11      C12
    00   00   0  0.8002536
  End

The //00 00 0// entries are angular momenta labels, all zero for static coefficients. Notice that the odd terms are zero for an isotropic dispersion model.

Important

All results are in atomic units.

Important

While isotropic dispersion models work reasonably well (see Misquitta and Stone, 2008), it may be better to use an anisotropic static polarizability model. Also, if you wish to use a simple $C_6$ dispersion model, consider scaling it to mimic the effects of the missing higher-ranking terms. This is an approximation, but it has been shown to give reasonable result in our 2008 Mol. Phys. paper.

Distributed molecular properties with CamCASP 6.0.x and later

Warning

This section is for CamCASP 6.0.x

Files and Quantities

We perform the properties calculation in two stages:

  • First perform an ISA calculation, then
  • perform an ISA-Pol calculation.

This can be done in one step, but as we are still learning about the ISA, I find it best to not treat the ISA as a black-box. Instead I always look at the ISA solution to see if it makes sense.

You will need:

  • Molecular geometry in either XYZ format or the Z-matrix format.
  • This should be a closed-shell system!!!
  • For calculations with an asymptotically corrected XC functional you will need the first vertical ionization potential (IP), and possibly the AC shift, $\Delta = {\rm IP} + \epsilon_{\rm HOMO}$. For the IP, either use an experimental value or calculate it using the PBE0 functional in a $\Delta\mbox{DFT}$ procedure. The energy of the HOMO level should be calculated with your functional of choice and without the AC!

Files needed:

  • Cluster file for the properties calculation.
  • Axis file for local axis definitions (if not used, all properties will be calculated in the global (lab) axis system).

We will illustrate this with the water molecule. Here is the cluster file for water (file name: H2O.clt):

H2O-ISA.clt

Title H2O : ISA calculation with CamCASP 6.0.x

Global
  Units Bohr Degrees
  Overwrite Yes
End

Molecule H2O
  I.P.       0.4638 a.u. ( NIST Chem Webbook )
  AC-Shift   0.1311 a.u. ( for PBE0 with an aTZ basis )
  ! Vibrationally averaged geometry from
  ! Mas and Szalewicz [J Chem Phys 104, 7606 (1996)].
  Units Bohr
  O   8.0    0.0000000000    0.0000000000    0.0000000000 Type O
  H1  1.0   -1.4536519600    0.0000000000   -1.1216873200 Type H
  H2  1.0    1.4536519600    0.0000000000   -1.1216873200 Type H
End

Show H2O in XYZ format  ( you could also use PDB as the format )

Run-Type
  Properties
  
  Molecule    H2O
  
  Main-Basis     aVTZ   Type  MC
  Aux-Basis      aVTZ   Type  MC   Cartesian   Use-ISA-Basis  ( this basis *could* be Spherical )
  AtomAux-Basis  aVQZ   Type  MC   Spherical   Use-ISA-Basis  ( this *must* be Spherical )
  ISA-Basis      set2   Min-S-exp-H = 0.2   
                               ( convergence tends to be better with this setting )

  Func         PBE0            ( other functionals possible, but this is probably the best )
  AC           CS00            ( [TH | LINEAR | GRAC] [CS00 | MULTIPOLE | LB94] : CS00 with NWChem)
  Kernel       ALDA+CHF        ( not needed for an ISA calculation )
  SCF-code     NWChem          ( DALTON2015, NWCHEM, etc )

  File-Prefix  H2O_aTZ         ( change this to something reasonable )

  #METHOD      isa-A           ( this will use the file isa-A from the methods/ directory )

End

Finish

The axis file takes the form (file name H2O.axes):

H2O.axes

Axes
  O  z between H1 and H2  y from H1 to H2
  H1 z from H1 to O       y from H1 to H2
  H2 z from H2 to O       y from H2 to H1
End

See the ORIENT manual for detail of the syntax.

Running CamCASP

Step (1): ISA

In this step we perform the DFT calculation, obtain the molecular orbitals, and calculate the ISA solution and ISA multipoles.

If the CamCASP code has been set up correctly, all you need to do is execute:

$ runcamcasp.py --clt H2O-ISA.clt -q bg --nproc 1 H2O_aTZ >& log &

where //H2O_aTZ// is the file prefix defined in the Cluster input file, //H2O-ISA.clt// is the name of the cluster input file (by default the scripts will look for a file with the name PREFIX.clt) and I have sent the job to the background using //-q bg//. For more details please see the User's Guide.

If all runs correctly, you will have a directory H2O_aTZ/OUT/ in which you should see the files

AIM-isosurface_H1.grid  H2O_atoms.ISA     H2O_ISA.mom             w-H2--fn1-CONV.dat
AIM-isosurface_H2.grid  H2O.cks           warnings.log            w-O--exp-analysis.dat
AIM-isosurface_O.grid   H2O.out           w-H1--exp-analysis.dat  w-O--fn1-CONV.dat
data-summary.data       H2O-ISA.clt       w-H1--fn1-CONV.dat
H2O_1.00E-03_iso.grid   H2O_ISA-GRID.mom  w-H2--exp-analysis.dat

There may be other files too. The CamCASP output is in H2O.out. Check and see what the ISA solution looks like. The ISA multipole moments are in H2O_ISA.mom and H2O_ISA-GRID.mom. We recommend you use the latter as they are generally more accurate.

Important

The ISA solution is in file H2O_atoms.ISA. This file will be needed for any subsequent restart of the ISA. We will need it for the next step.

Warning

What if the ISA does not converge?

See this page on possible solutions for help.

Step (2) : ISA-Pol

Now we need to calculate the ISA-Pol non-local polarizabilities. This will be done as follows:

In the future we will streamline these steps, but for now you do need to do some things by hand.

Here are some more details:

Warning

The above commands assume you have a directory structure like this:

  • H2O_aTZ/
    • /OUT/ /ISA-Pol/

If your directory structure differs from this, then correct the paths used in the link commands (ln -s) accordingly. For example, a common directory structure is

  • H2O/
    • /H2O_aTZ/
      • /OUT/
      /ISA-Pol/

In this case the link commands would be:

  • $ cd ISA-Pol $ ln -s ../H2O_aTZ/*movecs . $ ln -s ../H2O_aTZ/OUT/*ISA .

Cluster file for ISA-Pol

The Cluster file for this calculation is very similar to the one use for the ISA calculation. There are some changes to the RUN-TYPE block that you will need to make:

The complete Cluster file should look like this:

H2O-ISApol.clt

Title H2O : ISA-Pol calculation with CamCASP 6.0.x
 
Global
  Units Bohr Degrees
  Overwrite Yes
End
 
Molecule H2O
  I.P.       0.4638 a.u. ( NIST Chem Webbook )
  AC-Shift   0.1311 a.u. ( for PBE0 with an aTZ basis )
  ! Vibrationally averaged geometry from
  ! Mas and Szalewicz [J Chem Phys 104, 7606 (1996)].
  Units Bohr
  O   8.0    0.0000000000    0.0000000000    0.0000000000 Type O
  H1  1.0   -1.4536519600    0.0000000000   -1.1216873200 Type H
  H2  1.0    1.4536519600    0.0000000000   -1.1216873200 Type H
End
 
Show H2O in XYZ format  ( you could also use PDB as the format )
 
Run-Type
  Properties
 
  Molecule    H2O
 
  Main-Basis     aVTZ   Type  MC
  Aux-Basis      aVTZ   Type  MC   Cartesian    No-ISA-Basis  ( CHANGE MADE HERE )
  AtomAux-Basis  aVQZ   Type  MC   Spherical   Use-ISA-Basis  
  ISA-Basis      set2   Min-S-exp-H = 0.2   
                               
 
  Func         PBE0            ( other functionals possible, but this is probably the best )
  AC           CS00            ( [TH | LINEAR | GRAC] [CS00 | MULTIPOLE | LB94] : CS00 with NWChem)
  Kernel       ALDA+CHF        ( must be this for the PBE0 functional )
  SCF-code     NWChem          ( DALTON2015, NWCHEM, etc )
 
  File-Prefix  H2O_aTZ         ( change this to something reasonable )
 
  #METHOD      isa-pol-from-isa-A-restart    ( CHANGE MADE HERE )
 
End
 
Finish

The axis file we will use is

H2O.axes

Axes
  O  z between H1 and H2  y from H1 to H2
  H1 z from H1 to O       y from H1 to H2
  H2 z from H2 to O       y from H2 to H1
End

Warning

It is very important to make sure that the axis file correctly reflects the symmetry of the system. In the above example we have indicated that the sites H1 and H2 both have the same TYPE H. This means that we have to use a local axis system in which these types are indeed the same and the above one does this as this axis system obeys the $C_{2v}$ symmetry of H2O.

Running the ISA-Pol calculation

This time we do not need to run NWChem as we already have the molecular orbitals, so rather than use the runcamcasp.py - which cannot as yet skip the DFT part - we use an older Perl code which has the rather confusing name runcamcasp. But as this script does not run Cluster, we need to do this first:

$ cluster < H2O-ISApol.clt
$ runcamcasp -q bg H2O_aTZ

Warning

One version of the Cluster program needs the job name specified using --job as:

  • $ cluster --job H2O_aTZ < H2O-ISApol.clt

This will not take long. CamCASP will first restart the ISA from the solution in H2O_atoms.ISA. Then the ISA-Pol calculation will start. When done you should have an OUT/ directory containing the following files:

H2O_atoms_001.ISA      H2O_ISA-GRID_f11_NL4_fmtA.pol  warnings.log            w-H2--fn1-CONV.dat
H2O_atoms.ISA          H2O_ISA-GRID_f11_NL4_fmtB.pol  w-H1--exp-analysis.dat  w-O--exp-analysis.dat
H2O_aTZ-A-asc.movecs   H2O_ISA-GRID.mom               w-H1--fn1-CONV.dat      w-O--fn1-CONV.dat
H2O_aTZ.out            H2O_lim2.0_4.0_p2000_f11.p2p   w-H2--exp-analysis.dat

The important files are:

  • H2O_ISA-GRID_f11_NL4_fmtB.pol : The ISA-Pol non-local polarizabilites in the new format. The older, format A, will soon be deprecated.
  • H2O_lim2.0_4.0_p2000_f11.p2p : point-to-point polarizabilities.
  • H2O_aTZ.out (or H2O_aTZ_camcasp.out): CamCASP output file. Look at this to see if the results make sense. The total polarizabiities and some of the non-local polarizabilities are printed in this file.

Warning

I used a slightly different calculation to generate these examples. So file prefixes are sometimes H2O and other times H2O_aTZ. Keep this in mind!

You should have the following total polarizability tensor:

H2O total polarizabilities

    Spherical tensor form

    Origin is (   0.0000,   0.0000,   0.0000) BOHR

 Order:    0 by    0
   0.0000000E+00

 Order:    0 by    1
   0.3077651E-07   0.0000000E+00   0.0000000E+00

 Order:    0 by    2
   0.6482435E-07   0.0000000E+00   0.0000000E+00   0.3091289E-06   0.0000000E+00


 Order:    1 by    1
   0.9542880E+01   0.6063761E-05   0.0000000E+00
   0.6063696E-05   0.9985474E+01   0.0000000E+00
   0.0000000E+00   0.0000000E+00   0.9129133E+01
    Isotropic polarizability:  0.9552495E+01
  Anisotropic polarizability:  0.7417533E+00

These were calculated with a larger basis, but the aug-cc-pVTZ basis should result in something like this.

Localization of the ISA-Pol polarizabilities

The non-local polarizabilities can be localized using the localize.py code. This is a multi-step process which has been described above. Before performing the localization we need to make some decisions:

Let's have a look at both isotropic and anisotropic models.

Isotropic local models

Here we will create a rank 3 isotropic local model for water. Make sure you are in the ISA-Pol/ directory and issue the following command:

localize.py        --limit 3 --wsmlimit 3 --hlimit 3 H2O_aTZ

This will work, but I prefer to place the commands in a file (easy to keep a record and make modifications later), and also to perform the localization in a sub-directory as that way things are kept in order. So instead of the above use the following:

loc-iso-command.bash

localize.py        --limit 3 --wsmlimit 3 --hlimit 3          \
                   --isotropic  --subdir L3iso-wt3-coeff1e-3  \
                   --loc LW  --weight 3 --weightcoeff 1e-3    \
                   H2O_aTZ

Important

Place the above commands in a file. Here we have called it loc-iso-command.bash, but you could use any name. Make this file executable using the chmod command:

  • $ chmod +x loc-iso-command.bash

Now you can execute the commands in this file using:

  • $ ./loc-iso-command.bash

Here are the options used in the above:

When executed you will get the following output:

$ ./loc-iso-command.bash 
Localisation settings for H2O_aTZ 
Sub-directory    L3iso-wt3-coeff1e-3-new 
Axes file:       H2O.axes 
Pol file format: NEW 
Limit:           3 
WSM-Limit:       3 
H-Limit:         3 
Isotropic?:      isotropic 
Model file:      H2O.pdef 
Pol Cutoff:      0.0001 
Loc algorithm:   LW 
Weight:          3 
Weight coeff:    0.001 
SVD threshold:   0.0 
NoRefine?:       False 
Using polarizability file OUT/H2O_ISA-GRID_f11_NL4_fmtB.pol
Splitting OUT/H2O_ISA-GRID_f11_NL4_fmtB.pol into individual frequencies ...  done
Localizing polarizabilities using LW procedure ...
000 001 002 003 004 005 006 007 008 009 010  ... done
Preparing to refine the local polarizabilities
Using axis definition file H2O.axes
Refining the local polarizabilities
Creating new polarizability model definition file H2O.pdef
000 001 002 003 004 005 006 007 008 009 010 
Refinement finished
Refined localized polarizabilities, static and at imaginary frequency,
are in H2O_ref_wt3_L3iso_0f10.pol
Calculating the dispersion coefficients ...   done
Dispersion coefficients are in H2O_ref_wt3_L3iso_casimir.out.
The dispersion potential, in Orient form, is in H2O_ref_wt3_L3iso_C12iso.pot.

Warning

The localization script will print out the localization settings for your reference. These are also printed out in the dispersion model file (below). Notice that there seems to be an axis file present. There will be one listed even if you didn't provide any as localize.py will create an empty file if it finds none present.

And here is the content of the potential file:

H2O_ref_wt3_L3iso_C12iso.pot

!  Localisation settings for H2O 
!  Axes file:       H2O.axes  
!  Pol file format: NEW  
!  Limit:           3  
!  WSM-Limit:       3  
!  H-Limit:         3  
!  Isotropic?:      True  
!  Model file:      H2O.pdef  
!  Pol Cutoff:      0.0001  
!  Loc algorithm:   LW  
!  Weight:          3  
!  Weight coeff:    0.001  
!  SVD threshold:   0.0  
!  NoRefine?:       False  
!                       
! Dispersion coefficients for H2O (hartree bohr^n)
  O  O               C6             C7             C8             C9             C10            C11            C12
    00   00   0   24.56813         0.0          487.0036         0.0          12225.23         0.0          229547.1
  End
  H  O               C6             C7             C8             C9             C10            C11            C12
    00   00   0   4.297156         0.0          65.58407         0.0          1317.060         0.0          17518.63
  End
  H  H               C6             C7             C8             C9             C10            C11            C12
    00   00   0  0.7630528         0.0          8.072538         0.0          112.2174         0.0          1262.786
  End

One way of assessing if this is a good model is to have a look at the output files of the form H2O_ref_*_xxx.out. The one with xxx=000 contains the output of the refinement of the static polarizabilities. It's quite a large file, but the important parts are:

...
...
Parameter values:
  1            6.78769194     O_1_iso_A
  2           24.82298591     O_2_iso_A
  3          217.75336965     O_3_iso_A
  4            1.35712844     H1_1_iso_A
  5            2.20851916     H1_2_iso_A
  6            6.64427990     H1_3_iso_A
Residuals               per cent of range
Maximum      0.00130009       6.951
Minimum     -0.00122429      -6.545
R.m.s.       0.00006801       0.364
Sum of squared residuals =         0.00926

Parameter                Fitted      Anchor    Difference   Penalty
O_1_iso_A                6.78769     7.08513    -0.29744    1.72795E-06
O_2_iso_A               24.82299    29.95894    -5.13595    2.93566E-05
O_3_iso_A              217.75337   217.73979     0.01358    3.88949E-12
H1_1_iso_A               1.35713     1.23371     0.12342    6.03978E-06
H1_2_iso_A               2.20852     2.33266    -0.12415    2.39268E-06
H1_3_iso_A               6.64428     6.63712     0.00716    1.13848E-09

...
...
Total molecular polarizability
      00          z           x           y           20          21c         21s         22c         22s
     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000
     0.00000     9.50195     0.00000     0.00000    -6.08910     0.00000     0.00000     0.00000     0.00000
     0.00000     0.00000     9.50195     0.00000     0.00000    -5.27331     0.00000     0.00000     0.00000
     0.00000     0.00000     0.00000     9.50195     0.00000     0.00000    -5.27331     0.00000     0.00000
     0.00000    -6.08910     0.00000     0.00000    48.63565     0.00000     0.00000    -9.93419     0.00000
     0.00000     0.00000    -5.27331     0.00000     0.00000    56.69164     0.00000     0.00000     0.00000
     0.00000     0.00000     0.00000    -5.27331     0.00000     0.00000    39.48512     0.00000     0.00000
     0.00000     0.00000     0.00000     0.00000    -9.93419     0.00000     0.00000    46.44655     0.00000
     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000    46.44655
Mean polarizability       =   9.50195
Polarizability anisotropy =   0.00000
Finished

Look for the following:

Anisotropic local models

The only difference with anisotropic models is that you need to have the axis file present in the work directory (ISA-Pol/), and that the refinement takes a lot longer as there are more parameters to optimise.

Warning

Make sure you have copied or linked the H2O.axes file to the ISA-Pol/ directory!

Here is the localization command:

loc-aniso-command.bash

localize.py        --limit 3 --wsmlimit 3 --hlimit 3             \
                   --axes H2O.axes --subdir L3-wt3-coeff1e-3     \
                   --loc LW  --weight 3 --weightcoeff 1e-3       \
                   H2O

Execute this as before...

$ ./loc-aniso-command.bash

And the output will be in directory L3-wt3-coeff1e-3/.

The potential file now looks like this:

H2O_ref_wt3_L3_C12.pot

!  Localisation settings for H2O 
!  Axes file:       H2O.axes  
!  Pol file format: NEW  
!  Limit:           3  
!  WSM-Limit:       3  
!  H-Limit:         3  
!  Isotropic?:      False  
!  Model file:      H2O.pdef  
!  Pol Cutoff:      0.0001  
!  Loc algorithm:   LW  
!  Weight:          3  
!  Weight coeff:    0.001  
!  SVD threshold:   0.0  
!  NoRefine?:       False  
!                       
! Dispersion coefficients for H2O (hartree bohr^n)
  O  O               C6             C7             C8             C9             C10            C11            C12
    00   00   0   25.61645         0.0          570.7627         0.0          13821.21         0.0          263833.7
    00   10   1     0.0          12.23063         0.0          214.5696         0.0          4057.924
    00   11c  1     0.0         0.5356526E-02     0.0         0.8741571E-01     0.0          1.567216
    00   20   2 -0.1775111         0.0          2.132830         0.0         -165.2042         0.0         -5261.844
    00   21c  2     0.0            0.0         0.4491132E-02     0.0         0.2977902         0.0          8.375013
    00   22c  2 -0.3082694         0.0         -117.7108         0.0         -2372.298         0.0         -50530.35
    00   30   3     0.0         -3.660322         0.0         -70.86989         0.0         -1286.275
    00   31c  3     0.0        -0.9719078E-03     0.0        -0.2925074E-01     0.0        -0.5761495
    00   32c  3     0.0          6.174164         0.0          139.9797         0.0          2632.120
    00   33c  3     0.0        -0.3764183E-02     0.0        -0.7991419E-01     0.0         -1.477058
...
...

Once again we should assess the localization by looking at the file H2O_ref_wt3_L3_000.out:

H2O_ref_wt3_L3_000.out

...
...
Residuals               per cent of range
Maximum      0.00029275       1.565
Minimum     -0.00017147      -0.917
R.m.s.       0.00000528       0.028
Sum of squared residuals =         0.00006

Parameter                Fitted      Anchor    Difference   Penalty
O_10_10_A                6.93856     6.96757    -0.02901    1.69829E-08
O_10_20_A                0.28301     0.28126     0.00175    2.85077E-09
O_10_22c_A               1.29707     1.30912    -0.01205    5.35476E-08
O_10_30_A               -1.92469    -1.92309    -0.00160    5.47390E-10
O_10_31c_A               0.00117     0.00058     0.00059    3.45231E-10
O_10_32c_A              -5.06593    -5.12840     0.06247    1.42943E-07
O_10_33c_A               0.00030     0.00025     0.00005    2.40869E-12
O_11c_11c_A              6.75861     6.77750    -0.01889    7.60090E-09
O_11c_21c_A              2.60901     2.60524     0.00378    1.83090E-09
O_11c_31c_A              5.86913     5.86684     0.00229    1.47921E-10
...
...
there are lots of parameters this time...
...
...
Total molecular polarizability
      00          z           x           y           20          21c         21s         22c         22s
     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000
     0.00000     9.52357     0.00000     0.00000    -4.87982     0.00000     0.00000    -2.84111     0.00000
     0.00000     0.00000     9.97845     0.00000     0.00000   -10.02021     0.00000     0.00000     0.00000
     0.00000     0.00000     0.00000     9.09659     0.00000     0.00000    -4.55987     0.00000     0.00236
     0.00000    -4.87982     0.00000     0.00000    45.20779    -0.00055     0.00000     2.24896     0.00000
     0.00000     0.00000   -10.02021     0.00000    -0.00055    57.65207     0.00000     0.00034     0.00000
     0.00000     0.00000     0.00000    -4.55987     0.00000     0.00000    46.54062     0.00000    -0.00051
     0.00000    -2.84111     0.00000     0.00000     2.24896     0.00034     0.00000    47.39755     0.00000
     0.00000     0.00000     0.00000     0.00236     0.00000     0.00000    -0.00051     0.00000    44.94860
Mean polarizability       =   9.53287
Polarizability anisotropy =   0.76384
Finished

Notice that the residuals are much smaller this time. Also, notice that the mean and anisotropic dipole-dipole polarizabilities are very close to those from CamCASP (above).

Warning

GOTCHAs

  • If your axis system is inconsistent with the TYPE symmetry declared in the Cluster file, then the residuals will be quite a lot larger and the solution will be garbage.

  • Your axis file may not be in the work directory.

Hetero-dimers : Dispersion models

The localize.py script works for homogeneous dimers. To get dispersion coefficients for a heterogeneous dimer, say water and Cl$^-$, you will need to first follow the above procedure to calculate the WSM models for the homogeneous case for both molecules. We then need to create the Casimir command file for the dimer. For this we will use the Process code. Then run Casimir and get our dispersion model. In a sense, all we need to do is give Casimir the WSM frequency-dependent polarizabilities in the right format.

Let's say we have calculated the WSM polarizabilities and now have isotropic rank 3 static polarizabilities and dispersion coefficients for H$_2$O and Cl$^-$. We will make use of a few files from these calculations:

  • WSM polarizabilities at all frequencies for each of the monomers: These are in the files //H2O_aTZ_ref_wt4_L3_0f10.pol// and //Clminus_ref_wt4_L3_0f10.pol// (I used the file prefix //Clminus// for the Cl$^-$ calculation).

  • The Process command files that create the Casimir input files: These are //Clminus_casimir.prss// and //H2O_aTZ_casimir.prss//.

What we need to do is use the data in these files as the input to the Casimir code which uses the localized polarizabilities to calculate the dispersion coefficients. Unfortunately, at present, this requires creating an input file to the Process program by hand. Create a file //Clminus_H2O_casimir.prss// that contains the following:

TITLE  PROCESS file to write the CASIMIR input for Cl- and H2O

Set Global-data
   ! This should be the PATH to your CamCASP installation:
   CamCASP /home/alston/CamCASP/current
   Units  BOHR  DEGREES
   Overwrite
   Echo Off
End


Molecule Cl-
  Units BOHR
  Cl         17.00    0.00000000    0.00000000    0.00000000   Type Cl
End

Read local pols for Cl-
   Use ascii file Clminus_ref_wt4_L3_0f10.pol
   Maximum rank  3
   Limit rank to 3
   Frequencies STATIC + 10
End

Molecule H2O
  Units BOHR
  O           8.00    0.00000000    0.00000000    0.00000000   Type O
  H1          1.00   -1.45365196    0.00000000   -1.12168732   Type H
  H2          1.00    1.45365196    0.00000000   -1.12168732   Type H
End

Read local pols for H2O
   Use ascii file H2O_aTZ_ref_wt4_L3_0f10.pol
   Maximum rank  3
   Limit rank to 3
   Limit rank to 1 for sites +++
      H1       H2      
   Frequencies STATIC + 10
End


Write
   File-prefix Clminus-H2O
   CASIMIR file for Cl- +++
        and H2O +++
        with cutoff 0.0001
End

Finish

All the information in this file is contained in the files //Clminus_casimir.prss// and //H2O_aTZ_casimir.prss//. All we have done here is to merge them and add one additional line.

Now create the Casimir command file using:

$ process < Clminus_H2O_casimir.prss > Clminus_H2O_WSM_L3iso.casimir

This file looks like:

Title Cl- ... H2O
Frequencies   0.5    10
Skip  0
Print nonzero

Molecule  Cl-
Site  Cl  type  Cl
10 10
      27.59015300       27.31627800       25.80430000       21.76549300       15.39099900    
      8.552807300       3.515195700       1.006707300      0.1707998800      0.6797508900E-02
11c 11c
      27.59015300       27.31627800       25.80430000       21.76549300       15.39099900    
      8.552807300       3.515195700       1.006707300      0.1707998800      0.6797508900E-02
11s 11s
...
...
33s 33s
      1094.410600       1088.040400       1050.364500       927.8203500       676.7228300    
      373.9877100       156.0034600       46.76507900       7.607165400      0.2612380000    
End

Molecule  H2O
Site  O  type  O
10 10
      6.483185500       6.461184100       6.331261300       5.904744500       4.952804400    
      3.465239500       1.847456100      0.6679903900      0.1282095700      0.5212038600E-02
11c 11c
 ...
 ...
Site  H2  type  H
10 10
      1.427775400       1.419258000       1.370174900       1.221677900      0.9394025700    
     0.5798742900      0.2672248900      0.8300164000E-01  0.1333205500E-01  0.4515054600E-03
11c 11c
      1.427775400       1.419258000       1.370174900       1.221677900      0.9394025700    
     0.5798742900      0.2672248900      0.8300164000E-01  0.1333205500E-01  0.4515054600E-03
11s 11s
      1.427775400       1.419258000       1.370174900       1.221677900      0.9394025700    
     0.5798742900      0.2672248900      0.8300164000E-01  0.1333205500E-01  0.4515054600E-03
End

CGdir  /home/alston/CamCASP/current/data/realcg
Dispersion   12  Cl-  H2O

Finish

Important

The path to CamCASP is used in the last line where it points to the location of the real Clebsch-Gordon coefficients which are in file //realcg//.

If this path is not correct or if you do not have read privileges to the data in file //readcg// the next step will fail.

Now run the Casimir program:

$ casimir < Clminus_H2O_WSM_L3iso.casimir > Clminus_H2O_WSM_L3iso.casimir.out

The output contains the dispersion model. It will be at the end of this file. For this example we get:

Dispersion coefficients for Cl- and H2O
  Cl  O   C6    C7     C8     C9    C10     C11     C12
    00   00   0   72.95569   0.0    1653.591   0.0    42806.65   0.0  858440.7
  End
  Cl  H   C6    C7     C8     C9    C10     C11     C12
    00   00   0   14.14859   0.0    190.8377   0.0     2799.988
  End

Important

If anisotropic polarizabilities are used, the resulting dispersion model will be in the axis system used for those polarizabilities!

Basis Sets

A quick word about these: molecular properties like polarizabilities and dispersion coefficients require quite large and diffuse basis sets before they are sufficiently converged. The //aug-cc-pVTZ// basis set used in the above calculations is not good enough, particularly for the anion Cl$^-$. It is better to use a much more diffuse basis. Only bear in mind that there may not be an optimized auxiliary basis available for some of the very diffuse basis sets like the //d-aug-cc-pVTZ//. Some degree of experimentation with basis sets may be required.

Theory level

We will normally use linear-response time-dependent (LR-TDDFT) DFT with PBE0 and the hybrid ALDA+CHF response kernel. But for some systems, for example anions like Cl$^-$, PBE0 may not be adequate because of the self-interaction error. The asymptotic correction can correct for this in some cases, but when it fails, you may need to use LC-PBE0 and use IP-tuning to determine the range-separation parameter. However this leads to a problem: the response kernel in CamCASP cannot yet be used consistently for range-separated functionals. So in this case please contact us to see what your options are.

Energy evaluations with ORIENT

The dispersion models should be tested against calculated SAPT(DFT) dispersion energies (use the total dispersion energy $E^{(2)}_{\rm disp,tot} = E^{(2)}_{\rm disp,pol} + E^{(2)}_{\rm disp,exch}$ in the comparisons). To do this we need the damping parameter. A good choice for the dispersion models is the expression given above with $\beta$ a function of the vertical ionization energies of the interacting molecules. Alternatively, this coefficient can be obtained by fitting to the SAPT(DFT) total dispersion energies.

Model energies can be evaluated using the ORIENT program for dimer geometries specified in the angle-axis format. See the Orient: Energy Calculations section for more details.

AJMPublic/camcasp/pols-cn (last edited 2021-04-14 11:41:15 by apw109)