**2.2 FMO**

4 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Fig. 1. Schematics of the FMO-MD method exemplified by an ion solvation with four water molecules. The atomic nuclei are represented by black circles (the large one for the ion, medium ones for Oxygens, and small ones for Hydrogens) and the electron cloud by a grey shadow. The electronic structure is calculated by FMO to give force (**F**) and energy (*E*), which are then used to update the 3D coordinates of nuclei (**r**) by MD, i.e., by solving the

performed by the PEACH/ABINIT-MP software system composed of the PEACH MD program (Komeiji *et al*., 1997) and the ABINIT-MP 1 (F)MO program (Nakano *et al*., 2000). We have revised the system several times (Komeiji *et al*., 2004, 2009a), but here we describe the latest system, which has not yet been published. In the latest system, the PEACH program prepares the ABINIT-MP input file containing the list of fragments and 3D atomic coordinates, executes an intermediate shell script to run ABINIT-MP, receives the resultant FMO energy and force, and updates the coordinates by the velocity-Verlet integration algorithm. This procedure is repeated for a given number of

The above implementation of FMO-MD, referred to as the PEACH/ABINIT-MP system, has both advantages and disadvantages. The most important advantage is the convenience for the software developers; both FMO and MD programmers can modify their programs independently from each other. Also, if one wants to add a new function of MD, one can first write and debug the MD program against an inexpensive classical force field simulation and then transfer the function to FMO-MD, a costly *ab initio* MD. Nonetheless, the PEACH/ABINIT-MP system has several practical disadvantages as well, mostly related to the use of the systemcall command to connect the two programs. For example, frequent invoking of ABINIT-MP from PEACH sometimes causes a system error that leads to an abrupt end of simulations. Furthermore, use of the systemcall command is prohibited in many supercomputing facilities. To overcome these disadvantages, we are currently

1 Our developers' version of ABINIT-MP is named ABINIT-MPX, but it is referred to as ABINIT-MP

classical equation of motion.

time steps.

throughout this article.

FMO, the essential constituent of FMO-MD, is an approximate *ab initio* MO method (Kitaura *et al.*, 1999). FMO scales to *N*1-2, is easy to parallelize, and retains chemical accuracy during these processes. A vast number of papers have been published on the FMO methodology, but here we review mainly those closely related to FMO-MD. To be more specific, those on the FMO energy gradient, Energy Minimization (EM, or geometry optimization), and MD are preferentially selected in the reference list. Thus, those readers interested in FMO itself are referred to Fedorov & Kitaura (2007b, 2009) for comprehensive reviews of FMO. Also, one can find an extensive review of fragment methods in Gordon *et al.* (2011), where FMO is re-evaluated in the context of its place in the history of the general fragment methods.

### **2.2.1 Hartree-Fock (HF)**

We describe the formulation and algorithm for the HF level calculation with 2-body expansion (FMO2), the very fundamental of the FMO methodology (Kitaura *et al.*, 1999). Below, subscripts *I, J, K*... denote fragments, while *i*, *j*, *k*,... denote atomic nuclei.

First, the molecular system of interest is divided into *Nf* fragments. Second, the initial electron density, *ρI*(*r*), is estimated with a lower-level MO method, e.g., extended Hückel, for all the fragments. Third, self-consistent field (SCF) energy, *EI,* is calculated for each fragment monomer while considering the electrostatic environment. The SCF calculation is repeated until all *ρI*(*r*)'s are mutually converged. This procedure is called the self-consistent charge (SCC) loop. At the end of the SCC loop, monomer electron density *ρI*(**r**) and energy *EI* are obtained. Finally, an SCF calculation is performed once for each fragment pair to obtain dimer electron density *ρIJ*(**r**) and energy *EIJ*. Total electron density *ρ*(**r**) and energy *E* are calculated using the following formulae:

$$\rho\left(\mathbf{r}\right) = \sum\_{l>l} \rho\_{l|l}\left(\mathbf{r}\right) - \left(\mathbf{N}\_f - 2\right) \sum\_{l} \rho\_l\left(\mathbf{r}\right) \tag{1}$$

$$E = \sum\_{l>l} E\_{l\parallel} - (\mathcal{N}\_f - \mathcal{2}) \sum\_l E\_l \,. \tag{2}$$

In calculation of the dimer terms, electrostatic interactions between distant pairs are approximiated by simple Coulombic interactions (dimer-ES approximation, Nakano *et al.*, 2002). This approximation is mandatory to reduce the computation cost from *O*(*N*4) to *O*(*N*2).

The total energy of the molecular system, *U,* is obtained by adding the electrostatic interaction energy between nuclei to *E*, namely,

$$\text{CLI} = \sum\_{l>j} E\_{lj} - (\text{N}\_f - 2) \sum\_l E\_l + \sum\_{i>j} \frac{Z\_i Z\_j}{r\_{ij}} \tag{3}$$

where *rij* denotes the distance between nuclei *i* and *j* and *Zi* and *Zi* their charges, respectively.

Force (**F***i*) acting on atomic nucleus *i* can be obtained by differenciation of eq. (3) by **r**i as follows:

$$\mathbf{F}\_i = -\nabla\_i \mathcal{U} \tag{4}$$

Recent Advances in Fragment Molecular

energy incorporated by MP2 to describe a condensed phase.

**2.2.4 Configuration Interaction Singles (CIS)**

**2.2.5 Unrestricted Hartree-Fock (UHF)** 

**2.2.6 Model Core Potential (MCP)** 

been developed and made available (Fujiwara *et al*., 2011).

**2.2.7 Periodic Boundary Condition (PBC)** 

Orbital-Based Molecular Dynamics (FMO-MD) Simulations 7

value than that of HF/FMO-MD was. This result indicated the importance of the correlation

CIS is a useful tool to model low-lying excited states caused by transitions among near HOMO-LUMO levels in a semi-quantitative fashion (Foresman *et al.*, 1992). A tendency of CIS to overestimate excitation energies is compensated for by CIS(D) in which the orbital relaxation energy for an excited state of interest as well as the differential correlation energy from the ground state correlated at the MP2 level (Head-Gordon *et al.*, 1994). Both CIS and CIS(D) have been introduced to multilayer FMO (MFMO; Fedorov *et al.*, 2005) in ABINIT-MP (Mochizuki *et al.*, 2005a, 2007a). Very recently, Mochizuki implemented the parallelized FMO3-CIS gradient calculation, based on the efficient formulations with Fock-like contractions (Foresman *et al*., 1992). The dynamics of excited states is now traceable as long as the CIS approximation is qualitatively correct enough. The influence of hydration on the excited state induced proton-transfer (ESIPT) has been attracting considerable interest, and

UHF is the simplest method for handling open-shell molecular systems, as long as care for the associated spin contamination is taken. The UHF gradient was implemented by preparing - and β-density matrices. Simulation of hydrated Cu(II) has been underway at the FMO3-UHF level, and the Jahn-Teller distortion of hexa-hydration has been reasonably reproduced (Kato *et al.*, in preparation). The extension to a UMP2 gradient is planned as a future subject, where the computational cost may triple the MP2 gradient because of the three types of transformed integrals, (,), (,), and (,) (Aikens *et al.*, 2003).

Heavy metal ions play major roles in various biological systems and functional materials. Therefore, it is important to understand the fundamental chemical nature and dynamics of the metal ions under physiological or experimental conditions. Each heavy metal element has a large number of electrons to which relativistic effects must be taken into account, however. Hence, the heavy metal ions increase the computation cost of high-level electronic structure theories. A way to reduce the computation is the Model Core Potential (MCP; Sakai *et al.*, 1987; Miyoshi *et al.*, 2005; Osanai *et al.*, 2008ab; Mori *et al.*, 2009), where the proper nodal structures of valence shell orbitals can be maintained by the projection operator technique. In the MCP scheme, only valence electrons are considered, and core electrons are replaced with 1-electron relativistic pseudo-potentials to decrease computational costs. The MCP method has been combined with FMO and implemented in ABINIT-MP (Ishikawa *et al*., 2006), which has been used in the comparative MCP/FMO-MD simulations of hydrated *cis*-platin and *trans*-platin (see subsection 3.6). Very recently, the 4f-in-core type MCP set for trilvalent lanthanides has

PBC was finally introduced to FMO-MD in the TINKER/ABINIT-MP system by Fujita *et al.* (2011). PBC is a standard protocol for both classical and *ab intio* MD simulations.

we have started related simulations for several pet systems such as toropolone.

Analytical formulation of eq. (4) was originally derived for the HF level by Kitaura *et al.* (2001) and used in several EM calculations (for example, Fedorov *et al.*, 2007a) and in the first FMO-MD simulation (Komeiji *et al.,* 2003). Later on, the HF gradient was made fully analytic by Nagata *et al.* (2009, 2010, 2011a).

#### **2.2.2 FMO***n*

The procedure described in the previous subsection is called FMO2, with "2" indicating that the energy is expanded up to 2-body terms of fragments. It is possible to improve the precision of FMO by adding 3-body, 4-body, ..., and *n*-body terms (FMO*n*) at the expense of the computation cost of *O*(1). FMO3 has been implemented in both GAMESS and ABINIT-MP. The improvement by FMO3 is especially apparent in FMO-MD, as exemplified by a simulation of proton transfer in water (Komeiji *et al.*, 2010). Recently, FMO4 was implemented in ABINIT-MP (Nakano *et al.,* 2012), which will presumably make it possible to regard even a metal ion as a fragment.

#### **2.2.3 Second-order Moeller-Plesset perturbation (MP2)**

The HF calculation neglects the electron correlation effect, which is necessary to incorporate the so-called dispersion term. The electron correlation can be calculated fairly easily by the second-order Moeller-Plesset perturbation (MP2). Though the MP2/FMO energy formula was published as early as 2004 (Fedorov *et al.,* 2004b; Mochizuki *et al.*, 2004ab), the energy gradient formula for MP2/FMO was first published in 2011 by Mochizuki *et al.* (2011) and then by Nagata *et al.* (2011). In Mochizuki's implementation of MP2 to ABINIT-MP, an integral-direct MP2 gradient program module with distributed parallelism was developed for both FMO2 and FMO3 levels, and a new option called "FMO(3)" was added, in which FMO3 is applied to HF but FMO2 is applied to MP2 to reduce computation time, based on the relatively short-range nature of the electron correlation compared to the range of the Coulomb or electrostatic interactions.

The MP2/FMO gradient was soon applied to FMO-MD of a droplet of water molecules (Mochizuki *et al.*, 2011). The water was simulated with the 6-31G\* basis set with and without MP2, and the resultant trajectories were subjected to calculations of radial distribution functions (RDF). The RDF peak position of MP2/FMO-MD was closer to the experimental

The total energy of the molecular system, *U,* is obtained by adding the electrostatic

*IJ f I I J I ij ij*

*r*

where *rij* denotes the distance between nuclei *i* and *j* and *Zi* and *Zi* their charges,

Force (**F***i*) acting on atomic nucleus *i* can be obtained by differenciation of eq. (3) by **r**i as

Analytical formulation of eq. (4) was originally derived for the HF level by Kitaura *et al.* (2001) and used in several EM calculations (for example, Fedorov *et al.*, 2007a) and in the first FMO-MD simulation (Komeiji *et al.,* 2003). Later on, the HF gradient was made fully

The procedure described in the previous subsection is called FMO2, with "2" indicating that the energy is expanded up to 2-body terms of fragments. It is possible to improve the precision of FMO by adding 3-body, 4-body, ..., and *n*-body terms (FMO*n*) at the expense of the computation cost of *O*(1). FMO3 has been implemented in both GAMESS and ABINIT-MP. The improvement by FMO3 is especially apparent in FMO-MD, as exemplified by a simulation of proton transfer in water (Komeiji *et al.*, 2010). Recently, FMO4 was implemented in ABINIT-MP (Nakano *et al.,* 2012), which will presumably make it possible

The HF calculation neglects the electron correlation effect, which is necessary to incorporate the so-called dispersion term. The electron correlation can be calculated fairly easily by the second-order Moeller-Plesset perturbation (MP2). Though the MP2/FMO energy formula was published as early as 2004 (Fedorov *et al.,* 2004b; Mochizuki *et al.*, 2004ab), the energy gradient formula for MP2/FMO was first published in 2011 by Mochizuki *et al.* (2011) and then by Nagata *et al.* (2011). In Mochizuki's implementation of MP2 to ABINIT-MP, an integral-direct MP2 gradient program module with distributed parallelism was developed for both FMO2 and FMO3 levels, and a new option called "FMO(3)" was added, in which FMO3 is applied to HF but FMO2 is applied to MP2 to reduce computation time, based on the relatively short-range nature of the electron correlation compared to the range of the

The MP2/FMO gradient was soon applied to FMO-MD of a droplet of water molecules (Mochizuki *et al.*, 2011). The water was simulated with the 6-31G\* basis set with and without MP2, and the resultant trajectories were subjected to calculations of radial distribution functions (RDF). The RDF peak position of MP2/FMO-MD was closer to the experimental

*U EN E*

( 2) *<sup>i</sup> <sup>j</sup>*

*Z Z*

(3)

**F***i i U* (4)

interaction energy between nuclei to *E*, namely,

analytic by Nagata *et al.* (2009, 2010, 2011a).

to regard even a metal ion as a fragment.

Coulomb or electrostatic interactions.

**2.2.3 Second-order Moeller-Plesset perturbation (MP2)** 

respectively.

follows:

**2.2.2 FMO***n*

value than that of HF/FMO-MD was. This result indicated the importance of the correlation energy incorporated by MP2 to describe a condensed phase.
