Navigation:

• CamCASP main page

The ENERGY-SCAN module

In CamCASP it is possible to compute the first-order energies $E^{(1)}_{\rm elst}$ and $E^{(1)}_{\rm exch}$ for a large number of dimer configurations without repeating the wavefunction calculation on the monomers. This is possible only because the first-order energies can be accurately evaluated in a monomer basis set, that is, without the presence of the mid-bond and far-bond functions that are necessary for converging the second-order interaction energy components. In this tutorial we will cover how these energies can be computed using the ENERGY-SCAN module, and additionally, how the distributed density-overlap integrals can be computed.

Commands in ENERGY-SCAN

ENERGY-SCAN commands

BEGIN SCAN/ENERGY-SCAN
  {PROBE || SCAN} [molecule || atom] <molname> +++
      WITH [<molname> [[and] CHARGE <charge>]]
  UNITS [A || ANG || ANGSTROM || BOHR || AU || A.U.]  +++
      [DEG || DEGREE || DEGREES || RAD || RADIAN || RADIANS]
      NOTE: UNITS must come before POINTS/RANDOM
  CENTRE/CENTER PROBE [on] [COM]/[XYZ <x> <y> <z>]
  REGULARIZE [E1] [[and] E2] [with] ETA [=] <reg_eta>
  {ENERGIES || ENERGY} [DISP/DISPERSION/E2DISP] +++
   [&/and/,] [IND/INDUCTION/E2IND] [&/and/,] [ELST/E1ELST/ELECTROSTATICS] +++
   [&/and/,] [E1EXCH] [&/and/,] [E2EXIND] & [E2EXDISP] & [OVERLAP] [NOTHING]
  LATTICE
  POINTS
    [TRANSLATIONS || TRANSLATIONS-ONLY]
    [SKIP-FIRST-COLUMN]
    x1 y1 z1    [angle1 nx1 ny1 nz1]
    x2 y2 z2    [angle2 nx2 ny2 nz2]
    etc
  ---  <----end POINTS block with three hyphens.
  RANDOM
    POINTS [=] <number of points>
    DRMIN  [=] <dRmin>
    DRMAX  [=] <dRmax>
    RADIAL [points] [=] <number of radial points per angular configuration>
    SKIP   [=] <number to skip> [angular] [configurations]
    WRITE  <filename> [APPEND]
  ---
  ENERGY-FILE [PREFIX] <file prefix> [[in || using] STYLE <style index>]
  OVERLAP-FILE [PREFIX] <file prefix> [[in || using] STYLE <style index>]
  DEBUGGING T/TRUE/ON or F/FALSE/OFF
  {QUIET || VERBOSE}
END

Files needed for the examples

• Pyridine dimer geometry file: pyr2-scan.geom

• Pyridine ISA file: pyridine_daTZ_aQZset2_A+DFz0.9.ISA

First-order energies: $E^{(1)}_{\rm elst}$ and $E^{(1)}_{\rm exch}$

MC main basis and MC auxiliary basis

Important points to get right:

• REDO-DF-ON-ROTATION False : Needed to avoid re-calculation of the monomer density-fitting as the monomers are rotated.

• Auxiliary basis sets should be Spherical in this case as CamCASP can only rotate integrals that use spherical angular momenta.

• Integral switch = 0: force 1-electron integrals to use density-fitting. Results in far more accurate $E^{(1)}_{\rm elst}$ energies. But with an MC-MC type of basis this energy will always be inaccurate!

• DF-TYPE-MONOMER OO and DF-TYPE-DIMER OO: For first-order energies only integrals in the OO space are needed. This causes the DF-INTEGRAL module to change its defaults and use the smaller space. Results in much faster execution.

Important command for the MC-MC case

! This is important as it save a lot of computer time
! As both the main and auxiliary basis sets are type MC, 
! all objects, including the DF solution matrices can be 
! rotated with the molecules. So no re-calculations are needed.
 
SET DF
  REDO-DF-ON-ROTATION   False
END
 
! Various DF parameters. 
 
SET DF
  Molecule pyridine1
  Type OO
  Eta = 0.0
  Lambda = 0.0
  Gamma = 0.0
  Print only normalization constraints
  Solver LU
END
SET DF
  Molecule pyridine2
  Type OO
  Eta = 0.0
  Lambda = 0.0
  Gamma = 0.0
  Print only normalization constraints
  Solver LU
END
 
! Only occupied-occupied (OO) type integrals are needed.
! This saves computer time. The default is to compute 
! integrals in the OV and VV spaces too. Not needed for
! first-order energies.
 
Set DF-INTS
  DF-TYPE-MONOMER  OO
  DF-TYPE-DIMER    OO
End
 
! Set the integral switches in the energy modules. 
! Switch = 1 causes the code to compute the overlap and nuclear
! integrals without density-fitting (i.e. they are more accurate), but this introduces
! a large error in E(1)elst as there is no error cancellation. 
! For some reason this error cancellation is needed when we use an MC-type
! aux basis. So here we use switch = 0 to force CamCASP to use the 1-electron
! integrals with density-fitting. This results in more accurate E(1)elst, but
! bear in mind that E(1)elst with an MC-type main basis will never be 
! accurate enough.
 
Set E1elst
  Integral switch = 0
End
Set E1exch
  Integral switch = 0
End
 
Begin Energy-Scan
  Probe pyridine1 with pyridine2
  Energy  E1elst & E1exch
  Units Bohr
  ! Debug
  #include pyr2-scan.geom
End

MC main basis and DC auxiliary basis

This is the more accurate way of calculating the $E^{(1)}_{\rm elst}$ energy. For reasons we still do not understand, the DC-type auxiliary basis results in far more accurate electrostatic energies. The $E^{(1)}_{\rm exch}$ energy seems more insensitive to this. As the density-fitting will be repeated for each dimer configuration, this scan will be much more computationally expensive than the MC/MC case. However it will also be more accurate.

Distributed density-overlap integrals

The distributed density-overlap integrals needed for the density-overlap model can be calculated using the ENERGY-SCAN module using the keyword:

  Energy Overlap

This can be done using two algorithms:

• The density-fitting-based partitioning, and
• the ISA-based partitioning.

The latter is much more stable and may be used with extensive basis sets. The former is much faster but must be used with relatively small auxiliary basis sets.

DF-based distribution

Basis restrictions:

• Small auxiliary basis sets: TZVPP is OK but does lead to artifacts.
• The auxiliary basis must used spherical GTOs.
• Both the main and auxiliary basis sets must be type MC.

ISA-based distribution

Basis restrictions:

• At present, the s-functions in the auxiliary basis must be those required for the ISA. This will soon change (it may have already changed, but I need to check the code!)
• The auxiliary basis must used spherical GTOs.
• Both the main and auxiliary basis sets must be type MC.

Other important points to get right:

• REDO-DF-ON-ROTATION False : Needed to avoid re-calculation of the monomer density-fitting as the monomers are rotated.

• Auxiliary basis sets should be Spherical in this case as CamCASP can only rotate integrals that use spherical angular momenta.

• Integral switch = 0: force 1-electron integrals to use density-fitting. Results in far more accurate $E^{(1)}_{\rm elst}$ energies. But with an MC-MC type of basis this energy will always be inaccurate!

• DF-TYPE-MONOMER OO and DF-TYPE-DIMER OO: For first-order energies only integrals in the OO space are needed. This causes the DF-INTEGRAL module to change its defaults and use the smaller space. Results in much faster execution.

Accuracy : MC-DC versus MC-MC and the choice of integral switch

Here is some data for the thiophene dimer in ten random dimer orientations:

Thiophene2_MC_DC_aux_comparison.dat

#        Reference energies    ------------------MC main basis-------------------------
#      ---augA-Sadlej-MC+----  -----aTZ DC-aux sw=1-------   -----aTZ DC-aux sw=0------
# INDX    E1elst     E1exch    E1elst        E1exch          E1elst         E1exch
1       -0.82874     5.35489 -6.6592046E-01  4.9352442E+00  -2.1537427E-01  4.9184553E+00
2       -9.69716    29.56556 -9.6704016E+00  3.0054717E+01  -8.9839291E+00  3.0025835E+01
3       -6.73208    21.89199 -6.6285621E+00  2.1964070E+01  -5.9417571E+00  2.1931287E+01
4       -1.21089     6.55433 -1.0410530E+00  6.1161696E+00  -5.2299977E-01  6.0962712E+00
5       -2.48791     6.39493 -2.3995603E+00  6.2131489E+00  -2.2424074E+00  6.2058767E+00
6       -1.28705     3.16231 -1.2235851E+00  3.0334854E+00  -1.2353517E+00  3.0317530E+00
7      -23.70028    68.83931 -2.3378553E+01  6.9401212E+01  -2.5423431E+01  6.9406683E+01
8       -1.24095     3.04211 -1.1786025E+00  2.9168411E+00  -1.1988017E+00  2.9152978E+00
9        0.36216     1.18635  4.3647476E-01  1.2357284E+00   5.9249559E-02  1.2331486E+00
10      -1.36398    10.34951 -1.2889322E+00  1.0329523E+01  -4.4949354E-01  1.0303064E+01

Here's the same for the MC-MC case:

Thiophene2_MC_MC_aux_comparison.dat

#        Reference energies    ------------------MC main basis-------------------------
#      ---augA-Sadlej-MC+----  -----aTZ MC-aux sw=1-------   -----aTZ MC-aux sw=0------
# INDX    E1elst     E1exch    E1elst        E1exch          E1elst         E1exch
1       -0.82874     5.35489  4.3197305E+00  4.9326346E+00  -9.0835363E-01  4.9198501E+00
2       -9.69716    29.56556 -2.6562133E+00  3.0045234E+01  -1.1141111E+01  3.0029574E+01
3       -6.73208    21.89199  4.0413476E-02  2.1955698E+01  -7.8247516E+00  2.1935417E+01
4       -1.21089     6.55433  4.1474269E+00  6.1129764E+00  -1.3724563E+00  6.0980295E+00
5       -2.48791     6.39493  2.8214742E+01  6.2149486E+00  -5.0185247E+00  6.2141753E+00
6       -1.28705     3.16231  2.6948757E+01  3.0337061E+00  -3.0910482E+00  3.0362936E+00
7      -23.70028    68.83931  1.1345649E+01  6.9423383E+01  -2.7840501E+01  6.9372187E+01
8       -1.24095     3.04211  2.6857644E+01  2.9170139E+00  -3.0089173E+00  2.9196679E+00
9        0.36216     1.18635  8.4263178E+00  1.2338500E+00   1.8745135E-01  1.2334557E+00
10      -1.36398    10.34951  5.0912047E+00  1.0308759E+01  -3.4663106E-01  1.0303721E+01

Here's the above data analysed:

Conclusions:

• $E^{(1)}_{\rm exch}$ is relatively insensitive to the type of auxiliary basis. The small relative rms errors suggest that these energies can be used fairly reliably with either of the integral switches and either of the MC-DC or MC-MC auxiliary basis sets.

• $E^{(1)}_{\rm elst}$ is very sensitive to the switch. For the MC-DC type of basis, the switch = 1 should be used. Even so, the errors made in this energy are relatively large. However for the MC-MC type of basis, the switch = 0 must be used: errors are enormous with switch = 1.