Analytic Gradients for the Effective Fragment Molecular Orbital Method

Article pubs.acs.org/JCTC

Analytic Gradients for the Effective Fragment Molecular Orbital Method Colleen Bertoni and Mark S. Gordon* Department of Chemistry, Iowa State University, Ames, Iowa 50014, United States S Supporting Information *

ABSTRACT: The analytic gradient for the Coulomb, polarization, exchange-repulsion, and dispersion terms of the fully integrated effective fragment molecular orbital (EFMO) method is derived and the implementation is discussed. The derivation of the EFMO analytic gradient is more complicated than that for the effective fragment potential (EFP) gradient, because the geometry of each EFP fragment is flexible (not rigid) in the EFMO approach. The accuracy of the gradient is demonstrated by comparing the EFMO analytic gradient with the numeric gradient for several systems, and by assessing the energy conservation during an EFMO NVE ensemble molecular dynamics simulation of water molecules. In addition to facilitating accurate EFMO geometry optimizations, this allows calculations with flexible EFP fragments to be performed.

1. INTRODUCTION Many interesting chemical systems involve large molecules (such as protein−ligand complexes and enzyme catalysis) or many molecules (such as chemical reactions in solution). However, it is computationally expensive to perform ab initio calculations on large systems. Several methods have been developed to make such calculations feasible. These include using parametrized classical force fields to model interactions between molecules, hybrid quantum mechanics (QM)/molecular mechanics (MM) methods, and fragmentation schemes that perform ab initio calculations on fragments of a system and then combine the fragment results.1 The effective fragment molecular orbital (EFMO) method was developed to combine the sophisticated semi-classical effective fragment potential (EFP) method2 with the fragment molecular orbital (FMO) method,3 in order to take advantage of the computational efficiency of both.4 The FMO method is a fragmentation method based on a many-body expansion of the energy that has been applied extensively to molecular clusters and biological systems.5 The EFP method is a sophisticated model potential method that is derived from first principles, with no empirically fitted parameters. Fragment geometries in the EFP method are rigid. The EFP method decomposes the interaction energy into five terms: Coulomb, polarization, exchangerepulsion, dispersion, and charge transfer. It has enabled many studies of intermolecular interactions, including solvent effects on chemical processes.6 The original EFMO method combined the fragmentation scheme of the FMO method with just the Coulomb and polarization interaction energy terms of the EFP method. An approximate gradient for the original EFMO method was reported.4 © 2016 American Chemical Society

The EFMO method has recently been greatly improved, by incorporating the EFP dispersion, exchange-repulsion, and charge-transfer interaction terms.7 This improved method was called the fully integrated effective fragment molecular orbital (FIEFMO) method in ref 7. The gradient for the additional terms was not derived or implemented. Hereinafter the FIEFMO method will be referred to simply as the EFMO method, and the original method with just the Coulomb and polarization terms will be referred to as the original EFMO method. There are many motivations for the development of fully analytic gradients. Geometry optimizations of molecules are typically much more accurate, numerically stable, and less timeconsuming with analytic, rather than numeric gradients. Transition-state searches and reaction path following are enabled by analytic gradients, and fully analytic gradients are essential for molecular dynamics (MD) simulations.8 This work presents the derivation and implementation of the gradient terms that are needed to make the original EFMO gradient fully analytic, and the Coulomb, exchange-repulsion, polarization, and dispersion terms that are needed to make the fully integrated EFMO gradient fully analytic. The gradient of the charge-transfer term, usually the least important and most computationally demanding component of the EFP interaction energy,42 has not been derived or implemented, as discussed further in Section 3.3. Since the EFMO analytic gradient involves EFP interaction energy derivatives without assuming the fragments are rigid, an added benefit of the derivation presented here Received: April 4, 2016 Published: July 27, 2016 4743

DOI: 10.1021/acs.jctc.6b00337 J. Chem. Theory Comput. 2016, 12, 4743−4767

Article

Journal of Chemical Theory and Computation

2.5. Localized Molecular Orbital Notation and Definitions. LMOs |l⟩ are related to CMOs by a unitary transformation L

is that it provides insight regarding which EFP interaction energy terms are most important with regard to fragment flexibility. This paper is organized as follows: Section 2 introduces the notation used; Section 3 gives a brief overview of the EFMO energy expression; Section 4 presents the derivation of the EFMO gradient while noting the differences with the EFP gradient; Section 5 discusses the implementation of the EFMO analytic gradient; Section 6 presents test calculations on a variety of systems (a cluster of water molecules, a cluster of water molecules, methanol molecules, dimethyl sulfoxide molecules, and an ionic liquid pair) and discusses the consequent potential energy surfaces; Section 7 presents timing comparisons to FMO gradients. The final section concludes.

occ CMO

∂x

∑ q

=

∑

Llicμi

i

(2.2)

∑mLMO

L is a unitary transformation matrix, LmnLml = δln, calculated by a localization method, such as the Boys method,10 which was originally proposed by Edmiston and Ruedenberg.11,12 2.6. Derivative of a Localized Molecular Orbital Coefficient with Respect to a Nuclear Perturbation, Written in Terms of the Canonical Response Matrix and Localization Response Matrix. Following previous studies that considered perturbed LMOs,13−16 the nuclear derivative of the LMO coefficient is split into a term that includes a localization response matrix (which describes how the localization transform changes with geometry) and a term that includes the canonical response matrix (which describes how the CMOs change with geometry) ∂cμLl ∂x

occ CMO

=

∑ i

∂ (Llicμi) ∂x

=

∑

cμLnvnlx

+

n

where

vxnl

occ+vir CMO

occ CMO

LMO

∑ i

Lli

∑

Uqixcμq

q

(2.3)

is the localization transform response matrix.

3. EFMO METHOD The EFMO method is an integration of the FMO and EFP methods, designed to take advantage of the speed and accuracy of the two methods. The FMO and EFP methods are described briefly in Sections 3.1 and 3.2, respectively. Since the EFMO method has been described previously,7,4 only a brief overview is given in Section 3.3. 3.1. Fragment Molecular Orbital (FMO) Method. Reference 5 provides an excellent review of the FMO method. In general, the system is divided into fragments (monomers) in a chemically sensible way, for example, using common functional groups. Then, the energy of each monomer is calculated in a Coulomb field due to the other monomers. Since the field depends on the electron density of the monomers, the Coulomb field is converged self-consistently. This level of theory is called FMO1. After it has converged, the dimer (pair of fragments) and trimer (set of three fragments) energy may be computed in the self-consistently converged monomer Coulomb field as well. The monomer, dimer, and possibly trimer energies are added together to obtain the total energy for the system. The computational expense increases when one adds all dimers (FMO2) and (especially) trimers (FMO3) to the monomer calculations. The total FMO2 energy can be written as

CMO

=

Lli|i⟩

occ CMO

i

2. NOTATION AND DEFINITIONS Much of the notation and definitions are adopted from Yamaguchi et al.9 This work assumes a basis set that contains both contracted and uncontracted Gaussian functions. 2.1. Indices. • i,j,k denote occupied canonical molecular orbital (occ CMO) indices • l,m,n,o denote localized molecular orbital (LMO) indices • a,b,c denote virtual molecular orbital (vir) indices • p,q,r,s denote any canonical molecular orbital (occ or vir) indices • t,u denote primitive Gaussian (PG) indices • μ,ν,ξ,σ denote atomic orbital (AO) indices • A,B,C denote fragment indices • I,J,K denote nuclei indices or multipole expansion points • α,β,γ,κ denote directions x, y, or z 2.2. Definitions. • ZI is the nuclear charge on atom I • Sps is the overlap integral between orbitals p and s • cμp is the canonical or virtual MO coefficient of AO μ in MO p • cLμL is the LMO coefficient of AO μ in LMO l • Pμν (= 2∑occ i cμicνi is the restricted Hartree−Fock (RHF) density matrix element for AOs μ and ν 2.3. Superscript Notation. A variable with a superscript in parentheses, e.g., S(x) pq , denotes that the derivative with respect to x is taken only of the AO terms, and any MO coefficients are considered to be constant. A variable with a fragment index as a superscript denotes the variable for that fragment. However, if the appropriate fragment is clear by context, the superscript might be omitted. 2.4. Derivative of a Canonical Molecular Orbital Coefficient with Respect to a Perturbation. The derivative of a MO coefficient can be written in terms of the orbital response matrix9 Ux ∂cμp

∑

|l ⟩ =

cμLl

x Uqp cμq

(2.1)

fragments

x

U is the orbital response matrix to a perturbation x. In this work, there are nuclear perturbations, which will be denoted by an x, and field perturbations, which will be denoted by a Greek letter. Ux is an NMO × NMO matrix, where NMO is the number of MOs. It is convenient to think about the response matrix in terms of sub-matrices; i.e., the occupied orbital−occupied orbital (occ-occ) block, the virtual orbital−virtual orbital (vir-vir) block, and the virtual orbital−occupied orbital (vir-occ) block.

∑ A

fragments

EA +

∑ A>B

(EAB − EA − EB)

(3.1)

EA and EAB are the energies of the monomers and dimers, respectively. Approximations to the dimer energies can be used to decrease the computational cost. An additional advantage of the FMO method is that it is naturally parallelizable. The FMO method is parallelized with the general distributed data 4744


Article

Journal of Chemical Theory and Computation interface (GDDI).17 Because the FMO method is a generally applicable approach to dividing a system into smaller pieces, it can be combined with any electronic structure method. The usual notation is FMO/A, where A is a specific quantum chemistry method, such as second-order perturbation theory (MP2). 3.2. Effective Fragment Potential (EFP) Method. The EFP method was initially developed to model aqueous solvent effects.2 In the EFP method, the system is split into solute and solvent molecules. In this context, the solute molecules are typically calculated using an ab initio electronic structure method. The one-electron term in the solute Hamiltonian is modified by an explicit EFP solvent model potential. An EFP is generated by performing a single ab initio calculation on a solvent molecule, and then using the wave function to generate the input for the potential. Thus, it contains no empirical or fitted parameters. EFP internal geometries are rigid. More broadly, the EFP method can be used to explore intermolecular (non-covalent) interactions, without the need for an ab initio component. In this case, the system is divided into fragments that are modeled with EFPs. The EFP only (no ab initio solute) method is considered in this work. The EFP method decomposes the interaction energy of a system into the Coulomb energy, exchange repulsion energy, dispersion energy, charge-transfer energy, and many-body polarization energy terms. All terms are pairwise additive except for the polarization energy. The energy can be written as

is the multipole interaction tensor for sites I and J, and RIJ is the distance between expansion points I and J. The multipole moments on each site are calculated using the DMA. This is described in more detail in the Supporting Information. The multipole moments can be expressed as

+

+

1 3 1 3

B

⎡

I x ,y,z

J

⎢⎣

x ,y,z

∑

−2

α ,β ,γ

l

m

∑

∑

l

m

∑

∑

l

m

LMO ∈ A LMO ∈ B

+2

∑

∑

l

m

nuclei ∈ A

−

∑ I

2 ⎛ ⎞ −2Slm 1 ⎜ ⎟ R lm ⎝ −2 ln|Slm| ⎠

(3.5)

2

2 − 2lnSlm Slm π R lm

⎡ LMO ∈ A Slm⎢ ∑ FlnASnm + ⎢⎣ n

LMO ∈ B

∑ n

⎤ B Fmn Snl − 2Tlm⎥ ⎥⎦

LMO ∈ B ⎡ nuclei ∈ B Z 1 J 2⎢ − ∑ +2 ∑ Slm ⎢⎣ R R lJ ln J n

ZI +2 R lm

LMO ∈ A

∑ n

⎤ 1 1 ⎥ − R nm R lm ⎥⎦

(3.6)

RlJ is the distance between MO centroid l and atom J, Tlm is the kinetic energy integral between l and m, and FAln is the Fock matrix element between l and n on fragment A. 3.2.3. Polarization Term. The polarization energy (sometimes referred to as the induction interaction since it arises from multipole−induced multipole interactions) can be thought of as the interaction energy that occurs due to the change in the charge distribution of one molecule by the electric field due to the charge distribution of the other molecule. In the EFP method, the polarization energy is calculated by placing LMO dipole polarizability tensors, αβγ, on the LMO centroids of each fragment. The electric field of the other fragments (due to both the static multipole field and the induced dipoles on the other fragments) acts on the polarizability tensors and self-consistently generates induced dipoles, p, on the LMO centroids of the fragment.20 The induced dipole on LMO centroid l in the β-direction in fragment A, pAl,β, is

α

⎤ IJ + ...⎥ μαJ Θ Iβγ Tαβγ ⎥⎦

∑

LMO ∈ A LMO ∈ B

x ,y,z

α

∑

LMO ∈ A LMO ∈ B rep = −2 EAB

x ,y,z

I IJ IJ Tαβ + ∑ μαJ q I TαIJ − ∑ μαJ μβI Tαβ ∑ q J Θαβ α ,β

(3.4)

where Rlm is the distance between the LMO centroids of l and m, and ⟨l|x|l⟩ is the centroid in the x-direction for LMO l. 3.2.2. Exchange Repulsion Term. The exchange repulsion energy is a quantum mechanical contribution to the interaction energy that arises due to the Pauli exclusion principle. It is derived from approximations to the overlap of the wave functions of two isolated molecules.19 The exchange repulsion interaction energy between fragments A and B and can be expressed as

∑ ∑ ⎢q J q I T IJ − ∑ q J μαI TαIJ x ,y,z

Put′ ⟨u|(x − xI )|t ⟩

ut nearest I

LMO ∈ A LMO ∈ B chg‐pen = EAB

(3.2)

A

∑

Pμν

and so forth, where xI is the location of expansion center I and Put′ is the PG cross term that contains the product of the contraction coefficients for PG u and t. To account for charge penetration between interacting fragments A and B, an overlap-based damping term is computed, and added to the Coulomb interaction energy term.23 The expression for this term is

The EFMO method uses the interaction energy calculations, so they are considered in more detail below. Since the gradient involves taking the derivative of the energy terms, it is important to first consider the details of the energy expressions. The chargetransfer term is not considered here. 3.2.1. Coulomb Term. The Coulomb energy can be thought of as the energy produced from the interaction of the static charge density of two molecules. In the EFP method, the Coulomb energy is based on a Taylor series expansion of Coulomb’s law and a distributed multipole moment expansion using Stone’s distributed multipole analysis (DMA).18 Multipole moment expansion sites are distributed across each fragment in the system. The Coulomb contribution to the interaction energy between two EFP fragments is the sum of the interaction energy between all pairs of multipole moments Coul EAB =

PG u ∈ μ PG t ∈ ν

μν

pol EFP EAB + Etotal

A>B

∑

=−

Put′ ⟨u|t ⟩

ut nearest I

AO ∈ A

μxI

∑

Pμν

μν

fragments

∑

∑

q = ZI −

EFP Coul rep disp ct EAB = EAB + EAB + EAB + EAB EFP Etotal =

PG u ∈ μ PG t ∈ ν

AO ∈ A I

α ,β

(3.3)

In eq 3.3, ECoul AB is the Coulomb interaction energy between fragments A and B, I (J) runs over all multipole moment expansion points in A (B), qI is the monopole on site I, μI is the 1 dipole on site I, ΘI is the quadrupole on site I, TIJαβ...ν = ∇α ∇β ...∇ν R IJ

4745


Article

Journal of Chemical Theory and Computation {x , y , z}

plA, β =

The multipole interaction tensors are multiplied by a damping 23 2 2 function, Fpol damp,lI = 1 − exp(−RlI fg )(1 + RlI fg ). (The terms f and g are constants usually set to 0.6.) The damped multilI, damped = pole interaction tensors can be written as Tαβ...ν lI Fpol T . Substituting the damped multipole interaction damp,lI αβ...ν tensors into the static electric field, the damped static electric fields become

A αl , βγEltot, ,γ

∑

(3.7)

γ

where αl,βγ is the dipole polarizability tensor on LMO l, γ, κ are tot,A field directions, and El,γ is the total field in the γ-direction at tot,A LMO l. Since El,γ can be written in terms of a static electric field, 0,A El,γ , and a field due to induced dipoles, eq 3.7 can be rewritten as {x , y , z}

plA, β =

∑ γ

⎛ αl , βγ ⎜⎜El0,, γA + ⎝

fragments LMO ∈ B {x , y , z}

∑

∑ ∑

B≠A

m

κ

⎞ lm B ⎟ Tγκ pm , κ ⎟ ⎠

B≠A

∑ ∑ ElI0, γ B

⎛

{x , y , z}

∑ ∑ ⎜⎜q I TγlI + ∑

B≠A

⎝

I

lI μαI Tγα +

α

1 3

{x , y , z}

∑ αβ

⎞ I lI ⎟ Θαβ Tγαβ ⎟ ⎠

B

m

(D−1)lm , βα Em0,,Bα

α

(3.10)

Dlm , βγ = 0 (when l and m are on the same fragment) lm = −T βγ (when l and m are on different fragments)

disp EAB =

After generating the converged induced dipoles, the polarization energy can be calculated as

∑

Epol =

A

⎡ 1 = ∑ ⎢− ⎢ fragments ⎣ 2 A

LMO ∈ A {x , y , z}

∑ ∑

∑ B

n

α

∑ ∑ n

β

⎤ (D−1)nm , αβ Em0,,Bβ ⎥ ⎥⎦

occ CMO ∈ A vir ∈ A

∑ ∑ jk

∑ ∑ κ

(3.14)

C6, AB 6 RAB

C7, AB

+

7 RAB

+

C8, AB 8 RAB

+ ... (3.15)

C6, AB 6 RAB

+

1 C6, AB 6 , where all even terms with higher order than 3 RAB 1 C6, AB −7 6 RAB

. The R term has

25,26

but is not used in recently been derived and implemented, this work. The dispersion energy between fragments A and B can then be written in atomic units as

(3.11)

disp = EAB

4 ⎡⎢ 3 − 3 ⎢⎣ π 1

LMO ∈ A LMO ∈ B

∑

∑

l

m

1 ⎜⎛ 6 ⎝ R lm

∫0

∞

⎞⎤ α̅ l(iω)α̅ m(iω) dω⎟⎥ ⎠⎥⎦ (3.16)

{x , y , z}

l αββ (iω), and αlβγ(iω) is the distributed where α̅ l = 3 ∑β dynamic polarizability at LMO l for a frequency iω. Using a 12-point Gauss−Legendre quadrature and substitution of variables, the integral can be rewritten as a sum

LnjLnk UajγA⟨a|β|k⟩

a

m

the R−6 term have been approximated as 3

α

∑ ∑ m

Edisp AB =

En0,, αA

The dipole polarizability tensors on the LMO centroids of each fragment are calculated by decomposing the total dipole polarizability tensor for each fragment into contributions from each LMO.21,22 The dipole polarizability tensor on LMO centroid n on fragment A is αn , βγ = −4

(3.13)

In the EFP method, the dispersion energy is calculated by distributing isotropic dynamic polarizabilty tensors on the LMO centroids of each fragment. For this work, the total dispersion energy between fragments A and B is approximated as24

⎤ En0,, αApnA, α ⎥ ⎥⎦

LMO ∈ A {x , y , z}


×

⎞ lI ,damped ⎟ I Tγαβ Θαβ ⎟ ⎠

3.2.4. Dispersion Term. The dispersion energy can be thought of as the energy that arises from the interaction between induced multipoles on two molecules. The dispersion energy can be derived from Rayleigh−Schrodinger perturbation theory, starting from the sum of the Hamiltonians for two non-interacting molecules.18 The second-order correction to the energy contains the dispersion energy. The dispersion energy between fragments A and B can be written in terms of inverse powers of the distance between the molecules

Dll , βγ = (αl)−1 βγ

⎡ ⎢− 1 ⎢⎣ 2

∑ B≠A

where

fragments

lI ,damped μαI Tγα

α

⎞ lm ,damped B ,damped ⎟ Tγκ pm , κ ⎟ ⎠


+


∑ ∑

∑ γ

As in the Coulomb term, in this work, the expansion sites are only on the nuclei. Collecting the terms containing the induced dipoles, eq 3.7 can be written as

∑

αβ

∑

⎛ αl , βγ ⎜⎜El0,, γA ,damped ⎝

{x , y , z}

plA, β,damped =

(3.9)

plA, β =

∑

{x , y , z}

A,damped and the damped induced dipoles, pl,β , can be written as

I

fragments

⎛ ∑ ⎜⎜q I TγlI ,damped + I ⎝ B

{x , y , z}

1 3

+

B

B≠A

=

∑

=

where is the dipole moment interaction tensor for sites l and m. E0,A l,γ is the electric field at site l on fragment A due to the static DMA-calculated multipole moments, qI, μIβ, ΘIβγ on all multipole expansion points I on fragments other than A in the system fragments

I

fragments

Tlm γκ

=

∑ ∑ ElI0,damped ,γ B≠A

(3.8)

El0,, γA

B

fragments

El0,, γA ,damped =

(3.12)

In eq 3.12, γ is a field perturbation. 4746


Article

Journal of Chemical Theory and Computation ⎡ 4⎢ 3 − 3 ⎢⎣ π

disp EAB =

LMO ∈ A LMO ∈ B

∑

∑

l

m

1 6 R lm

12

⎞⎤ l m ⎟⎥ i i α ω α ω ( ) ( ) ̅ ̅ f f ⎟⎥ 2 ⎝ (1 − t f ) ⎠⎦ ⎛

∑ ⎜⎜wf f

The general EFMO energy expression is

2v0

(3.17)

where wf, v0, and tf are constants used in the numerical quadrature. The distributed dynamic polarizability on a fragment at LMO l for a frequency iω, αlβγ(iω), can be calculated as follows

b

(Hai(2), bjHbj(1), ckZckγ (iω)) − (iω)2 Zaiγ (iω) Hai(2), bj⟨b|γ |j⟩

j

Hai(2), bj = (aj|bi) − (ab|ij) + (ϵa − ϵi)δabδij Hai(1), bj = 4(ai|bj) − (aj|bi) − (ab|ji) + (ϵa − ϵi)δabδij (3.19)

where ϵa and ϵi are virtual and occupied orbital energies, respectively. The EFP method contains a multiplicative damping factor for the dispersion term.23 Incorporating the damping term, the dispersion energy becomes disp, damped EAB

4⎡ 3 = ⎢− 3 ⎢⎣ π ×

⎛ ⎝

⎜

∫0

∞

LMO ∈ A LMO ∈ B

∑

∑

l

m

⎛ −2 ln|Slm| ⎞ ⎟ = 1 − |Slm| ∑ ⎜ ⎝ ⎠ n! n=0 6

EFP pol (EAB ) + Etot

(3.22)

|rI − rJ| VI + VJ

is the relative

(3.23)

The charge-transfer term is not included in this work. As noted above, the charge-transfer term is the most computationally expensive component of the EFP energy, and it is usually the smallest term in the EFP interaction energy. Charged systems are an exception.23 Additionally, since charge transfer is a shortrange interaction,18 most of the charge-transfer interaction energy will be captured by the ab initio dimer interaction. Therefore, it is not necessary to have charge transfer in the long-range EFP interaction energy. Substituting eq 3.23 into eq 3.22, the energy expression becomes

(3.20)

In the EFMO method, the damping function is an overlap-based formula disp,damp Flm

∑

EFP Coul rep disp EAB = EAB + EAB + EAB

disp,damp 1 Flm 6 R lm

⎤ ⎞ α̅ l(iω)α̅ m(iω) dω⎟⎥ ⎠⎥⎦

A>B

minimum interatomic distance between atoms I on fragment A and atoms J on fragment B, weighted by the sum of the van der Waals radii, VI and VJ. RA,B is compared to Rcut to determine if the EFP method is to be used to calculate the dimer energy. The EFMO energy is calculated by summing the gas-phase ab initio energy of each monomer (fragment). Then, one loops over all pairs of monomers, and the dimer energy is added to the monomer energy. If the distance between two monomers is less than Rcut, the dimer energy is calculated with the chosen gasphase ab initio method (subtracting out the EFP polarization energy of the dimer to avoid double counting). If the distance is greater than Rcut, the dimer energy is approximated by the EFP interaction energy. The EFP polarization energy of the entire system is then added to the dimer and monomer energies. For this work, the long-range EFP energy is

c k occ vir CMO

∑

A RA , B> R cut

The inter-fragment distance RA,B = minI∈A,J∈B

occ occ vir CMO vir CMO

= −2 ∑

0 pol (ΔEAB − EAB )

polarization energy for fragments A and B. EFMO dimer calculations are performed with the chosen ab initio method (e.g., MP2) unless the two fragments in the dimer are farther apart than a predetermined cutoff Rcut. In the latter case, the dimer calculation is done using the EFP method.

where Zaiγ (iω) is the response vector that is calculated from solving the dynamic analog of the time-dependent coupled perturbed Hartree−Fock theory (TD CPHF) equations,27 and γ is a field perturbation. The TD CPHF equations are

j

∑

A>B

(3.18)

b

EA0 +

where E0A is the gas-phase energy of fragment A, ΔE0AB = E0AB − E0A − E0B (the dimer 2-body interaction energy), EEFP AB is the longrange EFP energy between fragments A and B, Epol tot is the EFP polarization energy for the entire system, and Epol AB is the EFP

occ occ ⎛ CMO ⎞⎛ CMO ⎞ γ ⎜ ⎟ ⎜ = −2 ∑ ∑ ⟨a|β|j⟩Llj ∑ Zai(iω)Lli⎟ ⎜ ⎟⎜ ⎟ a ⎝ j ⎠⎝ i ⎠

∑∑∑∑

∑ +

vir

l αβγ (iω)

RA , B≤ R cut

fragments

E EFMO =

n /2

2

RA , B≤ R cut

fragments

∑

E EFMO =

(3.21)

EA0 +

A RA , B> R cut

The damping function in eq 3.21 is a recent improvement on the original EFP damping function.28 3.3. General EFMO Energy Expression. The EFMO energy expression is a many-body expansion, similar to the FMO energy expression. As with the FMO method, the EFMO method begins by dividing the system into fragments. However, the monomer and dimer energy calculations differ: the EFMO method contains a many-body EFP polarization term, generated from all of the fragments and does not require the self-consistent convergence of the monomer Coulomb field. Importantly, the EFMO method inherits the GDDI parallelization of the FMO method.

∑

+

∑

0 pol (ΔEAB − EAB )

A>B Coul rep disp pol (EAB + EAB + EAB ) + Etot

A>B

(3.24)

This can be written as EFMO EFMO E EFMO = Eab initio + E EFP

(3.25)

where RA , B≤ R cut

fragments EFMO Eab initio =

∑ A

4747

EA0 +

∑

0 ΔEAB

A>B DOI: 10.1021/acs.jctc.6b00337 J. Chem. Theory Comput. 2016, 12, 4743−4767

Article

Journal of Chemical Theory and Computation EFP ∂(EAB ) = ∂xA

and RA , B≤ R cut EFMO E EFP =

RA , B> R cut pol (−EAB ) +

∑ A>B

∑

Coul rep disp pol (EAB + EAB + EAB ) + Etot

A>B

fragments

RA , B> R cut

+

∑ A>B

∑ A

∂EA0 + ∂xK

RAB≤ R cut

∑ A>B

⎛ ∂ΔE 0 ∂Epol ⎞ AB ⎜⎜ − AB ⎟⎟ ∂xK ⎠ ⎝ ∂xK

pol ⎛ ∂E Coul ∂Etot ∂Erep ∂E disp ⎞ ⎜⎜ AB + AB + AB ⎟⎟ + ∂xK ∂xK ⎠ ∂xK ⎝ ∂xK

A

∑ K

(4.1)

(

0 ∂EA0 ∂ Δ EAB , ∂x ∂xK K

∂⟨l|β|l⟩ = δβx , ∂xK

∂ ∂xK ∈ A

β = x, y, z (4.4)

occ CMO ∈ A vir ∈ A X = NRXAB , xK + EAB

∑ ∑ i

(4.2)

vir ∈ A

UijxKALXA,,ijB +

∑

+

UaixKALXA,,aiB

a

occ CMO ∈ A

4.1. Gas-Phase Gradient Terms. Two ab initio gas-phase terms

(4.3)

If bond midpoints are used as multipole expansion points, a similar expression applies for the derivative of the position of the bond midpoints. (3) Since the EFMO fragments are not rigid, the gradient with respect to each atom is calculated. Each EFP interaction energy gradient term between fragments A and B in the EFMO method can be written in the form

Each term in eq 3.24 is differentiated with respect to the x-coordinate of atom K (xK). The EFMO energy expression is a combination of gas-phase ab initio energy terms (E0A, ΔE0AB) and EFP interaction energy Coul rep disp pol terms (Epol AB , EAB , EAB , EAB , Etot ). Thus, the EFMO gradient is derived from ab initio gradient terms and EFP interaction energy gradient terms. To make the different types of terms clear, eq 3.25 can be used to write eq 4.1 as EFMO ∂E EFMO ∂E EFP ∂E EFMO = ab initio + ∂xK ∂xK ∂xK

K

EFP ∂(EAB ) ∂xK ∈ A

where xA is the translational motion of fragment A in the x-direction. (2) The derivative of an LMO centroid appears in the exchange-repulsion, polarization, and dispersion gradient terms. The derivative of an LMO centroid with respect to the translational motion of a rigid fragment is a delta function. That is, when a fragment translates, the LMO centroids move with it

4. ANALYTIC EFMO GRADIENT The expression for the EFMO gradient is ∂E EFMO = ∂xK

A

∑

∑

ij

) are computed using standard method-

xK A A , B Uab LX , ab

ab

LMO ∈ A LMO ∈ A

+

ology,29 so they are not discussed here. Note that if the gas-phase ab initio method chosen has response terms (e.g., MP2), then response equations for the monomers and dimers must be solved. For the monomer terms, the responses can be added to the response equations that arise from the EFP interaction energy gradient terms (formulated in later sections) and solved without additional cost. For the dimer terms, the response equations are solved separately and added into the gradient. 4.2. EFP Interaction Energy Gradient Terms. The gradient terms for the EFP method were derived and implemented previously.2 However, the EFP gradient terms cannot be used in EFMO directly, because the EFP method has rigid fragments while the EFMO method has flexible fragments. In the EFMO method, the internal geometry can change during a geometry optimization or molecular dynamics simulation, so the gradient must take this flexibility into account. For each term in the EFP interaction energy, a general formula for the nuclear gradient is presented below. The EFP terms in the EFMO gradient and the EFP translational gradient in the EFP method can both be derived from the general formula for the EFP nuclear gradient. After presenting the general formula for the nuclear gradient, the terms needed for the EFP method will be briefly discussed, since they are already implemented and can be reused in the EFMO method. Then, the remaining terms needed for the EFMO method will be discussed. To compare the EFP and EFMO gradient terms, it is useful to note the following points: (1) The translational gradient of the EFP interaction energy between fragments A and B with respect to the coordinates of fragment A can be derived by summing over the nuclear gradient of the EFP interaction energy with respect to the coordinates of each atom on fragment A

∑

∑

xK A A , B vml MX , ml

l

m occ {x , y , z} CMO ∈ A vir ∈ A

+

∑

∂UaiβA β , A , B NX , ai ∂xK

∑ ∑ i

β

a

occ 12 {x , y , z} vir ∈ A CMO ∈ A

+

∑ ∑ ∑ ∑ f

β

a

i

∂ZaiβA(iωf ) ∂xK ∈ A

β ,ω ,A ,B

NX , aif

(4.5)

In eq 4.5, the superscript/subscript X represents one of the EFP components Coul, rep, pol, or disp, corresponding to Coulomb, exchange-repulsion, polarizationn or dispersion. NRXAB,xK contains all “non-response” terms that do not contain a first- or second-order CMO response or a localization response. The response matrices UxKA/UβA, vxKA, and ZβA(iω)) are defined in Sections 2.4, 2.6, and 3.2.4, respectively. The superscript A indicates response matrices for fragment A. Using the Z-vector method (see Appendix A), the last three terms in eq 4.5 can be replaced with non-response terms and terms involving the CMO response matrix. The CMO response term can then be obtained using the Z-vector method. Throughout the following, the gradient for each term will be written in a manner that is consistent with eq 4.5. Since the EFP terms are based on MOs obtained from a separate gas-phase ab initio calculation on a particular monomer, the response equations for each monomer depend only on that monomer. Thus, in contrast to the FMO method,30 there is no response equation with the dimension of the entire system. 4.2.1. Coulomb Gradient Term. The gradient of the Coulomb interaction energy between fragments A and B can be written as 4748


Article

Journal of Chemical Theory and Computation A

Coul ∂EAB = ∂xK ∈ A

1 + 3 −

B

⎡

∑ ∑ ⎢q J I

J

x ,y,z

∑q

J

α ,β

⎢⎣

∂(q I T IJ ) − ∂xK ∈ A

I IJ ) ∂(Θαβ Tαβ

∑

∑ qJ

∑ μαJ

+

α

+

1 3

x ,y,z

μαJ

∑ α ,β ,γ

For the EFP method, the gradient of the EFP Coulomb term was derived by Day et al.2 In the EFP method, only the first term in eq 4.7 is included in the translational gradient, since the multipole moments depend only on the internal geometry of the fragment, and do not change as the fragment translates. The net Coulomb translational gradient on the fragment is calculated by summing the derivatives of the Coulomb energy with respect to each atom center on the fragment. The EFP implementation of the first term in eq 4.7 can be reused for the EFMO gradient, with the gradient stored separately for each atom. As shown in eq 3.4, the multipole moments are a sum of the product of a density matrix and a Gaussian function integral. The gradient of a multipole moment can therefore be calculated using the product rule. Consequently, each multipole moment derivative gives rise to a term involving AO derivatives and a term involving the CMO response matrix. The final EFMO Coulomb gradient can be written as

∂(μαI TαIJ ) ∂xK ∈ A

α

x ,y,z

∂xK ∈ A

IJ x ,y,z ) ∂(μβI Tαβ J μα ∂xK ∈ A α ,β

x ,y,z

∂(q I TαIJ ) ∂xK ∈ A

IJ ) ∂(Θ Iβγ Tαβγ

∂xK ∈ A

⎤ + ...⎥ ⎥⎦ (4.6)

Each term in eq 3.3 is differentiated with respect to the x-coordinate of atom K in fragment A. The multipole moments on fragment B are constant with respect to atoms on fragment A, so those terms are not included in the derivative. For this work, the expansion in eq 4.6 is terminated at the quadrupole− quadrupole term, and multipole expansion points are only on atomic centers. The gradient terms are derivatives of products, so the product rule can be used. Then, eq 4.6 can be written as

Coul ∂EAB = NRCoul AB , xK + ∂xK ∈ A

⎛⎧ ∂T IJ ⎫ ⎞ Coul ∂EAB αβ ...γ ⎪ Coul ⎜⎪ I J ⎟ ⎨ ⎬ = F AB , { m }, { m } ⎜⎪ ∂x ⎟ ⎪ ∂xK ∈ A ⎝⎩ K ∈ A ⎭ ⎠ ⎞ ⎛ ⎧ ∂m I ⎫ Coul IJ ⎜⎜{Tαβ ⎬, {m J }⎟⎟ + F AB ...γ }, ⎨ ⎩ ∂xK ∈ A ⎭ ⎠ ⎝ I

(4.7)

where m is an arbitrary multipole moment, is a multipole moment interaction tensor of the appropriate rank, ⎞

⎝⎩

⎠

αβ ...γ ⎬, {m I }, {m J }⎟ is the sum of all terms involving derivFCoul AB ⎜⎨ ∂x K∈A

⎭

∂m I ∂xK ∈ A

{ }, {m })

(

IJ atives of interaction tensors, and FCoul AB {Tαβ ... γ },

J

is the sum of all terms involving the derivative of multipole moments. These terms are expanded in Appendix B. rep ∂EAB = −2 ∂xK ∈ A

LMO ∈ A LMO ∈ B

∑

∑

l

m

⎡ ∂Slm ⎢ Slm ⎛ −2 ln Slm 2 ⎜⎜ − +4 ∂xK ⎢⎣ R lm ⎝ −π ln Slm π

LMO ∈ B ⎛ nuclei ∈ B Z 1 J − 2Slm⎜⎜ − ∑ +2 ∑ − R R lJ ln J n ⎝ LMO ∈ A LMO ∈ B

−2

∑

∑

l

m

LMO ∈ A nuclei ∈ B

+2

∑

∑

l

J

∂Tlm [−2Slm] + 2 ∂xK

nuclei ∈ A

∑ I

ZI +2 RIm

LMO ∈ A LMO ∈ B

∑

∑

l

m

nuclei ∈ A ∂R lJ ⎡ LMO ∈ B 2 ZJ ⎤ ⎢ ∑ Slm ⎥ + 2 ∑ ∂xK ⎢⎣ m R lJ2 ⎥⎦ I

∑ ∑ i

A ,B UaixKALCoul, ai

(4.8)

a

The details of the non-response term, NRCoul AB,xK, and the A,B coefficient of the CMO response matrix, LCoul,ai, are presented in Appendix B. In Section 4.3, eq 4.8 will be combined with the other EFP interaction energy gradient terms, and the Z-vector method will be applied to give the form of the EFMO gradient that was implemented in GAMESS. Coulomb Damping Term. The derivative of eq 3.5 can be easily added to the exchange repulsion gradient, so it is briefly discussed in the following subsection. 4.2.2. Exchange Repulsion Energy Term. The gradient of the EFMO exchange repulsion term can be expressed by taking the nuclear derivative of each term in eq 3.6, as follows

TIJαβ...γ

⎛⎧ ∂T IJ ⎫


LMO ∈ A

∑ n

LMO ∈ A LMO ∈ B ⎞ B ⎟⎟ + 2 ∑ FnlASnm + 2 ∑ Fnm Snl − 2Tlm ⎠ n n

LMO ∈ A ⎞⎤ 1 1 ⎟⎥ ∑ 2 − − R nm R lm ⎟⎠⎥⎦ l

2 2 ∂R lm ⎡ −2 ln Slm Slm Slm ⎢2 + + 2 2 π ∂xK ⎢⎣ R lm R lm LMO ∈ B

∑ m

LMO ∈ A ⎤ ∂RIm ⎡ 2 ZI ⎥ ⎢ ∑ Snm 2 ⎥⎦ ∂xK ⎢⎣ n RIm

Each term in eq 4.9 is differentiated with respect to the x-coordinate of atom K in fragment A. This means that gradient terms that depend only on the geometry of fragment B will be zero. In the EFP method, the internal geometries of the fragments do not change when the fragments translate, so the MO coefficients and Fock matrices do not change, and the LMO centroids move with the fragments. Thus, the kinetic energy and overlap integral derivatives (the first and third terms in eq 4.9) are computed by only taking the derivative of the AO integrals (ignoring the MO coefficient derivatives), and the Fock matrix

LMO ∈ A

∑ m

LMO ∈ B

∑ n

−2

∂FlmA ⎡ ⎢ ∂xK ⎢⎣

Sln2 2 R lm

LMO ∈ B

∑ n

⎤ SlnSmn⎥ ⎥⎦

LMO ∈ A

+

∑ n

−2

2 ⎤ Snm ⎥ 2 R lm ⎦⎥

(4.9)

derivative is not computed. Since the net translational gradient is calculated by summing over the nuclear gradients of the exchange repulsion energy with respect to each atom on the fragment, by eq 4.4 there is no need to explicitly calculate the derivative of each LMO centroid. The implementation of the first, third, fourth, fifth, and sixth terms in eq 4.9 for the EFP method can be reused for the EFMO gradient, but with additional terms added for the derivative of the LMO centroids and the CMO coefficients, and with the gradient stored separately for each atom. 4749


Article

Journal of Chemical Theory and Computation For the EFMO gradient, the LMO centroid derivatives in the fourth and fifth terms of eq 4.9 can be collected. The explicit expressions are shown in Appendix C. Equation 4.9 can then be written as rep ∂EAB = −2 ∂xK ∈ A

LMO ∈ A LMO ∈ B

∑

∑

l

m

∂Slm S [W lm ] ∂xK ∂FlmA ⎡ ⎢ ∂xK ⎢⎣

LMO ∈ A LMO ∈ A

−2

∑

∑

l

m

LMO ∈ A LMO ∈ B

∑

−2

∑

l m LMO ∈ A {x , y , z} l

∑ l

∑ n

⎤ SlnSmn⎥ ⎥⎦

pol ∂Etot =− ∂xK ∈ A

∂l α l [WlRα ] ∂xK

α

⎡ LMO ∈ A

LMO ∈ B

+2

LMO ∈ B

∂Tlm T [Wlm] ∂xK

∑ ∑

+2

In Section 4.3, eq 4.11 will be combined with the other EFP terms, and the Z-vector method will be applied to give the form of the EFMO gradient that was implemented. Coulomb Damping Function. The derivative of eq 3.5 can easily be added to the exchange repulsion term, since it only involves LMO dipole and overlap integrals. Since the derivative is derived in a similar manner to the exchange repulsion gradient, it is not shown here. 4.2.3. Polarization Energy Term. The EFMO polarization energy gradient can be derived beginning with eq 26 in ref 20 (here written in the notation of this paper)

(xK − ⟨l|x|l⟩)⎢ ⎢⎣

∑

2 Sml

m

+

ZK ⎤ ⎥ 3 ⎥⎦ RKl (4.10)

∂xK

∂xK

αl .

∑ ∑

+

∑

∑

∑

l

m

i,j

{x , y , z}

∑ β

⎡ {x , y , z} ∂⟨n|β|n⟩ ⎢ ∑ ∂xK ∈ A ⎢⎣ α

LMO ∈ A {x , y , z}

∑ ∑ n

m

n

α ,β

⎛ ∂Dnm , αβ ⎞ B (pnC̃ , α )⎜ ⎟p ⎝ ∂xK ∈ A ⎠ m , β

γ

∂⟨n|γ |n⟩ ∂xK ∈ A

pol ∂EAB , ∂xK ∈ A

can be derived in a

B

fragments

nI αβ ... γ

∂m I ∂xK ∈ A

n

The second term in eq 4.12 can be expanded using the definition in eq 3.10. The expansion of eq 4.12 is shown in more detail in Appendix D. Then, eq 4.12 can be written as

xK A A , B vml M rep, ml

nI ⎞⎤ ⎫ ⎛ ∂E 0 ⎞⎛ p A + p ̃ A ⎞⎤⎤ ⎡ ⎛⎧ ⎪ ∂Tαβ ... γ ⎪ nI , α n , α ⎟⎥⎥ I ⎢F P ⎜⎨ ⎟⎥ ⎟⎟⎜ n , α ⎬ − , { m }, { p } n ⎟⎥ ⎪ ⎜ ⎟⎥⎥ ⎢ A ⎜⎪ ∂x 2 ⎝ ∂⟨n|β|n⟩ ⎠⎝ ⎠⎦ ⎠⎦⎦ ⎣ ⎝⎩ K ∈ A ⎭

∑ ∑ ⎜⎜

B≠A

I

⎡ ⎛ ⎞⎤ ⎡ 1 ⎧ ∂m I ⎫ nI ⎬, {pn }⎟⎟⎥ + ⎢ − ⎢FAP ⎜⎜{Tαβ ... γ }, ⎨ ⎢⎣ ⎝ ⎩ ∂xK ∈ A ⎭ ⎠⎥⎦ ⎢⎣ 2 ⎡ 1 +⎢ ⎢⎣ 2

BC

∑ ∑

( ({T }, { }, {p })).

(4.11)

⎡ LMO ∈ A pol ∂Etot ⎢ = ⎢− ∑ ∂xK ∈ A n ⎣

∑

multipole moment derivatives FAP

LMO ∈ A LMO ∈ A A ,B U jixKAOrep, ji +

∑

{ }

A ,B UaixKALrep, ai

a

occ CMO ∈ A

α

similar fashion, so the details are only shown for the total polarization energy term. The static electric field in the first term in eq 4.12 is represented by the multipole moment expansion as in the Coulomb term. The derivative is handled in a similar manner here as in the Coulomb term. That is, it is split into a term with the sum of all multipole moment interaction tensors ⎛ P ⎛ ∂TαβnI ... γ ⎞⎞ ⎜F ⎜ , {m I }, {pn }⎟⎟ and a term with the sum of all ⎝ A ⎝ ∂xK ∈A ⎠⎠

occ CMO ∈ A vir ∈ A i

n

fragments LMO ∈ B LMO ∈ C {x , y , z}

The dimer polarization energy term,

It is important to note that the MO coefficient derivatives are derivatives of LMO coefficients, so the derivative results in a term with a CMO response matrix and a term with a localization response matrix, as shown in eq 2.3. Appendix C provides the details that lead from eq 4.10 to eq 4.11. Combing all non-response terms into NRrep AB , and writing out the response terms, one obtains rep ∂EAB = NRrep AB , xK + ∂xK ∈ A

∑ ∑

B

The derivative in eq 4.12 is taken with respect to atom K on ∑LMO∈B ∑{x,y,z} fragment A. As in ref 20, p̃Al,β = ∑fragments B m α −1 T 0,B (D )lm,βαEm,α. All other terms are defined in Section 3.2.3.

all the terms in the coefficient of ∂Tlm , and WRlα holds all the terms in the coefficient of

∑

⎛ ∂E 0, B ⎞⎛ p B + p ̃B ⎞ n,α ⎟ ⎜⎜ n , α ⎟⎟⎜ n , α ⎜ ⎟ ∂ x 2 ⎝ K ∈ A ⎠⎝ ⎠

(4.12)

where WSlm holds all the terms in the coefficient of ∂Slm , WTlm holds ∂l ∂xK

1 2


LMO ∈ A {x , y , z}

∑ ∑ n


∑ B≠A

∑ ∑ m

α ,β

The EFP polarization gradient for the EFP method was derived by Li et al.20 and Day et al.2 As in the Coulomb and exchange-repulsion gradient, the net polarization translational gradient with respect to fragment A can be calculated by summing the nuclear derivatives with respect to each atom on fragment A. Only the first, second, and fifth terms in eq 4.13 are needed for the EFP translational gradient. The third and fourth terms have derivatives of the multipole moments and the dipole polarizability tensor, respectively, which depend

β ,γ

⎛ ∂(α −1 ) ⎞ ⎤ n , βγ ⎟⎟p A ⎥ n,γ ⎝ ∂xK ∈ A ⎠ ⎥⎦

(pñ A, β )⎜⎜

⎡ ⎛ ∂Dmn , αβ ⎞ A ⎤⎤ ⎛ ∂Dnm , αβ ⎞ B B ⎢(pñ A, α )⎜ ⎟pm , β + (pm̃ , α )⎜ ⎟p ⎥⎥ ⎝ ∂⟨n|γ |n⟩ ⎠ n , β ⎦⎥⎦ ⎝ ∂⟨n|γ |n⟩ ⎠ ⎣

(4.13)

only on the internal geometry. As with the exchange repulsion term, the terms in eq 4.13 that contain derivatives of the LMO centroids can be expressed without explicitly calculating ∂⟨l | β | l⟩ , by using eq 4.4 instead. The EFP implementation of ∂xK

the first, second, and fifth terms can be used for the EFMO method, with additional terms added for the derivative of the LMO centroids, and the gradient stored separately for each atom. 4750


Article


The CMO response terms and the localization response terms can also be removed using the Z-vector method. This will be done in Section 4.3 for all EFP interaction energy terms. Polarization Damping Function. The polarization energy derivative can be modified to include damping. The damping term is a function of the distance between two LMO centroids or an LMO centroid and an atom center, so the derivative is straightforward. Using the expression for polarization damping in Section 3.2.3

The LMO centroid derivatives in eq 4.13 can be combined. The third term can be replaced with two terms arising from the derivative of the multipole moments, as in the Coulomb term. Since this involves the derivative of the CMO density matrix, a CMO response matrix term is necessary. The fourth term in eq 4.13 can be manipulated using matrix derivative operations20 and the defintion in eq 3.7 1 2

⎛ ∂(α −1 ) ⎞ n , βγ ⎟⎟p A n,γ ⎝ ∂xK ∈ A ⎠

LMO ∈ A {x , y , z}

(pñ A, β )⎜⎜

∑ ∑ n

β ,γ

=−

1 2

LMO ∈ A {x , y , z}

∑ ∑ n

β ,γ

lm ,damped ∂Tαβ ... γ

∂x

⎛ ⎞ tot, A ∂(αn , βγ ) A ⎟Entot, Eñ , β ⎜ ,γ ∂ x ⎝ K∈A ⎠

tot, A

⎡ disp ∂EAB ∂ ⎢ 4 − = ∂xK ∈ A ⎢⎣ π ∂xK ∈ A

{x , y , z}

=

∑

αl−, κβ1 pl̃ ,Aκ

κ

The LMO dipole polarizability tensor in eq 4.14 is expanded as21 −

1 2

LMO ∈ A {x , y , z}

∑ ∑ n

1 =− 2

β ,γ

n

f

∑

=

jk

a

⎤ γA A LnjLnk Uaj ⟨a|β|k⟩)⎥Entot, ⎥ ,γ ⎦

(

∂xK ∈ A

∂xK ∈ A

+

∑

disp ∂EAB = ∂xK ∈ A

A ,tot UaixcALpol, ai

⎡ 4 × ⎢− ⎢⎣ π

a vir ∈ A

A ,tot UijxKAOpol, ij +

ij

∑

xK A A ,tot Uab Vpol, ab

∑

ab

∑

+

∑ β

∑ ∑ i

a

∂UaiβA β , A ,tot Npol, ai ∂xK

2v0

∂α̅ l(iωf ) ⎡ 4 ⎛ 2v0 ⎞ ⎢ − ⎜w ⎟ f ⎜ ∂xK ∈ A ⎢⎣ π ⎝ (1 − t f )2 ⎟⎠

(4.18)

l

β

∂⟨l|β|l⟩ ∂xK ∈ A

⎞⎤⎤ l m ⎟⎥⎥ i i α ω α ω ( ) ( ) ̅ f f ⎟ ⎥ 2 ̅ t (1 ) − ⎠⎥⎦⎦ ⎝ f ⎛

12

∑ ⎜⎜wf f

f

2v0

∂α̅ l(iωf ) ∂xK ∈ A

LMO ∈ B ⎤ ⎡ ⎛ 2v0 ⎞ 4 ⎟ ∑ 1 α̅ m(iωf )⎥ × ⎢− ⎜⎜wf ⎟ 2 6 ⎥⎦ ⎢⎣ π ⎝ (1 − t f ) ⎠ m R lm

m

+

∑ ∑

∑ ∑ l

xK A A ,tot vml M pol, ml

l occ {x , y , z} CMO ∈ A vir ∈ A

−6 ∂(R lm ) ∂xK ∈ A

⎤ 1 m ⎥ ( i ) α ω ̅ f 6 ⎥⎦ R lm

LMO ∈ A {x , y , z}

LMO ∈ A 12

LMO ∈ A LMO ∈ A

+

m

⎡ LMO ∈ B − 6(⟨l|β|l⟩−⟨m|β|m⟩) ×⎢ ∑ 8 ⎢ R lm ⎣ m

occ CMO ∈ A vir ∈ A i

l

1 6 R lm

The first term in eq 4.18 can be written in terms of an LMO centroid derivative as in the polarization energy gradient. Equation 4.18 can then be expressed as

Once eq 4.15 has been expanded in terms of non-response terms, localization transform derivative terms, and second-order CMO field response terms, the polarization energy gradient can be rewritten as

occ CMO ∈ A

∑ m

UajγA

∂

∑ ∑

f

LMO ∈ B

×

); and one ⟨a|β|k⟩). term with the derivative of the dipole, ( the derivative of the CMO field response,

(4.17)

⎞⎤ l m ⎟⎥ ( ) ( ) i i α ω α ω ̅ f f ⎟ 2 ̅ − (1 t ) ⎝ ⎠⎥⎦ f

f

l

The right-hand side (RHS) of eq 4.15 results in three terms: one term with the derivative of the LMO transforms; one term with

∑

∑ ⎜⎜wf

∑ ∑

+

∑

⎛

12

LMO ∈ A 12

(4.15)

∂

m

⎡ 4 × ⎢− ⎢⎣ π

β ,γ

pol ∂Etot = NRpol A ,tot, xK + ∂xK ∈ A

∑

l

∑ ∑

∂x

2v0

LMO ∈ A LMO ∈ B

tot, A Eñ , β


lm ∂Tαβ ... γ

⎞⎤ l m ⎟⎥ ( i ) ( i ) α ω α ω ̅ f f ⎟ 2 ̅ ⎝ (1 − t f ) ⎠⎥⎦ ⎛

12

LMO ∈ A {x , y , z}

∑ ∑

LMO ∈ A LMO ∈ B

∑ ⎜⎜wf

×

⎛ ⎞ tot, A ∂(αn , βγ ) A ⎟Entot, Eñ , β ⎜ ,γ ⎝ ∂xK ∈ A ⎠

⎡ ∂ (−4 ×⎢ ⎢ ∂xK ∈ A ⎣

∂x

pol lm Tαβ ... γ + Fdamp, lm

All multipole interaction tensor derivatives can be replaced with the above, and the gradient can be evaluated in the same way. 4.2.4. Dispersion Energy Term. The EFMO dispersion gradient can be expressed by taking the derivative of eq 3.17

(4.14)

where El̃ , β

pol ∂Fdamp, lm

=

(4.19)

In the EFP method, only the first term in eq 4.19 is needed for the translational gradient. As in the exchange repulsion and polarization terms, there is no need to calculate the nuclear derivative of the LMO centroids explicitly to get the EFP translational dispersion gradient. The second term contains the derivative of the dynamic polarizability tensor, which only depends on the internal geometry of the fragment, and does not change as the fragment translates. For the EFMO gradient, the first term in eq 4.19 can be calculated

(4.16)

A,tot A,tot A,tot The terms in eq 4.16 (such as NRpol A,tot,xK, Lpol,ai, Opol,ij, Vpol,ab, β,A,tot MA,tot pol,ml, and Npol,ai ) are similar to those in eq 4.5, but with a tot superscript/subscript instead of a B superscript/subscript, to denote that this is a gradient contribution from the total polarization energy instead of a gradient contribution from a dimer interaction energy between fragments A and B.

4751


Article


Dispersion Damping Function. The damping function, shown in eq 3.21 adds a factor that depends only on the overlap. The energy and gradient then become

using the implementation from the EFP method, but with additional terms added for the derivative of the LMO centroids, and the gradient stored separately for each atom. The first term in eq 4.19 is a LMO centroid derivative. The LMO centroid derivative has been discussed in the subsections on the exchange-repulsion and polarization gradient previously. The second term in eq 4.19 contains the derivative of the dynamic polarizability tensor. This is derived in a similar manner as the derivative of the static polarizability tensor, and is discussed in Appendix E2. The dispersion energy gradient can then be written as

LMO ∈ A LMO ∈ B disp,damped EAB =

∂xK ∈ A

disp,damped ∂EAB = ∂x

A ,B UaixKALdisp, ai

∑ ∑

= NRdisp AB , xK +

i

vir ∈ A A ,B UijxKAOdisp, ij +

∑

+

∑

ji

l

m

∑

l

m

+

∑ ∑ ∑ ∑ f

a

β

∂xK ∈ A

i

β ,ω ,A ,B Ndisp,fai

RA , B≤ R cut

∑

+

i occ CMO ∈ A

+

∑ ij

⎛ RA , B≤ R cut

UijxK A ⎜⎜ ⎝

LMO ∈ A LMO ∈ A

+

∑

∑

l

∑

∑

i

A ,B (−Opol, ij) +

xK A ⎜ vml ⎜ ⎝

∑ A>B

∑

⎞ A ,B A ,B A ,tot ⎟ (Orep, ij + Odisp, ij) + Opol, ij ⎟ + ⎠

A ,B (−M pol, ml ) +

∑

A>B

∑ A>B

∑

∂xK ⎜⎝

A>B

EFMO ∂E EFP = NRtot A , xK + ∂xK

i

vir ∈ A

UijxKA(LijA,tot) +

∑

∑

xK A A Uab (Lab ,tot)

ab

LMO ∈ A LMO ∈ A

∑

∑

l

m

xK A (MmlA ,tot) vml


∑ ∑ ∑ ∑ f

a i occ {x , y , z} CMO ∈ A vir ∈ A

+

β

∑ β

∞ ⎡ 4 1 ⎤ ∂⎢ − π 6 ∫ α̅ l(iω)α̅ m(iω) dω ⎥ ⎣ Rlm 0 ⎦

(

)

(4.21)

∂x

∑ ∑ i

a

∑ ab

⎛ RA , B≤ R cut xK A ⎜ A ,B Uab ⎜ ∑ (−Vpol, ab) + ⎝ A>B occ 12 {x , y , z} vir ∈ A CMO ∈ A

∑ ∑ ∑ ∑ f

β

a

i

RA , B> R cut

∑ A>B

⎞ A ,B A ,tot ⎟ (Vdisp, ab) + Vpol, ab⎟ ⎠

∂ZaiβA(iωf ) ⎛ ⎜ ∂xK ∈ A ⎜⎝

RA , B> R cut

∑

∂ZaiβA(iωf ) ∂xK ∈ A

∂UaiβA β , A Nai ,tot ∂xK

β ,ω ,A ,B

Ndisp,fai

A>B

⎞ ⎟ ⎟ ⎠

The Z-vector method can be used to replace the term that involves the derivative of the CMO response (the last term in eq 4.23), the derivative of the time-dependent response term (the second to last term in eq 4.23), and the localization response term (the third to last term in eq 4.23). After solving the Z-vector equations, the second-order canonical response term, secondorder time-dependent response term, and localization response term are replaced with first-order canonical response terms and non-response terms. The first-order canonical response terms can then be collected with the other first-order canonical response terms. Using details given in Appendix D3, Appendix E1, and Appendix A8, eq 4.23 can be written as

UaixKALaiA,tot

ij

+

disp,damp Flm

(4.22)

a

occ CMO ∈ A

+

m

⎞ β ,A ,B β , A ,tot ⎟ (−Npol, ai ) + Npol, ai ⎟ ⎠


∑ ∑

vir ∈ A

⎞ A ,B A ,B A ,tot ⎟ (M rep, ml + Mdisp, ml ) + M pol, ml ⎟ + ⎠

where the derivative is taken with respect to atom K on fragment A. Then, the non-response terms and the coefficients of the response terms can be collected and combined into terms with the superscript/subscript tot, as follows

+

∑

l

⎞ A ,B A ,B A ,B A ,tot ⎟ (LCoul, ai + Lrep, ai + Ldisp, ai) + Lpol, ai ⎟ ⎠

RA , B> R cut

RA , B≤ R cut ∂UaiβA ⎛ ⎜

a

RA , B> R cut

A>B

⎛ RA , B≤ R cut

∑ ∑

β

∑

⎞⎤ α̅ l(iω)α̅ m(iω) dω⎟⎥ ⎠⎦

rep disp pol (NRCoul AB , xK + NRAB , xK + NRAB , xK ) + NRA ,tot, xK

RA , B> R cut

A>B

m occ {x , y , z} CMO ∈ A vir ∈ A

+

∑

⎛ RA , B≤ R cut A ,B UaixK A ⎜⎜ ∑ (−Lpol, ai) + ⎝ A>B

a

∞

A>B

A>B

∑ ∑

m

disp,damp ∂Flm ∂x

RA , B> R cut

(−NRpol AB , xK ) +


l

⎞⎤ α̅ l(iω)α̅ m(iω) dω⎟⎥ ⎠⎦

The damping function depends only on the overlap, and the damping function gradient can be computed in a similar manner to the LMO overlap derivatives in the exchange-repulsion energy term. 4.3. Combined Gradient. The terms in eq 4.1 (and likewise, eq 4.2) that are EFP interaction energy derivatives can be expressed using eqs 4.8, (4.11), (4.16), and (4.20)

(4.20)

As shown in Section 4.3, eq 4.20 can be combined with the other EFMO gradient terms, and the Z-vector method can be used to calculate the CMO response terms and the localization response terms. EFMO ∂E EFP = ∂xK

∑

⎡ 4 1 ⎛ ⎜ × ⎢− 6 ⎝ ⎣ π R lm

×

∂ZaiβA(iωf )

∑

+

xK A A , B vml Mdisp, ml


∞

LMO ∈ A LMO ∈ B

xK A A , B Uab Vdisp, ab

ab

∑

∫0

LMO ∈ A LMO ∈ B

LMO ∈ A LMO ∈ A

+

disp,damp Flm

∫0

a

occ CMO ∈ A

∑

⎡ 4 1 ⎛ ⎜ × ⎢− 6 ⎝ ⎣ π R lm


disp ∂EAB

∑

EFMO ∂E EFP = NRtot,3 A , xK + ∂xK

β ,ω ,A

Nai ,totf


∑ ∑ i

occ CMO ∈ A

+ (4.23)

∑ ij

4752

,3 UaixKA(LaiA,tot )

a vir ∈ A

,3 UijxKA(LijA,tot )+

∑ ab

xK A A ,2 Uab (Lab ,tot)

(4.24)


Article

Journal of Chemical Theory and Computation NRtot,3 A,xK is a non-response term resulting from the Z-vector A,3 methods, and LA,3 ai,tot, Lij,tot are the coefficients of the occ-occ and vir-occ CMO response matrices after the terms from the Z-vector methods have been added in. Next, the terms that involve the vir-vir and occ-occ parts of the CMO response matrix are considered. More details are provided in Appendix A6. Equation 4.24 becomes EFMO ∂E EFP = NRtot,3 A , xK + ∂xK core act CMO ∈ A CMO ∈ A

+

∑

∑

k

j

core CMO ∈ A

+

ij

+

∑ ab

∑ ∑ i

a

⎡ ⎢ ,3 UaixKA ⎢LaiA,tot ⎣

⎤ ⎛ A ′ L A ,3 A′ L A ,3 ⎞ ⎜ kj , ai kj ,tot + jk , ai jk ,tot ⎟⎥ ⎜ (ϵ −ϵ ) (ϵk −ϵj) ⎟⎠⎥ ⎝ j k ⎦

⎛ 1 (xK )A ⎞ A ,3 ⎜− ⎟(L S )+ ⎝ 2 ij ⎠ ij ,tot

∑ vir ∈ A


second-order CMO response, and second-order time-dependent CMO response are solved. The Z-vectors that result from solving the Z-vector equations are summed with a non-response term and a term that involves the CMO response matrix. Since the application of the Z-vector method to the localization response, second-order CMO response, and second-order time-dependent CMO response terms contributes to the coefficient of the CMO response matrix, these Z-vector equations must be solved first. Then the Z-vector equation for the CMO response matrix is solved.

⎛ 1 (xK )A ⎞ A ,2 ⎜− ⎟(L S )+ ⎝ 2 ab ⎠ ab ,tot

act CMO ∈ A

∑ ij

6. TEST CALCULATIONS To evaluate the accuracy of the gradient, two methods were used. First, the analytic gradient was compared to the numeric gradient for several systems (Section 6.1). Second, the EFMO method and analytic gradient were used in MD simulations to test energy conservation in a Velocity-Verlet NVE ensemble33,8 (Section 6.2). Again, note that neither the energy expression used in the numeric gradient nor the analytic gradient contain the EFP charge-transfer term. Although the charge-transfer term might make a significant contribution to the ionic liquid dimer, it is stressed that this section is meant only to assess the accuracy of the analytic gradient. 6.1. Analytic to Numeric Comparison. For the comparisons, the analytic gradient was computed for several systems and compared to the numeric gradient. A 6-31++G(d,p) basis set was used for all calculations, and Rcut was set to 0.3, forcing all dimer interaction energies to be evaluated as EFP interaction energies. The multipole moments are only expanded through the quadrupole−quadrupole term, and all multipole moment expansion points are exclusively on atomic centers, in contrast to the EFP method in which bond midpoints are also expansion centers. The numeric gradient was calculated using a two-point formula. For the three systems, the maximum absolute difference and root-mean-square deviation (RMSD) are presented. The RMSD, for N gradient elements, is calculated as

⎛ 1 (xK )A ⎞ A ,3 ⎜− ⎟(L S ) ⎝ 2 ij ⎠ ij ,tot

core act CMO ∈ A CMO ∈ A

∑

∑

i

j

⎡⎛ B xKA L A ,3 ⎞ ⎢⎜ ij ij ,tot ⎟ ⎢⎣⎜⎝ (ϵj −ϵi) ⎟⎠

⎛ B xKA L A ,3 ⎞⎤ ji ji ,tot ⎥ ⎟ + ⎜⎜ ⎟ ( ϵ −ϵ ⎝ i j) ⎠⎥⎦

(4.25)

A′kj,ai, BxjiKA,

where and ϵq are defined in eq A.1. The non-response terms in eq 4.25 can be combined, and then eq 4.25 becomes EFMO ∂E EFP = NRtot,4 A , xK + ∂xK

⎡ ⎢ A ,3 ⎢Lai ,tot + ⎣


∑ ∑

i core act CMO ∈ A CMO ∈ A

∑

∑

k

j

UaixKA

a

⎤ ⎛ A ′ L A ,3 A′ L A ,3 ⎞ ⎜ kj , ai kj ,tot + jk , ai kj ,tot ⎟⎥ ⎜ (ϵ − ϵ ) (ϵk − ϵj) ⎟⎠⎥ k ⎝ j ⎦ (4.26)

NRtot,4 A,xK

where is the sum of the non-response terms in eq 4.25. Finally, all CMO response terms are collected, and the Z-vector method can be used to replace the CMO response matrices. This gives the EFP interaction energy part of the EFMO gradient EFMO ∂E EFP = NRtot,3 A , xK + ∂xK


∑ ∑ i

a

BaixKA ZaiA ,cphf

N

∑i (analytic gradient element i − numeric gradient element i)2 N

(6.1)

The max interaction gradient value is the maximum contribution to the analytic gradient from the EFP interaction energy gradient. The interaction gradient is calculated by subtracting the one-body (ab initio) gradient from the total gradient for each gradient element. 6.1.1. 64 Water Molecules. The numeric gradient for the system of 64 water molecules (shown in Figure 1) was calculated using a 0.005 Å step size. Table 1 shows the maximum interaction analytic gradient value, the maximum analytic gradient value, the RMSD, and the maximum absolute difference between the numeric and analytic gradients. The RMSD and maximum absolute difference values are small and comparable to the RMSDs for other analytic gradients.30,33 The values in Table 1 demonstrate that the gradient is accurate for the system of 64 water molecules. 6.1.2. 5 Dimethyl Sulfoxide (DMSO) Molecules, 5 Methanol Molecules, and 10 Water Molecules. The numeric gradient for the system of 5 DMSO molecules, 5 methanol molecules, and 10 water molecules (shown in Figure 2) was calculated using a 0.005 Å step size and a 0.001 Å step size.

(4.27)

where vir ∈ A occ ∈ A

∑ ∑ b

Abj , ai ZbjA ,cphf =

j

⎡ ⎢ A ,3 ⎢Lai ,tot + ⎣

core act CMO ∈ A CMO ∈ A

∑

∑

k

j

⎤ ⎛ A ′ L A ,3 A′ L A ,3 ⎞ ⎜ kj , ai kj ,tot + jk , ai jk ,tot ⎟⎥ ⎜ (ϵ − ϵ ) (ϵk − ϵj) ⎟⎠⎥ k ⎝ j ⎦

5. IMPLEMENTATION The EFMO gradient has been implemented in the GAMESS quantum chemistry software package.32 The coefficients of the LMO centroid derivative term, Fock matrix derivative term, CMO response matrix term, localization response matrix term, and second-order response matrix terms are collected separately. The Z-vector equation for the localization response, 4753


Article


slightly different allocations of charge density in the multipole moment calculation. That is, some charge density components were allocated to different expansion centers. In some cases, this led to large enough energy differences in the calculation of the numeric gradient that the numeric and analytic gradient elements differed by ∼10−3 Hartree/Bohr. These differences, which are not observed for the water cluster discussed in the previous subsection, are too large to be considered accurate. Reducing the step size to 0.001 Å removed all of the instances in which the forward and backward energy calculations allocated charge density to different expansion points. The analytic and numeric gradients match well, as shown in Table 2. These results show that the analytic gradient is accurate. Table 2. Comparison of Analytic and Numeric Gradient (Hartree/Bohr) for a System of 5 Dimethyl Sulfoxide Molecules, 5 Methanol Molecules, and 10 Water Molecules (See Text for Details)

Figure 1. Geometry of 64 water molecules used in the numeric and analytic gradient comparison. Hydrogen atoms are light gray and oxygen atoms are red.

max absolute analytic gradient value

RMSD

max absolute difference

0.026037

0.025670

8.1 × 10−6

3.7 × 10−5

max absolute analytic gradient value

RMSD

max absolute difference

0.213522

0.208892

2.8 × 10−7

1.2 × 10−6

6.1.3. Ionic Liquid Dimer. The numeric gradient for the system of two hexafluorophosphate (PF6)− anions and two 1-Nbutyl-3-methylimidazolium (bmim)+ cations (shown in Figure 3) was calculated using a 0.005, 0.001, 0.0005, 0.0001, and 0.00005 Å step size.

Table 1. Comparison of Analytic and Numeric Gradient (Hartree/Bohr) for a System of 64 Water Molecules (See Text for Details) max interaction analytic gradient value

max interaction analytic gradient value

Figure 2. Geometry of a cluster of 5 DMSO, 5 methanol, 10 water molecules used in the numeric and analytic gradient comparison. Hydrogen atoms are light gray, carbon atoms are dark gray, oxygen atoms are red, and sulfur atoms are yellow.

Figure 3. Geometry of 2[bmim]PF6 used in the numeric and analytic gradient comparison. Hydrogen atoms are light gray, carbon atoms are dark gray, nitrogen atoms are blue, and fluorine atoms are green.

Using a 0.005 Å step size for the numeric gradient resulted in instances for which the forward and backward steps in the two energy calculations for each numeric gradient element had

The differences between elements of the numeric and analytic gradient are large (10−2 to 10−4) until the step size is decreased to 4754


Article


randomly reassigns the velocities to a Maxwell−Boltzman distribution every N fs, denoted Nosé−Hoover (N), was used. The details of the simulations are as follows: 1. The initial configuration was generated by randomly placing 32 water molecules in a box with a volume that matches the density of water at 300 K. 2. A 6 ps NVT classical MD simulation of the water molecules was performed with the EFP method. The temperature was set to 300 K and the time step size to 0.5 fs. A Nosé−Hoover (500) thermostat was used to regulate the temperature. 3. The last configuration of the previous run was used as the initial configuration for a 500 fs NVT equilibration run performed with the EFMO method, with a 1.0 fs time step size and a Nosé−Hoover (100) thermostat, at 300 K. The 6-31++G(d,p) basis set and Rcut = 0.3 were used. This set of equilibration runs was done to match previous MD simulations used to check energy conservation.8 Periodic boundary conditions were not used, since this work is a test of the gradient, not a production simulation. As discussed by Nakata et al.33 and Brorsen et al.,8 the energy conservation in an NVE simulation using the Velocity-Verlet algorithm can be tested by comparing the RMSD(E) to the time step size. The RMSD(E), for M steps, is calculated as

0.0001 Å or below. Using a 0.0001 or 0.00005 Å step size results in fewer instances in which the numeric gradient forward and backward steps have density components allocated to different expansion points. Once the step size decreases to 0.0001 Å, the RMSD and maximum gradient difference between the numeric and analytic gradient were both on the order of 10−6 or 10−7 Hartree/Bohr. Similar to the previous case, this suggests that there may be discontinuities in the potential energy surface. Small step sizes in numeric gradients can be suspect, so in this case, the there is not an accurate numeric to analytic gradient comparison. However, when the multipole moment allocation algorithm is modified to always allocate to the same expansion points (this was done by modifying the nearest-site allocation algorithm to choose not the nearest atom, but rather the nearest of the two atoms upon which the two Gaussian basis functions that comprise the piece of charge density are centered), the RMSD between the numeric and analytic gradient is on the order of 10−6 Hartree/Bohr and the maximum gradient difference between the numeric and analytic gradient is on the order of 10−5 Hartree/Bohr, both in the acceptable range. 6.1.4. Discussion of Potential Energy Surface. In two of the above cases, the numeric and analytic gradients did not match until the step size was decreased or when the allocation algorithm for the multipole moments was changed to always allocate density components to the same expansion point for forward and backward displacements. This can be understood by considering how the multipole moments are calculated. To calculate the multipole moments, the nearest-site allocation algorithm is used to place multipole moments on expansion centers. The nearestsite allocation algorithm involves evaluating multipole moments at every Gaussian basis function overlap center (that is, at each piece of charge density), and then shifting the multipole moments to the nearest expansion center. In the EFMO method, all expansion centers are atom centers. If, during a MD simulation or geometry optimization, atoms in a single fragment move in such a way that the multipole moments at a Gaussian basis function overlap center are suddenly closer to a different atom center, then the multipole moments on the atoms are calculated differently, and thus the final energy is different. If the energy is significantly different, the PES will not be smooth. As noted above, one way to solve this problem is to require the energy calculations to use the same set of expansion centers. Although this is useful for testing the gradient, this is not an ideal solution since it changes the energy calculation. Alternatively, it is possible that including bond midpoints as expansion centers in the EFMO multipole expansion (as is done in the EFP method) might decrease or eliminate the problem. This possibility will be explored. Currently, for the systems studied, the maximum difference between the numeric and analytic gradient is on the order of 10−5 Hartree/Bohr or less, and the RMSD is on the order of 10−6 Hartree/Bohr or less once the numeric gradient step size is small enough or if the allocation algorithm is modified so that the same expansion points are always used for the forward and backward steps in the numerical gradient procedure. The small differences between the analytic and numerical gradients imply that the analytic gradient is accurate. In addition, for small molecules such as water, displacing the atoms by 0.005 Å in the forward and backward directions generally uses the same expansion points. The problem discussed here is most likely to arise for larger molecules. 6.2. MD Simulations. A system of 32 water molecules was equilibrated before each production run, as summarized below. For several of the steps, a Nosé−Hoover thermostat that

M

∑ j (Energy at step j − Average energy of all steps)2 M (6.2)

For the Velocity-Verlet algorithm the relationship between the MD simulation time step size and the RMSD(E) for NVE ensembles should be RMSD(E) ∝ (time step size)2

(6.3)

Equation 6.3 can be rewritten as log(RMSD(E)) ∝ 2 log(time step size)

(6.4)

Thus, a log−log plot should show a straight line with a slope of about 2. To check that the EFMO MD simulation using the analytic energy gradient closely follows eq 6.3, seven NVE EFMO MD simulations were run for 50 fs each. The initial configuration and velocity were taken from the last step of the equilibration runs. The seven runs had time step sizes of 0.1, 0.2, 0.25, 0.35, 0.5, 0.6, and 0.75 fs and Rcut = 0.3. Figure 4 shows a log−log plot of the time step size vs the RMSD(E) for the 7 runs. The plot shows a straight line with a slope of about 2.03, close to that which is expected when the energy is conserved. This suggests that for the seven time step sizes, EFMO MD simulations using the analytic gradient properly conserve energy in NVE ensembles.

7. TIMINGS Timing comparisons between EFMO/MP2 and FMO2/MP2 gradient calculations, with the 6-31++G(d,p) basis set, are presented in Table 3. All calculations were done on four compute nodes. Each compute node has two quad-core 3.0 GHz Intel Xenon E5450 CPUs connected by Mellanox 4X DDR Infiniband. Multi-level parallelism with GDDI was used to split each calculation into four groups. The timings were done for Rcut = 1 and Rcut = 2. As can be seen in the table, the EFMO/MP2 method gives a speed up ranging from 1.58× to 2.06× compared to 4755


Article


As demonstrated in Section 7, the EFMO gradient is up to 2.06× faster than the FMO gradient. In testing the analytic gradient, it was discovered that the allocation of charge density in the multipole moment calculation during a numeric gradient calculation for a large molecule can differ for the forward and backward steps, thereby causing the numeric and analytic gradients to differ by too much. It is anticipated that adding bond midpoints as expansion points should decrease the impact of the allocation difference. This will be explored in a future paper. Future work will also include the derivation and implementation of the gradient of the chargetransfer term.

■

Figure 4. A log−log plot of seven EFMO MD simulations of a 32-water cluster in the NVE ensemble using six different time step sizes vs the RMSD(E) of the energy.

A1. Coupled Perturbed Hartree−Fock Equation

Since the molecular orbital coefficients are calculated using a variational minimization, the derivative of the MO coefficients can be calculated using the derivative of the variational conditional. For the RHF self consistent field (SCF) equation, the variational condition is that Fia = Fai = 0, where i is an occupied orbital and a is virtual orbital. For canonical molecular orbitals, the variational condition becomes Fpq = Fqp = 0 where p ≠ q. This was described in detail previously,9 and results in the CPHF equation

Table 3. Timing Comparison for EFMO/MP2 and FMO2/ MP2 Gradient Calculations on Water Clusters Rcut = 1

20 water molecules 30 water molecules 40 water molecules 64 water molecules

Rcut = 2

EFMO wall clock time (s)

FMO wall clock time (s)

FMO time/ EFMO time

EFMO wall clock time (s)

FMO wall clock time (s)

FMO time/ EFMO time

20.00

33.00

1.65

35.90

61.30

1.71

30.20

52.20

1.73

62.60

122.50

1.96

42.30

69.80

1.65

74.20

138.30

1.86

64.40

101.50

1.58

134.90

277.60

2.06

APPENDIX A: RESPONSE EQUATIONS, DISCUSSION OF RESPONSES, AND THE Z-VECTOR METHOD

vir occ

∑ ∑ A pq ,ck Uckx = Bpqx c

(A.1)

k

where A′pq , ck = 4(pq|ck) − (pc|qk) − (pk|qc) A pq , ck = (ϵq − ϵp)δpcδqk − A′pq , ck

FMO2/MP2. Recall that EFMO includes explicit many-body interactions via the self-consistent EFP polarizability, whereas FMO2 does not. As seen in ref 7, EFMO can attain the same level of accuracy as FMO but with a smaller Rcut value. Thus, it is possible that the speed-up might be greater. As mentioned above, the dimension of the largest response equation is the dimension of the largest monomer. In this work, the timings were obtained for systems in which the largest monomer is a water molecule. Since the dimension of the largest response equation in one iteration of the self-consistent Z-vector method in the FMO gradient30 is also the dimension of the largest monomer, the comparison should hold for larger molecules.

occ x Bpq = F (pqx) − S(pqx)ϵq −

∑ Skj(x)(2(pq|kj) − (pk|qj)) kj

and ϵq is the orbital energy of MO q. Only the virtual−occupied block of the response matrix is uniquely defined. For the RHF energy, unitary transformations between occupied−occupied (occ-occ) and virtual−virtual (vir-vir) MOs do not change the energy (but must still follow the orthonormality constraint on orbitals), and are not uniquely defined. They are often referred to as “non-independent”, while the vir-occ block is “independent”.9 If it is assumed that the CMOs are used (Fij = Fji = δijϵi, where i is an occupied orbital and j is an occupied orbital), then the occ-occ part of the response matrix can be written as

8. CONCLUSIONS As shown in Section 6, the current implementation of the EFMO gradient is fully analytic for the Coulomb, exchange-repulsion, polarization, and dispersion terms. For the EFP interaction energy part of the gradient, response equations must be solved, but the response equations are separable, so there is no response equation with the dimension of the full system. That is, for the EFP interaction energy part of the gradient, the dimension of the largest response equation is the dimension of the largest monomer. If the gradient of the chosen ab initio method has response terms, then the gas-phase monomer and dimer energy gradients will have response terms. The monomer response terms can be combined with the EFP interaction energy response terms, but the dimer response terms must be solved separately.

vir occ

Uijx =

1 (∑ ∑ Aij′ , ck Uckx + Bijx ) (ϵj − ϵi) c k

(A.2)

although it is undefined when there is a degeneracy in orbital energy. Alternately, if the energy expression is invariant to unitary transformations among the CMO occupied orbitals or the virtual orbitals then34,35 1 Uijx = − Sij(x) 2 4756

1 (x) x Uab = − Sab 2

(A.3)


Article


(or similar localization methods) in detail,13,15,16,36 so the result is presented here. An overview of the derivation used in the code implemented for this study is below. Equation A.5 is expanded and rearranged to form an equation for the localization response matrix (as is similarly done in the CPHF equation derivation). The derivative of the LMO coefficient as shown in eq 2.3 is used throughout. Note also that vxnl is antisymmetric. Using

for the occ-occ or vir-vir part of the response matrix. This can only be used in certain formulations of gradients. For example, eq A.3 is valid if the gradient is directly derived from an energy formula that does not assume a particular unitary transform of the MOs. The CPHF equation can be solved for each canonical response matrix element. In practice, this is avoided, as described in the Z-vector method. A2. Derivative of the Orthonormality Constraint with Respect to a Nuclear Perturbation

Spq = δpq implies Uxpq + Uxqp = −S(x) pq when x is a nuclear perturbation (and the basis set depends on nuclear coordinates). For a field perturbation, α, Uαpq + Uαqp = 0 since the basis set does not depend on the field.

∂l αm ∂x

=

∂x

occ CMO

∑ i occ CMO

+

=

∑ i occ CMO

=

∑ i occ CMO

=

∑ i LMO

=

∑ n

∂Lli ∂x

∑

⎣

i

LMO

=

∑

+

cμLnvnlx +

n

∑ ∑ j

∑

−

Uqixcμq

+

∑

Lli

occ+vir CMO

∑

Lli

∑ q

occ+vir CMO

∑

Lli

i

∑

− Lmj rqm) +

∑ cμLl μ

∂μ αl ∂x

(A.6)

After rearranging, and separating out the response matrices, eq A.6 can be written as

Uqixcμq

occ CMO

∂μ αm ∂x

⎞ ∂μ ⎟ αm ⎟=0 ⎟ ∂x ⎠

cμLm

μ

i

occ+vir CMO

occ CMO

⎛ AO ∂μ αm (rll − rmm)⎜⎜∑ cμLl ∂x ⎝ μ

Uqixcμq

q

⎛ AO ∂μ αl + 2rlm⎜⎜∑ cμLl ∂x ⎝ μ

Uqixcμq

q

AO

+

∑ cμLm μ

AO

−

∑ cμLm μ

⎞ ∂μ α l ⎟⎟ ∂x ⎠

⎞ ∂μ α m ⎟⎟ ∂x ⎠

LMO

As previously known,16 the localization condition for Boys LMOs is

+

∑

vnox(δol rnm(roo − rmm) + δom rnl(rll − roo)

o

Analytic Gradients for the Effective Fragment Molecular Orbital Method

Recommend Documents