
<<TableOfContents()>>

Navigation:
  * [[AJMPublic/camcasp|CamCASP]]

= 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:
 1. Perform an ISA-A calculation with CamCASP. Save the ISA solution that will be in a *.ISA file.
 1. Do this for both molecules. So you will have, say, files molA.ISA and molB.ISA
 1. Now we will use the ''dens-ovr-integrals'' method (in ''$CAMCASP/methods/'') to do the following:
   1. Define a CamCASP file with the two monomers using MC-type basis sets. 
     1. 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.
     1. 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.
     1. 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. 
   1. Restart the molA ISA from molA.ISA
   1. Test the solution by computing the ISA multipoles. These should agree with previously calculated ISA multipoles or else something may be wrong.
   1. Do the same for molB.
   1. 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)
   1. 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.
 1. You will get two output files: 
   1. 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. 
   1. 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:

[[attachment:dens-ovr-integrals]]

{{{
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:

[[attachment:dens-ovr-integrals.clt_tmpl]]

{{{
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
}}}