Contents
- Distributed polarizabilities using the ISA-Pol algorithm
- ISA-Pol with Cartesian auxiliary basis set
Navigation:
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:
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
FinishComments:
- 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.
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
EndHow 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
FinishedCompare 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.5629The 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.