Contents
Navigation:
Distributed molecular properties
Introduction
CamCASP makes it very easy to calculate molecular properties like:
- Distributed multipole moments
- Distributed (local) polarizabilities, and
- Distributed (two-centred) dispersion coefficients
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:
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.
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.
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
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.
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):
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):
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.11142748The 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:
- 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).
- Run Cluster to generate the ORIENT input files for localization.
- Run ORIENT to localize the non-local polarizabilities (using, in this example and by default, the Lillestolen-Wheatley localization scheme).
- Run the Process code to write out the PFIT commands for the refinement stage.
- Run PFIT to refine the localized polarizabilities.
- Run Process to write out the Casimir commands.
- 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
EndThe 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):
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
FinishThe axis file takes the form (file name 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.
Step (2) : ISA-Pol
Now we need to calculate the ISA-Pol non-local polarizabilities. This will be done as follows:
- Create a directory for the ISA-Pol calculation.
- Create a new Cluster file for the ISA-Pol calculation.
- Execute Cluster to create a new CamCASP command file.
- Link the MO file and the ISA solution file to the ISA-Pol directory.
Run CamCASP only using the 'runcamcasp' command (note, there is no .py)
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:
- Go to directory H2O_aTZ/. This is the work directory for this job. - Make directory ISA-Pol and go there. - Copy the H2O.clt: cp ../H2O.clt ./H2O-ISApol.clt. We will edit it as shown below. - Make links: ln -s ../*.movecs . and ln -s ../OUT/*.ISA . - Check: Do the links actually point to the files? Use ls -al to check. - Edit the H2O-ISApol.clt file as shown next.
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/
- /H2O_aTZ/
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 RUN-TYPE will be PROPERTIES
- The AUX-BASIS will be the standard auxiliary basis with CARTESIAN GTOs. In the ISA calculation we used the ISA s-function set as part of this basis; here will will not do so. Why not? Because it does not seem to be needed.
Change the METHOD to isa-pol-from-isa-A-restart
The complete Cluster file should look like this:
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
FinishThe axis file we will use is
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:
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+00These 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:
- The maximum rank of the localized polarzabilities.
- Is this to be an isotropic or anisotropic model? At present we do not have the ability to generate hybrid models, but this will be done soon (it can be done with a manual editing of the intermediate files created in the localization process).
- The local axis frame for the anisotropic model. We have already defined the axis frame for H2O (see above).
- Some settings used in the refinement stage.
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:
localize.py --limit 3 --wsmlimit 3 --hlimit 3 \
--isotropic --subdir L3iso-wt3-coeff1e-3 \
--loc LW --weight 3 --weightcoeff 1e-3 \
H2O_aTZImportant
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:
- --limit 3 : perform the initial localization to rank 3.
- --wsmlimit 3 : Set the rank of the refined polarizabilities to 3. This applies to all but hydrogen atoms, for these use
- --hlimit 3 : to set the rank of the polarizabilities on hyrdogen atoms also to 3.
- --isotropic : force isotropic models
--subdir L3iso-wt3-coeff1e-3 : do the calculations in sub directory L3iso-wt3-coeff1e-3/
- --loc LW : use the Lillestolen-Wheatley localization for the initial step.
- --weight 3 : use weight scheme 3 for the refinement step
- --weightcoeff 1e-3 : and associated with this weight is a strength, set this to 0.001 in arbitrary units. This is fairly strong and it keeps the polarizabilities from wandering too far from the initial values.
- H2O : the file prefix.
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:
! 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
EndOne 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
FinishedLook for the following:
- The residuals: they should be small. 6% is OK for an isotropic model.
- The total polarizability: This should be close to what was calculated by CamCASP.
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:
localize.py --limit 3 --wsmlimit 3 --hlimit 3 \
--axes H2O.axes --subdir L3-wt3-coeff1e-3 \
--loc LW --weight 3 --weightcoeff 1e-3 \
H2OExecute 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:
! 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:
...
...
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
FinishedNotice 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
FinishAll 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
FinishImportant
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
EndImportant
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.