Navigation:

# Distributed density-overlap

The distributed density overlap integrals are used as a means to fit the short-range energies to a distributed form. If the density is partitioned as $$

\rho

^{A(r) = \sum_a \rho}a(r),

$$ where $a$ is a site label (this need not be an atomic site) in the molecule A, and likewise for molecule B, then the density overlap integral for A and B is defined as $$

S

^{{AB}({\bf R}) = \int \rho}A(r) \rho^B(r) dr,

$$ where $\bf R$ describes the relative position and orientation of the molecules. Using the partitioned densities we can write the above as $$

S

^{{AB}({\bf R}) = \sum_a \sum_b \int \rho}a(r) \rho^{b(r) dr = \sum_a \sum_b S}{ab}({\bf R}),

$$ where the $S^{ab}({\bf R})$ are the //distributed// density overlap integrals. Clearly the values of these integrals depends on the type of distribution scheme used.

In the past we used to use a density-fitting based distribution scheme. This algorithm was primitive but worked (see our papers from 2007 and 2008), but necessitated the use of relatively small auxiliary basis sets or else the partitioning could be nonsensical. In fact, in our 2016 JCTC on pyridine we show just how nonsensical the domains could be. The GEM approach of Cisneros may be a better method and could be relatively easily coded in CamCASP, but this has not yet been done. In any case we have a much better approach in the BS-ISA algorithm present in CamCASP.

## ISA-based distribution

The BS-ISA algorithm is a basis-space implementation of Lillestolen and Wheatley's iterated stockholder atoms (ISA) algorithm for defining an atom-in-a-molecule (AIM). The unique features of this method that make it ideal for the distribution we need are:

- It is unique.
- Large basis sets can be used.
- The atoms are maximally spherical.
- They reflect chemical charge movement.
- The domains are consistent with the ISA multipoles, and ISA polarizabilities. So we get a consistent short- and long-range distribution.

That's the sale's pitch. How can the distribution be done?

In principle it is easy: get ISA solutions for the monomers. You also get the partitioned densities. So all that needs to be done is to calculate the overlap integrals of the partitioned densities.

However, this gives you the integrals for one dimer orientation only. For other orientations you need to rotate the ISA solutions appropriately. CamCASP can do this rotation for SPHERICAL GTOs only, so this forces the ISA basis sets to use SPHERICAL GTOs.

Additionaly, we know that using SPHERICAL GTOs in the standard ISA basis sets leads to atomic densities with small negative terms in the density tails. So the smaller of the overlap integrals are non-sensical.

We get around this problem using a GRID-based partitioning. This works but adds to the computational cost, though it is not too computationally expensive.

Finally, as you may also compute energies along side the density overlap integrals, and as these energies will use the molecular auxiliary basis set, and since quantities associated with this basis will need to be rotated, the molecular AUX basis will also need to use SPHERICAL GTOs only.

One day I will remove this restriction and allow rotation of CARTESIAN GTOs in CamCASP.

Here are the steps we will use:

- Perform an ISA-A calculation with CamCASP. Save the ISA solution that will be in a *.ISA file.
- Do this for both molecules. So you will have, say, files molA.ISA and molB.ISA
Now we will use the

*dens-ovr-integrals*method (in*$CAMCASP/methods/*) to do the following:- Define a CamCASP file with the two monomers using MC-type basis sets.
- The main basis is the one used to get the molecular orbitals. We assume this has been done.This will likely be a basis of SPHERICAL GTOs and will be a monomer (MC) basis.
- The auxiliary basis is a SPHERICAL GTO basis that will also be a monomer basis. This basis needs to use spherical GTOs as we will rotate objects computed in this basis using Wigner D-matrices. CamCASP cannot as yet rotate Cartesian GTOs.
The atom-aux basis will have the same $s$-functions as the basis used for the ISA solutions. The higher angular momenta functions can differ. This basis too needs to use SPHERICAL GTOs as here too we will perform rotations.

- Restart the molA ISA from molA.ISA
- Test the solution by computing the ISA multipoles. These should agree with previously calculated ISA multipoles or else something may be wrong.
- Do the same for molB.
Use the

*ENERGY-SCAN*module to calculate the distributed density overlap integrals for specified dimer geometries that you need to specify. (CamCASP can generate these)We optionally calculate $E^{(1)}_{\rm elst}$ and $E^{(1)}_{\rm exch}$. The former is calculated almost for free but the latter will cost more as the dimer density-fitting will be done.

- Define a CamCASP file with the two monomers using MC-type basis sets.
- You will get two output files:
energy_file.dat : This contains the energies $E^{(1)}_{\rm elst}$ and $E^{(1)}_{\rm exch}$ (if computed) and the total overlap integral. This integral is not an energy and so does not have units, but CamCASP treats it as an energy and converts (so to speak) from a.u. to the defined units.

- overlap_file.dat : This contains the distributed density overlap integrals. The structure of the file is easy to figure out when you see the file.

### The ''dens-ovr-integrals'' method

The method file for the above is *$CAMCASP/methods/dens-ovr-integrals*. The contents of this file may change but at present it contains the following:

CamCASP-commands ! The grid matters here. These parameters should be sufficient. BEGIN GRID Molecule A Angular 200 Radial 100 END BEGIN GRID Molecule B Angular 200 Radial 100 END ! Set this to FALSE as the auxiliary basis sets are type-MC SET DF REDO-DF-ON-ROTATION False END ! Perform the OO-type density fitting for each of the monomers BEGIN DF Molecule A Type OO Eta = 0.0 Lambda = 0.0 Gamma = 0.0 Print only normalization constraints Solver LU END BEGIN DF Molecule B Type OO Eta = 0.0 Lambda = 0.0 Gamma = 0.0 Print only normalization constraints Solver LU END ! Setting for the ISA calculation. But we do not perform the ISA here: ! We read a pre-computed solution from file using the RESTART block. Begin Stockholder Molecule A DF = Drho W-INIT = ONE-GTO ALPHA0 = 1.0 ISA-Algorithm A Zeta = 1.0 DF-PARAMETERS Lambda = 1000.0 Solver LU Convergence Convergence-Type W EPS-Norm = 1.0e-10 EPS-Q = 1.0e-4 Max-Iter = 60 W-Damping = 0.0 W-Mix-Fraction = 0.0 Skip-Iterations = 20 W-Eps = 0.17 S-Block-Only Couple Switch-On-Eps-Norm = 1.0e-04 Positive-W Lambda = 0.0001 for Max-Alpha = 2.0 AUTO Tail-Iterations = 30 End W-TAILS Func = 1 R1-Multiplier = 1.5 R2-Multiplier = 2.5 Fit-Type = 3 W-Tests END Restart File DEFAULT Iterate No Test Yes End End ! Useful to calculate multipole to make sure it was all OK ! Note that without iterations these may not come out accurate. Begin Multipoles Molecule A DF Type ISA Rank 4 End Begin Multipoles Molecule A DF Type ISA-GRID Rank 4 End ! Same for the second molecule Begin Stockholder Molecule B DF = Drho W-INIT = ONE-GTO ALPHA0 = 1.0 ISA-Algorithm A Zeta = 1.0 DF-PARAMETERS Lambda = 1000.0 Solver LU Convergence Convergence-Type W EPS-Norm = 1.0e-10 EPS-Q = 1.0e-4 Max-Iter = 60 W-Damping = 0.0 W-Mix-Fraction = 0.0 Skip-Iterations = 20 W-Eps = 0.17 S-Block-Only Couple Switch-On-Eps-Norm = 1.0e-04 Positive-W Lambda = 0.0001 for Max-Alpha = 2.0 AUTO Tail-Iterations = 30 End W-TAILS Func = 1 R1-Multiplier = 1.5 R2-Multiplier = 2.5 Fit-Type = 3 W-Tests END Restart File DEFAULT Iterate No Test Yes End End Begin Multipoles Molecule B DF Type ISA Rank 4 End Begin Multipoles Molecule B DF Type ISA-GRID Rank 4 End ! Integral switches in the energy modules Set E1elst Integral switch = 0 End Set E1exch ! Overlap integrals Algorithm DF without constraints Integral switch = 0 End ! This is where the distributed density-overlap parameters are set. ! We use the ISA-based partitioning on a density obtained using DF without ! constraints. The ISA-based partitioning algorithm can be done ! in basis-space (not a good idea) or in real-space. The latter is ! much more accurate, but takes much longer. ! We set the grid parameters for the integration. These parameters ! are known to work well. Set Dist-Dens-Overlap Type ISA DF-Type OO without constraints ISA-DIST-TYPE GRID-ALG1 DF-INTEGRAL-SWITCH = 0 Integration-Grid Radial 40 Angular 200 ISA-Options TAIL-FIX = False Print-Results No End ! Energy-scan for E(1)elst and overlap only. The former cost nothing ! more and server as a check on the results. ! ! If you also include E1exch then you can make a plot of ! E1exch vs the total overlap 'energy' in the file energy_file.dat ! and see if there is a correlation. There should be! ! ! E1elst will not be believable as we have used monomer basis sets ! in this calculation, and we need at least a DC auxiliary basis to ! get this energy correctly. Begin Energy-Scan Energy E1elst & Overlap ( E1exch ) Units Bohr ! The dimer geometries should be in a file called dimers.geom: #include dimers.geom End End-CamCASP-commands

This method implements the list of steps described in the previous section.

The associated Cluster file is:

Title Template Cluster file for a distributed density overlap integral calculation ! We assume the following: ! (1) You have got MO files for the two molecules in their reference ! geometries. These are the geometries defined below. ! (2) You have for the ISA solution files (*.ISA) for the two molecules. ! (3) You have a list of dimer geometries in a file called: ! dimers.geom ! (4) At present, the AUX and ATOMAUX basis sets both need to have ! SPHERICAL GTOs. ! (5) ATOMAUX must have the same s-functions as used to get the ISA solutions. ! But other functions may differ. ! (6) AUX need not be the same as ATOMAUX. ! ! Procedure: ! * Edit this file as needed. ! * Have the items listed above ready. ! * Run: ! cluster < this_file.clt ! * If it ran correctly you will get many files, but we need only the *.cks file. ! * See if the *.cks file is correct. Edit it as needed and then run: ! camcasp < *.cks ! Global Units Bohr Degrees Overwrite Yes End Molecule A ... End Molecule B ... End Run-Type SAPT(DFT) Molecules A and B ! Change these basis sets as appropriate ! But all must be type MC Main-Basis aVDZ Type MC Aux-Basis aVTZ Type MC Spherical No-ISA-Basis ! Make sure that the s-functions in this basis are the same as ! those used to get the ISA solutions. AtomAux-Basis aVTZ Type MC Spherical Use-ISA-Basis ISA-Basis set2 Min-S-exp-H = 0.2 ( convergence tends to be better with this setting ) Func PBE0 ( other functionals possible, but this is probably the best ) AC MULTIPOLE ( [TH | LINEAR | GRAC] [CS00 | MULTIPOLE | LB94] : CS00 wit h NWChem) Kernel ALDA+CHF ( not needed for an ISA calculation ) SCF-code <name of code> ( DALTON2015, NWCHEM, etc ) MO-FILE-A <name of *.movecs> Format ASCII ( Set this to the acutal MO file ) MO-FILE-B <name of *.movecs> Format ASCII ( Set this to the acutal MO file ) File-Prefix A-B-ovr ( change this to something reasonable ) #METHOD dens-ovr-integrals ( this is the method file ) End Finish