**3. Applications of FMO-MD**

FMO-MD has been extensively applied to hydrated small molecules to simulate their solvation and chemical reactions. Some benchmark FMO-MD simulations were described briefly in the previous section. In this section, we review genuine applications of FMO-MD in detail.

### **3.1 Excitation energy of hydrated formaldehyde**

FMO-MD and MFMO-CIS(D) were combined to evaluate the lowest n\* excitation energy of hydrated formaldehyde (H2CO) molecules (Mochizuki *et al.*, 2007b). The shift of excitation energy of a solute by the presence of a solvent, known as solvatochronism, has drawn attention of both experimentalists and theorists and has been studied by various computational methods, mostly by the quantum mechanics and molecular mechanics (QM/MM) method. Alternatively, Mochizuki *et al.* (2007b) tried a fully *ab initio* approach, in which FMO-MD sampled molecular configurations for excited calculations.

The DF algorithm gracefully handles molecular systems consisting of small solute and solvent molecules, but not those containing large molecules such as proteins and DNA, which should be fragmented at covalent bonds. Currently, Mode 0 is the only choice of fragmentation for these large molecules, in which the initial fragmentation should be used throughout and no fragment rearrangement is allowed (Nakano *et al.*, 2000; Komeiji *et al.*, 2004). This limitation of the DF algorithm will be abolished soon by the introduction of a

The blue moon ensemble method (Sprik & Ciccotti, 1998) is a way to calculate the free energy profile along a reaction coordinate (RC) while constraining RC to a specified value. The method was implemented in FMO-MD (Komeiji, 2007) and was successfully applied to

The nuclei were handled by the classical mechanics in most of the FMO-MD simulations performed to date (Fig. 1), but PIMD (Marx & Parrinello, 1996) has been introduced into FMO-MD to incorporate the nucleic quantum effect (Fujita *et al.*, 2009). FMO-PIMD consumes tens of times more computational resource than the classical FMO-MD does but is

Miscellaneous MD methods implemented in the PEACH/ABINIT-MP system include the Nosé-Hoover (chain) thermostat, RATTLE bond constraint, RC constraint, spherical solvent boundary, and so on (Komeiji *et al.*, 2009a). Another research group has implemented the Hamiltonian Algorithm (HA) to FMO-MD to enhance conformation sampling of, for

FMO-MD has been extensively applied to hydrated small molecules to simulate their solvation and chemical reactions. Some benchmark FMO-MD simulations were described briefly in the previous section. In this section, we review genuine applications of FMO-MD

FMO-MD and MFMO-CIS(D) were combined to evaluate the lowest n\* excitation energy of hydrated formaldehyde (H2CO) molecules (Mochizuki *et al.*, 2007b). The shift of excitation energy of a solute by the presence of a solvent, known as solvatochronism, has drawn attention of both experimentalists and theorists and has been studied by various computational methods, mostly by the quantum mechanics and molecular mechanics (QM/MM) method. Alternatively, Mochizuki *et al.* (2007b) tried a fully *ab initio* approach, in which FMO-MD sampled molecular configurations for excited

drawing a free energy profile of the Menschutkin reaction (Komeiji *et al.*, 2009a).

necessary for a better description of, for example, a proton transfer reaction.

example, polypeptides (Ishimoto *et al.*, 2004, 2005; Tamura *et al.*, 2008).

mixed algorithm of DF and a static fragmentation.

**2.3.3 Path Integral Molecular Dynamics (PIMD)** 

**2.3.2 Blue moon ensemble**

**2.3.4 Miscellaneous** 

in detail.

calculations.

**3. Applications of FMO-MD** 

**3.1 Excitation energy of hydrated formaldehyde** 

Fig. 3. An FMO-MD snapshot of the solvated H2CO (left). Histogram of excitation energies for CIS and CIS(D) calculations (right). Reproduced from Mochizuki *et al.* (2007b) with permission.

In the configuration sampling, H2CO was solvated within a droplet of 128 water molecules (Fig. 3 left), and the molecular system was simulated by FMO-MD at the FMO2-HF/6-31G level to generate a 2.62-ps trajectory at 300 K. From the last 2-ps portion of the trajectory, 400 conformations were chosen and were subjected to MFMO-CIS(D) calculations at the FMO2/HF/6-31G\* level. In MFMO, the chromophore region contained H2CO and several water molecules and was the target of CIS(D) calculation. The calculated excitation energy was averaged over the 400 configurations (Fig. 3 right). A similar protocol was applied to an isolated H2CO molecule to calculate the excitation energy in a vacuum. The blue-shift by solvatochromism thus estimated was 0.14 eV, in agreement with preceding calculations.

The solvatochromism of H2CO is frequently challenged by various computational methods, but this study distinguishes itself from preceding studies in that all the calculations were fully quantum, without classical force field parameters.

### **3.2 Hydrolysis of a methyl diazonium ion**

The hydrolysis of the methyl-diazonium ion (CH3N2+) is an SN2-type substitution reaction that proceeds as follows:

$$\text{H}\_{2}\text{O} + \text{CH}\_{3}\text{N}\_{2}\text{\*} \rightarrow [\text{H}\_{2}\text{O}...\text{CH}\_{3}\text{\*}...\text{N}\_{2}]\text{\*} \rightarrow \text{\*} + \text{H}\_{2}\text{OCH}\_{3}\text{\*} + \text{N}\_{2}\text{!} \tag{5}$$

Traditionally, this reaction is believed to occur in an enforced concerted mechanism in which a productive methyl cation after N2 leaving is too reactive to have a finite lifetime, and consequently the attack by H2O and the bond cleavage occur simultaneously. This traditional view was challenged by Sato *et al.* (2008) using FMO-MD. The FMO-MD simulations exhibited diverse paths, showing that the chemical reaction does not always proceed through the lowest energy paths.

This reaction was simulated as follows. FMO-MD simulations were conducted at the FMO2/HF/6-31G level. CH3-N2+ was optimized in the gas phase and then hydrated in a sphere of 156 water molecules. The water was optimized at 300 K for 0.5 ps with the RATTLE bond constraint. The temperature of the molecular system was raised to 1000 K,

Recent Advances in Fragment Molecular

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

Fig. 5. Charge transfer interaction energy between attacking H2O and CH3N2+ as functions of *R*O-C (left) and *R*C-N (right). The open circles show trajectory A, and the filled triangles

Fig. 6. *R*C-N-*R*O-C plot of the ten trajectories that resulted in product formation. Those

Reproduced from Sato *et al.* (2008) by permission.

stepwise mechanism with ZW as a stable intermediate.

**3.3 Amination of formaldehyde** 

trajectories that proceeded along the diagonal line are regarded as tight SN2, in which attack by water and the exit of N2 occurred simultaneously, while a trajectory that deviated from the diagonal line is regarded as loose SN2, in which N2 left before the attack by water.

Sato *et al.* (2010) tackled the reaction mechanism of the amination of H2CO by FMO-MD simulations. In particular, they focused on whether the reaction proceeds via a zwitterion (ZW) intermediate (Fig. 7). The results indicated that the reaction proceeds through a

show trajectory B. Reproduced from Sato *et al*. (2008) by permission.

and the simulation was continued for 5 ps. From the 1000 K trajectory, 15 configurations were taken and subjected to a further run at 700 K without any constraint. Ten trajectories out of fifteen produced the final products (CH3-OH2 ++N2). The ten productive trajectories were classified into three groups: tight SN2, loose SN2, and intermediate.

Fig. 4. Initial droplet structure and structures of substrate and nearby water molecules along type A and B trajectories. Numbers are atomic distances in Å. Reproduced from Sato *et al.* (2008) by permission.

Trajectory A in Fig. 4 is of the tight SN2 type, in which the attack by H2O and C-N bond cleavage, i.e. release of N2, occur concertedly. Trajectory B is of the loose SN2 type, which shows a two-stage process in which C-N bond cleavage precedes the attack by H2O.

The difference between trajectories A and B was further analyzed by the configuration analysis for fragment interaction (CAFI; Mochizuki *et al.*, 2005b), and the results are plotted in Fig. 5. Charge-transfer (CT) interaction between the two fragments increases rapidly when the C-N distance increases to 1.6 Å for trajectory A, but for trajectory B the CT increased only when *R*CN was 2.4 Å or longer. In trajectory B, the C-N bond cleavage and O-C bond formation events take place in a two-stage fashion. The CT interaction energy is larger for trajectory B than for A at *R*C-O = 2.6 Å, because at the same C-O distance the C-N bond is cleaved to a larger extent, and hence the CH3 moiety has more positive charge for trajectory B than for trajectory A.

Most of the other productive trajectories exhibited intermediate characteristics between those of trajectories A and B. The diversity of the reaction path can be illustrated by the twodimensional *R*C-N-*R*O-C plot (Fig. 6). The existence of different paths indicates that the reaction does not always proceed through the lowest energy pathway with optimal solvation.

In summary, this series of simulations illustrated for the first time how the atoms in reacting molecules, from reactant to product, behave in solution at the molecular level. This was made possible by the advent of the full *ab initio* FMO-MD method.

and the simulation was continued for 5 ps. From the 1000 K trajectory, 15 configurations were taken and subjected to a further run at 700 K without any constraint. Ten trajectories

Fig. 4. Initial droplet structure and structures of substrate and nearby water molecules along type A and B trajectories. Numbers are atomic distances in Å. Reproduced from Sato *et al.*

Trajectory A in Fig. 4 is of the tight SN2 type, in which the attack by H2O and C-N bond cleavage, i.e. release of N2, occur concertedly. Trajectory B is of the loose SN2 type, which shows a two-stage process in which C-N bond cleavage precedes the attack

The difference between trajectories A and B was further analyzed by the configuration analysis for fragment interaction (CAFI; Mochizuki *et al.*, 2005b), and the results are plotted in Fig. 5. Charge-transfer (CT) interaction between the two fragments increases rapidly when the C-N distance increases to 1.6 Å for trajectory A, but for trajectory B the CT increased only when *R*CN was 2.4 Å or longer. In trajectory B, the C-N bond cleavage and O-C bond formation events take place in a two-stage fashion. The CT interaction energy is larger for trajectory B than for A at *R*C-O = 2.6 Å, because at the same C-O distance the C-N bond is cleaved to a larger extent, and hence the CH3 moiety has more positive charge for

Most of the other productive trajectories exhibited intermediate characteristics between those of trajectories A and B. The diversity of the reaction path can be illustrated by the twodimensional *R*C-N-*R*O-C plot (Fig. 6). The existence of different paths indicates that the reaction does not always proceed through the lowest energy pathway with optimal

In summary, this series of simulations illustrated for the first time how the atoms in reacting molecules, from reactant to product, behave in solution at the molecular level. This was

made possible by the advent of the full *ab initio* FMO-MD method.

++N2). The ten productive trajectories

out of fifteen produced the final products (CH3-OH2

(2008) by permission.

trajectory B than for trajectory A.

by H2O.

solvation.

were classified into three groups: tight SN2, loose SN2, and intermediate.

Fig. 5. Charge transfer interaction energy between attacking H2O and CH3N2 + as functions of *R*O-C (left) and *R*C-N (right). The open circles show trajectory A, and the filled triangles show trajectory B. Reproduced from Sato *et al*. (2008) by permission.

Fig. 6. *R*C-N-*R*O-C plot of the ten trajectories that resulted in product formation. Those trajectories that proceeded along the diagonal line are regarded as tight SN2, in which attack by water and the exit of N2 occurred simultaneously, while a trajectory that deviated from the diagonal line is regarded as loose SN2, in which N2 left before the attack by water. Reproduced from Sato *et al.* (2008) by permission.

#### **3.3 Amination of formaldehyde**

Sato *et al.* (2010) tackled the reaction mechanism of the amination of H2CO by FMO-MD simulations. In particular, they focused on whether the reaction proceeds via a zwitterion (ZW) intermediate (Fig. 7). The results indicated that the reaction proceeds through a stepwise mechanism with ZW as a stable intermediate.

Recent Advances in Fragment Molecular

was equilibrated for 0.3 ps and sampled for a further 0.3 ps.

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

for RC = 0.9, 1.0,..., 1.8 Å starting from configuration B. For each RC value the configuration

The diagram thus obtained clearly favored the stepwise mechanism over the concerted mechanism (Fig. 9). Nevertheless, there remained a possibility of the MD trajectory being trapped in a local minimum. To investigate the possibility, we conducted additional FMO-MD simulations starting from configuration C, the concerted TS-like one. These additional trajectories all diverted from the TS-like structure toward the trajectory of the stepwise path (see Sato *et al.*, 2010, for details), thus confirming the validity of the stepwise mechanism.

Fig. 9. Reaction profile obtained by FMO-MD simulations (left). The concerted TS-like

In summary, the constraint FMO-MD simulations indicated that the H2CO amination in

The divalent zinc ion, Zn(II), plays bio-chemically relevant roles, e.g., as the reaction center of superoxide dismutase. By using a droplet model of the Zn(II) ion with 64 water molecules, FMO2- and FMO3-MD simulations were performed at the HF/6-31G level, supposing that the electrostatic and coordination interactions are dominant in this system (Fujiwara *et al.*, 2010b). The Zn-O peak positions at the first hydration shell were investigated, and a better accuracy of FMO3-MD than that of FMO2-MD was demonstrated, where the FMO3 value of 2.05 Å agreed well with the experimental value of 2.06±0.02 Å (Fig. 10). The coordination number of the first hydration shell was 6 consistently. Additionally, the charge fluctuations on the Zn atom were evaluated by the natural population analysis (NPA) as well as the conventional Mulliken population analysis (MPA). The NPA results showed a consistent picture with the coordination bond with reasonable fluctuation (around a net charge of 1.8), while MPA yielded an artificially enhanced

structure (right). Reproduced from Sato *et al.* (2010) by permission.

**3.4 Hydration of Zn(II)** 

water solvent occurs by the stepwise mechanism, not by the concerted one.

Fig. 7. Two contradictory schemes of H2CO amination. RT: reactant; ZW: zwitterion; PD: product.

The FMO-MD simulations were designed as follows. RC was defined as *R*N-C-*R*N-H. With RC constrained, structural changes of the reactant (RT) molecules in MD simulations are confined to the line that has the slope=1 and intercept=RC in a More O'Ferrall–Jencks-type diagram (Fig. 8). This diagram allows the reader to distinguish between the stepwise process and the concerted one.

Fig. 8. Schematic representation of the More O'Farrall-Jencks-type diagram of carbinolamine formation of formaldehyde and ammonia (left). Three optimized initial configurations (right). Reproduced from Sato *et al.* (2010) by permission.

By FMO-MD, a More O'Ferrall–Jencks-type diagram was drawn for the H2CO amination. Three initial configurations were prepared, (A) zwitterion-like, (B) reactant-like, and (C) concerted TS-like (Fig. 8), each solvated with ca. 200 water molecules. After appropriate optimization and equilibration by classical and FMO-EM/MD methods, average *R*NH and *R*NC were calculated at 300 K for RC = -0.4, -0.3,..., 0.9 Å starting from configuration A and

Fig. 7. Two contradictory schemes of H2CO amination. RT: reactant; ZW: zwitterion; PD:

The FMO-MD simulations were designed as follows. RC was defined as *R*N-C-*R*N-H. With RC constrained, structural changes of the reactant (RT) molecules in MD simulations are confined to the line that has the slope=1 and intercept=RC in a More O'Ferrall–Jencks-type diagram (Fig. 8). This diagram allows the reader to distinguish between the stepwise

Fig. 8. Schematic representation of the More O'Farrall-Jencks-type diagram of carbinolamine formation of formaldehyde and ammonia (left). Three optimized initial configurations

By FMO-MD, a More O'Ferrall–Jencks-type diagram was drawn for the H2CO amination. Three initial configurations were prepared, (A) zwitterion-like, (B) reactant-like, and (C) concerted TS-like (Fig. 8), each solvated with ca. 200 water molecules. After appropriate optimization and equilibration by classical and FMO-EM/MD methods, average *R*NH and *R*NC were calculated at 300 K for RC = -0.4, -0.3,..., 0.9 Å starting from configuration A and

(right). Reproduced from Sato *et al.* (2010) by permission.

product.

process and the concerted one.

for RC = 0.9, 1.0,..., 1.8 Å starting from configuration B. For each RC value the configuration was equilibrated for 0.3 ps and sampled for a further 0.3 ps.

The diagram thus obtained clearly favored the stepwise mechanism over the concerted mechanism (Fig. 9). Nevertheless, there remained a possibility of the MD trajectory being trapped in a local minimum. To investigate the possibility, we conducted additional FMO-MD simulations starting from configuration C, the concerted TS-like one. These additional trajectories all diverted from the TS-like structure toward the trajectory of the stepwise path (see Sato *et al.*, 2010, for details), thus confirming the validity of the stepwise mechanism.

Fig. 9. Reaction profile obtained by FMO-MD simulations (left). The concerted TS-like structure (right). Reproduced from Sato *et al.* (2010) by permission.

In summary, the constraint FMO-MD simulations indicated that the H2CO amination in water solvent occurs by the stepwise mechanism, not by the concerted one.

#### **3.4 Hydration of Zn(II)**

The divalent zinc ion, Zn(II), plays bio-chemically relevant roles, e.g., as the reaction center of superoxide dismutase. By using a droplet model of the Zn(II) ion with 64 water molecules, FMO2- and FMO3-MD simulations were performed at the HF/6-31G level, supposing that the electrostatic and coordination interactions are dominant in this system (Fujiwara *et al.*, 2010b). The Zn-O peak positions at the first hydration shell were investigated, and a better accuracy of FMO3-MD than that of FMO2-MD was demonstrated, where the FMO3 value of 2.05 Å agreed well with the experimental value of 2.06±0.02 Å (Fig. 10). The coordination number of the first hydration shell was 6 consistently. Additionally, the charge fluctuations on the Zn atom were evaluated by the natural population analysis (NPA) as well as the conventional Mulliken population analysis (MPA). The NPA results showed a consistent picture with the coordination bond with reasonable fluctuation (around a net charge of 1.8), while MPA yielded an artificially enhanced

Recent Advances in Fragment Molecular

application between *cis*- and *trans*-platins.

Reproduced from Mori *et al.* (2012) by permission.

solvent water molecule.

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

Their hydration should be investigated to understand the difference in the medical

Fig. 11. Structures of *cis*- and *trans*-platins and schematic representations of DNA adducts.

FMO-MD simulations were performed for hydrated *cis*- and *trans*-platins. The simulation conditions were set as described below. Each platin complex was hydrated with a spherical droplet of water centred at the Pt atom with a diameter of 10.5 Å. This diameter was determined to include up to the second solvation shell, so that the physicochemical properties of the first shell should be reproduced. In the FMO-MD simulations, the electronic states of the hydrated platin complexes were described by FMO(3)-MP2. The basis sets were MCPdz for Pt, MCPdzp for Cl, and 6-31G(d) for the others, respectively. The MCP basis sets were applied for heavy elements (see subsection 2.2.6). The central platin and each of the water molecules were regarded as independent fragments. DF was applied to allow for the generation of proton-transferred species during the production MD runs. For each *cis*- and *trans*-platin system, a 1-ps equilibration and a subsequent 2-ps production MD run were performed using the Nose-Hoover Chains NVT ensemble at 300 K. NPA was also performed during the FMO-MD run to analyze the differences in charge fluctuations between *cis*- and *trans*-platin, illuminating the differences in the hydration environment around polarized Pt+-Cl- bonds, which should be cleaved by the nucleophilic attack of a

The time evolution of the natural charge on each ligand in *cis*- and *trans*-platin, and that of Pt-Cl bond lengths are shown in Fig. 12. Relatively larger charge fluctuations were observed on the Pt/Cl sites than on the NH3 sites in both platins. This difference among the sites was attributable to the fact that NH3 has no amplitude in the highest occupied molecular orbital. A close comparison of the left and right graphs in Fig. 12 revealed a correlation between fluctuation of the Pt/Cl sites and that of the Pt-Cl bond. By applying the Fourier transform technique to the charge fluctuation, we calculated the frequency of the fluctuation to be 334 cm-1. This frequency can also be assigned to the Pt-Cl stretching mode coupled with intermolecular vibrations between the solute platin and solvent water molecules. The correlation observed in charge fluctuation on Pt and Cl sites means that there is a CT interaction between them. Since the frontier MO that participates in the CT process is a Pt-Cl antibonding orbital, the CT interaction coupled with the fluctuation of the solvent water should induce a Pt-Cl bonds fluctuation. Since *trans*-platin has inversion symmetry, the

fluctuation with a larger extent of electron donation (net charge of 1.3-1.4). Discussion with NPA was found to be preferable for hydrated metal ions.

Fig. 10. Zn-O RDFs and coordination numbers (CN) calculated by FMO2/3-MD simulations. Reproduced from Fujiwara *et al.* (2010b) by permission.
