Navigation:

• CamCASP

Distributed polarizabilities using the ISA-Pol algorithm

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

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.

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

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

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:

1. Rename the ISA solution in unsat_m2_x1_psi4/OUT directory to unsat_m2_x1_psi4/OUT_0_ISA.

2. Create a new Cluster file for the polarizabilities calculation.
3. Run this to create a new CamCASP (*.cks) command file.
4. Replace the existing *.cks file with the one from step (2).

5. Link the ISA solution file unsat_m2_x1_atoms.ISA from unsat_m2_x1_psi4/OUT/ to unsat_m2_x1_psi4.

6. Re-run the runcamcasp.py command with the --restart option to tell it to use existing files.

7. The new output will be placed in unsat_m2_x1_psi4/OUT which is why we renamed the old OUT directory in step (0).

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

1. The ISA-Pol *.cks file in the work directory (unsat_m2_x1_psi4)

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

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

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

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.

unsat_m2_x1.axes

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

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.

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:

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

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

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.00063031       2.004
  2 Minimum     -0.00051790      -1.647
  3 R.m.s.       0.00002702       0.086
  4 Sum of squared residuals =         0.00146
  5 
  6 Parameter                Fitted      Anchor    Difference   Penalty  
  7 H1_10_10_A              -0.59202     3.37695    -3.96897    1.26999E-03
  8 H1_10_20_A              -1.65285    -5.10261     3.44976    4.40174E-04
  9 H1_10_30_A               6.43701     6.98081    -0.54380    5.94631E-06
 10 H1_11c_11c_A             0.77522     0.67085     0.10437    7.51241E-06
 11 H1_11c_21c_A             0.05638     0.14368    -0.08730    7.46749E-06
 12 H1_11c_31c_A            -1.02213    -1.03108     0.00895    3.87950E-08
 13 H1_11s_11s_A             0.76875     0.67085     0.09791    6.61107E-06
...
...

The R.M.S. error is computed against the reference point-to-point polarizabilities. A value of 0.086% of the range is good.

BUT, the (10,10) component of the polarizability tensor for H1 starts of as positive (Anchor value 3.37 a.u.) but ends up negative (Fitted value -0.59 a.u.). This indicates a problem. Let's see if we can resolve it by changinf the default weights.

Changing the refinement weights

$ localize.py --axes unsat_m2_x1.axes --limit 3 --wsmlimit 3 --hlimit 3 -d L3_wt3-c0.01  --weight 3 --weightcoeff 0.01  unsat_m2_x1

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              -0.59202     3.37695    -3.96897    1.26999E-03
  2 H1_10_20_A              -1.65285    -5.10261     3.44976    4.40174E-04
  3 H1_10_30_A               6.43701     6.98081    -0.54380    5.94631E-06
  4 H1_11c_11c_A             0.77522     0.67085     0.10437    7.51241E-06. <---- (A)
  5 H1_11c_21c_A             0.05638     0.14368    -0.08730    7.46749E-06
  6 H1_11c_31c_A            -1.02213    -1.03108     0.00895    3.87950E-08
  7 H1_11s_11s_A             0.76875     0.67085     0.09791    6.61107E-06. <---- (B)
  8 H1_11s_21s_A             0.06380     0.14368    -0.07988    6.25149E-06
  9 H1_11s_31s_A            -1.02159    -1.03108     0.00948    4.36026E-08

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.

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

ISA-Pol with Cartesian auxiliary basis set

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.