Contents

Navigation:

# Distributed polarizabilities using the ISA-Pol algorithm

Important

This tutorial applies to CamCASP 6.0.17 and higher. At the time of writing this version of the code is private. If you are interested in these calculations, please contact AJM.

Papers:

- //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. - //ISA-Pol: Distributed polarizabilities and dispersion models from a basis-space implementation of the iterated stockholder atoms procedure//, A. J. Misquitta & A. J. Stone (2018)

**submitted**Version on ArXiv

In an earlier tutorial we have described how the CamCASP code can be used to calculate distributed non-local polarizabilities. There we have described how these non-local polarizabilities can be localised, and from these (frequency-dependent) local polarizabilities, we can derive dispersion models (either isotropic or anisotropic) with terms from $C_6^{ab}$ to $C_{12}^{ab}$. We have used this procedure quite extensively, and while it is satisfactory, right from the start (way back in 2006) we felt that the algorithm had shortcomings that led to un-physical terms, poor convergence with basis set, and perhaps a sub-optimal convergence with rank.

The problem lay with the way we partitioned the transition-densities: we used a constrained version of density-fitting (Misquitta & Stone 2006) which gave us a numerically stable way of partitioning the FDDS, but there never was any good reason for the results to be physical. Further, while the dominant terms in the resulting partitioned polarizability were usually well-behaved, it was common to get negative terms that then led to negative $C_n^{ab}$ dispersion coefficients. This could sometimes be avoided by altering the weighting scheme used in the refinement process, but it was never satisfactory.

An alternative is to use the ISA for the partitioning. The theoretical basis of this partitioning will be presented in a forthcoming paper (Misquitta 2016). Here we describe the way these calculations proceed.

## Methods

In CamCASP 6.0.017 and later we have used the concept of METHODs that can be accessed from the Cluster file. The methods are contained in $CAMCASP/methods/ and usually consist of two file:

- method
- method.clt_tmpl

For example, an ISA calculation using the 'A' algorithm is a method called *isa-A* (case sensitive!) and its associated sample Cluster command file is *isa-A.clt_tmpl*. The method should be altered only if you know what you are doing. The sample Cluster command file gives you guidance on how to construct the Cluster file for the method.

Important

The files in $CAMCASP/methods/ will usually be consistent with the CamCASP version, but this wiki page may not be. So the contents of this directory should over-rule any statements made on this page.

## The ISA-Pol method

In principle the ISA-Pol algorithm can be used to calculate distributed polarizabilities in a single calculation, but it may be useful to perform the calculation in two steps, at least until the code gets more robust. These steps are:

- ISA calculation: See the tutorial on the ISA. - ISA-Pol calculation: This will use the ISA solution from the previous step.

If you wish to compare the ISA-Pol results with the earlier constrained-DF (cDF) polarizabilities, then first use the tutorial on distributed properties to calculate these polarizabilities.

Let's walk through a typical calculation using HC$_4$H as an example.

### HC$_4$H ISA calculation

Here is the Cluster file for this calculation:

Title H-[-C=C-]2-H properties : basis = aug-cc-pVDZ : func = PBE0 Global Units Bohr Degrees Overwrite Yes End Molecule unsat_m2_x1 Units Bohr H1 1 0.0 0.0 5.716421 Type H1 C1 6 0.0 0.0 3.713312 Type C2 C2 6 0.0 0.0 1.237771 Type C1 C3 6 0.0 0.0 -1.237771 Type C1 C4 6 0.0 0.0 -3.713312 Type C2 H2 1 0.0 0.0 -5.716421 Type H1 End Show unsat_m2_x1 in XYZ format Run-Type Properties ( If you already have the *-asc.movecs file then replace this ) ( with CamCASP-Only as you will need only a new CamCASP command file ) ( for this calculation. ) Molecule unsat_m2_x1 File-prefix unsat_m2_x1 Basis aug-cc-pVDZ Type MC Aux-Basis aug-cc-pVDZ Type MC Cartesian Use-ISA-Basis ( this could be spherical ) AtomAux-Basis aug-cc-pVTZ Type MC Spherical Use-ISA-Basis ( this MUST be spherical ) ISA-Basis set2 Min-S-exp-H = 0.2 ( This limit leads to better convergence ) Functional PBE0 ! Type of method for polarizabilities Kernel ALDA+CHF ! Program used for the DFT part: SCFcode NWChem ! If you already have a *-asc.movecs file then uncomment this line and set the name of ! the file. Use this only if the RUN-TYPE is CamCASP-Only. ! MO-File <name of MO file> Format ASCII ( set this to the name of the *-asc.movecs file ) Memory 4000 MB #method isa-A ( this will use the isa-A method from the $CAMCASP/methods/ directory ) End Finish

Comments:

- The METHOD is read from the $CAMCASP/methods/ directory. This name is case-sensitive.
- In the METHODs directory is a sample Cluster template file associated with this method. In this case it is called isa-A.clt_tmpl. Make sure that you read this file and ensure that your Cluster file (the one above) is consistent with the sample.
- There are three basis sets defined.
- The MAIN basis is the one used for the DFT calculation.
- The AUX basis is the one used for the density-fitting.
- And the ATOMAUX basis is the one used for the ISA atomic expansions. This one must use Spherical GTOs.
- Convergence tends to be best if the AUX and ATOMAUX basis sets share the s-functions. This is done using ISA-BASIS. It is also possible to use different s-function basis sets for the two auxiliary basis sets.

This file **includes** the *isa-display* template file from the CamCASP library. You will need to copy this from *$CAMCASP/data/camcasp-cmnds/* to your working directory. Edit it as you see fit, but only if you know what you are doing! For some guidance on the contents see the ISA tutorial.

Run this calculation the usual way:

$ runcamcasp.py --scfcode nwchem --nproc 2 --q batch --clt unsat_m2_x1.clt --direct unsat_m2_x1

Important

On some systems the *batch* queue does not work. If this is the case, use the *bg* (background) queue. Since this is not a queue at all, but just runs jobs directly, make sure the machine is not being used when you do this or you may risk overloading it.

The output will be in *./OUT/* and should contain the files (and a few more besides):

$ ls AIM-isosurface_C.grid unsat_m2_x1_1.00E-03_iso.grid w-C--fn1-CONV.dat AIM-isosurface_H.grid unsat_m2_x1_ISA.mom w-H--exp-analysis.dat DMA_Z0_L5.mom unsat_m2_x1_atoms.ISA w-H--fn1-CONV-TAIL-ITER.dat DMA_Z4_L5.mom unsat_m2_x1_camcasp.out w-H--fn1-CONV.dat unsat_m2_x1-data-summary.data w-C--exp-analysis.dat unsat_m2_x1.ornt w-C--fn1-CONV-TAIL-ITER.dat

As usual the CamCASP output is in file *unsat_m2_x1_camcasp.out*. Check it to see if the ISA calculation has converged appropriately. It is very important that this is the case or else the results of the ISA-Pol calculation will be meaningless.

Important

TIPS: Things to look for that are hallmarks of a converged ISA calculation are:

- A relatively small (0.01e or so) charge violation in the first ISA step, before the density-fitting constraint is included. If you find a larger charge violation then something has probably gone wrong.
Ideally the shape-functions (tabulated in the

*w-*-CONV-TAIL-ITER.dat*files) should be positive everywhere and exponentially (or nearly so) decaying. Plot them to see if this is the case. Sometimes the shape-functions for the hydrogen atoms may get negative past $6.0$ Bohr. This is usually OK. If you see appreciable wiggles (in a semilog plot of $w(r)$) or the on-set of negative values at small $r$ then something may have gone wrong.Things usually go wrong because the basis set used for the ISA was not flexible enough, or if the

**W-EPS**term was too large. See the ISA tutorial for details.

The ISA solution for this calculation is in file *unsat_m2_x1_atoms.ISA*. We will need this file in the next stage.

### HC$_4$H ISA-Pol calculation

Files needed from the previous step:

*unsat_m2_x1_atoms.ISA*: Contains the ISA solution.*unsat_m2_x1-A-asc.movecs*: The molecular orbital energies and coefficients.

The Cluster file for this calculation is essentially identical to the previous one, but we include a few different commands:

RUN-TYPE CamCASP-Only ( use this as you already have the *-asc.movecs file ) ... ...commands as above... ... MO-file unsat_m2_x1-A-asc.movecs Format ASCII #method isa-pol-from-isa-A-restart END

Call this new cluster file *unsat_m2_x1_pol.clt*. The method *isa-pol-from-isa-A-restart* is also located in *$CAMCASP/methods/* and this method contains commands that allow the ISA to restart from a previous ISA solution with the ISA-Pol calculation following this. There is no need to edit this file unless you have edited the *isa-A* method used in the previous step. Just make sure that the ISA settings in the two methods are the same.

To run the job perform the following steps. Here I have assumed the following directory structure:

./ISA/ ISA/OUT/ ./ISA-Pol/ The ISA directory is the one that contains the ISA solution from the previous step. It is currently called unsat_m2_x1, so rename it: $ mv unsat_m2_x1 ISA In this directory you will find the *.movecs file and in ISA/OUT/ you will find the *.ISA file. We need these for the ISA-Pol calculation. There is no need to copy these files, as symbolic links will suffice: $ cd ISA-Pol $ ln -s ../ISA/unsat_m2_x1-A-asc.movecs . $ ln -s ../ISA/OUT/unsat_m2_x1_atoms.ISA . Now create the CamCASP command file (and other files - not all needed): $ cluster < unsat_m2_x1_pol.clt And run the job using a fairly old bash script that was set up to re-run just the CamCASP code: $ runcamcasp unsat_m2_x1 -q batch And wait for a while...

The ISA-Pol algorithm is not fast, but is much faster than it was when it was first coded up one night after a CECAM conference dinner at Nancy, so just the fact that it worked first time round is a wonder. When the calculation is done you will obtain the files:

$ cd OUT/ $ ls unsat_m2_x1-data-summary.data unsat_m2_x1_atoms.ISA w-C--fn1-CONV.dat unsat_m2_x1_ISA-GRID.mom unsat_m2_x1_atoms_001.ISA w-H--exp-analysis.dat unsat_m2_x1_ISA-GRID_f11_NL4_fmtA.pol unsat_m2_x1_camcasp.out w-H--fn1-CONV.dat unsat_m2_x1_ISA-GRID_f11_NL4_fmtB.pol w-C--exp-analysis.dat

The output is in file *unsat_m2_x1_camcasp.out* and the polarizabilities should be present in either one or two formats. If you get two formats, then you will see two **.pol* files as shown above. The **_fmtB.pol* is a new format that I am still working on. These polarizability files contain the non-local polarizabilities at 11 imaginary frequencies. To localize these polarizabilities we use the techniques described in the Distributed Properties tutorial and outlined in the next step.

### Localizing the non-local polarizabilities

Details of how this is done can be found in the Distributed Properties tutorial. Here I will be brief.

We use the *localize.py* code to drive the localization process. This code uses a number of other Fortran programs and accepts a number of options that you can see using

$ localize.py --help

Important

Read the localize.py help page as it will probably be more up-to-date than this wiki!

Here's how we will perform the localization (see the water example for more details) for the HC$_4$H molecule. We will work in the *ISA-Pol* directory and perform the localization in a sub-directory named *L3* (I've also got a sub-directory named *L2* in the example below):

$ ls unsat_m2_x1.axes unsat_m2_x1_ISA-Pol.bash L2 OUT unsat_m2_x1_atoms.ISA unsat_m2_x1_ISA-Pol.cks L3

The *OUT/* directory contains the output of the ISA-Pol calculation. unsat_m2_x1.axes contain the local-axis definition:

Axes H1 z from C1 to H1 Global x C1 z from C1 to H1 Global x C2 z from C1 to H1 Global x C3 z from C4 to H2 Global x C4 z from C4 to H2 Global x H2 z from C4 to H2 Global x End

Besides this file we will need a *unsat_m2_x1.clt* file with the molecule definition and basic commands for Cluster to generate the relevant input files. The *unsat_m2_x1.clt* file should be placed in L3:

Title unsat_m2_x1 : Localization of polarizabilities Global Units kJ/mol Bohr Degrees Overwrite Yes End MOLECULE unsat_m2_x1 Units Bohr H1 1 0.0 0.0 5.716421 Type H1 C1 6 0.0 0.0 3.713312 Type C2 C2 6 0.0 0.0 1.237771 Type C1 C3 6 0.0 0.0 -1.237771 Type C1 C4 6 0.0 0.0 -3.713312 Type C2 H2 1 0.0 0.0 -5.716421 Type H1 END Run-Type Localization Molecule unsat_m2_x1 File-prefix unsat_m2_x1 ! If you intend to localize non-local polarizabilities created with an older ! version of CamCASP then also include the following line: ! Orient Pol-Format OLD End Finish

Run Cluster on this file and make a link from the *OUT/* directory to this one:

$ cluster < unsat_m2_x1-loc.clt $ ln -s ../OUT .

This will generate a bunch of files. Normally these will not need editing.

Now run the *localize.py* code with the following options:

$ localize.py --axes ../unsat_m2_x1.axes --limit 3 --wsmlimit 3 --hlimit 3 unsat_m2_x1

Here's a brief explanation of the options:

*--axes*: location and name of the axis file.*--limit 3*: localize non-local polarizabilities to rank 3. This rank can be treated independently from the other rank limits. Best to leave it 3.*--wsmlimit 3*: limit of WSM (this is a misnomer for ISA-Pol polarizabilities!) refined polarizabilities. This limit is applied to heavy atoms only. By default, for H-atoms the limit is kept to 1, so we need the next command...*--hlimit 3*: set the rank of polarizabilities on the H-atoms to 3.

Important

If you are localizing non-local polarizabilities written in the old format then include a *--format OLD* in the localize.py commands. In this case you should also have included the line **Orient Pol-Format OLD** in the Cluster file.

If all goes well, the localized polarizabilities will be in file *unsat_m2_x1_ref_L3_xxx.pol*, where

- xxx=000 : static pols
- xxx=001..010 : pols at imaginary frequencies.
The plain-text output files for the refinement step are

*unsat_m2_x1_L3_xxx.out*

And the dispersion model will be in file *unsat_m2_x1_ref_wt3_L3_C12.pot*

Important

To obtain an isotropic dispersion model include the option *--isotropic* in the *localize.py* command:

- $ localize.py --isotropic --limit 3 --wsmlimit 3 --hlimit 3 unsat_m2_x1

You do not need to specify an axis file in this case.

Important

It is advisable to perform different localization calculations in different directories.

# ISA-Pol with Cartesian auxiliary basis set

Warning

This part of the tutorial can only be used with CamCASP version 6.0.010 and newer. This is a developer's version of the code and it not available for download as yet.

The above calculations used auxiliary basis sets with spherical GTOs and this leads to a slight loss in accuracy as the density-fitting step does benefit from an auxiliary basis with Cartesian GTOs.

## ISA-Pol with ISA A

This calculation can be done in one step, but we strongly recommend doing it in two stages as follows:

Step 1: Perform a standard ISA A calculation with Cartesian GTOs in the AUX basis and spherical GTOs in the ATOM-AUX basis. Use method

**isa-A**. Here it is advisable to have the spherical parts of AUX and ATOM-AUX the same.Step 2: Perform an ISA-Pol calculation by restarting from this solution. Use method

**isa-pol-from-isa-A-restart**.

If you do this in two stages, then you can trap errors made in the ISA stage before going on to the more computationally demanding ISA-Pol calculation.

## ISA-Pol with A+DF ISA

This one must be done in two steps:

Step 1: Perform a standard ISA A+DF calculation with spherical GTOs in the AUX and ATOM-AUX basis sets. This is done using METHOD

**isa-A+DF**. In this step we must have AUX == ATOM-AUX!- Step 2a: Restart from the ISA solution obtained from the above, this time using the ISA A-algorithm (no DF!) and a Cartesian AUX basis and spherical ATOM-AUX basis. There must be no ISA iterations performed in this restart.
- Step 2b: Perform an ISA-Pol calculation.

Steps 2a and 2b can be performed using method **isa-pol-from-isa-A-restart**.

For more details on Step 1 and the restart part of Step 2 see the tutorial on using the ISA with Cartesian GTOs in the auxiliary basis.