**2.3 MD**

The MD portion of FMO-MD resembles the conventional classical MD method, but several algorithms have been introduced to facilitate FMO-MD.

#### **2.3.1 Dynamic Fragmentation (DF)**

DF refers to the redefinition of fragments depending on the molecular configuration during FMO-MD. For example, in an H+-transfer reaction (AH+ + B → AHB+ → A + BH+), AH+ and B can be separate fragments before the reaction but should be unified in the transition state AHB+, and A and BH+ may be separated after the reaction. The DF algorithm handles this fragment rearrangement by observing the relative position and nuclear species of the constituent atoms at each time step of a simulation run.

The need for DF arose for the first time in an FMO-MD simulation of solvated H2CO (Mochizuki *et al.* 2007b; see subsection 3.1). During the equilibration stage of the simulation, an artifactual H+-transport frequently brought about an abrupt halt of the simulation. To avoid the halt by the H+-transport, T. Ishikawa developed a program to unite the donor and acceptor of H+ by looking up the spatial formation of the water molecules. This program was executed at each time step of the simulation. This was the first implementation of the DF algorithm (see Komeiji *et al.*, 2009a, for details). A similar *ad hoc* DF program was written for a simulation of hydrolysis methyl-diazonium (Sato *et al.*, 2008; see subsection 3.2). Thus, at the original stage, different DF programs were needed for different molecular systems.

Nonetheless, partly due to the complexity of PBC in formulation but mostly due to its computation cost, FMO-MD simulations reported in the literature had been performed under a free boundary condition, usually with a cluster solvent model restrained by a harmonic spherical potential. This spherical boundary has the disadvantage of exposing the simulated molecular system to a vacuum condition and altering the electronic structure of the outer surface (Komeiji *et al.*, 2007). Hence, PBC is expected to avoid the disadvantage and to extend FMO-MD to simulations of bulk solvent and crystals. For PBC simulations to be practical, efficient approximations in evaluating the ESP matrix elements will need to be developed. A technique of multipole expansion may be worth

Analytic gradient formulae have been derived for several FMO methods and implemented in the GAMESS software, including those for the adaptive frozen orbital bond detachment scheme (AFO; Fedorov *et al*., 2009), polarizable continuum model method (PCM; Li *et al.*, 2010), time-dependent density functional theory (TD-DFT; Chiba *et al.*, 2009), MFMO with active, polarisable, and frozen sites (Fedorov *et al.*, 2011), and effective fragment potential (EFP; Nagata *et al.*, 2011c). Also, Ishikawa *et al.* (2010) implemented partial energy gradient (PEG) in their software PACIS. These gradients have been used for FMO-EM calculations of appropriate molecules. Among them, the EFP gradient has already been applied successfully to FMO-MD (Nagata *et al.*, 2011c), and the others will be combined with FMO-

The MD portion of FMO-MD resembles the conventional classical MD method, but several

DF refers to the redefinition of fragments depending on the molecular configuration during FMO-MD. For example, in an H+-transfer reaction (AH+ + B → AHB+ → A + BH+), AH+ and B can be separate fragments before the reaction but should be unified in the transition state AHB+, and A and BH+ may be separated after the reaction. The DF algorithm handles this fragment rearrangement by observing the relative position and nuclear species of the

The need for DF arose for the first time in an FMO-MD simulation of solvated H2CO (Mochizuki *et al.* 2007b; see subsection 3.1). During the equilibration stage of the simulation, an artifactual H+-transport frequently brought about an abrupt halt of the simulation. To avoid the halt by the H+-transport, T. Ishikawa developed a program to unite the donor and acceptor of H+ by looking up the spatial formation of the water molecules. This program was executed at each time step of the simulation. This was the first implementation of the DF algorithm (see Komeiji *et al.*, 2009a, for details). A similar *ad hoc* DF program was written for a simulation of hydrolysis methyl-diazonium (Sato *et al.*, 2008; see subsection 3.2). Thus, at the original stage, different DF programs were

algorithms have been introduced to facilitate FMO-MD.

constituent atoms at each time step of a simulation run.

**2.3.1 Dynamic Fragmentation (DF)** 

needed for different molecular systems.

considering.

**2.2.8 Miscellaneous** 

MD in the near future.

**2.3 MD** 

The DF algorithm was generalized later to handle arbitrary molecular systems (Komeiji *et al.*, 2010). The algorithm requires each atom's van der Waals radius and instantaneous 3D coordinate, atomic composition and net charge of possible fragment species, and certain threshold parameters.

Presently, PEACH has four fragmentation modes, as follows:

Mode 0: Use the fragmentation data in the input file throughout the simulation.

Mode 1: Merge covalently connected atoms, namely, those constituting a molecule, into a fragment.

Mode 2: Fragments produced by Mode 1 are unified into a larger fragment if they are forming an H-bond.

Mode 3: Fragments produced by Mode 2 are unified if they are an ion and coordinating solvent molecules.

The modes are further explained as follows. Heavy atoms located significantly close to each other are united as a fragment, and each H atom is assigned to its closest heavy atom (Mode 1). Then, two fragments sharing an H atom are unified (Mode 2). Finally, an ion and surrounding molecules are united (Mode 3). See Figure 2 for typical examples of DF. Usually, Mode 1 is enough, but Mode 2 or 3 sometimes become necessary.

Fig. 2. Typical examples of fragment species generated by the generalized DF scheme. Expected fragmentation patterns are drawn for three solute molecules, A–C. Reproduced from Komeiji *et al.* (2010) with permission.

Recent Advances in Fragment Molecular

permission.

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

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

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

The hydrolysis of the methyl-diazonium ion (CH3N2+) is an SN2-type substitution reaction

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

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,

+...N2] → + H2OCH3 + N2. (5)

fully quantum, without classical force field parameters.

**3.2 Hydrolysis of a methyl diazonium ion** 

proceed through the lowest energy paths.

H2O + CH3N2+ → [H2O...CH3

that proceeds as follows:

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 mixed algorithm of DF and a static fragmentation.
