Acceleration of Environmental Electrostatic Potential Using Cholesky Decomposition with Adaptive Metric (CDAM) for Fragment Molecular Orbital-based Molecular Dynamics (FMO-MD) Simulation

Abstract

The fragment molecular orbital (FMO) method is an efficient quantum chemical method suitable for calculating the electronic structures of large molecular systems. FMO can be accelerated by several approximations, an important one being the approximation to the environmental electrostatic potential (ESP) exerted on the fragment monomers or dimers. The environmental ESP is often approximated using the Mulliken atomic orbital charge (AOC) for proximal fragment dimers (ESP-AOC approximation) and the Mulliken point charge (PTC) for distant dimers (ESP-PTC approximation). Recently, another approximation method based on Cholesky decomposition with adaptive metric has been proposed for environmental ESP (ESP-CDAM, Okiyama et al., Bull. Chem. Soc. Jpn., 2021, 94, 91). In the current article, the energy gradient is derived under the ESP-CDAM approximation and implemented in FMO-based molecular dynamics (FMO-MD). Several test FMO-MD simulations are performed to compare the accuracy of the ESP-CDAM approximation with that of the conventional methods. The results show that ESP-CDAM is more accurate than ESP-AOC.

1 INTRODUCTION

The ab initio fragment molecular orbital (FMO) method, first published in 1999 [1], is a highly efficient and accurate tool for electronic structure calculations of large molecular systems, especially biomolecules [2,3]. In FMO, a molecular system is divided into fragments; self-consistent field (SCF) calculations are performed for fragment monomers, dimers (FMO2) [1], trimers (FMO3) [4], tetramers (FMO4) [5], and so on; and finally the total energy of the system is reconstructed. FMO2 is usually used as the default option, but the accuracy of FMO can be systematically improved by including higher terms such as FMO3 and FMO4. However, to achieve the full improvement in FMO performance, we need to accelerate FMO2 because FMO3 and FMO4 are performed by adding correction terms to FMO2.

FMO2 can be accelerated by several approximations. The most important is the dimer-ES approximation [6], which reduces the overall cost by replacing the SCF calculation of a distant fragment pair with a simple electrostatic calculation. In addition, we can accelerate the environmental electrostatic potential (ESP) terms to further improve the throughput performance of FMO2. In the conventional FMO2 protocol [7,8], the environmental ESP is often calculated by the ESP-AOC approximation for the proximal region of a monomer using the atomic orbital (AO) charge (AOC) population analysis, and also by the ESP-PTC approximation for the distant region using the Mulliken atomic point charge (PTC). However, in another FMO2 protocol for environmental ESP [9], the ESP from the proximal region is treated by the ESP-CDAM approximation based on the Cholesky decomposition [10] with adaptive metric (CDAM) [11].

In this article, we report on the implementation of the latter protocol, namely the approximation of the environmental ESP using CDAM, in the ABINIT-MP program [3, 12, 13]. We present the semi-analytic formula of the potential energy gradient (negative force) at the Hartree-Fock (HF) and second-order Møller-Plesset perturbation theory (MP2) levels under the approximation. We show single-point energy and gradient calculations [14] as well as FMO-based molecular dynamics (FMO-MD) simulations for several molecular systems [15,16,17,18] to demonstrate the power of the ESP-CDAM/FMO implemented here.

2 METHODS

In the subsequent mathematical discussion, a fragment monomer is referred to simply as a “monomer” and a pair of fragments as a “dimer.” Indices I, J, and K specify monomers; i, j, and k, orbitals; and μ, ν, λ and σ, AOs. Atomic nuclei are represented by indices A, B, and C; their charges by Z_A, Z_B, and Z_C; and their 3D coordinates by A, B, and C.

2.1 2.1 HF/FMO2.

We review the outline of the FMO2 energy and gradient calculation at the HF level [1,14] to facilitate description of the implementation of ESP-CDAM to FMO2. See Figure 1 for a schematic diagram of the FMO2 algorithm.

Figure 1.

 Algorithm of FMO2.

In FMO2, the total energy of the molecular system, E, is approximated by the following equation, where E_I and E_IJ refer to the total energies of monomer I and dimer IJ, respectively, and N_f is the number of fragments.

  

$\begin{array}{c}E\approx \sum _{I>J}{E}_{IJ}-\left(N_f-2\right)\sum _I{E}_I\end{array}$(2.1.1)

Similarly, the electron density of the molecular system, ρ(r), is approximated as follows.

  

$\begin{array}{c}\rho \left(r\right)\approx \sum _{I>J}{\rho }_{IJ}\left(r\right)-\left(N_f-2\right)\sum _I{\rho }_I\left(r\right)\end{array}$(2.1.2)

     

The electronic structure of each monomer is calculated with an SCF loop, which is iterated until all the monomer densities become self-consistent with each other by the self-consistent charge (SCC) loop. In contrast, the SCF loop for the dimers is performed only once with the charges of monomers obtained in the previous SCC loop.

Below, we describe the monomer and dimer SCF calculations introduced in equations (2.1.3–11), namely, how to calculate equations (2.1.1) and (2.1.2) under the HF approximation. We postulate that each of the monomers and dimers forms a closed shell. Then, the Hartree‒Fock‒Roothaan equation [19] of a monomer or dimer X( = I or IJ) becomes

  

$\begin{array}{c}{F}^X{C}^X=S^X{C}^X{\epsilon }^X\end{array}$(2.1.3)

using the variables defined below.

  

$\begin{array}{c}{F}^X=H^X+G^X\end{array}$(2.1.4)

  

$\begin{array}{c}{H}_{\mu \nu }^X=H_{\mu \nu }^{\left(\text{core}\right)X}+V_{\mu \nu }^X+\sum _k{B}_k\mu \text{|}{\theta }_k{\theta }_k\text{|}\nu \end{array}$(2.1.5)

     

  

$\begin{array}{c}{V}_{\mu \nu }^X=\sum _{K\ne X}\left(u_{\mu \nu }^K+v_{\mu \nu }^K\right)\end{array}$(2.1.6)

  

$\begin{array}{c}{u}_{\mu \nu }^K=\sum _{A\in K}\mu \left|\left(-Z_A/\left|r-A\right|\right)\right|\nu \end{array}$(2.1.7)

  

$\begin{array}{c}{v}_{\mu \nu }^K=\sum _{\lambda \sigma \in K}{D}_{\lambda \sigma }^K\left(\mu \nu |\sigma \lambda \right)\end{array}$(2.1.8)

  

$\begin{array}{c}{D}_{\mu \nu }^X=2\sum _{i=1}^{occ}{C}_{\mu i}^{\text{*}X}{C}_{\nu i}^X\end{array}$(2.1.9)

  

$\begin{array}{c}{G}_{\mu \nu }^X=\sum _{\lambda \sigma \in X}{D}_{\lambda \sigma }^X\left[\left(\mu \nu |\sigma \lambda \right)-\frac{1}{2}\left(\mu \lambda |\sigma \nu \right)\right]\end{array}$(2.1.10)

  

$\begin{array}{c}{E}_X=\frac{1}{2}\sum _{\mu \nu \in X}{D}_{\nu \mu }^X\left(H_{\mu \nu }^X+F_{\mu \nu }^X\right)+\sum _{A\in X}\sum _{A>B\in X}\frac{Z_A{Z}_B}{\left|A-B\right|}\end{array}$(2.1.11)

$V_{\mu \nu }^X$(2.1.11)

appearing in equation (2.1.5), an element of matrix V^X, is the environmental ESP comprised of ESPs of nuclei (u^K_µν) and electrons (v^K_µν) of surrounding monomers (Eq. 2.1.6). The third term of the right-hand-side of equation (2.1.5) is the projection operator that removes localized molecular orbital θ_k from the variational space; in this term, parameter B_k should be set to 10⁶ or larger. A molecule is usually fragmented at an sp³ carbon, and the atom located at the fragmentation point is called the bond detached atom (BDA). Z_A in Eqs (2.1.7) and (2.1.11) corresponds to the nuclear charge of BDAs after the partitioning of the charge [14].

The naive algorithm of FMO2 described above can be accelerated by several approximations, two of which are especially effective [7]. One is the approximation of the environmental ESP contained in the monomer and dimer, equation (2.1.3), and the other is the dimer-ES approximation for distant dimers [6].

To introduce approximations to the environmental ESP, we need to reformulate the FMO2 energy, equation (2.1.1), as

  

$\begin{array}{c}E\approx \sum _{I>J}\left(E{\text{'}}_{IJ}-E{\text{'}}_I-E{\text{'}}_J\right)+\sum _{I>J}\text{Tr}\left(\Delta D^{IJ}{V}^{IJ}\right)\\ +\sum _{I}E{\text{'}}_I.\end{array}$(2.1.12)

Here, E´_X is the energy of monomer or dimer X from which the environmental ESP has been deleted. Defining the difference density matrix (ΔD) by

  

$\begin{array}{c} \Delta D^{IJ}=D^{IJ}-D^I\oplus D^J,\end{array}$(2.1.13)

we can apply different ESP approximations for monomers and dimers using equation (2.1.12), which contains ESP for dimers only [7]. $\Delta {\tilde{E}}_{IJ}$ defined by equation (2.1.14) can be interpreted as the interaction energy between fragments I and J(inter-fragment interaction energy, IFIE; or pair interaction energy, PIE) [7, 20, 21, 22].

  

$\begin{array}{c}\Delta {\tilde{E}}_{IJ}=\left(E{\text{'}}_{IJ}-E{\text{'}}_I-E{\text{'}}_J\right)+Tr\left(\Delta D^{IJ}{V}^{IJ}\right)\end{array}$(2.1.14)

Under the ESP-CDAM approximation [9], the environmental ESP acting on a monomer or dimer, $v_{\mu \nu }^K$, reduces to

  

$\begin{array}{c} \left(\mu \nu |\lambda \sigma \right)\approx \sum _{MN}\left(\mu \nu \text{|}M\right)G_{MN}^{-1}\left(N\text{|}\lambda \sigma \right).\end{array}$(2.1.15)

When matrix G is symmetric positive semidefinite, the decomposition of its inverse can be given by

  

$\begin{array}{c}{G}_{MN}^{-1}\approx \sum _R{Z}_{MR}{Z}_{NR}.\end{array}$(2.1.16)

In the CDAM method, the Cholesky basis, $L_{R, \mu \nu }$, is calculated as,

  

$\begin{array}{c}{L}_{R,\mu \nu }=\sum _M{Z}_{MR}\left(M\text{|}\mu \nu \right),\end{array}$(2.1.17)

where M is an AO pair from which the linear dependence has been deleted, and G_MN = (M|N). Let K be the index of the monomer exerting the environmental ESP, whose Cholesky basis is $L_{R, \mu \nu }^K$. Defining $Q_R^K$ and $Q{\text{'}}_M^K$ by equations (2.1.18) and (2.1.19), respectively, and calculating $Q{\text{'}}_M^K$ we obtain the environmental ESP acting on X from the electrons of monomer K, $v_{\mu \nu }^K$, simply by calculating a pseudo-three-center integral (μν | M).

  

$\begin{array}{c}{Q}_R^K=\sum _{\lambda \sigma \in K}{L}_{R,\lambda \sigma }^K{D}_{\lambda \sigma }^K\end{array}$(2.1.18)

  

$\begin{array}{c}Q{\text{'}}_M^K=\sum _R{Z}_{MR}^K{Q}_R^K\end{array}$(2.1.19)

  

$\begin{array}{c}{v}_{\mu \nu }^K\approx \sum _M\left(\mu \nu \text{|}M\right)Q{\text{'}}_M^K\end{array}$(2.1.20)

The dimer-ES approximation of ESP for distant dimers is expressed as

  

$\begin{array}{c}{E^{\prime }}_{IJ}\approx {E^{\prime }}_I+{E^{\prime }}_J+Tr\left(D^I{u}^J\right)+Tr\left(D^J{u}^I\right)\\ +\sum _{\mu \nu \in I}\sum _{\sigma \lambda \in J}{D}_{\mu \nu }^I{D}_{\sigma \lambda }^J\left(\mu \nu |\sigma \lambda \right)+\sum _{A\in I}\sum _{B\in J}\frac{Z_A{Z}_B}{\left|A-B\right|},\end{array}$(2.1.21)

which avoids SCF calculations for many dimers.

We have so far described only the energy calculations. We describe additional formula necessary to calculate the energy gradient of FMO2 [14]. The contribution of the third term of equation (2.1.5), i.e., the projection operator, is neglected because it would amount to only $O\left(B_k^{-1}\right)$. Note that this contribution can be incorporated easily if needed [23]. To obtain the energy gradient via FMO2-HF, we need the reverse environmental ESP from the electrons of a dimer to a monomer, which is calculated by

  

$\begin{array}{c}{v}_{\lambda \sigma }^{IJ}=\sum _{\mu \nu \in IJ}\Delta D_{\mu \nu }^{IJ}\left(\mu \nu |\sigma \lambda \right)\end{array}$(2.1.22)

if no ESP approximations are applied. In the ESP-CDAM approximation, $v_{\lambda \sigma }^{IJ}$ is approximated as

  

$\begin{array}{c}{v}_{\lambda \sigma }^{IJ}\approx \sum _M\left(\mu \nu \text{|}M\right)Q{\text{'}}_M^{IJ}\end{array}$(2.1.23)

$Q{\text{'}}_M^{IJ}$ in equation (2.1.23) can be obtained by replacing monomer index K in equations (2.1.18) and (2.1.19) with dimer index IJ. The one-center CD (1C-CD) approximation [24] is introduced to make the potential energy surface smoother than the full-CD does while accelerating ESP-CDAM [25].

In the differentiation of ESP from the surrounding N_f − 2 fragments in the dimer calculation, from a strict viewpoint, it is necessary to solve the coupled perturbed Hartree‒Fock (CPHF) equation because the MO of the monomer has not been determined variationally [26, 27]. Nonetheless, the utilization of equations in Ref. [14] was shown to give sufficient precision based on a numerical analysis.

The FMO2 energy is thus calculated at the HF level by the above formulas.

Note that we have chosen not to use the fully analytical energy gradient by the Self-Consistent Z-vector (SCZV) method [23, 26, 27] for CDAM, despite its accuracy superiority over the original semi-analytical gradient [14], due to a particular drawback. The SCZV approach requires that ESP is not switched by distance, which makes only the ESP-PTC approximation valid and invalidates more practical ESP approximations, such as ESP-AOC-PTC. Unfortunately, the ESP-PTC approximation is too coarse to be applicable in FMO-MD, while acceleration algorithms such as Ewald and PPPC are required for FMO-MD simulations of large molecules. Since these methods divide the interactions into proximal and distant ones, they cannot be implemented with the SCZV method. We have therefore decided not to include SCZV in CDAM.

2.2 MP2/FMO2.

In ABINIT-MP, the implementation of the MP2 first derivative (or nuclear gradient) [29,30,31,32] was based on a series of expressions given by Ref. [28], where an AO‒MO hybrid strategy was adopted for an efficient parallelized integral-direct processing of one-body density matrix generations with MO indices. These generated D and W are to be back-transformed to the corresponding AO-based matrices, as was just done for the non-separable two-body density term (Г^NS) evaluated as the MP2 amplitudes. The construction of the energy-weighted density matrices W implemented is somewhat lengthy, as shown in equations 181 to 187 of Gordon’s paper [28]; therefore, the listing is omitted here. For simplicity of description, the MP2 derivative equations including the HF contributions are hereafter given AO label. Namely, the gradient with respect to nuclear coordinate a( = x, y, or z element of A) is expressed as

  

$\begin{array}{c}\frac{\partial E^{\left(2\right)}}{\partial a}=\sum _{\mu \nu }{D}_{\mu \nu }^{\left(2\right)}\frac{\partial }{\partial a}{H}_{\mu \nu }^{\text{core}}+\sum _{\mu \nu }{W}_{\mu \nu }^{\left(2\right)}\frac{\partial }{\partial a}{S}_{\mu \nu }\\ +\sum _{\mu \nu \lambda \sigma {\text{Γ}}_{\mu \nu \lambda \sigma }^{\left(2\right)}}\\ \partial \end{array}\partial a\mu \nu |\lambda \sigma ;$(2.2.1)

where

  

$\begin{array}{c}{\text{Γ}}_{\mu \nu \lambda \sigma }^{\left(2\right)}={\text{Γ}}_{\mu \nu \lambda \sigma }^{\text{(2)NS}}+{\text{Γ}}_{\mu \nu \lambda \sigma }^{\left(2\right)\text{S}}\end{array}$(2.2.2)

and

  

$\begin{array}{c}{\text{Γ}}_{\mu \nu \lambda \sigma }^{\left(2\right)S}=D_{\mu \nu }^{\left(2\right)}{D}_{\lambda \sigma }^{\text{HF}}-\frac{1}{2}{D}_{\mu \sigma }^{\left(2\right)}{D}_{\lambda \nu }^{\text{HF}}.\end{array}$(2.2.3)

The total FMO2-MP2 energy, $E^{\text{FMO2-MP2}}$, is expressed as

  

$\begin{array}{c}{E}^{\text{FMO2-MP2}}=\sum _I\left({E^{\prime }}_I+E_I^{\left(2\right)}\right)\\ +\sum _{I>J}\left[\Delta {\tilde{E}}_{IJ}+\left(E_{IJ}^{\left(2\right)}-E_I^{\left(2\right)}-E_J^{\left(2\right)}\right)\right].\end{array}$(2.2.4)

The first derivative of $E^{\left(2\right)X}$ by a has a form similar to that of the conventional MP2 method,

  

$\begin{array}{c}\frac{\partial E^{\left(2\right)X}}{\partial a}=\sum _{\mu \nu }{D}_{\mu \nu }^{\left(2\right)X}\frac{\partial }{\partial a}{H}_{\mu \nu }^X+\sum _{\mu \nu }{W}_{\mu \nu }^{\left(2\right)X}\frac{\partial }{\partial a}{S}_{\mu \nu }^X\\ +\sum _{\mu \nu \lambda \sigma }{ﾃ}_{\mu \nu \lambda \sigma }^{\left(2\right)X}\frac{\partial }{\partial a}\left(\mu \nu | \lambda \sigma \right)X,\\ \end{array}$(2.2.5)

where X = I for a monomer and X = IJ for a dimer. The following holds if we denote the environmental ESP calculated by the electron density obtained by the HF level monomer SCC as $V_{\mu \nu }^X$,

  

$\begin{array}{c}{H}_{\mu \nu }^X=H_{\mu \nu }^{\left(\text{core}\right)X}+V_{\mu \nu }^X\end{array}$(2.2.6)

If we define equation (2.2.7) as

  

$\begin{array}{c}\frac{\partial {E^{\prime }}^{\left(2\right)X}}{\partial a}=\sum _{\mu \nu }{D}_{\mu \nu \left(2\right)X}^{}\\ \partial \end{array}\partial aH\mu \nu \text{core}X+\sum \text{undOvron}\mu \nu W\mu \nu 2X\partial \partial aS\mu \nu X+\sum \text{undOvron}\mu \nu \lambda \sigma ﾃ\mu \nu \lambda \sigma 2X\partial \partial a\mu \nu |\lambda \sigma X,$(2.2.7)

then it follows that

  

$\frac{\partial }{\partial a}{E}^{\text{FMO2-MP2}}=\sum _I\frac{\partial }{\partial a}\left(E{\text{'}}_I+E{\text{'}}_I^{\left(2\right)}\right) \begin{array}{c}+\sum _{I>J}\left[\begin{array}{c}\frac{\partial }{\partial a}\Delta {\tilde{E}}_{IJ}+\frac{\partial }{\partial a}\left(E{\text{'}}_{IJ}^{\left(2\right)}-{E^{\prime }}_I^{\left(2\right)}-E{\text{'}}_J^{\left(2\right)}\right)\\ +\sum _{\mu \nu }\Delta D_{\mu \nu }^{\left(2\right)IJ}\frac{\partial }{\partial a}{V}_{\mu \nu }^{IJ}\end{array}\right],\end{array}$(2.2.8)

where

  

$\begin{array}{c}\Delta D^{\left(2\right)IJ}=D^{\left(2\right)IJ}-D^{\left(2\right)I}\oplus D^{\left(2\right)J}\end{array}$(2.2.9)

and

  

$\begin{array}{c}\sum _{\mu \nu \in IJ}\Delta D_{\mu \nu }^{\left(2\right)IJ}\frac{\partial V_{\mu \nu }^{IJ}}{\partial a}\\ =\sum _{\mu \nu \in IJ}\Delta D_{\mu \nu }^{\left(2\right)IJ}\sum _{K\ne IJ}\left[\frac{\partial u_{\mu \nu }^K}{\partial a}+\sum _{\lambda \sigma \in K}{D}_{\lambda \sigma }^K\frac{\partial \left(\mu \nu |\lambda \sigma \right)}{\partial a}\right].\end{array}$(2.2.10)

The MP2/FMO2 energy and gradient thus calculated are used in FMO-MD at the MP2 level [32].

2.3 Test calculations.

The new approximations were implemented in a local version of ABINIT-MP Open series Ver. 1.0 Rev. 15 [12]. Three molecular systems were used for test calculations: a cluster of 400 water molecules W400, the Trp-cage polypeptide (1L2Y) [33], and the GESKG droplet. The GESKG droplet is our in-house benchmark system consisting of a polarizable neutral peptide (NH₃⁺-Gly-Glu-Ser-Lys-Gly-CO₂⁻) and a surrounding sphere of 206 water molecules (Figure 2).

Figure 2.

 The GESKG droplet visualized with BioStation Viewer graphical software [41].

FMO2 calculations were performed at MP2/6-31G*. Each water molecule was regarded as a fragment. The polypeptides were fragmented so that an amino acid residue should form a fragment unit, but at C_α to avoid breaking peptide bonds. The dimer-ES approximation was applied only to those dimers separated within a preset threshold (L_dimer) measured in van der Waals units [6]. Test calculations were performed on a 32-core Intel Xenon computer at Rikkyo University and a 40-core computer of the same type at Hoshi University.

The GESKG droplet was subjected to FMO-MD simulation as the final test of the ESP-CDAM energy gradient with MP2. The droplet was energy-minimized by the Broyden–Fletcher–Goldfarb–Shanno method [35] for 50 steps, given a kinetic energy of 200 K, heated to 300 K with the Nosé‒Hoover chains thermostat [34] with a time step of 0.5 fs for 50 fs, and then subjected to several test runs with various time steps each for 50 fs without any thermostat or constraint to draw a double logarithmic plot of [time step] vs. [root mean square deviation of total energy (RMSD (E_t))]. The simulations were performed with the FMO-MD option [36] of the ABINIT-MP program.

3 RESULTS AND DISCUSSION

Before investigating ESP-CDAM, we determined the optimal value of L_dimer, the threshold for the dimer-ES approximation. The errors in energy and gradient were examined for W400, 1L2Y, and the GEKSG droplet at FMO2-MP2/6-31G* (Figures 3 and 4). The energy error dropped sharply at of L_dimer = 2.0. Therefore, we concluded that of L_dimer should be appropriate and used this value for subsequent calculations.

Figure 3.

  A. Dependence on L_dimer of the computation error in the energy at FMO2-MP2/6-31G*. The error was calculated as the difference in energy between the presence and absence of the dimer-ES approximation without the environmental ESP approximation. B. An enlarged view of the region where (1.5 ≤ L_dimer ≤ 3.0) in Figure A.

Figure 4.

  A. Dependence on L_dimer of the computation error in the gradient (RMSD) at FMO2-MP2/6-31G*. The error was calculated as the difference in gradient between the presence and absence of the dimer-ES approximation without the environmental ESP approximation. B. An enlarged view of the region where (1.5 ≤ L_dimer ≤ 3.0) in Figure A.

Next, we investigated the accuracy of the ESP-CDAM for W400 at FMO2-MP2/6-31G*. We calculated the energy and gradient with no ESP approximation (ESP-NON) as a reference and also with the ESP-CDAM and ESP-AOC approximations to evaluate the accuracy of the energy and gradient of the proximal fragment dimers. From the results, we calculated the simple difference for the total energy and the RMSD for the gradient between the reference and the ESP-CDAM and ESP-AOC data to represent the errors due to the approximations (Table 1). The error due to ESP-CDAM in the total energy was as small as −0.5578 mhartree, while the error due to ESP-AOC was 10 times larger. Similarly, the error in the energy gradient was only 0.4574 mhartree/bohr for ESP-CDAM, but as large as approximately 38 mhartree/bohr for ESP-AOC. On the other hand, the calculation speed of ESP-AOC was more than six times faster than that of ESP-CDAM. This decrease in speed is due to the overhead of Cholesky decomposition for fragments as small as water molecules.

Table 1. Errors in total energy and energy gradient due to the ESP-CDAM and ESP-AOC approximations calculated for W400, along with the wall-clock time measured on a 32-core computer.

Environmental ESPapproximation	ESP-NON	ESP-CDAM	ESP-AOC
Total energy(mhartree)	0.0	-0.5578	-11.2070
Gradient(mhartree/bohr)	0.0	0.4574	38.1551
Time (s)	6602.0	5870.3	995.2

We then tested the ESP approximations to the polypeptides, i.e., Trp-cage (Table 2) and the GESKG droplet (Table 3). We found that the acceptable errors in total energy and gradient for Trp-cage were −0.0851 mhartree and 0.7578 mhartree/bohr, respectively, and that those for the GESKG droplet were 0.7516 mhartree and 1.8183 mhartree/bohr.

Table 2. Errors in total energy and energy gradient due to the ESP-CDAM and ESP-AOC approximations calculated for Trp-cage (1L2Y), along with the elapsed time measured on a 32-core computer.

Environmental ESPapproximation	ESP-NON	ESP-CDAM	ESP-AOC
Total energy(mhartree)	0.0	-0.0851	12.2303
Gradient(mhartree/bohr)	0.0	0.7578	2.6267
Time (s)	27897.2	17462.7	17001.6

Table 3. Errors in total energy and energy gradient due to the ESP-CDAM and ESP-AOC approximations calculated for the GESKG droplet, along with the elapsed time measured on a 40-core computer.

Environmental ESPapproximation	ESP-NON	ESP-CDAM	ESP-AOC
Total energy(mhartree)	0.0	0.7516	-59.2469
Gradient(mhartree/bohr)	0.0	1.8183	1.3140
Time (s)	10345.5	5259.6	3469.7

We also compared the accuracy and speed of the ESP-CDAM approximation with that of the conventional approximation, ESP-AOC, using the GESKG droplet system as the target (Table 3). The error in total energy due to ESP-CDAM (approximately 0.8 mhartree) was significantly lower than that due to ESP-AOC (approximately −59 mhartree). The error in energy gradient due to ESP-CDAM (~1.8 mhartree/bohr) was similar to that of the other two (~1.3 mhartree/bohr), but not as pronounced as the error in total energy. The reason for the relatively low gradient accuracy compared to the energy accuracy is not clear and should be investigated in the future. The computational speed of ESP-CDAM was about 1.5 times lower than that of the ESP-AOC approximation. Thus, the ESP-CDAM approximation was more accurate than the other approximation, although slower. ESP-CDAM was only 3% slower in the case of the dehydrated polypeptides (Tables 2), but 50% slower in the case of the GESKG droplet vs. ESP-AOC. Again, we suspect that this difference in behavior is due to the overhead of Cholesky decomposition in calculating fragments as small as water molecules, the dominant species in the GESKG droplet. To overcome this overhead, we plan to speed up ESP-CDAM by reducing the number of small fragments using the tree method [37].

To test the ESP-CDAM approximation, we ran several short FMO-MP2-MD simulations of the GESKG droplet and compared the accuracies and computational speeds of ESP-CDAM and ESP-AOC with those of ESP-NON. We excluded ESP-AOC-PTC from the test because its accuracy is similar to that of ESP-AOC and also because we wanted to avoid the influence of the error in the environmental ESP due to distant fragments. We performed FMO-MD simulations with different time steps each for 50 fs to draw double logarithmic plots of [time step] vs. [RMSD (E_t)]. The plot should become a straight line with a slope of 2 in the case of the velocity‒Verlet time integration used this time; hence, deviation of the slope from 2 indicates degradation of the accuracy of the energy gradient [38].

The resulting plots showed that ESP-CDAM has an accuracy comparable to that of ESP-NON and better than that of ESP-AOC (Figure 5). However, all three plots deviated from the ideal slope of 2. The deviation could have been mitigated, if not avoided, by using the analytical gradient [8, 39] and/or by including higher terms [40]. Nevertheless, the current result is sufficient for the evaluation of ESP-CDAM, since we focus only on the difference between ESP-CDAM and the reference ESP-NON or the conventionally used ESP-AOC. Taken together, the present results indicate that ESP-CDAM may be a promising alternative to ESP-AOC, especially for gradient calculation in FMO-MD.

Figure 5.

 Double logarithmic plots of [Time step] vs [RMSD (Et)] of three ESP approximations.

4 CONCLUSION

We investigated the accuracy of the newly implemented ESP-CDAM approximation for environmental ESP. The errors in energy and gradient due to ESP-CDAM against the reference values (ESP-NON) were less than 1 mhartree and about 0.5–2 mhartree/bohr, respectively, for W400, Trp cage, and the GESKG droplet. The errors were similar to or smaller than those obtained using the conventional ESP-AOC approximation. Additional FMO-MD testing of the GESKG droplet showed that the ESP-CDAM gradient was practically as accurate as that of ESP-NON and significantly better than that of ESP-AOC. However, further improvement of the accuracy of the energy gradient in the ESP-CDAM approximation is a future challenge.

The computational speed of ESP-CDAM was twice that of ESP-NON, but significantly slower than that of ESP-AOC. We attribute the slower speed of ESP-CDAM to the overhead of Cholesky decomposition for fragments as small as water molecules. We nonetheless believe that the ESP-CDAM option can be accelerated by using the tree method, which approximates multiple fragments as a multipole [37].

ACKNOWLEDGMENTS

We thank Prof. Kaori Fukuzawa of Osaka University and Dr. Chiduru Watanabe of RIKEN for their support and cooperation. Some of the FMO calculations were performed on the Oakforest-PACS supercomputer (project ID: hp200101) and the Fugaku supercomputer (project ID: hp210130). This study was supported by JSPS Kakenhi (JP19K12010), the Platform Project for Supporting Drug Discovery and Life Science Research (Basis for Supporting Innovative Drug Discovery and Life Science Research, BINDS) from the Japan Agency for Medical Research and Development (AMED, grant number JP21am0101113), and the FMO Drug Design Consortium (FMODD). Y. M would thank Rikkyo University Special Fund for Research (SFR).

REFERENCES

• [1]   K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett., 313, 701 (1999). doi:10.1016/S0009-2614(99)00874-X
• [2]   D. G. Fedorov, T. Nagata, K. Kitaura, Phys. Chem. Chem. Phys., 14, 7562 (2012). , doi:10.1039/c2cp23784a PMID:22410762
• [3]   S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, K. Fukuzawa, Phys. Chem. Chem. Phys., 16, 10310 (2014). , doi:10.1039/C4CP00316K PMID:24740821
• [4]   D. G. Fedorov, K. Kitaura, J. Chem. Phys., 120, 6832 (2004). , doi:10.1063/1.1687334 PMID:15267582
• [5]   T. Nakano, Y. Mochizuki, K. Yamashita, C. Watanabe, K. Fukuzawa, K. Segawa, et al., Chem. Phys. Lett., 523, 128 (2012). doi:10.1016/j.cplett.2011.12.004
• [6]   T. Nakano, T. Kaminuma, T. Sato, Y. Akiyama, M. Uebayasi, K. Kitaura, Chem. Phys. Lett., 318, 614 (2000). doi:10.1016/S0009-2614(00)00070-1
• [7]   T. Nakano, T. Kaminuma, T. Sato, K. Fukuzawa, Y. Akiyama, M. Uebayasi, et al., Chem. Phys. Lett., 351, 475 (2002). doi:10.1016/S0009-2614(01)01416-6
• [8]   T. Nagata, D. G. Fedorov, K. Kitaura, Chem. Phys. Lett., 475, 124 (2009). doi:10.1016/j.cplett.2009.05.004
• [9]   Y. Okiyama, T. Nakano, C. Watanabe, K. Fukuzawa, Y. Komeiji, K. Segawa, et al., Bull. Chem. Soc. Jpn., 94, 91 (2021). doi:10.1246/bcsj.20200227
• [10]   T. B. Pedersen, F. Aquilante, R. Lindh, Theor. Chem. Acc., 124, 1 (2009). doi:10.1007/s00214-009-0608-y
• [11]   Y. Okiyama, T. Nakano, K. Yamashita, Y. Mochizuki, N. Taguchi, S. Tanaka, Chem. Phys. Lett., 490, 84 (2010). doi:10.1016/j.cplett.2010.03.001
• [12]   ABINIT-MP-Open-Series, http://www.cenav.org/abinitmpopen1/.
• [13]   Y. Mochizuki, T. Nakano, K. Sakakura, Y. Okiyama, H. Watanabe, K. Kato, et al., in Recent Advances of the Fragment Molecular Orbital Method, ed. by Yuji Mochizuki, Shigenori Tanaka, Kaori Fukuzawa, Springer Singapore, 2021, p. 53.
• [14]   K. Kitaura, S. I. Sugiki, T. Nakano, Y. Komeiji, M. Uebayasi, Chem. Phys. Lett., 336, 163 (2001). doi:10.1016/S0009-2614(01)00099-9
• [15]   Y. Komeiji, T. Nakano, K. Fukuzawa, Y. Ueno, Y. Inadomi, T. Nemoto, et al., Chem. Phys. Lett., 372, 342 (2003). doi:10.1016/S0009-2614(03)00430-5
• [16]   Y. Komeiji, Y. Mochizuki, T. Nakano, D. G. Fedorov, J. Mol. Struct. THEOCHEM, 898, 2 (2009). doi:10.1016/j.theochem.2008.07.001
• [17]   Y. Komeiji, Y. Mochizuki, T. Nakano, H. Mori, in Molecular Dynamics: Theoretical Developments and Applications in Nanotechnology and Energy, ed. by Lichang Wang, InTech, Goatia, 2012, p. 3
• [18]   T. Nagata, D. G. Fedorov, K. Kitaura, Theor. Chem. Acc., 131, 1136 (2012). doi:10.1007/s00214-012-1136-8
• [19]   A. Szabo, N. S. Ostlund, Modern Quantum Chemistry, Dover, New York 1996.
• [20]   Y. Mochizuki, K. Fukuzawa, A. Kato, S. Tanaka, K. Kitaura, T. Nakano, Chem. Phys. Lett., 410, 247 (2005). doi:10.1016/j.cplett.2005.05.079
• [21]   S. Amari, M. Aizawa, J. Zhang, K. Fukuzawa, Y. Mochizuki, Y. Iwasawa, et al., J. Chem. Inf. Model., 46, 221 (2006). , doi:10.1021/ci050262q PMID:16426058
• [22]   D. G. Fedorov, K. Kitaura, J. Comput. Chem., 28, 222 (2007). , doi:10.1002/jcc.20496 PMID:17109433
• [23]   T. Nagata, D. G. Fedorov, K. Kitaura, Chem. Phys. Lett., 492, 302 (2010). doi:10.1016/j.cplett.2010.04.043
• [24]   F. Aquilante, R. Lindh, T. Bondo Pedersen, J. Chem. Phys., 127, 114107 (2007). , doi:10.1063/1.2777146 PMID:17887828
• [25]   F. Aquilante, M. G. Delcey, T. B. Pedersen, I. F. Galván, R. Lindh, Mol. Phys., 115, 2052 (2017). doi:10.1080/00268976.2017.1284354
• [26]   T. Nagata, K. Brorsen, D. G. Fedorov, K. Kitaura, M. S. Gordon, J. Chem. Phys., 134, 124115 (2011). , doi:10.1063/1.3568010 PMID:21456653
• [27]   T. Nagata, D. G. Fedorov, K. Ishimura, K. Kitaura, J. Chem. Phys., 135, 044110 (2011). , doi:10.1063/1.3611020 PMID:21806093
• [28]   C. M. Aikens, S. P. Webb, R. L. Bell, G. D. Fletcher, M. W. Schmidt, M. S. Gordon, Theor. Chem. Acc., 110, 233 (2003). doi:10.1007/s00214-003-0453-3
• [29]   Y. Mochizuki, S. Koikegami, T. Nakano, S. Amari, K. Kitaura, Chem. Phys. Lett., 396, 473 (2004). doi:10.1016/j.cplett.2004.08.082
• [30]   Y. Mochizuki, T. Nakano, S. Koikegami, S. Tanimori, Y. Abe, U. Nagashima, et al., Theor. Chem. Acc., 112, 442 (2004). doi:10.1007/s00214-004-0602-3
• [31]   D. G. Fedorov, K. Kitaura, J. Chem. Phys., 121, 2483 (2004). , doi:10.1063/1.1769362 PMID:15281845
• [32]   Y. Mochizuki, T. Nakano, Y. Komeiji, K. Yamashita, Y. Okiyama, H. Yoshikawa, et al., Chem. Phys. Lett., 504, 95 (2011). doi:10.1016/j.cplett.2011.01.039
• [33]   J. W. Neidigh, R. M. Fesinmeyer, N. H. Andersen, Nat. Struct. Biol., 9, 425 (2002). , doi:10.1038/nsb798 PMID:11979279
• [34]   G. J. Martyna, M. L. Klein, M. Tuckerman, J. Chem. Phys., 97, 2635 (1992). doi:10.1063/1.463940
• [35]   W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, Numerical Recipes in Fortran 77, 2nd ed., Cambridge Univ. Press, Cambridge, 1992.
• [36]   Y. Komeiji, T. Fujiwara, Y. Okiyama, Y. Mochizuki, Chem-Bio Inform. J., 13, 45 (2013). doi:10.1273/cbij.13.45
• [37]   S. Pfalzner, P. Gibbon, Many-body tree methods in physics, Cambridge Univ. Press, Cambridge, 1996.
• [38]   M. P. Allen, D. J. Tildesley, Computer simulation of liquids, 2nd ed., Oxford Univ. Press, Oxford, 2017.
• [39]   K. R. Brorsen, N. Minezawa, F. Xu, T. L. Windus, M. S. Gordon, J. Chem. Theory Comput., 8, 5008 (2012). , doi:10.1021/ct3007869 PMID:26583213
• [40]   Y. Komeiji, Y. Mochizuki, T. Nakano, Chem. Phys. Lett., 484, 380 (2010). doi:10.1016/j.cplett.2009.11.045
• [41]   BioStation Viewer, https://fmodd.jp/biostationviewer-dl