Navigation:

Distributed polarizabilities using the ISA-Pol algorithm

Important

This tutorial applies to CamCASP 6.0.17 and higher. At the time of writing this version of the code is private. If you are interested in these calculations, please contact AJM.

Papers:

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:

For example, an ISA calculation using the 'A' algorithm is a method called isa-A (case sensitive!) and its associated sample Cluster command file is isa-A.clt_tmpl. The method should be altered only if you know what you are doing. The sample Cluster command file gives you guidance on how to construct the Cluster file for the method.

Important

The files in $CAMCASP/methods/ will usually be consistent with the CamCASP version, but this wiki page may not be. So the contents of this directory should over-rule any statements made on this page.

The ISA-Pol method

In principle the ISA-Pol algorithm can be used to calculate distributed polarizabilities in a single calculation, but it may be useful to perform the calculation in two steps, at least until the code gets more robust. These steps are:

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  Cartesian  Use-ISA-Basis  ( this could be spherical )
  AtomAux-Basis aug-cc-pVTZ  Type MC  Spherical  Use-ISA-Basis  ( this MUST be spherical  )
  ISA-Basis     set2  Min-S-exp-H = 0.2  ( This limit leads to better convergence )
  
  Functional PBE0
  ! Type of method for polarizabilities
  Kernel ALDA+CHF
  ! Program used for the DFT part:
  SCFcode NWChem
  
  ! If you already have a *-asc.movecs file then uncomment this line and set the name of 
  ! the file. Use this only if the RUN-TYPE is CamCASP-Only.
  ! MO-File  <name of MO file>  Format ASCII  ( set this to the name of the *-asc.movecs file )
  
  Memory 4000 MB
  
  #method isa-A  ( this will use the isa-A method from the $CAMCASP/methods/ directory )
  
End

Finish

Comments:

This file includes the isa-display template file from the CamCASP library. You will need to copy this from $CAMCASP/data/camcasp-cmnds/ to your working directory. Edit it as you see fit, but only if you know what you are doing! For some guidance on the contents see the ISA tutorial.

Run this calculation the usual way:

  $ runcamcasp.py --scfcode nwchem --nproc 2 --q batch --clt unsat_m2_x1.clt --direct unsat_m2_x1

Important

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

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

$ ls
AIM-isosurface_C.grid         unsat_m2_x1_1.00E-03_iso.grid w-C--fn1-CONV.dat
AIM-isosurface_H.grid         unsat_m2_x1_ISA.mom           w-H--exp-analysis.dat
DMA_Z0_L5.mom                 unsat_m2_x1_atoms.ISA         w-H--fn1-CONV-TAIL-ITER.dat
DMA_Z4_L5.mom                 unsat_m2_x1_camcasp.out       w-H--fn1-CONV.dat
unsat_m2_x1-data-summary.data w-C--exp-analysis.dat
unsat_m2_x1.ornt              w-C--fn1-CONV-TAIL-ITER.dat

As usual the CamCASP output is in file unsat_m2_x1_camcasp.out. Check it to see if the ISA calculation has converged appropriately. It is very important that this is the case or else the results of the ISA-Pol calculation will be meaningless.

Important

TIPS: Things to look for that are hallmarks of a converged ISA calculation are:

  • A relatively small (0.01e or so) charge violation in the first ISA step, before the density-fitting constraint is included. If you find a larger charge violation then something has probably gone wrong.
  • Ideally the shape-functions (tabulated in the w-*-CONV-TAIL-ITER.dat files) should be positive everywhere and exponentially (or nearly so) decaying. Plot them to see if this is the case. Sometimes the shape-functions for the hydrogen atoms may get negative past $6.0$ Bohr. This is usually OK. If you see appreciable wiggles (in a semilog plot of $w(r)$) or the on-set of negative values at small $r$ then something may have gone wrong.

  • Things usually go wrong because the basis set used for the ISA was not flexible enough, or if the W-EPS term was too large. See the ISA tutorial for details.

The ISA solution for this calculation is in file unsat_m2_x1_atoms.ISA. We will need this file in the next stage.

HC$_4$H ISA-Pol calculation

Files needed from the previous step:

The Cluster file for this calculation is essentially identical to the previous one, but we include a few different commands:

RUN-TYPE
   CamCASP-Only  ( use this as you already have the *-asc.movecs file )
   ...
   ...commands as above...
   ... 
   MO-file  unsat_m2_x1-A-asc.movecs  Format ASCII
   
   #method isa-pol-from-isa-A-restart
   
END

Call this new cluster file unsat_m2_x1_pol.clt. The method isa-pol-from-isa-A-restart is also located in $CAMCASP/methods/ and this method contains commands that allow the ISA to restart from a previous ISA solution with the ISA-Pol calculation following this. There is no need to edit this file unless you have edited the isa-A method used in the previous step. Just make sure that the ISA settings in the two methods are the same.

To run the job perform the following steps. Here I have assumed the following directory structure:

./ISA/
  ISA/OUT/
./ISA-Pol/

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

$ mv unsat_m2_x1 ISA  

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

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

Now create the CamCASP command file (and other files - not all needed):

$ cluster < unsat_m2_x1_pol.clt

And run the job using a fairly old bash script that was
set up to re-run just the CamCASP code:

$ runcamcasp unsat_m2_x1 -q batch  


And wait for a while...

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

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

The output is in file unsat_m2_x1_camcasp.out and the polarizabilities should be present in either one or two formats. If you get two formats, then you will see two *.pol files as shown above. The *_fmtB.pol is a new format that I am still working on. These polarizability files contain the non-local polarizabilities at 11 imaginary frequencies. To localize these polarizabilities we use the techniques described in the Distributed Properties tutorial and outlined in the next step.

Localizing the non-local polarizabilities

Details of how this is done can be found in the Distributed Properties tutorial. Here I will be brief.

We use the localize.py code to drive the localization process. This code uses a number of other Fortran programs and accepts a number of options that you can see using

$ localize.py --help

Important

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

Here's how we will perform the localization (see the water example for more details) for the HC$_4$H molecule. We will work in the ISA-Pol directory and perform the localization in a sub-directory named L3 (I've also got a sub-directory named L2 in the example below):

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:

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:

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 and make a link from the OUT/ directory to this one:

$ cluster < unsat_m2_x1-loc.clt
$ ln -s ../OUT .

This will generate a bunch of files. Normally these will not need editing.

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

$ localize.py --axes ../unsat_m2_x1.axes --limit 3 --wsmlimit 3 --hlimit 3 unsat_m2_x1

Here's a brief explanation of the options:

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

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

<note 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:

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:

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.