Contents

Distributed polarizabilities using the ISA-Pol algorithm
ISA-Pol with Cartesian auxiliary basis set
1. ISA-Pol with ISA A
2. ISA-Pol with A+DF ISA

Navigation:

CamCASP

Distributed polarizabilities using the ISA-Pol algorithm

Important

This tutorial applies to CamCASP 7.x and higher.

Updated: 10th November 2024.

Papers:

Distributed polarizabilities obtained using a constrained density-fitting algorithm, A. J. Misquitta and A. J. 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) Version on ArXiv Published version at Theoretical Chemistry Accounts

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.

Codes needed

Orient : needs to be on your PATH.
CamCASP : All codes contained in the distribution. These include: camcasp, cluster, pfit, casimir.

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 that particular CamCASP version. If you see differences in the commands used here and those in the method files in your distribution, it is likely that this wiki has got outdated.

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:

unsat_m2_x1.clt

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  Spherical  Use-ISA-Basis  ( this could be spherical or Cartesian )
  AtomAux-Basis aug-cc-pVTZ  Type MC  Spherical  Use-ISA-Basis  ( this MUST be spherical  )
  ISA-Basis     set2  Min-S-exp-H = 0.0  ( This sets the smallest s-exponent for H-atom basis sets )
  
  Functional PBE0
  ! Type of method for polarizabilities
  Kernel ALDA+CHF
  ! Program used for the DFT part: You could use NWChem, Psi4, DALTON, etc. We will use Psi4.
  SCFcode Psi4
  
  ! 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 8 GB
  
  #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 using the runcamcasp.py code. Here we specify which Cluster input file (--clt) to use, and also which directory (-d) to run in. As usual, the last argument (a positional argument) is the FILE-PREFIX used in the Cluster command file (unsat_m2_x1).

  $ runcamcasp.py --clt unsat_m2_x1.clt -d unsat_m2_x1_psi4 unsat_m2_x1

Important

On some computing systems you need to submit to a queuing system using an appropriate queue. The default is to use the bg or background to run the job: i.e., not to queue it at all. Information about queues will be included elsewhere.

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

$ ls
AIM-isosurface_C1_1.00E-03.grid data-summary.data               unsat_m2_x1_A.out               w-C3--exp-analysis.dat
AIM-isosurface_C2_1.00E-03.grid unsat_m2_x1-data-summary.data   unsat_m2_x1_ISA-GRID.mom        w-C4--exp-analysis.dat
AIM-isosurface_C3_1.00E-03.grid unsat_m2_x1.cks                 unsat_m2_x1_ISA.mom             w-H1--exp-analysis.dat
AIM-isosurface_C4_1.00E-03.grid unsat_m2_x1.log                 unsat_m2_x1_atoms.ISA           w-H2--exp-analysis.dat
AIM-isosurface_H1_1.00E-03.grid unsat_m2_x1.out                 w-C1--exp-analysis.dat          warnings.log
AIM-isosurface_H2_1.00E-03.grid unsat_m2_x1_1.00E-03_iso.grid   w-C2--exp-analysis.dat

The CamCASP output is in file unsat_m2_x1.out. Check it to see if the ISA calculation has converged appropriately (search for CONVERG and the second hit should be the one you want. It is very important that this is the case or else the results of the ISA-Pol calculation will be meaningless. Here's the relevant piece of the output file unsat_m2_x1.out (vim line numbers are in the first column):

   2  CHECKStockholderConverged T
   1  Stockholder convergence on main_iteration :           36
1746  BEGIN TEST SHAPE-FUNCTIONS: CONVERGED NORMAL-ITERATIONS
   1 Index Label     W0-norm    W-Norm    Charge    Delta  Conv?  KLa-BASIS KLa-GRID 
   2 =================================================================================
   3      1 H1         0.705     0.705     0.295  0.99829E-09 Y    0.0131488 0.0107864
   4      2 C1         6.358     6.358    -0.358  0.39663E-11 Y    0.0457978 0.0359672
   5      3 C2         5.922     5.922     0.078  0.96834E-12 Y    0.0336793 0.0255089
   6      4 C3         5.922     5.922     0.078  0.96845E-12 Y    0.0336793 0.0255089
   7      5 C4         6.358     6.358    -0.358  0.39664E-11 Y    0.0457978 0.0359672
   8      6 H2         0.705     0.705     0.295  0.99829E-09 Y    0.0131488 0.0107864
   9 =================================================================================
  10  Residual charge =    0.031409151
  11  KL(BASIS) =    0.185251722     KL(GRID) =    0.144524932
  12 =================================================================================
  13  END TEST SHAPE-FUNCTIONS: CONVERGED NORMAL-ITERATIONS

The calculation has converged in 36 iterations. The large residual charge of 0.031 indicates we have an inadequate basis set for the atomic bases. This should ideally be less than 0.001, but for this example it is OK. Another important metric is the Kullback-Liebler information loss KL(GRID). But this is something you can monitor *between* runs: the lower it is the better the solution (usually).

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. Here we have a 0.03 charge error which is a little large, but which is typical of aug-cc-pVDZ main basis sets. In this case, simply using the aug-cc-pVTZ main basis reduces the charge conservation error to 0.0009.
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.
unsat_m2_x1-A.basis : Contains the molecular basis set from the Psi4 calculation.

Here's what we are going to do:

Rename the ISA solution in unsat_m2_x1_psi4/OUT directory to unsat_m2_x1_psi4/OUT_0_ISA.
Create a new Cluster file for the polarizabilities calculation.
Run this to create a new CamCASP (*.cks) command file.
Replace the existing *.cks file with the one from step (2).
Link the ISA solution file unsat_m2_x1_atoms.ISA from unsat_m2_x1_psi4/OUT/ to unsat_m2_x1_psi4.
Re-run the runcamcasp.py command with the --restart option to tell it to use existing files.
The new output will be placed in unsat_m2_x1_psi4/OUT which is why we renamed the old OUT directory in step (0).
Rename the ISA-Pol output directory unsat_m2_x1_psi4/OUT as unsat_m2_x1_psi4/OUT_1_ISApol.

All of these steps can (and should) be automated.

Step (0) is trivial. Let's look at steps (1) to (4) now.

Creating the ISA-Pol input files

To create the CamCASP command files we start with a Cluster file. I find it best to simply create a directory unsat_m2_x1_psi4/clt_files/ where I do my work on the cluster files. So in this directory copy the old cluster file and make a couple of changes:

Title ...
...
...
...

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. ( this tells CamCASP to use the pre-computed MO vectors )
   
   #method isa-pol-from-isa-A-restart ( this is the ISA-Pol method file we need ) 
   
End

Finish

Name this new cluster file unsat_m2_x1_pol.clt and place it in unsat_m2_x1_psi4/clt_files/. 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.

Now run this through cluster and copy the *cks file to the work directory:

$ cd unsat_m2_x1_psi4/clt_files/
$ ls
unsat_m2_x1_pol.clt
$ cluster < unsat_m2_x1_pol.clt
$ ls
unsat_m2_x1.cks     unsat_m2_x1.xyz
$ cp unsat_m2_x1.cks ..
$ cd ..

Now link the ISA solution to the work directory:

$ pwd
unsat_m2_x1_psi4
$ ln -s OUT_0_ISA/unsat_m2_x1_atoms.ISA .

And get back to the top-level directory:

$ cd ..
$ ls
unsat_m2_x1_psi4
unsat_m2_x1_psi4.clt

Now you should have

The ISA-Pol *.cks file in the work directory (unsat_m2_x1_psi4)
The ISA solution *.ISA file linked to the work directory

Running the ISA-Pol calculation

Now all you need to do is run runcamcasp.py once more, but this time with the --restart option:

$ pwd
unsat_m2_x1_psi4
$ runcamcasp.py --clt unsat_m2_x1_psi4.clt -d unsat_m2_x1_psi4 --restart unsat_m2_x1

The calculation will take about 2 minutes and you will then have these files in unsat_m2_x1_psi4/OUT:

$ cd unsat_m2_x1_psi4/OUT
$ ls
data-summary.data                     unsat_m2_x1_ISA-GRID.mom              w-C1--exp-analysis.dat                w-H2--exp-analysis.dat
unsat_m2_x1-data-summary.data         unsat_m2_x1_ISA-GRID_f11_NL4_fmtA.pol w-C2--exp-analysis.dat                warnings.log
unsat_m2_x1.cks                       unsat_m2_x1_ISA-GRID_f11_NL4_fmtB.pol w-C3--exp-analysis.dat
unsat_m2_x1.log                       unsat_m2_x1_atoms_001.ISA             w-C4--exp-analysis.dat
unsat_m2_x1.out                       unsat_m2_x1_lim2.0_4.0_p2000_f11.p2p  w-H1--exp-analysis.dat

Before doing anything else, move the OUT/ directory to something more descriptive in name:

$ pwd
unsat_m2_x1_psi4
$ mv OUT OUT_1_ISApol_from_0

Here I've include the string 'ISApol_from_0' to remind me that this ISA-Pol calculation has been restarted from the ISA solution in OUT_0_ISA.

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

The ISA-Pol non-local polarizabilities

The output is in file unsat_m2_x1.out. If all has gone well you should see these lines:

...
...
 26 BEGIN Polarizability      
  25   DIST-ALG ISA-GRID   ( this sets the ISA-GRID-based partitioning )
  24   ! Sq-freq  0.0      ( use this if you want only static pols )
  23   Quad   10           ( and this is you want 1+10 static+freq-dependent pols )
  22                       ( the settings in QUAD determine the frequencies )
  21 
  20   Invert No
  19   Spherical
  18   Rank 4
  17   Calculate only total and distributed polarizabilities and perturbations
  16   Print only static total pols upto rank 2
  15   Pert-file DEFAULT
  14   Pol-file  DEFAULT       
  13 
  12 
  11   ISA-OPTIONS  TAIL-FIX  OFF
  10   INTEGRATION-OPTIONS  GRID-ALGORITHM = 1
   9    Numerical Integrals: local atomic grid only
   8 END
   7  Finished reading in the POLARIZABILITY block
   6  Now for the calculations...
   5       The perturbation matrix is written in file:unsat_m2_x1_lim2.0_4.0_p2000_f11.p2p
   4       The polarization matrix in format A is written in file:unsat_m2_x1_ISA-GRID_f11_NL4_fmtA.pol
   3       The polarization matrix in format B is written in file:unsat_m2_x1_ISA-GRID_f11_NL4_fmtB.pol
   2  options  RANDOM  4.0000    2.0000    1.0000   seed =   1   
   1  Lattice created for unsat_m2_x1 
987   Type of lattice = RANDOM
   1  Number of lattice points =         2000
 28 
  27  End of Module LATTICE
  26  =====================
  25 
  24        Polarizability calculations for w^2 =   0.0000000000000000
  23        i.e., w =    0.0000000000000000
  22 
  21  Multipole moment tensor form: SPHERICAL
  20  Here preface all terms by Q_
  19  See The Theory of Intermolecular Forces by Stone
  18  l=0: 00
  17  l=1: 10  11c  11s
  16  l=2: 20  21c  21s  22c  22s
  15  l=3: 30  31c  31s  32c  32s  33c  33s
  14  l=4: 40  41c  41s  42c  42s  43c  43s  44c  44s
  13 
  12  Polarizability calculation
  11     Spherical tensor form
  10 
   9     Origin is (   0.0000,   0.0000,   0.0000) BOHR
   8 
   7  Order:    0 by    0
   6    0.5590315E-06
   5 
   4  Order:    0 by    1
   3    0.0000000E+00   0.0000000E+00   0.0000000E+00
   2 
   1  Order:    0 by    2
1017   -0.1831507E-02   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  22 
  21  Order:    1 by    1
  20    0.1065772E+03   0.0000000E+00   0.0000000E+00
  19    0.0000000E+00   0.3143528E+02   0.0000000E+00
  18    0.0000000E+00   0.0000000E+00   0.3143528E+02
  17     Isotropic polarizability:  0.5648258E+02
  16   Anisotropic polarizability:  0.7514191E+02
  15 
  14  Order:    1 by    2
  13    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  12    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  11    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  10 
   9  Order:    2 by    2
   8    0.1990004E+04   0.0000000E+00   0.0000000E+00   0.2890542E-06   0.0000000E+00
   7    0.0000000E+00   0.1211361E+04   0.0000000E+00   0.0000000E+00   0.0000000E+00
   6    0.0000000E+00   0.0000000E+00   0.1211361E+04   0.0000000E+00   0.0000000E+00
   5    0.2890543E-06   0.0000000E+00   0.0000000E+00   0.1391085E+03   0.0000000E+00
   4    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.1390965E+03
...
...

The above listing shows the total (molecular) static polarisability tensors at various orders of rank x rank. Units are always atomic units. The order of the terms is given in the comment lines above the polarisability tensors.

After the molecular static polarisability tensors are printed out the static distributed polarisability tensors are printed:

...
...
 28  Multipole moment tensor form: SPHERICAL
  27  Here preface all terms by Q_
  26  See The Theory of Intermolecular Forces by Stone
  25  l=0: 00
  24  l=1: 10  11c  11s
  23  l=2: 20  21c  21s  22c  22s
  22  l=3: 30  31c  31s  32c  32s  33c  33s
  21  l=4: 40  41c  41s  42c  42s  43c  43s  44c  44s
  20 
  19  Distributed polarizability calculation
  18  ALGORITHM: ISA-GRID : ISA partitioning using a numerical grid
  17     Spherical tensor form
  16 
  15  Order:    0 by    0 for sites   H1 (  1) and   H1 (  1)
  14    0.3639433E+00
  13 
  12  Order:    0 by    1 for sites   H1 (  1) and   H1 (  1)
  11    0.1787329E+00   0.0000000E+00   0.0000000E+00
  10 
   9  Order:    0 by    2 for sites   H1 (  1) and   H1 (  1)
   8    0.3686968E-01   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
   7 
   6  Order:    1 by    1 for sites   H1 (  1) and   H1 (  1)
   5    0.2568351E+00   0.0000000E+00   0.0000000E+00
   4    0.0000000E+00   0.2235048E+00   0.0000000E+00
   3    0.0000000E+00   0.0000000E+00   0.2235048E+00
   2     Isotropic polarizability:  0.2346149E+00
   1   Anisotropic polarizability:  0.3333030E-01
...
...

For a molecule with $n$ sites there are $n^2$ pairs of sites, and so you have a lot of polarisability tensors printed out.

Here the format is:

Order:  rank1 by rank2 for sites site1 (index) and site2 (index)
   <polarizability tensor follows>

The site indices are provided just in case you did not label all sites with unique labels.

Subsequently the frequency-dependent distributed polarizabilities are computed, but these are not printed out here: they are instead saved to files *.pol in two formats:

Format A: unsat_m2_x1_ISA-GRID_f11_NL4_fmtA.pol
Format B: unsat_m2_x1_ISA-GRID_f11_NL4_fmtB.pol

Both as ASCII files.

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.

Why localise?: Because most codes cannot use non-local polarisability tensors. That's the short answer. And for most systems, non-locality is not really needed, that is, unless you are interested in charge flow in a system, or unless (as we found in 2010 and later) you have long wavelength fluctuations that can cause dramatic changes in the dispersion energy.

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
usage: localize.py [-h] [--verbosity VERBOSITY | --verbose | --quiet] [--axes AXES] [--sites SITES] [--format {A,OLD,old,B,NEW,new}] [--polfile POLFILE]
                   [--limit LIMIT] [--wsmlimit WSMLIMIT] [--hlimit HLIMIT] [--isotropic] [--subdir SUBDIR] [--model MODEL] [--loc {LW,LS}]
                   [--cutoff CUTOFF] [--weight {0,1,2,3,4,5,6}] [--weightcoeff WEIGHTCOEFF] [--svd SVD] [--force {loc,refine,disp} [{loc,refine,disp} ...]]
                   [--clean] [--cleanall] [--norefine] [--nodisp] [--debug] [-i INFILE] [--polfile-prefix POLFILE_PREFIX] [-o OUTFILE]
                   name

...
...
many more lines of output!
...
...

Important

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

Localization steps

Here's how we will perform the localization (see the water example for more details) for the HC$_4$H molecule.

Make a directory in which the localisation will be performed, and link the ISA-Pol output directory here as OUT.

$ pwd
unsat_m2_x1_psi4
$ mkdir loc_pol
$ cd loc_pol
$ ln -s OUT_1_ISApol_from_0 OUT
$ ls
OUT

Create the input files

Now we need a Cluster file for the localisation step. This is created from the Cluster file/s you already have: The Run-Type is Localization.

unsat_m2_x1_loc.clt

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 which I have names unsat_m2_x1_loc.clt"

$ cluster < unsat_m2_x1_loc.clt
$ ls
OUT                      unsat_m2_x1.prss         unsat_m2_x1.xyz          unsat_m2_x1_loc.clt
unsat_m2_x1.ornt         unsat_m2_x1.sites        unsat_m2_x1_casimir.prss

Local Axes

For all but isotropic local polarisability models you will need to have a local axis file ready. This is a file that contains the definitions of the local axes for each atom in the molecule. Here you need unique names for the atoms or it will not work. The whole idea of local axes is so that symmetries are enforced. You may even make approximate symmetries exact, or make atoms that are manifestly not symmetric be treated as if they were - usually with non-sensical results. So make this file with care.

For isotropic models you do not need a local axis frame.

Using local axes can make the localisation process more data-efficient and the solution can have fewer artefacts.

For more information about these axes see the Orient and CamCASP user guides.

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

Define the $z$-axis using a direction vector generated from site positions. H1, C1 and C2 have their $z$-axes pointing one way, and C3, C4 and H2 have it pointing the opposite way.
Define the $x$-axis using the global $X$-axis. But again, three sites have it one way and the other three, the other way. This enforces the mirror-symmetry of the molecule.

$ more unsat_m2_x1.axes
Axes
  H1 z from C1 to H1 x  Global  X
  C1 z from C1 to H1 x  Global  X
  C2 z from C1 to H1 x  Global  X
  C3 z from C4 to H2 x  Global -X
  C4 z from C4 to H2 x  Global -X
  H2 z from C4 to H2 x  Global -X
End

Run the Localize code

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

$ localize.py --axes unsat_m2_x1.axes --limit 3 --wsmlimit 3 --hlimit 3 -d L3 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.
-d L3 : Do the localisation in the L3 directory. This keeps things tidy.

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.

You should get an output to screen that looks like this:

$ localize.py --axes unsat_m2_x1.axes --limit 3 --wsmlimit 3 --hlimit 3 -d L3 unsat_m2_x1
Splitting the point-to-point polarizability file. Please wait ... 
Using p2p file OUT/unsat_m2_x1_lim2.0_4.0_p2000_f11.p2p
Splitting point-to-point polarizabilities for unsat_m2_x1
Frequencies: 0 to 10
Number of points in grid: 2000
Total p2p pols per frequency: 2001000
Split p2p pols in files with prefix unsat_m2_x1_010
done

Localisation settings for unsat_m2_x1
Sub-directory    L3
Axes file:       unsat_m2_x1.axes
Pol file format: NEW
Limit:           3
WSM-Limit:       3
H-Limit:         3
Isotropic?:      
Model file:      unsat_m2_x1.pdef
Pol Cutoff:      0.0001
Loc algorithm:   LW
Weight:          3
Weight coeff:    0.001
SVD threshold:   0.0
NoRefine?:       False

Using polarizability file OUT/unsat_m2_x1_ISA-GRID_f11_NL4_fmtB.pol
Splitting OUT/unsat_m2_x1_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 unsat_m2_x1.axes
Refining the local polarizabilities
Creating new polarizability model definition file unsat_m2_x1.pdef
000 001 002 003 004 005 006 007 008 009 010 
Refinement finished
Refined localized polarizabilities, static and at imaginary frequency,
are in unsat_m2_x1_ref_wt3_L3_0f10.pol
Calculating the dispersion coefficients ...   done
Dispersion coefficients are in unsat_m2_x1_ref_wt3_L3_casimir.out.
The dispersion potential, in Orient form, is in unsat_m2_x1_ref_wt3_L3_C12.pot.

Let's see what we have:

$ ls
L3                       unsat_m2_x1.prss         unsat_m2_x1_001.p2p      unsat_m2_x1_005.p2p      unsat_m2_x1_009.p2p
OUT                      unsat_m2_x1.sites        unsat_m2_x1_002.p2p      unsat_m2_x1_006.p2p      unsat_m2_x1_010.p2p
unsat_m2_x1.axes         unsat_m2_x1.xyz          unsat_m2_x1_003.p2p      unsat_m2_x1_007.p2p      unsat_m2_x1_casimir.prss
unsat_m2_x1.ornt         unsat_m2_x1_000.p2p      unsat_m2_x1_004.p2p      unsat_m2_x1_008.p2p      unsat_m2_x1_loc.clt
$ cd L3
OUT                                 unsat_m2_x1_L3_0f10.pol             unsat_m2_x1_ref_wt3_L3_001.out      unsat_m2_x1_ref_wt3_L3_007.data
unsat_m2_x1.axes                    unsat_m2_x1_NL4_000.pol             unsat_m2_x1_ref_wt3_L3_001.pol      unsat_m2_x1_ref_wt3_L3_007.out
unsat_m2_x1.ornt                    unsat_m2_x1_NL4_001.pol             unsat_m2_x1_ref_wt3_L3_002.data     unsat_m2_x1_ref_wt3_L3_007.pol
unsat_m2_x1.pdef                    unsat_m2_x1_NL4_002.pol             unsat_m2_x1_ref_wt3_L3_002.out      unsat_m2_x1_ref_wt3_L3_008.data
unsat_m2_x1.prss                    unsat_m2_x1_NL4_003.pol             unsat_m2_x1_ref_wt3_L3_002.pol      unsat_m2_x1_ref_wt3_L3_008.out
unsat_m2_x1.sites                   unsat_m2_x1_NL4_004.pol             unsat_m2_x1_ref_wt3_L3_003.data     unsat_m2_x1_ref_wt3_L3_008.pol
unsat_m2_x1_L3_000.out              unsat_m2_x1_NL4_005.pol             unsat_m2_x1_ref_wt3_L3_003.out      unsat_m2_x1_ref_wt3_L3_009.data
unsat_m2_x1_L3_001.out              unsat_m2_x1_NL4_006.pol             unsat_m2_x1_ref_wt3_L3_003.pol      unsat_m2_x1_ref_wt3_L3_009.out
unsat_m2_x1_L3_002.out              unsat_m2_x1_NL4_007.pol             unsat_m2_x1_ref_wt3_L3_004.data     unsat_m2_x1_ref_wt3_L3_009.pol
unsat_m2_x1_L3_003.out              unsat_m2_x1_NL4_008.pol             unsat_m2_x1_ref_wt3_L3_004.out      unsat_m2_x1_ref_wt3_L3_010.data
unsat_m2_x1_L3_004.out              unsat_m2_x1_NL4_009.pol             unsat_m2_x1_ref_wt3_L3_004.pol      unsat_m2_x1_ref_wt3_L3_010.out
unsat_m2_x1_L3_005.out              unsat_m2_x1_NL4_010.pol             unsat_m2_x1_ref_wt3_L3_005.data     unsat_m2_x1_ref_wt3_L3_010.pol
unsat_m2_x1_L3_006.out              unsat_m2_x1_casimir.prss            unsat_m2_x1_ref_wt3_L3_005.out      unsat_m2_x1_ref_wt3_L3_0f10.pol
unsat_m2_x1_L3_007.out              unsat_m2_x1_ref_wt3_L3_000.data     unsat_m2_x1_ref_wt3_L3_005.pol      unsat_m2_x1_ref_wt3_L3_C12.pot
unsat_m2_x1_L3_008.out              unsat_m2_x1_ref_wt3_L3_000.out      unsat_m2_x1_ref_wt3_L3_006.data     unsat_m2_x1_ref_wt3_L3_casimir.data
unsat_m2_x1_L3_009.out              unsat_m2_x1_ref_wt3_L3_000.pol      unsat_m2_x1_ref_wt3_L3_006.out      unsat_m2_x1_ref_wt3_L3_casimir.out
unsat_m2_x1_L3_010.out              unsat_m2_x1_ref_wt3_L3_001.data     unsat_m2_x1_ref_wt3_L3_006.pol

An explanation:

The *.p2p files in the loc_pol directory (from where we issued the localize.py command) are the point-to-point polarizabilities taken from OUT/unsat_m2_x1_lim2.0_4.0_p2000_f11.p2p and split into individual files for each frequency. These are large files and when you are all done, delete them.
The actual localisation was performed in L3/. And there you will find:
- Outputs of the LS or LW localisation done using Orient. These are in files of the form unsat_m2_x1_L3_nnn.out. Where nnn = 000 (static polarizabilities), 001, 002, etc. (for the non-static polarizabilities).
- Non-Local polarisability tensors in files of the form unsat_m2_x1_NL4_nnn.pol
- PFIT Input files for the refinement step which are of the form: unsat_m2_x1_ref_wt3_L3_nnn.data. In these you will find the local polarizabilities taken from the Orient localisation output.
- PFIT outputs in the files unsat_m2_x1_ref_wt3_L3_nnn.out. More on this below.
- Local polarisability tensors in the files unsat_m2_x1_ref_wt3_L3_nnn.pol.
- Casimir input file: unsat_m2_x1_ref_wt3_L3_casimir.data
- Casimir output file: unsat_m2_x1_ref_wt3_L3_casimir.out
- Dispersion Model for the dimer (the default is to assume a dimer, but this can be changed) in: unsat_m2_x1_ref_wt3_L3_C12.pot

Trouble Shooting

Check to see if the localisation was sensible

There is no fool-proof way to do this, but large errors in the refinement (particularly for L3 models) are usually an indication of a problem. Additionally, always check for the total molecular polarizabilities: you have a reference calculation in the ISA-Pol output. Compare this with the total molecular polarisability tensor printed out in the refinement step. Use the static polarizabilities here.

First: did the refinement work well? Here's the relevant part of file unsat_m2_x1_ref_wt3_L3_nnn.out at line 375:

375 Residuals               per cent of range               
  1 Maximum      0.00023428       0.745
  2 Minimum     -0.00037284      -1.186
  3 R.m.s.       0.00001300       0.041
  4 Sum of squared residuals =         0.00034
  5 
  6 Parameter                Fitted      Anchor    Difference   Penalty  
  7 H1_10_10_A               3.03480     3.37695    -0.34214    9.43765E-06
  8 H1_10_20_A              -5.57960    -5.10261    -0.47699    8.41530E-06
  9 H1_10_30_A               6.78994     6.98081    -0.19087    7.32551E-07
 10 H1_11c_11c_A             0.67375     0.67085     0.00291    5.82854E-09
 11 H1_11c_21c_A             0.14666     0.14368     0.00298    8.71495E-09
 12 H1_11c_31c_A            -1.02830    -1.03108     0.00277    3.72334E-09
 13 H1_11s_11s_A             0.67298     0.67085     0.00213    3.12973E-09
 14 H1_11s_21s_A             0.14664     0.14368     0.00296    8.61171E-09
 15 H1_11s_31s_A            -1.02822    -1.03108     0.00285    3.94894E-09
 16 H1_20_20_A              10.74893    10.77427    -0.02534    5.48535E-09
 17 H1_20_30_A              28.73506    28.50924     0.22581    6.26599E-08
 ...
 ...

The R.M.S. error is computed against the reference point-to-point polarizabilities. A value of 0.041% of the range (the range of point-to-point polarisability values) is good.

Usually we do not see the Fitted values differ by much from the Anchors: this is the case here. The Anchors have been obtained from localisation of the ISA-Pol non-local polarizabilities, in this case, using LW localization.

Also check the total molecular (static) polarisability tensor from the refinement and ISA-Pol. Here is the reference ISA-Pol data from

1005  Polarizability calculation
   1     Spherical tensor form
   2 
   3     Origin is (   0.0000,   0.0000,   0.0000) BOHR
   4 
   5  Order:    0 by    0
   6    0.5590315E-06
   7 
   8  Order:    0 by    1
   9    0.0000000E+00   0.0000000E+00   0.0000000E+00
  10 
  11  Order:    0 by    2
  12   -0.1831507E-02   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  13 
  14  Order:    1 by    1
  15    0.1065772E+03   0.0000000E+00   0.0000000E+00
  16    0.0000000E+00   0.3143528E+02   0.0000000E+00
  17    0.0000000E+00   0.0000000E+00   0.3143528E+02
  18     Isotropic polarizability:  0.5648258E+02
  19   Anisotropic polarizability:  0.7514191E+02
  20 
  21  Order:    1 by    2
  22    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  23    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  24    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00
  25 
  26  Order:    2 by    2
  27    0.1990004E+04   0.0000000E+00   0.0000000E+00   0.2890542E-06   0.0000000E+00
  28    0.0000000E+00   0.1211361E+04   0.0000000E+00   0.0000000E+00   0.0000000E+00
  29    0.0000000E+00   0.0000000E+00   0.1211361E+04   0.0000000E+00   0.0000000E+00
  30    0.2890543E-06   0.0000000E+00   0.0000000E+00   0.1391085E+03   0.0000000E+00
  31    0.0000000E+00   0.0000000E+00   0.0000000E+00   0.0000000E+00   0.1390965E+03
  32

And here is the tensor from file L3/unsat_m2_x1_ref_wt3_L3_000.out:

503 Total molecular polarizability
  1       00          z           x           y           20          21c         21s         22c         22s
  2      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000
  3      0.00000   106.48910     0.00000     0.00000    -0.00000     0.00000     0.00000     0.00000     0.00000
  4      0.00000     0.00000    31.40565     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000
  5      0.00000     0.00000     0.00000    31.40425     0.00000     0.00000     0.00000     0.00000     0.00000
  6      0.00000    -0.00000     0.00000     0.00000  2035.58242     0.00000     0.00000     0.00000     0.00000
  7      0.00000     0.00000     0.00000     0.00000     0.00000  1211.48110     0.00000     0.00000     0.00000
  8      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000  1211.04284     0.00000     0.00000
  9      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000   137.72288     0.00000
 10      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000   137.78565
 11 Mean polarizability       =  56.43300
 12 Polarizability anisotropy =  75.08415

The mean polarisability and polarisability anisotropy are very close. This is good.
The differences in the quadrupole-quadrupole polarisability tensors are larger ( $\alpha_{20,20}$ is 1990 a.u. in the reference, but 2035 a.u. in the refinement.). But is something that I would usually look into, possibly by altering the weights to see if we can get a better match with different weights.

Changing the refinement weights

Weights and anchors are used to constrain the refinement process to remain close to the computed localised polarizabilities, which define the anchor values. In the absence of any weights/anchors, the refinement usually strays into unphysical values and leads to useless models.

By weakening the weights, the refinement is allowed to make bigger changes to the polarizabilities. But this may not be a good idea.

The ambiguity in the weights exposes a weakness in this approach to obtaining local polarizabilities: they are not well-defined in this approach. The problem is that they are not well defined in any approach as the only well-defined quantities are the non-local polarizabilities. And even these depend on the partitioning method used.

So here, we aim for simplicity: choose the simplest model that works.

The refinement proceeds by using the localised polarizabilities as anchor points and assigns weights to each anchor. The higher the weight, the closer the refined solution will be to the anchor value.

The default weight is type 3 with weight coefficient 0.001: $w(t,u) = \frac{0.001}{1 + \alpha_{tu}^2}$

For example, to suppress refinement we could change the weight coefficient to 1.0 using --weightcoeff 1.0:

$ localize.py --axes unsat_m2_x1.axes --limit 3 --wsmlimit 3 --hlimit 3 -d L3_wt3-c1.0  --weight 3 --weightcoeff 1.0  unsat_m2_x1

This turns out to be a bad idea for this linear system. Instead --weightcoeff 0.01 might be a better compromise. For details of the weights implemented, see the CamCASP User's Guide.

PDef File

The definition of the polarisability model is in the *.pdef file, here this is unsat_m2_x1.pdef. This file is created for you, but you can modify it and if you do, specify this file in the localisation using the --pdef <file name> (or --model <file name>). Here's how this file can look:

Polarizabilities
     H1  H1    10   10 =      H1_10_10_A
     H1  H1    10   20 =      H1_10_20_A
     H1  H1    10   30 =      H1_10_30_A
     H1  H1   11c  11c =    H1_11c_11c_A
     H1  H1   11c  21c =    H1_11c_21c_A
     H1  H1   11c  31c =    H1_11c_31c_A
     H1  H1   11s  11s =    H1_11s_11s_A
     H1  H1   11s  21s =    H1_11s_21s_A
     H1  H1   11s  31s =    H1_11s_31s_A
     ...
     ...
     H2  H2   COPY   H1  H1

     C1  C1    10   10 =      C1_10_10_A
     C1  C1    10   20 =      C1_10_20_A
     C1  C1    10   30 =      C1_10_30_A
     C1  C1   11c  11c =    C1_11c_11c_A
     C1  C1   11c  21c =    C1_11c_21c_A
     ...
     ...

The format for entries is:

     Site1  Site2   t1   t2 =  variable name

where the sites are indicated by Site1/2 which would be the same for localised polarizabilities. The rank of the particular term is given by (t1,t2), and the variable name in which the term is stored is given after the equals sign. If sites are equavilant, the polarizabilities from one site pair are COPYed onto another site pair.

The automatically generated definition file does not always capture all symmetries. For example, in this example, because of cylindrical symmetry, we should expect the (11c 11c) term to be the same as the (11s 11s) term. But if you look at the file unsat_m2_x1_ref_wt3_L3_000.out you see that the c terms are only approximately the same as the corresponding s terms. See (A) and (B) for example.

381 Parameter                Fitted      Anchor    Difference   Penalty
  1 H1_10_10_A               3.03480     3.37695    -0.34214    9.43765E-06
  2 H1_10_20_A              -5.57960    -5.10261    -0.47699    8.41530E-06
  3 H1_10_30_A               6.78994     6.98081    -0.19087    7.32551E-07
  4 H1_11c_11c_A             0.67375     0.67085     0.00291    5.82854E-09 <--- (A)
  5 H1_11c_21c_A             0.14666     0.14368     0.00298    8.71495E-09
  6 H1_11c_31c_A            -1.02830    -1.03108     0.00277    3.72334E-09
  7 H1_11s_11s_A             0.67298     0.67085     0.00213    3.12973E-09 <--- (B)
  8 H1_11s_21s_A             0.14664     0.14368     0.00296    8.61171E-09
  9 H1_11s_31s_A            -1.02822    -1.03108     0.00285    3.94894E-09

We can impose an exact symmetry using this using the modification:

Polarizabilities
     H1  H1    10   10 =      H1_10_10_A
     H1  H1    10   20 =      H1_10_20_A
     H1  H1    10   30 =      H1_10_30_A
     H1  H1   11c  11c =    H1_11c_11c_A
     H1  H1   11c  21c =    H1_11c_21c_A
     H1  H1   11c  31c =    H1_11c_31c_A
     H1  H1   11s  11s =    H1_11c_11c_A
     H1  H1   11s  21s =    H1_11c_21c_A
     H1  H1   11s  31s =    H1_11c_31c_A
     ...

Where I have replaced all variable names with s in them with their equivalents with c. In VIM, %s/s_/c_/ will do the substitution. Let's fix this as save the new PDEF file in loc_pol/unsat_m2_x1_L3.pdef.

Now, re-run the localisation using

$ localize.py --axes unsat_m2_x1.axes --limit 3 --wsmlimit 3 --hlimit 3 -d L3_pdef --pdef unsat_m2_x1_L3.pdef unsat_m2_x1

The result should be quite close to the previous one (with the default PDEF file) in this case, but the number of free variablesin the polarizability tensors are lower. And this is always a good thing. In some cases imposing proper symmetries can lead to much improved refinement of the polarisability tensors.

One important result of using a good *.pdef file is that the dispersion model can be dramatically simplified. For example, the dispersion model with the standard (automatically generated) PDEF file has more than 4000 lines in it, but the same with the above PDEF file has less than 1500 lines.

Of course, if simplicity in the dispersion model is your goal, then you may be better off with isotropic models.

Reducing the Rank

You can reduce the rank of the model in two ways:

Compute a high-ranking model and then modify it to keep only terms less than or equal to the rank you have chosen. So, in this case, to get a rank 1 (L1 == dipole-dipole) model, we could retain only localised dipole-dipole terms for the localised polarizability tensors and only terms for $C_6$ in the dispersion model that come from $lk = 00, 20$ and $j = 0, 2, 4$.
Alternatively, you could re-do the localisation using a lower maximum rank.

For example, to get a L1 or dipole-dipole model you could localise using this command:

$ localize.py --axes unsat_m2_x1.axes --limit 1 --wsmlimit 1 --hlimit 1 -d L1 --pdef unsat_m2_x1_L1.pdef unsat_m2_x1

where I have created a model (PDEF) file unsat_m2_x1_L1.pdef that includes only dipole-dipole terms of the right symmetries:

$ cat unsat_m2_x1_L1.pdef
Polarizabilities
     H1  H1    10   10 =      H1_10_10_A
     H1  H1   11c  11c =    H1_11c_11c_A
     H1  H1   11s  11s =    H1_11c_11c_A
     H2  H2   COPY   H1  H1

     C1  C1    10   10 =      C1_10_10_A
     C1  C1   11c  11c =    C1_11c_11c_A
     C1  C1   11s  11s =    C1_11c_11c_A
     C4  C4   COPY   C1  C1

     C2  C2    10   10 =      C2_10_10_A
     C2  C2   11c  11c =    C2_11c_11c_A
     C2  C2   11s  11s =    C2_11c_11c_A
     C3  C3   COPY   C2  C2

End

How does this compare with the L3 model?

Look at the output of the static polarizability refinement in unsat_m2_x1_ref_wt3_L1_000.out:

Parameter values:
  1            1.62676935     H1_10_10_A
  2            0.91823304     H1_11c_11c_A
  3           15.03146481     C1_10_10_A
  4           10.14612118     C1_11c_11c_A
  5           36.12743430     C2_10_10_A
  6            5.21366335     C2_11c_11c_A
Residuals               per cent of range
Maximum      0.00186568       5.877
Minimum     -0.00149335      -4.704
R.m.s.       0.00010859       0.342
Sum of squared residuals =         0.02360

Parameter                Fitted      Anchor    Difference   Penalty
H1_10_10_A               1.62677     3.25841    -1.63164    2.29164E-04
H1_11c_11c_A             0.91823     0.64484     0.27340    5.27932E-05
C1_10_10_A              15.03146    22.62315    -7.59169    1.12389E-04
C1_11c_11c_A            10.14612     9.60689     0.53923    3.11676E-06
C2_10_10_A              36.12743    27.65259     8.47484    9.38046E-05
C2_11c_11c_A             5.21366     5.87367    -0.66001    1.22708E-05

Sum of penalties =     5.03538E-04
Writing ASCII file unsat_m2_x1_ref_wt3_L1_000.pol
Done

Total molecular polarizability
      00          z           x           y
     0.00000     0.00000     0.00000     0.00000
     0.00000   105.57134     0.00000     0.00000
     0.00000     0.00000    32.55604     0.00000
     0.00000     0.00000     0.00000    32.55604
Mean polarizability       =  56.89447
Polarizability anisotropy =  73.01530

First of all, the RMS residual (as a % of range) is now nearly 10 times larger than the L3 case at 0.342%. This should not be a surprise as the model has only dipole-dipole terms and so cannot reproduce all of the point-to-point polarizabilities as well as it could when we include terms to octopole-octopole.
Then the difference between anchor and fitted values is much larger than what we saw in the L3 case. This is usually troubling and it happens because the model is now really messing around with the available terms to best fit the point-to-point polarizabilities. We could suppress this by increasing the weight coefficient from the default of 0.001 to, say, 0.01. Let's see what this gives:

$ localize.py --axes unsat_m2_x1.axes --limit 1 --wsmlimit 1 --hlimit 1 --weight 3 --weightcoeff 0.01 -d L1_wt3-c0.01 --pdef unsat_m2_x1_L1.pdef unsat_m2_x1

Localisation settings for unsat_m2_x1
Sub-directory    L1_wt3-c0.01
Axes file:       unsat_m2_x1.axes
Pol file format: NEW
Limit:           1
WSM-Limit:       1
H-Limit:         1
Isotropic?:      
Model file:      unsat_m2_x1_L1.pdef
Pol Cutoff:      0.0001
Loc algorithm:   LW
Weight:          3
Weight coeff:    0.01
SVD threshold:   0.0
NoRefine?:       False

Using polarizability file OUT/unsat_m2_x1_ISA-GRID_f11_NL4_fmtB.pol
Splitting OUT/unsat_m2_x1_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 unsat_m2_x1.axes
Refining the local polarizabilities
Using polarizability model definition file unsat_m2_x1_L1.pdef
000 001 002 003 004 005 006 007 008 009 010 
Refinement finished
Refined localized polarizabilities, static and at imaginary frequency,
are in unsat_m2_x1_ref_wt3_L1_0f10.pol
Calculating the dispersion coefficients ...   done
Dispersion coefficients are in unsat_m2_x1_ref_wt3_L1_casimir.out.
The dispersion potential, in Orient form, is in unsat_m2_x1_ref_wt3_L1_C6.pot.

$ cd L1_wt3-c0.01
$ vi unsat_m2_x1_ref_wt3_L1_000.out
...
...
Parameter values:
  1            2.33974282     H1_10_10_A
  2            0.72911359     H1_11c_11c_A
  3           14.03098317     C1_10_10_A
  4           10.29619174     C1_11c_11c_A
  5           36.44224257     C2_10_10_A
  6            5.23574069     C2_11c_11c_A
Residuals               per cent of range
Maximum      0.00187116       5.894
Minimum     -0.00149558      -4.711
R.m.s.       0.00010993       0.346
Sum of squared residuals =         0.02418

Parameter                Fitted      Anchor    Difference   Penalty
H1_10_10_A               2.33974     3.25841    -0.91867    7.26466E-04
H1_11c_11c_A             0.72911     0.64484     0.08428    5.01658E-05
C1_10_10_A              14.03098    22.62315    -8.59217    1.43963E-03
C1_11c_11c_A            10.29619     9.60689     0.68930    5.09298E-05
C2_10_10_A              36.44224    27.65259     8.78965    1.00903E-03
C2_11c_11c_A             5.23574     5.87367    -0.63793    1.14636E-04

Sum of penalties =     3.39086E-03
Writing ASCII file unsat_m2_x1_ref_wt3_L1_000.pol
Done

Total molecular polarizability
      00          z           x           y
     0.00000     0.00000     0.00000     0.00000
     0.00000   105.62594     0.00000     0.00000
     0.00000     0.00000    32.52209     0.00000
     0.00000     0.00000     0.00000    32.52209
Mean polarizability       =  56.89004
Polarizability anisotropy =  73.10385
Finished

Compare the new output with the previous one: The errors made are nearly the same, but now the difference in Anchor and Fitted values is smaller. Maybe we could use --weightcoeff 0.1 to see if we can get similar errors with less shift in the polarizabilies. It turns out 0.1 may be a little too large in this case. But this is the kind of investigation that needs to be done. Get the most accurate model with the least differences in Fitted and Anchor values, i.e., make the least amount of change.

Isotropic Polarizability Models

Important

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

$ localize.py --isotropic --limit 3 --hlimit 3 -d L3iso unsat_m2_x1

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

BUT: for a system as anisotropic as this one, forcing the model to be isotropic can result in not very convincing results.

Here as the static refined isotropic results from L3iso/unsat_m2_x1_ref_wt3_L3iso_000.out

139 Residuals               per cent of range               
  1 Maximum      0.00707370      22.494
  2 Minimum     -0.00794595     -25.267
  3 R.m.s.       0.00121722       3.871
  4 Sum of squared residuals =         2.96474
  5 
  6 Parameter                Fitted      Anchor    Difference   Penalty  
  7 H1_1_iso_A              -0.26777     1.57288    -1.84065    9.75255E-04
  8 H1_2_iso_A               4.26801     3.98743     0.28057    4.65810E-06
  9 H1_3_iso_A             -60.15342   -62.81350     2.66008    1.79297E-06
 10 C1_1_iso_A              15.28063    13.79531     1.48533    1.15320E-05
 11 C1_2_iso_A              72.74769    65.41175     7.33594    1.25747E-05
 12 C1_3_iso_A             133.76498   132.05362     1.71136    1.67940E-07
 13 C2_1_iso_A               8.47141    12.87356    -4.40215    1.16231E-04
 14 C2_2_iso_A             -25.83199   -34.57907     8.74708    6.39347E-05
 15 C2_3_iso_A             864.05001  1064.38780  -200.33779    3.54263E-05
  ...
  ...
258 Total molecular polarizability
  1       00          z           x           y           20          21c         21s         22c         22s     
  2      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000 
  3      0.00000    46.96854     0.00000     0.00000    -0.00000     0.00000     0.00000     0.00000     0.00000 
  4      0.00000     0.00000    46.96854     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000 
  5      0.00000     0.00000     0.00000    46.96854     0.00000     0.00000     0.00000     0.00000     0.00000 
  6      0.00000    -0.00000     0.00000     0.00000  1821.79692     0.00000     0.00000     0.00000     0.00000 
  7      0.00000     0.00000     0.00000     0.00000     0.00000  1391.93954     0.00000     0.00000     0.00000 
  8      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000  1391.93954     0.00000     0.00000
  9      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000   102.36740     0.00000 
 10      0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000     0.00000   102.36740 
 11 Mean polarizability       =  46.96854
 12 Polarizability anisotropy =   0.00000

Compare these total polarizabilities with those from the anisotropic model.

Likewise, the dispersion model too is unsatisfactory. It is often better to construct an isotropic dispersion model from an anisotropic one by retaining only the rank zero terms.

Dispersion Models

Here is what the dispersion model from the L3_pdef calculation looks like:

  15 !  Localisation settings for unsat_m2_x1
  14 !  Axes file:       unsat_m2_x1.axes
  13 !  Pol file format: NEW   
  12 !  Limit:           3
  11 !  WSM-Limit:       3
  10 !  H-Limit:         3
   9 !  Isotropic?:      False
   8 !  Model file:      unsat_m2_x1_L3.pdef
   7 !  Pol Cutoff:      0.0001
   6 !  Loc algorithm:   LW
   5 !  Weight:          3     
   4 !  Weight coeff:    0.001 
   3 !  SVD threshold:   0.0   
   2 !  NoRefine?:       False 
   1 ! 
16   ! Dispersion coefficients for unsat_m2_x1 (hartree bohr^n)
   1   H1  H1             C6             C7             C8             C9             C10            C11            C12
   2     00   00   0  0.8650165         0.0          10.67602         0.0         -76.88432         0.0         -2631.487
   3     00   10   1     0.0         -3.126294         0.0         -10.87203         0.0          695.5964
   4     00   20   2  0.4195764         0.0          9.648186         0.0         -42.33255         0.0         -3022.498
   5     00   30   3     0.0         -2.247714         0.0         -7.804549         0.0          298.5705
   6     00   40   4     0.0            0.0          6.388538         0.0         -26.95332         0.0         -1926.729
   7     00   50   5     0.0            0.0            0.0          2.927911         0.0          26.53243
   8     00   60   6     0.0            0.0            0.0            0.0         -72.12813         0.0         -723.5023
   9     10   00   1     0.0         -3.126294         0.0         -10.87203         0.0          695.5964
  10     10   10   0     0.0            0.0         -3.203316         0.0          49.95388         0.0         -647.3020
  11     10   10   2     0.0            0.0          10.25061         0.0         -142.7254         0.0          1726.139
   ...
   ...

There will be terms like these for every distinct pair of types: these are the symmetry-distinct sites listed in the Cluster input file. For this system we have only types H1, C1 and C2, so the dispersion model will contain terms for (H1, H1), (H1, C1), (H1, C2), (C1, C1), (C1, C2) and (C2, C2).

Each line has the format:

   lk(a) lk(b)  j   C6  C7  C8  C9  C10  C11  C12 (all in atomic units)

Very few codes can handle the anisotropic dispersion models (Orient can), so you may want to make it isotropic. There are two ways of doing this:

Take the above anisotropic model and keep only isotropic terms (with terms: 00 00 0)
Use the --isotropic flag and create an isotropic model.

These are not equivalent! Often the first approach is better as the model is less likely to distort to accommodate any strong anisotropies that may be present in the system.

How to get good models

Some tips:

Improve the basis sets: The aug-cc-pVDZ basis is too small to correctly represent responses beyond dipole-dipole on H-atoms, and quadrupole-quadrupole on heavier atoms. And these too will be underestimated.
Choose a good set of local axes.
Use a good PDEF file with all symmetries made manifest.
Make sure the ISA solution is a good one: look for a low KL-divergence and small charge violations.
Localise to a rank 3 anisotropic model (if you can) and then reduce rank by eliminating higher ranking terms.
For the dispersion model: construct one using a rank 3 anisotropic model, and then discard the anisotropies. This is often better than trying to construct an isotropic model directly, particularly if the system has strong anisotropies.

In the above example of a dispersion model you will have noticed that the C10 and C12 coefficients for (H1, H1) are negative. This is likely because the aug-cc-pVDZ basis is just too small. Here are the same dispersion coefficients computed with the aug-cc-pVTZ main basis set:

! Dispersion coefficients for unsat_m2_x1 (hartree bohr^n)
  H1  H1             C6             C7             C8             C9             C10            C11            C12
    00   00   0  0.8124537         0.0          10.95922         0.0          32.02953         0.0         -618.7762
    00   10   1     0.0         -2.838402         0.0         -24.27693         0.0          160.8397
    00   20   2  0.3970648         0.0          10.72614         0.0          4.612501         0.0         -1925.300
    00   30   3     0.0         -2.054181         0.0         -8.999696         0.0          127.6601
    00   40   4     0.0            0.0          5.451760         0.0         -39.96477         0.0         -1486.017
    00   50   5     0.0            0.0            0.0          5.144957         0.0          44.85830
    00   60   6     0.0            0.0            0.0            0.0         -57.50496         0.0         -600.5629

The C6 and C8 terms are similar to those from the aug-cc-pVDZ calculation, but the C10 term is now positive and the C12 term is, while still negative, much smaller in magnitude. Typically, as the basis sets get larger, the terms tend to be more well-behaved.

Dispersion Models for hetero-mers (A..B complexes)

See this tutorial for dispersion models for A..B complexes.

Some steps will need to be done by hand, but these are easy to automate if needed.

Damping the dispersion models

This takes us beyond the scope of this tutorial. In fact you might ask how to damp the polarisation model too. More on this later, but read our published papers on ISA-Pol and Reg-SAPT(DFT) for information.

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

Distributed polarizabilities using the ISA-Pol algorithm

Codes needed

Methods

The ISA-Pol method

HC$_4$H ISA calculation

HC$_4$H ISA-Pol calculation

Creating the ISA-Pol input files

Running the ISA-Pol calculation

The ISA-Pol non-local polarizabilities

Localizing the non-local polarizabilities

Localization steps

Create the input files

Local Axes

Run the Localize code

Trouble Shooting

Changing the refinement weights

PDef File

Reducing the Rank

Isotropic Polarizability Models

Dispersion Models

How to get good models

Dispersion Models for hetero-mers (A..B complexes)

Damping the dispersion models

ISA-Pol with Cartesian auxiliary basis set

ISA-Pol with ISA A

ISA-Pol with A+DF ISA