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

AJMPublic/camcasp/isa-pol (last edited 2021-04-14 11:46:41 by apw109)