Navigation:

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.

Important

Reference Geometries The energy scan with module ENERGY-SCAN works on the principle that we computer the molecular orbitals and energies for the molecule in their reference orientations, and then rotate/translate them as needed to compute the interaction energies for the interacting dimer. Consequently, it is important that this reference orientation must be used to define all molecules in the CamCASP command files.

Important

//Integral Accuracy//: When the DF is done accurately - here when the DC-aux basis is used - then it seems that we must use 1-electron integrals computed as accurately as is possible. So we use Integral switch = 1. But if we use a relatively poor DF (as is the case with the MC-MC case when we use the MC-aux basis type), then we need to rely on an error cancellation in the terms that appear in $E^{(1)}_{\rm elst}$. So we need to use Integral switch = 0.

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

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

MC main basis and MC auxiliary basis

Important points to get right:

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

ISA-based distribution

Basis restrictions:

Other important points to get right:

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

Warning

The basis sets used in these calculations are not equivalent: the augA-Sadlej* basis is nearly as good as the Dunning aug-cc-pVTZ basis, but it is not the same. Besides, it is used in the MC+ type, with mid-bonds and far-bonds. So treat the comparison with some caution.

* : augA-sadlej is a basis developed for intermolecular interactions by AJM. It is available in CamCASP 6.x.

Here's the above data analysed:

MC-DC

MC-MC

E(1)elst

E1exch

E(1)elst

E1exch

sw=1

sw=0

sw=1

sw=0

sw=1

sw=0

sw=1

sw=0

rmse

0.14

0.77

0.32

0.32

20.04

1.85

0.32

0.31

mad

0.11

0.61

0.25

0.25

16.09

1.42

0.25

0.25

rel rmse

0.10

0.45

0.04

0.04

13.14

0.77

0.04

0.04

Conclusions: