# 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:

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:

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

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.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):

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):

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.

### 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:

### 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:

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

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+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: • 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_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: • --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
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:

• 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       \
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)