Differences between revisions 6 and 23 (spanning 17 versions)

Contents

Distributed polarizabilities using the ISA-Pol algorithm
1. Methods
2. The ISA-Pol method
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.

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.

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 it is the basis set that's not great.
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, which, in this case has all atoms with the same $x$-axis, but H1, C1 and C2 have their $z$-axes pointing one way, and C3, C4 and H2 have it pointing the opposite 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 unsat_m2_x1

Here's a brief explanation of the options:

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

Important

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

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

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

And the dispersion model will be in file unsat_m2_x1_ref_wt3_L3_C12.pot

Important

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

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

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

Important

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

ISA-Pol with Cartesian auxiliary basis set

Warning

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

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

ISA-Pol with ISA A

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

Step 1: Perform a standard ISA A calculation with Cartesian GTOs in the AUX basis and spherical GTOs in the ATOM-AUX basis. Use method isa-A. Here it is advisable to have the spherical parts of AUX and ATOM-AUX the same.
Step 2: Perform an ISA-Pol calculation by restarting from this solution. Use method isa-pol-from-isa-A-restart.

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

ISA-Pol with A+DF ISA

This one must be done in two steps:

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

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

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

-  ⇤ ← Revision 6 as of 2021-03-26 17:13:54 → 
  Size: 18692
  Editor: bsw388
  Comment:
+   ← Revision 23 as of 2024-11-09 17:39:51 → ⇥
  Size: 30406
  Editor: apw185
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 1:
-## page was renamed from ajm/camcasp/isa-pol
#acl apw185:read,write,delete Known:read All:
-Line 7:
+Line 4:
-  * [[ajm/camcasp/start|CamCASP]]
+  * [[AJMPublic/camcasp|CamCASP]]
-Line 13:
+Line 10:
-This tutorial applies to CamCASP 6.0.17 and higher. At the time of writing this version of the code is private. If you are interested in these calculations, please contact AJM.
+This tutorial applies to CamCASP 7.x and higher.
-Line 18:
+Line 15:
-  - //ISA-Pol: Distributed polarizabilities and dispersion models from a basis-space implementation of the iterated stockholder atoms procedure//, A. J. Misquitta & A. J. Stone (2018) '''submitted''' [[https://arxiv.org/abs/1806.06737|Version on ArXiv]]

In an earlier tutorial we have described how the CamCASP code can be used to calculate [[ajm:camcasp:pols-cn|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.
+  - //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) [[https://arxiv.org/abs/1806.06737|Version on ArXiv]] [[https://doi.org/10.1007/s00214-018-2371-4|Published version at Theoretical Chemistry Accounts]]

In an earlier tutorial we have described how the CamCASP code can be used to calculate [[AJMPublic/camcasp/pols-cn|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.
-Line 35:
+Line 32:
-The files in $CAMCASP/methods/ will usually be consistent with the CamCASP version, but this wiki page may not be. So the contents of this directory should over-rule any statements made on this page.
+The 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.
-Line 41:
+Line 38:
-  - ISA calculation: See the [[ajm/camcasp/isa|tutorial on the ISA]].
+  - ISA calculation: See the [[AJMPublic/camcasp/isa|tutorial on the ISA]].
-Line 43:
+Line 40:
-If you wish to compare the ISA-Pol results with the earlier constrained-DF (cDF) polarizabilities, then first use the [[ajm:camcasp:pols-cn|tutorial on distributed properties]] to calculate these polarizabilities.
+If you wish to compare the ISA-Pol results with the earlier constrained-DF (cDF) polarizabilities, then first use the [[AJMPublic/camcasp/pols-cn|tutorial on distributed properties]] to calculate these polarizabilities.
-Line 81:
+Line 78:
-  Aux-Basis     aug-cc-pVDZ  Type MC  Cartesian  Use-ISA-Basis  ( this could be spherical )
+  Aux-Basis     aug-cc-pVDZ  Type MC  Spherical  Use-ISA-Basis  ( this could be spherical or Cartesian )
-Line 83:
+Line 80:
-  ISA-Basis     set2  Min-S-exp-H = 0.2  ( This limit leads to better convergence )
+  ISA-Basis     set2  Min-S-exp-H = 0.0  ( This sets the smallest s-exponent for H-atom basis sets )
-Line 88:
+Line 85:
-  ! Program used for the DFT part:
  SCFcode NWChem
+  ! Program used for the DFT part: You could use NWChem, Psi4, DALTON, etc. We will use Psi4.
  SCFcode Psi4
-Line 95:
+Line 92:
-  Memory 4000 MB
+  Memory 8 GB
-Line 113:
+Line 110:
-Run this calculation the usual way:
{{{
  $ runcamcasp.py --scfcode nwchem --nproc 2 --q batch --clt unsat_m2_x1.clt --direct unsat_m2_x1
}}}

{{{#!wiki  important
Important

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

The output will be in ''./OUT/'' and should contain the files (and a few more besides):
{{{
$ ls
AIM-isosurface_C.grid         unsat_m2_x1_1.00E-03_iso.grid w-C--fn1-CONV.dat
AIM-isosurface_H.grid         unsat_m2_x1_ISA.mom           w-H--exp-analysis.dat
DMA_Z0_L5.mom                 unsat_m2_x1_atoms.ISA         w-H--fn1-CONV-TAIL-ITER.dat
DMA_Z4_L5.mom                 unsat_m2_x1_camcasp.out       w-H--fn1-CONV.dat
unsat_m2_x1-data-summary.data w-C--exp-analysis.dat
unsat_m2_x1.ornt              w-C--fn1-CONV-TAIL-ITER.dat
}}}
As usual the CamCASP output is in file ''unsat_m2_x1_camcasp.out''. Check it to see if the ISA calculation has converged appropriately. It is very important that this is the case or else the results of the ISA-Pol calculation will be meaningless.
+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
}}}

{{{#!wiki  important
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).
-Line 140:
+Line 159:
-  * 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.
+  * 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 it is the basis set that's not great.
-Line 151:
+Line 170:
-The Cluster file for this calculation is essentially identical to the previous one, but we include a few different commands:
{{{
RUN-TYPE
+  * ''unsat_m2_x1-A.basis'' : Contains the molecular basis set from the Psi4 calculation.

Here's what we are going to do:
  0. Rename the ISA solution in ''unsat_m2_x1_psi4/OUT'' directory to ''unsat_m2_x1_psi4/OUT_0_ISA''.
  1. Create a new Cluster file for the polarizabilities calculation.
  2. Run this to create a new CamCASP (*.cks) command file.
  3. Replace the existing *.cks file with the one from step (2).
  4. Link the ISA solution file ''unsat_m2_x1_atoms.ISA'' from ''unsat_m2_x1_psi4/OUT/'' to ''unsat_m2_x1_psi4''. 
  5. Re-run the runcamcasp.py command with the ''--restart'' option to tell it to use existing files. 
  6. The new output will be placed in ''unsat_m2_x1_psi4/OUT'' which is why we renamed the old ''OUT'' directory in step (0).
  7. 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
-Line 159:
+Line 199:
-   MO-file  unsat_m2_x1-A-asc.movecs  Format ASCII
+   MO-file  unsat_m2_x1-A-asc.movecs  Format ASCII. ( this tells CamCASP to use the pre-computed MO vectors )
-Line 161:
+Line 201:
-   #method isa-pol-from-isa-A-restart
+   #method isa-pol-from-isa-A-restart ( this is the ISA-Pol method file we need )
-Line 163:
+Line 203:
-END
}}}
Call this new cluster file ''unsat_m2_x1_pol.clt''.
+End

Finish
}}}
Name this new cluster file ''unsat_m2_x1_pol.clt'' and place it in ''unsat_m2_x1_psi4/clt_files/''.
-Line 168:
+Line 210:
-To run the job perform the following steps. Here I have assumed the following directory structure:
{{{
./ISA/
  ISA/OUT/
./ISA-Pol/

The ISA directory is the one that contains the ISA solution from the previous step. 
It is currently called unsat_m2_x1, so rename it:

$ mv unsat_m2_x1 ISA  

In this directory you will find the *.movecs file and in ISA/OUT/ you will find the *.ISA file.
We need these for the ISA-Pol calculation.
There is no need to copy these files, as symbolic links will suffice: 

$ cd ISA-Pol
$ ln -s ../ISA/unsat_m2_x1-A-asc.movecs .
$ ln -s ../ISA/OUT/unsat_m2_x1_atoms.ISA .

Now create the CamCASP command file (and other files - not all needed):
+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
-Line 190:
+Line 216:
-And run the job using a fairly old bash script that was
set up to re-run just the CamCASP code:

$ runcamcasp unsat_m2_x1 -q batch  


And wait for a while...
}}}

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

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

The output is in file ''unsat_m2_x1_camcasp.out'' and the polarizabilities should be present in either one or two formats. If you get two formats, then you will see two ''*.pol'' files as shown above. The ''*_fmtB.pol'' is a new format that I am still working on.
These polarizability files contain the non-local polarizabilities at 11 imaginary frequencies. To localize these polarizabilities we use the techniques described in the [[ajm/camcasp/pols-cn|Distributed Properties tutorial]] and outlined in the next step.
+$ 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 [[ajm/camcasp/pols-cn|Distributed Properties tutorial]] and outlined in the next step.
-Line 217:
+Line 412:
+{{{#!wiki  important
'''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. 
}}}
-Line 220:
+Line 419:
+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!
...
...
-Line 228:
+Line 438:
-Here's how we will perform the localization (see the [[ajm/camcasp/pols-cn|water example for more details]]) for the HC$_4$H molecule. We will work in the ''ISA-Pol'' directory and perform the localization in a sub-directory named ''L3'' (I've also got a sub-directory named ''L2'' in the example below):

[[attachment:Localization steps]]

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

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

[[attachment: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
}}}

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

[[attachment:unsat_m2_x1-loc.clt]]
+==== Localization steps ====

Here's how we will perform the localization (see the [[ajm/camcasp/pols-cn|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''. 

[[attachment:unsat_m2_x1_loc.clt]]
-Line 291:
+Line 492:
-Run Cluster on this file and make a link from the ''OUT/'' directory to this one:
{{{
$ cluster < unsat_m2_x1-loc.clt
$ ln -s ../OUT .
}}}
This will generate a bunch of files. Normally these will not need editing.
+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 =====

{{{#!wiki important
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, which, in this case has all atoms with the same $x$-axis, but H1, C1 and C2 have their $z$-axes pointing one way, and C3, C4 and H2 have it pointing the opposite way. This enforces the mirror-symmetry of the molecule. 

[[attachment: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 =====

Diff for "AJMPublic/camcasp/isa-pol"

Distributed polarizabilities using the ISA-Pol algorithm

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

ISA-Pol with Cartesian auxiliary basis set

ISA-Pol with ISA A

ISA-Pol with A+DF ISA