Large-Scale Quantum-Mechanical Molecular Dynamics Simulations

Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Letter

Large-Scale Quantum-Mechanical Molecular Dynamics Simulations Using Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method Yoshio Nishimoto, Hiroya Nakata, Dmitri G. Fedorov, and Stephan Irle J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.5b02490 • Publication Date (Web): 01 Dec 2015 Downloaded from http://pubs.acs.org on December 2, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Large-Scale Quantum-Mechanical Molecular Dynamics Simulations Using Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method Yoshio Nishimoto,∗,†,‡ Hiroya Nakata,¶,§ Dmitri G. Fedorov,∗,∥ and Stephan Irle†,⊥ †Department of Chemistry, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan ‡Fukui Institute for Fundamental Chemistry, Kyoto University, 34-4 Takano Nishihiraki-cho, Sakyo-ku, Kyoto 606-8103, Japan ¶Center for Biological Resources and Informatics, Tokyo Institute of Technology, B-62 4259 Nagatsuta-cho, Midoriku, Yokohama 226-8501, Japan §RIKEN, Research Cluster for Innovation, Nakamura Lab, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan ∥Research Center for Computational Design of Advanced Functional Materials (CD-FMat), National Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan ⊥Institute of Transformative Bio-Molecules (WPI-ITbM), Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan E-mail: [email protected]; [email protected] Phone: +81 (0)75 711-7894; +81 (0)29 851-5426

1 ACS Paragon Plus Environment


Abstract The fully analytic gradient is developed for density-functional tight-binding (DFTB) combined with the fragment molecular orbital (FMO) method (FMO-DFTB). The response terms arising from the coupling of the electronic state to the embedding potential are derived, and the gradient accuracy is demonstrated on water clusters and a polypeptide. The radial distribution functions (RDFs) obtained with FMO-DFTB are found to be similar to those from conventional DFTB, while the computational cost is greatly reduced; for 256 water molecules one molecular dynamics (MD) step takes 73.26 and 0.68 s with full DFTB and FMO-DFTB, respectively, showing a speed-up factor of 108. FMO-DFTB/MD is applied to 100 ps MD simulations of liquid hydrogen halides and is found to reproduce experimental RDFs reasonably well.

Graphical TOC Entry For 2000 atoms (HX)1000, X=Halogen Full DFTB 1816 sec.

FMO-DFTB

g(r)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1.8 sec.

r


Page 2 of 20

Page 3 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Molecular dynamics (MD) simulations are very important in chemistry because they make it possible to study time evolution of chemical systems and describe various phenomena such as chemical reactions, diffusion, etc. Many of them occur on a time scale of nanoseconds or longer, and such MD simulations are very time consuming. Molecular mechanics (MM) has been successfully used to perform MD simulations, 1 and there is an ongoing effort to develop new kinds of force fields. 2 Although there are some ways to describe chemical reactions with force fields, 3 a more general approach is to use quantum-mechanical (QM) calculations. Hybrids of QM and MM (QM/MM) 4 have been very successful. Considerable progress has also been achieved with Car-Parinello MD 5 and with ab initio MD on modern computer hardware. 6 One can also use very fast parameterized QM approaches: semi-empirical 7 and density-functional tight-binding (DFTB). 8,9 Although DFTB/MD can be used for small and medium systems (containing hundreds of atoms), 10,11 the computational cost scales with the system size N as N 3 , and simulations become prohibitively expensive for larger systems. To reduce the scaling, fragment-based approaches 12–14 and other techniques 15,16 are efficient, and large scale MD simulations have been reported. 17 An approximate energy gradient has been developed for DFTB combined with the fragment molecular orbital (FMO) method, 18–21 which neglected 22,23 the contribution of the response terms arising from the coupling of the electronic state to an embedding electrostatic potential (ESP) determined by the electron density of fragments. Ab initio based FMO/MD using gradient without response terms has been applied to short MD simulations. 24,25 In this work, response terms are derived and implemented for FMO-DFTB, enabling fast large scale QM/MD simulations. FMO-DFTB/MD simulations (100 ps) of hydrogen halides containing 2,000 atoms are performed and the results are compared to experiment. An FMO-DFTB calculation is conducted by dividing a large system into fragments and performing fragment-based DFTB calculations in the presence of an embedding ESP. After fragment (monomer) calculations converge with respect to ESP, fragment pair (dimer) calculations are performed once to take into account higher order effects such as charge



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 20

transfer between fragments. More details can be found elsewhere. 19–23 The total energy of FMO-DFTB 22 is E=

N ∑

EI′ +

N ∑

I

′ (EIJ − EI′ − EJ′ ) +

I>J

N ∑

V ∆EIJ ,

(1)

I>J

′ where N is the number of fragments, EX is the internal energy of monomer X = I or V dimer X = IJ, and the ∆EIJ represents the coupling of the charge transfer to the embed-

ding potential for dimer IJ. The exact definitions of these terms are given elsewhere. 22,23 Differentiating the total energy E with respect to a nuclear coordinate a, one obtains the gradient. ′ To obtain the exactly analytic derivative of the internal energy ∂EX /∂a (eq 21 in ref 22

for FMO-DFTB2 and eq 14 in ref 23 for FMO-DFTB3 ), one should add previously neglected response contributions −U

a,X,Y

,

U

a,X,Y

=4

virt ∑ occ ∑

a,X Y Umi Vmi ,

(2)

m∈X i∈X

a,X and Umi arises from the derivatives of the molecular orbital (MO) expansion coefficients

Cµi . 26 X and Y can denote both monomers (I) and dimers (IJ). VijX is the ESP matrix of fragment X in MO basis. For FMO-DFTB3, it is

VijX

=

∑ µν∈X

X X X Cµi Cνj Sµν

N ∑ { ∑ 1 K̸=X C∈K

2

(γAC + γBC )

1 1 + (ΓAC ∆qAX + ΓBC ∆qBX ) + (ΓCA + ΓCB )∆qCK 3 6

} ∆qCK , (3)

while for FMO-DFTB2 all Γ containing terms are omitted. Atomic orbitals (AOs) µ and ν are centered on atoms A and B, respectively. γAC and ΓAC describe the distance dependence X are the AO-based overlaps, of the Coulomb interaction in the monopole approximation. 8,9 Sµν

and ∆qAX is the Mulliken charge of atom A in X. Next, to obtain the exactly analytic derivative of the embedding energy coupled to charge


Page 5 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


V transfer, ∂∆EIJ /∂a (eq 24 in ref 22 for FMO-DFTB2 and eq 17 in ref 23 for FMO-DFTB3

), one should add the response contributions

U

a,IJ,IJ

−U

a,IJ,I

−U

a,IJ,J

N virt ∑ occ ∑ ∑

+

a,K IJ,K Umi Lmi ,

(4)

K̸=I,J m∈K i∈K

where for FMO-DFTB3

LIJ,K mi

=2

∑ ∑{

γAC ∆∆qAIJ

A∈I,J C∈K

} 2 1 K IJ IJ + ΓAC ∆∆QA + ΓCA ∆∆qA ∆qC 3 3 ∑∑ K K K K K , (5) Cνi + Cµi Cνm )Sµν × (Cµm µ∈C ν∈K

and for FMO-DFTB2 all Γ containing terms are omitted. ∆∆qAIJ is the difference between Mulliken charges in dimer IJ and monomers I and J; ∆∆QIJ A is a similar difference, but for squares of charges. Collecting all U terms for the total gradient of E in eq 1 (note that there is a cancellation of terms, see eq 16 in ref 27), one obtains the total response contribution Ra , which should be added to the previously developed gradient 22,23 neglecting response terms, and then one obtains the fully analytic FMO-DFTB gradient.

R =− a

N ∑

U

a,I,I

+

I

+

N ( ∑

N ( ∑

−U

a,IJ,IJ

+U

a,I,I

a,J,J

)

I>J

U

a,IJ,IJ

−U

a,I,IJ

−U

a,J,IJ

) +

N ∑ virt ∑ occ ∑

I>J

=

+U

a,K K Umi Lmi

(6)

K m∈K i∈K


a,K K Umi Lmi ,

K m∈K i∈K

where the Lagrangian is LK mi

=

N ∑

LIJ,K mi .

I>J I̸=K,J̸=K


(7)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 20

Ua can be obtained by solving the coupled-perturbed DFTB (CP-DFTB) equation, 28

AUa = Ba0 .

(8)

Here, A, Ua , and Ba0 are supermatrices composed of blocks (matrices) corresponding to fragments. 29 The orbital Hessian related elements are I I AI,J ij,kl = δIJ δik δjl (εj − εi ) } ∑{ ∑ 2 2 I I I I I J (γAC + γBC ) + (ΓAC ∆qA + ΓBC ∆qB ) + (ΓCA + ΓBC )∆qC qCkl,J − Cµi Cνj Sµν 3 3 C∈J µν∈I N ∑ ( ) ∑ ∑ 2 I I I − δIJ Cµi Cνj Sµν qAkl,J ΓAC + qBkl,J ΓBC ∆qCK 3 K C∈K µν∈I

(9) and a,I B0,ij =

∑ µν∈I

I + Sµν

[ I I a,I Cµi Cνj Hµν − εIj N ∑ { ∑ 1 K C∈K

2

(γAC

N ∑ I ) ∑ ∂Sµν 1( I + Sµν ΓAC qAa,I + ΓBC qBa,I ∆qCK ∂a 3 K C∈K

1 1 + γBC ) + (ΓAC ∆qAI + ΓBC ∆qBI ) + (ΓCA + ΓCB )∆qCK 3 3

]

} qCa,K

,

(10) where the population-like contribution for a pair of orbitals i and j is

qAij,I =

∑∑

I I I I I (Cµi Cνj + Cµj Cνi )Sµν ,

(11)

µ∈A ν∈I

and the contribution arising from the derivatives of integrals in the Mulliken charges is

qAa,I =

∑∑

( a,I I ˜ µν Sµν W

µ∈A ν∈I

I ∂Sµν I + Dµν ∂a

) .

(12)

I a,I is the electron density is the derivative of integrals in the Hamiltonian matrix, 22,23 Dµν Hµν


Page 7 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


˜ a,I is the density weighted overlap derivative, defined as matrix, and W µν I 1 ∑ I ∂Sρσ a,I ˜ µν DI . W =− Dµρ 2 ρσ∈I ∂a σν

(13)

The response contribution (eq 6) is rearranged as

Ra =


a,K K Umi Lmi

K m∈K i∈K

= LT Ua

(14)

= LT A−1 Ba0 = ZT Ba0 , and instead of eq 8, one solves the Z-vector equations, 26,30,31

AT Z = L .

(15)

In FMO, eq 15 is solved using the self-consistent Z-vector (SCZV) method. 29 The fully analytic gradient for FMO-DFTB2 and FMO-DFTB3 was implemented in GAMESS-US 32 and parallelized with the generalized distributed data interface (GDDI), 33 similar to RHF and DFT formulations of SCZV. 29,34–36 Typically, in FMO-DFTB one CPU core is assigned to each GDDI group. The ES-DIM approximation 37 was used in FMODFTB with a unitless threshold of 2.0. Hybrid orbital projection (HOP) operators were used to treat fragment boundaries across covalent bonds. Polyalanine was divided as 1 residue per fragment and molecular clusters were divided as 1 molecule per fragment. The mio 8 and 3ob 9,38 sets of DFTB parameters were used for DFTB2 and DFTB3, respectively, for water and polyalanine, and another set was used for halogens. 39 For hydrogen halides, the UFF 40,41 dispersion was used. We also tested the recently published halogen parameter set 42 for hydrogen halides using FMO-DFTB3 combined with the D3(BJ) 43,44 dispersion



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

correction but 10 ps MD simulations did not show improvement in the RDF peaks, so for production MD runs of hydrogen halides we used FMO-DFTB2 with UFF dispersion. A summary of 10 ps FMO-DFTB/MD results using various parameters is found in the Supporting Information. The damping function was used in all DFTB3 and in hydrogen halides calculations for atomic pairs containing hydrogen atoms with the damping exponent of ζ = 4.0. 9 In MD simulations, the spherical solvent boundary potential (SSBP) 45 was used with the force constant of 2.0 kcal/(mol·Å2 ). The NVT ensemble was used for all MD simulations except for the evaluation of the slope (Figure 2) for which NVE was used. The bath temperature was set to 300 K, unless otherwise stated. The velocity Verlet algorithm was used to propagate nuclear coordinates. After an equilibration for 10 ps at the experimental temperature (HF: 296 K, 46 HCl: 313 K, 47 HBr: 216.7 K, 48 and HI: 253 K 49 ), FMO-DFTB2/MD simulations were performed for 100 ps with the Nose-Hoover chain thermostat and a time step of 0.1 fs (one million MD steps). Atomic coordinates from every 10 fs were used to calculate RDFs. The density of liquid hydrogen halides was taken from experiment and used in setting the SSBP radius (19.87, 25.86, 25.15, and 26.87 Å for HF, HCl, HBr, and HI, respectively). The accuracy of the implemented fully analytic FMO-DFTB gradient was evaluated for several representative systems computing the difference between the analytic and numerical gradients. The results are shown in Figure 1 and Table 1. Response terms make a contribution on the order of 10−4 a.u./Bohr, and therefore for reliable MD simulations they should be included in the gradient. Computational timings were obtained as an average of 10,000 NVT MD steps (SSBP was not used). The additional cost for FMO-DFTB3 relative to FMO-DFTB2 is less than 10%. The computation of the response contributions increases the timing by about 20%. In Table 1, the water cluster and the capped polyalanine (α-helix) have roughly the same numbers of atoms (192 and 212, respectively), and the latter timing is almost three times longer because the fragment size is larger (3 and 10 atoms, respectively).


Page 8 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Gradient Error / a.u.∙(Bohr)-1

Page 9 of 20

0.0006 0.0004 0.0002 0.0000 -0.0002 -0.0004 -0.0006

Gradient Element (a)

Figure 1: Differences between analytic and numerical FMO-DFTB3 gradients for (H2 O)64 . Red filled and black empty circles indicate the values (gradient errors) with and without the response terms (eq 14).

Table 1: Maximum (MAX) and root-mean-square deviations (RMSD) between the analytic and numerical gradients (a.u./Bohr) for FMO-DFTBm (m=2 or 3) with and without response terms (eq 14). system (H2 O)64 (H2 O)64 (H2 O)64 (H2 O)64 (Ala)20 (Ala)20 (Ala)20 (Ala)20 a

m 2 2 3 3 2 2 3 3

response yes no yes no yes no yes no

T , sa 0.162 0.130 0.176 0.136 0.460 0.369 0.484 0.396

RMSD 0.000 002 0.000 140 0.000 005 0.000 186 0.000 014 0.000 142 0.000 013 0.000 157

MAX 0.000 006 0.000 657 0.000 018 0.000 595 0.000 039 0.000 650 0.000 045 0.000 731

Computational time for a single point analytic gradient on a dual Xeon E5-2650 node (8×2 cores), an average of 10,000 MD steps.



Another way to check the gradient accuracy is to evaluate the slope of the energy fluctuations versus the time step. For (H2 O)64 , NVE FMO-DFTB/MD simulations were performed for 10 ps (time step 0.1 fs) using the Nose-Hoover chain thermostat at 300 K for equilibration, and the following 100 steps were used to compute the slope, shown in Figure 2. The theoretical slope should be 2, and the computed FMO-DFTB2 and FMO-DFTB3 slopes are 2.00 and 1.99, respectively, when the gradient with the response terms is used in MD simulations. For the calculation using the gradient without the response terms, the respective slopes are 1.40 and 0.87. This clearly emphasizes the necessity of including the response terms during MD simulations with FMO-DFTB. These terms are thus included in all calculations presented below. 1

RMSE / kcal∙(mol)-1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.1

0.01

0.001

FMO-DFTB2 without Response (slope: 1.40) FMO-DFTB3 without Response (slope: 0.87) FMO-DFTB2 with Response (slope: 2.00) FMO-DFTB3 with Response (slope: 1.99)

0.1

0.5

1.0

Time Step / fs

Figure 2: Double logarithmic plot of the root-mean-square deviations of total energy (RMSE; in kcal/mol) versus time step (in fs) from the simulations of (H2 O)64 with FMO-DFTB2 (red squares) and FMO-DFTB3 (black circles) with (filled) and without (empty) response terms. To elucidate the accuracy of FMO-DFTB3 in comparison to full DFTB3 (i.e., without fragmentation), RDFs from NVT MD simulations (100,000 steps with a time step of 0.1 fs) are compared for (H2 O)64 . Table 2 shows that FMO-DFTB reproduces RDFs from full DFTB with the deviation for peaks typically below 0.04 Å. The largest deviation (0.08 Å) is for the very broad peak at 5.2 Å . A comparison of timings is shown in Figure 3 (per MD step averaging timings for 10,000 steps). Even for the smallest system containing 192 atoms, FMO-DFTB is ten times faster 10 ACS Paragon Plus Environment

Page 10 of 20

Page 11 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Peaks in the RDFs in (H2 O)64 obtained with FMO-DFTB3 (full DFTB3 values are given in parentheses). Pair peak O-H 1 O-H 2 O-H 3 O-H 4 O-O 1 O-O 2

position / Å height 0.97 (0.97) 31.26 (31.34) 1.87 (1.84) 1.15 ( 1.14) 3.19 (3.19) 1.19 ( 1.18) 5.19 (5.27) 0.61 ( 0.60) 2.77 (2.73) 3.37 ( 3.23) 5.05 (5.02) 0.71 ( 0.68)

than the corresponding full DFTB calculation. For the medium system of 384 atoms, the speed-up is a factor of 30, and for the large system of 768 atoms, the speed-up is 108. Finally, FMO-DFTB2/MD is applied to liquid hydrogen halides HF, HCl, HBr, and HI. NVT MD simulations were performed at the experimental temperatures listed in Computational Methods for clusters containing 1000 hydrogen halides molecules (2000 atoms). For such a large system, full DFTB2 for (HF)1000 took 1816 seconds for a single point gradient calculation, whereas FMO-DFTB2 took only 1.8 seconds. Calculated RDFs are plotted in Figure 4, and a comparison with experiment is shown in Table 3. The first X-X peak (X = halogen) has a clear trend of shifting to the right (to longer distances) as one goes down the periodic table. This indicates that the intermolecular interactions are weaker for heavier atoms, which must be related to the weaker hydrogen bonding (and shielding of valence electrons by the increased core). On the other hand, the dispersion interaction in general increases with the atomic weight. Comparing FMO-DFTB2 peaks with experiment, one can observe a reasonable agreement, with the largest deviation found for HBr (0.2 Å). With the exception of HCl, FMO-DFTB2 appears to have a tendency to overestimate the peak positions, which might be related to some drawbacks in describing hydrogen bonding and dispersion (the UFF model). In this work, fully analytic FMO-DFTB2 and FMO-DFTB3 gradients have been derived, and it has been shown that the gradients are accurate and that FMO-DFTB MD simulations


80

Page 12 of 20

73.26 s (47.19)

60

40

56

O )2

(H

2

4

O )1 2

(H

O )6

(H

Conventional DFTB3 (DFTB2)

2

28

(H

O )2

4

O )1 2

(H

O )6 2

(H

56

9.25 s (8.48) 1.93 s (1.89)

28

0.68 s (0.52) 0.31 s (0.28) 0.19 s (0.16)

2

20

0

FMO-DFTB3 (FMO-DFTB2)

Figure 3: Wall-clock timings of 1 MD step for water clusters performed on a dual Xeon E5-2650 node (8×2 cores) with full DFTB3 and FMO-DFTB3 methods (the timings for full DFTB2 and FMO-DFTB2 are shown in parentheses). 3.0 Radial Distribution Function

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Wall-Clock Timing per MD Step / second


F

2.5

Br Cl

2.0

I

F-F Cl-Cl Br-Br I-I

1.5 1.0 0.5 0.0

1

2

3

4 5 Distance / Å

6

7

8

Figure 4: X-X radial distribution functions, shown with black solid (X = F), red dotted (Cl), blue dashed (Br), and purple dashed-dot (I) lines. 12 ACS Paragon Plus Environment

Page 13 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: First X-X peaks (in Å) in the calculated and experimental RDFs of liquid hydrogen halides HX. X FMO-DFTB2 F 2.60 Cl 3.86 Br 4.13 I 4.54

Exp. 2.51 46 3.88 47 3.93 48 4.42 49

reasonably reproduce both full DFTB RDF peaks as well as experiment. FMO-DFTB has been shown to greatly improve the efficiency of full DFTB MD simulations, with a speed-up of 108 obtained for a water cluster containing 256 molecules. 100 ps QM/MD simulations of clusters containing 2000 atoms have been reported. FMO-DFTB/MD looks like a promising tool for QM/MD simulations of large systems, considering its low computational scaling and high parallel efficiency. 22 In future, FMO-DFTB/MD may be combined with enhanced sampling techniques such as replica-exchange MD, 50,51 further improving the efficiency of sampling.

Acknowledgement Y.N. and H.N. were supported by a Research Fellowship of the Japan Society for Promotion of Science for Young Scientists. This work is partially supported by JSPS KAKENHI Grant Number 15H06316 (Grant-in-Aid for Research Activity Start-up). D.G.F. acknowledges support of the Next Generation Super Computing Project, Nanoscience Program (MEXT, Japan) and D.G.F. and S.I. thank Computational Materials Science Initiative (CMSI, Japan) for financial support. D.G.F. thanks Prof. Kazuo Kitaura for many fruitful discussions about the FMO method. Supporting Information Available: A comparison of RDF peaks obtained with FMODFTB/MD for hydrogen halides using several theoretical models is provided.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Zhao, G.; Perilla, J. R.; Yufenyuy, E. L.; Meng, X.; Chen, B.; Ning, J.; Ahn, J.; Gronenborn, A. M.; Schulten, K.; Aiken, C. et al. Mature HIV-1 Capsid Structure by Cryo-Electron Microscopy and All-Atom Molecular Dynamics. Nature 2013, 497, 643–646. (2) Gordon, M. S.; Smith, Q. A.; Xu, P.; Slipchenko, L. V. Accurate First Principles Model Potentials for Intermolecular Interactions. Ann. Rev. Phys. Chem. 2013, 64, 553–578. (3) Schmitt, U. W.; Voth, G. A. Multistate Empirical Valence Bond Model for Proton Transport in Water. J. Phys. Chem. B 1998, 102, 5547–5551. (4) Warshel, A.; Karplus, M. Calculation of Ground and Excited State Potential Surfaces of Conjugated Molecules. I. Formulation and Parametrization. J. Am. Chem. Soc. 1972, 94, 5612–5625. (5) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and DensityFunctional Theory. Phys. Rev. Lett. 1985, 55, 2471–2447. (6) Ufimtsev, I. S.; Martinez, T. J. Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients, Geometry Optimization, and First Principles Molecular Dynamics. J. Chem. Theory Comput. 2009, 5, 2619–2628. (7) Stewart, J. J. P. Application of the PM6 Method to Modeling Proteins. J. Mol. Model. 2009, 15, 765–805. (8) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260–7268. (9) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge


Page 14 of 20

Page 15 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Density-Functional Tight-Binding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931–948. (10) Irle, S.; Zheng, G.; Wang, Z.; Morokuma, K. The C60 Formation Puzzle “Solved”: QM/MD Simulations Reveal the Shrinking Hot Giant Road of the Dynamic Fullerene Self-Assembly Mechanism. J. Phys. Chem. B 2006, 110, 14531–14545. (11) Goyal, P.; Qian, H.-J.; Irle, S.; Lu, X.; Roston, D.; Mori, T.; Elstner, M.; Cui, Q. Molecular Simulation of Water and Hydration Effects in Different Environments: Challenges and Developments for DFTB Based Models. J. Phys. Chem. B 2014, 118, 11007–11027. (12) Otto, P.; Ladik, J. Investigation of the Interaction between Molecules at Medium Distances: I. SCF LCAO MO Supermolecule, Perturbational and Mutually Consistent Calculations for Two Interacting HF and CH2 O Molecules. Chem. Phys. 1975, 8, 192– 200. (13) Gao, J. Toward a Molecular Orbital Derived Empirical Potential for Liquid Simulations. J. Phys. Chem. B 1997, 101, 657–663. (14) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632–672. (15) Hu, H.; Lu, Z.; Elstner, M.; Hermans, J.; Yang, W. Simulating Water with the SelfConsistent-Charge Density Functional Tight Binding Method: From Molecular Clusters to the Liquid State. J. Phys. Chem. A 2007, 111, 5685–5691. (16) Giese, T. J.; Chen, H.; Dissanayake, T.; Giambau, G. M.; Heldenbrando, H.; Huang, M.; Kuechler, E. R.; Lee, T.-S.; Panteva, M. T.; Radak, B. K. et al. A Variational LinearScaling Framework to Build Practical, Efficient Next-Generation Orbital-Based Quantum Force Fields. J. Chem. Theory Comput. 2013, 9, 1417–1427.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(17) Xie, W.; Orozco, M.; Truhlar, D. G.; Gao, J. X-Pol Potential: An Electronic StructureBased Force Field for Molecular Dynamics Simulation of a Solvated Protein in Water. J. Chem. Theory Comput. 2009, 5, 459–467. (18) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment Molecular Orbital Method: an Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 1999, 313, 701–706. (19) Fedorov, D. G.; Kitaura, K. Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method. J. Phys. Chem. A 2007, 111, 6904–6914. (20) Fedorov, D. G.; Nagata, T.; Kitaura, K. Exploring Chemistry with the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2012, 14, 7562–7577. (21) Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. ElectronCorrelated Fragment-Molecular-Orbital Calculations for Biomolecular and Nano Systems. Phys. Chem. Chem. Phys. 2014, 16, 10310–10344. (22) Nishimoto, Y.; Fedorov, D. G.; Irle, S. Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2014, 10, 4801–4812. (23) Nishimoto, Y.; Fedorov, D. G.; Irle, S. Third-Order Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2015, 636, 90–96. (24) Komeiji, Y.; Nakano, T.; Fukuzawa, K.; Ueno, Y.; Inadomi, Y.; Nemoto, T.; Uebayasi, M.; Fedorov, D. G.; Kitaura, K. Fragment Molecular Orbital Method: Application to Molecular Dynamics Simulation, ‘Ab Initio FMO-MD’. Chem. Phys. Lett. 2003, 372, 342–347.


Page 16 of 20

Page 17 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(25) Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y.; T. Ishikawa, T. N. How Does an SN 2 Reaction Take Place in Solution? Full Ab Initio MD Simulations for the Hydrolysis of the Methyl Diazonium Ion. J. Am. Chem. Soc. 2008, 130, 2396–2397. (26) Yamaguchi, Y.; Schaefer III, H. F.; Osamura, Y.; Goddard, J. A New Dimension to Quantum Chemistry: Analytical Derivative Methods in Ab Initio Molecular Electronic Structure Theory; Oxford University Press: New York, 1994. (27) Nagata, T.; Fedorov, D. G.; Kitaura, K. Analytic Gradient for the Embedding Potential with Approximations in the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2012, 544, 87–93. (28) Witek, H. A.; Irle, S.; Morokuma, K. Analytical Second-Order Geometrical Derivatives of Energy for the Self-Consistent-Charge Density-Functional Tight-Binding Method. J. Chem. Phys. 2004, 121, 5163–5170. (29) Nagata, T.; Brorsen, K.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. Fully Analytic Energy Gradient in the Fragment Molecular Orbital Method. J. Chem. Phys. 2011, 134, 124115. (30) Handy, N. C.; Schaefer III, H. F. On the Evaluation of Analytic Energy Derivatives for Correlated Wave Functions. J. Chem. Phys. 1984, 81, 5031–5033. (31) Ochsenfeld, C.; Gordon, M. S. A Reformulation of the Coupled Perturbed SelfConsistent Field Equations Entirely within a Local Atomic Orbital Density MatrixBased Scheme. Chem. Phys. Lett. 1997, 270, 399–405. (32) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. J.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347–1363.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(33) Fedorov, D. G.; Olson, R. M.; Kitaura, K.; Gordon, M. S.; Koseki, S. A New Hierarchical Parallelization Scheme: Generalized Distributed Data Interface (GDDI), and an Application to the Fragment Molecular Orbital Method (FMO). J. Comput. Chem. 2004, 25, 872–880. (34) Nakata, H.; Fedorov, D. G.; Nagata, T.; Yokojima, S.; Ogata, K.; Kitaura, K.; Nakamura, S. Unrestricted Hartree-Fock Based on the Fragment Molecular Orbital Method: Energy and Its Analytic Gradient. J. Chem. Phys. 2012, 137, 044110. (35) Nakata, H.; Schmidt, M. W.; Fedorov, D. G.; Kitaura, K.; Nakamura, S.; Gordon, M. S. Efficient Molecular Dynamics Simulations of Multiple Radical Center Systems Based on the Fragment Molecular Orbital Method. J. Phys. Chem. A 2014, 118, 9762–9771. (36) Brorsen, K. R.; Zahariev, F.; Nakata, H.; Fedorov, D. G.; Gordon, M. S. Analytic Gradient for Density Functional Theory Based on the Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2014, 10, 5297–5307. (37) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Use of Approximate Electrostatic Potential. Chem. Phys. Lett. 2002, 351, 475–480. (38) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9, 338–354. (39) Kuba˘r, T.; Bodrog, Z.; Gaus, M.; Köhler, C.; Aradi, B.; Frauenheim, T.; Elstner, M. Parametrization of the SCC-DFTB Method for Halogens. J. Chem. Theory Comput. 2013, 9, 2939–2949. (40) Zhechkov, L.; Heine, T.; Patchkovskii, S.; Seifert, G.; Duarte, H. A. An Efficient a Posteriori Treatment for Dispersion Interaction in Density-Functional-Based Tight Binding. J. Chem. Theory Comput. 2005, 1, 841–847.


Page 18 of 20

Page 19 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(41) Rappé, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard III, W. A.; Skiff, W. M. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024–10035. (42) Kubillus, M.; Kubař, T.; Gaus, M.; Řezáč, J.; Elstner, M. Parameterization of the DFTB3 Method for Br, Ca, Cl, F, I, K, and Na in Organic and Biological Systems. J. Chem. Theory Comput. 2015, 11, 332–342. (43) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (44) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465. (45) Beglov, D.; Roux, B. Finite Representation of an Infinite Bulk System: Solvent Boundary Potential for Computer Simulations. J. Chem. Phys. 1994, 100, 9050–9063. (46) McLain, S. E.; Benmore, C. J.; Siewenie, J. E.; Urquidi, J.; Turner, J. F. C. On the Structure of Liquid Hydrogen Fluoride. Angew. Chem. Int. Ed. 2004, 43, 1952–1955. (47) Andreani, C.; Ricci, M. A.; Nardone, M.; Ricci, F. P.; Soper, A. K. Orientational Correlations and Hydrogen Bonding in Liquid Hydrogen Chloride. J. Chem. Phys. 1997, 107, 214–221. (48) Andreani, C.; Menzinger, F.; Ricci, M. A.; Soper, A. K.; Dreyer, J. Neutron Diffraction from Liquid Hydrogen Bromide: Study of the Orientational Correlations. Phys. Rev. B 1994, 49, 3811–3820. (49) Andreani, C.; Nardone, M.; Ricci, F. P.; Soper, A. K. Neutron-Diffraction Study of Liquid Hydrogen Iodide. Phys. Rev. A 1992, 46, 4709–4716.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(50) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141–151. (51) Fedorov, D. G.; Sugita, Y.; Choi, C. H. Efficient Parallel Implementations of QM/MMREMD (Quantum Mechanical/Molecular Mechanics-Replica-Exchange MD) and Umbrella Sampling: Isomerization of H2 O2 in Aqueous Solution. J. Phys. Chem. B 2013, 117, 7996–8002.


Page 20 of 20

Large-Scale Quantum-Mechanical Molecular Dynamics Simulations

Recommend Documents