12 Fragment-Based Methods

12.10 The Adiabatic ALMO-EDA Method and VFB Analysis

(May 16, 2021)

Despite the huge success and usefulness of today’s most popular EDA methods, they still face some limitations in their capabilities. For instance, EDAs are usually performed at complex geometries that are obtained from unconstrained electronic structure calculations (e.g., optimized equilibrium geometries). For strongly interacting systems, close intermolecular contacts driven by POL and particularly CT often result in largely unfavorable FRZ interaction, which offers little physical insights besides indicating obviously substantial intermolecular overlap. Another limitation is that the conventional EDA methods often partitions a “single-point" interaction energy evaluated at a given geometry. Therefore, the influence of FRZ, POL and CT on the structural and vibrational properties of an intermolecular complex cannot be directly characterized.

Recently Mao et al. reformulated the original ALMO-EDA method in an adiabatic picture, 715 Mao Y., Horn P. R., Head-Gordon M.
Phys. Chem. Chem. Phys.
(2017), 19, pp. 5944.
Link
where the term “adiabatic" is borrowed from spectroscopy and indicates that energy differences are evaluated at relaxed geometry on each potential energy surface (PES). In this scheme, the total binding energy (including monomer geometry distortions) is repartitioned into adiabatic FRZ, POL and CT terms:

ΔEbind=ΔEfrz(ad)+ΔEpol(ad)+ΔEct(ad). (12.23)

The adiabatic frozen interaction energy is given by the difference between the energy minimum of the frozen PES (on which the energy of each point is computed using the corresponding frozen wave function) and the sum of fully relaxed, non-interacting fragment energies:

ΔEfrz(ad)=E[𝐏frz(frz)]-AEA(0). (12.24)

Similarly, the adiabatic POL and CT terms can be obtained by performing geometry optimizations on the polarized (SCFMI) and fully relaxed (unconstrained SCF) PESs:

ΔEpol(ad) =E[𝐏pol(pol)]-E[𝐏frz(frz)], (12.25)
ΔEct(ad) =E[𝐏full(full)]-E[𝐏pol(pol)]. (12.26)

With this method, the changes in monomer structures and intermolecular coordinates due to FRZ, POL and CT and the accompanied energetics are provided. Moreover, at the energy minimum (or other stationary points) on each PES, the other properties such as multipole points, vibrational frequencies and intensities can also be computed, therefore the effect of different intermolecular interaction components on them can also be characterized.

The geometry optimization on the frozen PES is facilitated by the analytical gradient of the frozen wave function energy implemented in Q-Chem. As for the geometry optimization on the polarized PES, the nuclear gradient of the SCFMI energy has the same form as that of the full SCF energy if the original ALMO model is used. These analytical gradients can also be used for finite difference calculations of harmonic frequencies by setting IDERIV = 1. We note that the analytical gradients of SCFMI calculations that use FERFs are not available yet, and SCFMI_MODE = 0 is required for computing the forces on the frozen and polarized PESs. Also, the current implementation of this method requires users to perform geometry optimization on the three PESs separately (see the example below) and evaluate the energy components by taking several Q-Chem outputs (including geometry optimizations for the monomers) together, which is probably not so convenient. We look forward to extending the functionality of this method and improving its implementation in the near future.

As we mentioned in 12.7.4, for systems containing radicals of highly symmetric geometries, the frozen wavefunction obtained from concatenating the fragment MOs might be non-unique. In those cases, we recommend the user to set EDA_ALIGN_FRGM_SPIN = 1 or 2 when performing geometry optimization on the frozen PES. The job will then go through the fragment spin alignment procedure in each optimization cycle.

FRZ_GEOM
       Compute forces on the frozen PES.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not compute forces on the frozen PES. TRUE Compute forces on the frozen PES.
RECOMMENDATION:
       Set it to TRUE when optimized geometry or vibrational frequencies on the frozen PES are desired.

POL_GEOM
       Compute forces on the polarized (converged SCFMI) PES.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not compute forces on the polarized PES. TRUE Compute forces on the polarized PES.
RECOMMENDATION:
       Set it to TRUE when optimized geometry or vibrational frequencies on the polarized PES are desired.

Example 12.26  Geometry optimization of the ammonia-borane complex on the fully relaxed, polarized, and frozen potential energy surfaces successively.

$molecule
0 1
--
0 1
H       0.000000     0.000000     0.000000
H       0.000000     0.000000     1.629090
H       1.417687     0.000000     0.814543
N       0.473683    -0.370067     0.814542
--
0 1
N       3.494032    -1.531250     0.814538
H       3.967715    -1.901317    -0.000008
H       2.550028    -1.901319     0.814537
H       3.967715    -1.901317     1.629083
$end

$rem
   JOBTYPE          opt  !optimization on the fully relaxed PES
   GEN_SCFMAN       true
   METHOD           wb97x-d
   BASIS            6-31+g*
   XC_GRID          1
   THRESH           14
   SCF_CONVERGENCE  9
   SCF_GUESS        fragmo
   SYMMETRY         false
   SYM_IGNORE       true
$end

@@@

$molecule
   read
$end

$rem
   JOBTYPE          opt
   POL_GEOM         true !optimization on the polarized PES
   GEN_SCFMAN       true
   METHOD           wb97x-d
   BASIS            6-31+g*
   XC_GRID          1
   THRESH           14
   SCF_CONVERGENCE  9
   SYMMETRY         false
   SYM_IGNORE       true
   SCFMI_MODE       0
$end

@@@

$molecule
   read
$end

$rem
   JOBTYPE          opt
   FRZ_GEOM         true !optimization on the frozen PES
   GEN_SCFMAN       true
   METHOD           wb97x-d
   BASIS            6-31+g*
   XC_GRID          1
   THRESH           14
   SCF_CONVERGENCE  9
   SYMMETRY         false
   SYM_IGNORE       true
   SCFMI_MODE       0
$end

View output

Example 12.27  Geometry optimization of the [Cu(CO)]+ complex on the frozen PES, followed by a frequency calculation which is computed via finite differences.

$molecule
1 1
--
0 1
C         0.0000000000    0.0000000000    1.3792049588
O         0.0000000000    0.0000000000    2.4988670685
--
1 1
Cu        0.0000000000    0.0000000000   -0.9778656750
$end

$rem
   JOBTYPE            opt
   FRZ_GEOM           true
   METHOD             b3lyp
   BASIS              def2-svp
   UNRESTRICTED       false
   SYMMETRY           false
   SYM_IGNORE         false
   IDERIV             1
   FD_MAT_VEC_PROD    false
$end


@@@


$molecule
   read
$end

$rem
   JOBTYPE           freq
   FRZ_GEOM          true
   METHOD            b3lyp
   BASIS             def2-svp
   UNRESTRICTED      false
   SYMMETRY          false
   SYM_IGNORE        false
   IDERIV            1
   FD_MAT_VEC_PROD   false
$end

View output

To further understand the charge-transfer effects in dative complexes, in Q-Chem 5.2.2 and after, one is allowed to separate the overall CT into contributions from forward and backward donations using the variational forward-backward (VFB) approach. 695 Loipersberger M., Mao Y., Head-Gordon M.
J. Chem. Theory Comput.
(2020), 16, pp. 1073.
Link
Such a decomposition is achieved by introducing two additional constrained intermediate states in which only one direction of CT is permitted. These two “one-way” CT states are variationally relaxed such that the associated nuclear forces can be readily obtained. This allows for a facile integration into the adiabatic ALMO-EDA scheme introduced above:

ΔEctf(ad) =E[𝐏ctf(ctf)]-E[𝐏pol(pol)], (12.27)
ΔEctb(ad) =E[𝐏ctb(ctb)]-E[𝐏pol(pol)], (12.28)

and thus the molecular property changes arising from forward and backward donations can be separately assigned. Note that in its Q-Chem implementation, the evaluation of a VFB state always follows a polarization (standard SCFMI) calculation. Also, since the definition of VFB states is based on the generalized SCFMI technique (Sec. 12.7.1), SCFMI_MODE = 1 is required.

VFB_CTA
       Use the Variational Forward-Backward (VFB) approach to obtain “one-way” CT PESs.
TYPE:
       STRING
DEFAULT:
       NONE
OPTIONS:
       FORWARD Allow 12 CT only (1 and 2 are two fragments). BACKWARD Allow 21 CT only.
RECOMMENDATION:
       None

Example 12.28  Geometry optimization on one-side CT surface (2->1) using the variational forward-backward (VFB) approach

$molecule
   0 1
   --
   0 1
   O  -1.551007  -0.114520   0.000000
   H  -1.934259   0.762503   0.000000
   H  -0.599677   0.040712   0.000000
   --
   0 1
   O   1.350625   0.111469   0.000000
   H   1.680398  -0.373741  -0.758561
   H   1.680398  -0.373741   0.758561
$end

$rem
   JOBTYPE           opt
   METHOD            wb97x-d
   BASIS             6-31g
   VFB_CTA           backward
   THRESH            14
   SCF_CONVERGENCE   9
   SYMMETRY          FALSE
   SCF_ALGORITHM     DIIS
   IDERIV            1
   SCFMI_MODE        1
$end

View output

In Q-Chem 5.4 or later, analytical gradients for the polarized and two VFB “one-way” CT states with implicit solvent models PCM and SMD are supported so that one can perform part of the adiabatic ALMO-EDA steps (POL CTf/CTb Full) in solvation environments. To do this, one only needs to set the $rem variable SOLVENT_METHOD to PCM or SMD, which is similar to the usage of ALMO-EDA(solv) (see Sec. 12.7.5). The calculation of analytical forces on the frozen surface with implicit solvents is currently unavailable, and we look forward to enabling that in future releases of Q-Chem.

Example 12.29  Geometry optimization on the polarized surface with SMD solvent model

$molecule
   0 1
   --
   0 1
   H1
   O1 H1 0.95641
   H2 O1 0.96500  H1 104.77306
   --
   0 1
   O2 H2 dist     O1 171.85474 H1 180.000
   H3 O2 0.95822  H2 111.79807 O1 -58.587
   H4 O2 0.95822  H2 111.79807 O1 58.587

   dist = 2.0
$end

$rem
   JOBTYPE           OPT
   METHOD            wB97X-D
   BASIS             cc-pVDZ
   POL_GEOM          TRUE
   THRESH            14
   SCF_CONVERGENCE   9
   MEM_TOTAL         2000
   MEM_STATIC        500
   SCF_GUESS         FRAGMO
   SYMMETRY          FALSE
   SYM_IGNORE        TRUE
   IDERIV            1
   SCFMI_MODE        0
   SOLVENT_METHOD    SMD
$end

$smx
   solvent water
$end

View output