Can Systematic Molecular Fragmentation Be ... - ACS Publications

Subscriber access provided by RYERSON UNIVERSITY

Article

Can Systematic Molecular Fragmentation be Applied to Direct Ab Initio Molecular Dynamics Michael Anthony Collins J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b08739 • Publication Date (Web): 31 Oct 2016 Downloaded from http://pubs.acs.org on November 1, 2016

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 A 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 55

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

Can Systematic Molecular Fragmentation be Applied to Direct Ab Initio Molecular Dynamics

Michael A. Collins* Research School of Chemistry, Australian National University, Canberra, ACT, Australia

ACS Paragon Plus Environment

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

Page 2 of 55

Abstract

This paper demonstrates that systematic molecular fragmentation can be applied to direct ab initio molecular dynamics of solvated molecules under periodic boundary conditions. A method for rapidly updating the fragmentation of water at each time step in a simulation is presented and tested. This approach reduces the time required for implementation of systematic molecular fragmentation at each time step, in a highly connected system like water, from an excessively long time to a feasible value.


2

Page 3 of 55

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. Introduction Within the Born-Oppenheimer approximation, ab initio quantum chemistry calculates the total electronic energy of molecules, and using energy gradients, can evaluate equilibrium geometries. With the additional calculation of second derivatives of the energy with respect to the atomic positions, thermochemical and other properties can be evaluated, for comparison with experiment or prediction of property values. Similarly, so-called reaction paths can be identified on the molecular energy surface, including saddle points or variationally determined "transition states". Using various forms of transition state theory, chemical reaction rates can be estimated.

1

However, a complete, reliably accurate, estimation of dynamical processes such as chemical reactions requires an evaluation of the motion of the atomic nuclei on the "global" molecular potential energy surface (PES) which ab initio quantum chemistry can provide, in principle. For some reactions in the gas phase, involving up to about 6 atoms, fully quantum calculations

2,3

of the nuclear motion on an accurate ab initio PES have

been reported.4 At present, such fully quantum nuclear motion can only be evaluated in 8 dimensions, or less, due to the rapid increase in computation time and computer memory with increasing number of nuclear degrees of freedom.5 Classical mechanics can be used to approximate the nuclear motion, and so-called quasiclassical trajectory simulations of gas phase reactions have been pursued for several decades,6,7 and continue to provide insight into reaction mechanisms.8 However, even this approach is limited by the availability of accurate PESs. The configurations of a nonlinear molecule, of N atoms, occupy a space of dimension 3N-6. An ab initio calculation of the electronic energy at just one configuration may be very computationally expensive. A priori, one cannot know


3


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 55

which regions of this high dimensional space are relevant to any particular dynamical process. Nonetheless, some methods have been developed over the last twenty years to successfully construct highly accurate PES in the relevant regions of configuration space, using ab initio calculations at (only) several tens of thousands of molecular configurations.

9,10,11,12

It may be that these approaches to constructing PESs can be

extended to reactions involving tens of atoms, as many nuclear degrees of freedom can be expected to be "spectator modes". That is, many degrees of freedom will play little part in the chemical reaction, and only vary over a small range of configurations during the process. Time will tell. The progression of highly accurate dynamics on reliably accurate ab initio PESs for small gas phase reactions is very different from the path of progress followed for reactions of very large molecules in condensed phases; for example, the reactions of proteins in aqueous solution. Again over several decades, classical trajectory or molecular dynamics (MD) simulations have been used to study the motion of large biological systems in water.

13,14,15

Originally, these studies

employed so-called molecular dynamics force fields. This approach expresses the electronic energy of a molecule, and surrounding solvent, in terms of empirical functions of the bond lengths, valance angles and torsion (dihedral) angles in the molecule(s), plus functions to account for long range Coulomb interactions. Over the years, these empirical functions have been refined to improve the agreement between classical simulation and experiment. However, such methods cannot describe chemical reactions. To account for the making and breaking of bonds in such large systems, so-called QM/MM methods have been developed.16-20 There are a number of alternative approaches, but in general terms, the motion of most atoms in the large molecule-solvent system are governed by a molecular dynamics force field. In a


4

Page 5 of 55

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


relatively small part of the system, in the region where bonds are broken or formed, the motion of the atoms is governed by an ab initio or semi-empirical force field. Generally, due to computational cost, only a relatively modest level of ab initio theory or density function theory can be employed, and normally only equilibrium geometries and saddle points are investigated, not dynamics. However, semiempirical methods based on valence bond theory are used to describe chemical reaction dynamics. Global PES based solely on reliable ab initio calculations have not been attempted due to the exorbitant cost of such calculations at even one molecular configuration, and the absence of any method to combine such energy calculations into a global PES. However, over the last decade or so, methods have been developed to radically reduce the cost of ab initio calculations on large molecules. Based on the idea that chemical functionality is a local phenomena, one breaks the molecule into fragments, performs calculations on the fragments, then reconstitutes the value of the energy or property from the corresponding values for the fragments. Such methods include the effective fragment potential method

21,22

, the X-Pol method23-26, the

fragment molecular orbital method27-29, and energy-based fragmentation methods. The later involves breaking the molecule into fragments, evaluating the energy of each fragment, and recombining the fragment energies to estimate the total electronic energy. Several groups have developed such approximations to the energy, and other molecular properties, in recent times.30-50 A more complete description of this area was reported in a recent review.51 The ab initio quantum chemistry computational time for these fragmentation methods scales only linearly with the size of the molecule. The variation in computation time between different methods is mostly determined by the different sizes of the fragments. This paper employs the systematic


5


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 55

34-36,52,53

molecular fragmentation by annihilation method (SMFA)

, which will be

briefly reprised below. Breaking a large molecule, or molecule-solvent system, into small fragments reduces the cost of the ab initio calculation enormously because this cost is proportional to a high power (commonly from 4 to 7) of the number of basis functions in a calculation. A calculation for a system of 1000 atoms, compared to a calculation for 10 typical atoms in that system, requires a computation time that is longer by a factor of 100 raised to the power 4 - 7. For SMFA, and other related methods, the number of fragments is only linearly proportional to the total number of functional groups in the system. Hence, the total computational cost is only of the order of N times the cost of a calculation on a small fragment. Most importantly, each ab initio calculation on a fragment is completely independent of all other such calculations. Hence, if one has sufficient processors available, all these fragment calculations can be carried out simultaneously. Then, the "wall-time" or actual time required to carry out an ab initio calculation of a large system is independent of the size of the system; it is just the time for the ab initio calculation of the largest fragment. How large is this fragment? The size varies with fragmentation methods and the structure of the molecules. However, for SMFA the largest fragments of proteins contain of the order of 20 atoms. So, if the energy of a large molecule at some molecular configuration is given by a sum of the energies of many "small" fragments, then the global PES for the whole molecule is given by the same sum of PESs for the fragments. As noted above, the construction of an ab initio PES for a molecule of about 20 atoms is probably at or beyond current capability. To construct many hundreds of such PES is a daunting task. Moreover, as discussed below, a classical simulation requires that all the PES


6

Page 7 of 55

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


that contribute to the global PES must be evaluated in a short time, for the simulation to be feasible. Depending on the PES construction method, PES evaluation time may be a problem. However, this approach has been successfully carried out in a few cases. The minimum energy path for reaction of hydrogen with n-pentane was evaluated with such fragment PES.54 Classical simulations of reactive scattering of a hydrogen atom from a hydrogen covered Silicon crystal surface has been reported.

55

Adsorption of H2 in a metal-organic-framework material has also been studied with quantum diffusion Monte-Carlo using a fragment-based global PES of this type.

56

However, in this case the PES was only a function of six coordinates. In the case of hydrogen reacting with a hydrogen covered Silicon surface, the largest fragment contained 9 atoms. However, 9 is still much less than 20. Perhaps the solution of the PES construction problem is simply not to construct them at all. We could perform ab initio calculations of the forces on each atom in each fragment, whenever such forces are needed in a classical simulation. Indeed, this approach has also been developed over many years in the study of reactions in the gas phase.

57,58,59

This approach to reaction dynamics is usually

called "direct dynamics" or "direct ab initio molecular dynamics" or "ab initio molecular dynamics". Aside: Here we are not concerned with the methods for simultaneous evaluation of nuclear and electronic motion of the type developed by Car and Parinello,

60

which are also usually referred to as "ab initio molecular

dynamics". Direct ab initio molecular dynamics has thus far been applied to small molecule reaction dynamics.61,62,63 The question addressed in this paper, is whether fragmentation methods like SMFA can be used to extend direct ab initio dynamics to reactions of large (biological) molecules in solution. A question that has attracted recent interest.64


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 8 of 55

A classical simulation of a biological process is carried out by numerical integration of the classical equations of motion. This entails calculating the forces on all the atoms at the current molecular configuration. We can do this by simultaneously calculating the forces for each fragment and adding them up appropriately. The molecular configuration is then estimated at a short "time step" later. Typically this time step is of the order of 1 - 3 fsec. To follow a chemical or physical change in the system requires us to repeat this process for a nanosecond or, more probably, hundreds of nanoseconds. That means at least 106 - 108 time steps. Unfortunately, even using SMFA to avoid calculations of large molecules, a reliable ab initio calculation of the energy and forces for an organic molecule (fragment) of even just 20 atoms is still seriously time consuming in this context. If an ab initio calculation of the forces in a typical organic molecule of about 20 atoms took just 10 sec, then a simulation of a protein in water would take 107 - 109 sec, or about 116 - 11,600 days, to describe the molecular motion for 1 - 100 nsec. Chemically accurate ab initio methods might require much longer calculation times. Clearly, a direct ab initio molecular dynamics approach to biological dynamics will require substantial advances in reducing the computation time for accurate forces in moderate-sized organic molecules. This is not the subject of this paper. Here, we examine an essential auxiliary task: The methods that are necessary to implement fragmentation (here the SMFA approach) in direct ab initio dynamics in a computationally feasible way. Any fragmentation method applied to a given molecule/solvent configuration, results in a set of molecular fragments. In a MD simulation the molecule/solvent configuration changes at each time step. To retain the accuracy of the fragmentation approximation, the fragments must change with time. That is, the SMFA fragmentation procedure must be repeated, again and again.


8

Page 9 of 55

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


Unfortunately, the time required for SMFA to fragment a set of hundreds or thousands of atoms varies enormously depending on the degree of connectivity of the atoms, from a few seconds to as much as 101 - 103 secs for highly connected materials (such as bulk water). The purpose of this paper is to introduce a new algorithm for sequentially modifying the SMFA fragmentation, which reduces the computation time for fragmentation to a sufficiently short time to make SMFA feasible for direct ab initio MD. In SMFA the fragmentation is determined by the bonding in the molecule/solvent system. In a covalent molecule the bonding only changes during a chemical reaction. Hence, during a dynamical simulation, the fragments comprising a covalent molecule change only rarely, when a reaction occurs. However, a molecule in water, for example, is connected to surrounding water molecules only by relatively weak hydrogen bonds, or perhaps by Van der Waals interactions only. Similarly, the water molecules are only bonded to each other by hydrogen bonds. As we shall see in detail below, such bonds are formed and broken often. Hence, it is the solvent and the solvent-molecule connections that change most rapidly with time, and that require the fragmentation approximation to be re-evaluated again and again in a classical simulation. For this reason, this paper deals in most detail with the simulation of a solvent (water) which poses the most demands on SMFA, or any fragmentation scheme, in terms of the need to repeat the fragmentation frequently. The paper is arranged as follows. Section 2 presents some details of the test case employed herein to illustrate the method. This is a simulation "box" of water molecules under periodic boundary conditions. Section III describes the methods and algorithms, while results are presented in Section 4. The paper finishes with a discussion.


9


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 10 of 55

2. Test Case The molecular dynamics program Gromacs

15

was used to simulate the

classical dynamics of water at a temperature of 300 K and atmospheric pressure. The SPC model65 was used to describe the interaction between the water molecules, and 66

the Particle Mesh Ewald method

was employed to account for long range

interactions. The simulation "box" contains 270 water molecules and was cubic in shape, with a side of approximately 20 Å. Since the simulation was carried out at constant pressure the box dimensions varied slightly during the simulation. The simulation was carried out under periodic boundary conditions. Hence, whenever the centre of mass of a water molecule moved outside the box, this molecule was replaced with a molecule of the same structure, translated to the other side of the box. Under periodic boundary conditions, boundary effects are minimised because water molecules near the side of the box interact with water molecules outside the box; these later molecules are translated images of those just inside the opposite side of the box. Hence, all 270 water molecules experience similar interaction environments. The simulation was carried out for 1 nsec with a time step of 1 fsec. The atomic coordinates for the last 10,000 time steps (10 psec) were used here as a sample of typical atomic configurations. No doubt a more sophisticated model of water would produce a slightly different sample of atomic configurations. However, it seems reasonable to use this model to provide a sample of typical environments around each H2O molecule and typical rates of change of those environments under standard conditions. Figure 1 depicts the configuration of the water molecules at one time step.


10

Page 11 of 55

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


Figure 1. A figurative representation of the positions of all 270 water molecules at one time step in the simulation. 3. Methods 3.1. Review of SMFA The original algorithm for systematic molecular fragmentation was modified some years ago to a method entitled systematic molecular fragmentation by annihilation (SMFA). 52 The approach has been described in detail recently, 53 so only a brief summary of relevant aspects is presented herein. Bonding and groups Single bonds between atoms are assigned when the atom-atom distance falls below a tolerance which is set by the covalent radii of the atoms. Multiple bonds between atoms are assigned when the atom-atom distance falls below a tolerance which is set from experimentally determined values

67,68

for each pair of atoms. See

the Supporting Information for details. A hydrogen bond is taken to exist between a hydrogen atom and an atom (A) from groups 5, 6 or 7 in the periodic table if the atom-atom distance RHA obeys:

RHA ≤ VdW (H ) + VdW (A) − 0.32


(3.1)

11


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 12 of 55

where VdW denotes the Van der Waals radius of an atom (in Angstrom). Again, see the Supporting Information for details. For hydrogen bonds to oxygen atoms, RHO = 2.4 Å. The atoms are collected in functional groups; atoms connected by multiple bonds are in the same group, hydrogen atoms are in the group of the nearest "heavy atom". A molecule (or solvent) is then viewed as a set of functional groups connected by single bonds, or perhaps by hydrogen bonds. Fragmentation of molecules and clusters Molecules are broken into small molecular fragments using an algorithm wherein one or more groups are annihilated or removed in a sequence of steps, that preserves account of the proximity of groups in terms of bonding. For a given value of the parameter Level, for each group in the molecule, the final set of fragments contains at least one fragment containing that group and any group that is within Level bonds from it. As a simple example, n-pentane is a chain of 5 groups, and we would write the fragmentation pattern of n-pentane, at Level = 1, 2 and 3, respectively, as 12345 → 12 + 23 + 34 + 45 - 2 - 3 - 4; 12345 → 123 + 234 + 345 - 23 - 34; 12345 → 1234 + 2345 - 234.

(3.2)

Each fragment as an assigned sign. Both sides of these equations contain exactly the same groups and atoms. When covalent bonds are broken in forming these fragments, "capping" hydrogen atoms are appended to each group to restore the valence, as previously described.

69

When a hydrogen bond is broken, obviously a cap is not

required. For a general molecule-solvent system, M, the decomposition into a total of Nfrag fragments at some value of Level is written as


12

Page 13 of 55

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


N frag

M → ∑ fn Fn ,

(3.3)

n=1

where Fn denotes a fragment and fn is an integer coefficient (typically ±1). The electronic energy, E, of the original molecule or cluster is given by the rhs of these fragment equations as: N frag

E=

∑

fn E ( Fn ) ,

(3.4)

n=1

where E(Fn) denotes the electronic energy of a molecular fragment. Previous papers have detailed how additional terms can be added to the rhs of Eq. (3.4) to account for long range interactions between groups in molecules. In this paper, we shall simply consider the basic fragmentation approximation of Eqs (3.3) and (3.4), and discuss such additional terms in the last section. Several previous papers have reported extensive tests of the SMF (and SMFA) approximation to the molecular energy,

35,36

and for other properties including the

molecular structure, vibrational frequencies 70,71 and NMR chemical shifts 72. The test molecules include most types of organic molecules, proteins, water clusters, and systems containing formal charges. When covalent bonds are broken in the fragmentation, tests have shown that the SMF and SMFA approximations (including long range corrections) give energies that are accurate to a few kJ mol-1, when Level = 3 or Level = 4. Of particular relevance to this paper, we note that extensive tests were carried out for many water clusters of up to 64 monomers, with various basis sets and levels of ab initio theory.

73

For such clusters, only hydrogen bonds are broken in the

fragmentation procedure. Those tests demonstrated that for Level = 2, the total binding energy of the clusters was given on average to within 1.3%, compared to an ab initio calculation for the whole cluster.

73

Moreover, both the structure and

vibrational frequencies of water clusters were accurately reproduced using SMFA. 71


13


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 14 of 55

Fragmentation of periodic systems SMFA has been applied to crystals and crystal surfaces.

69,74,55,70

A simple

example illustrates the method. Consider a periodic lattice of atoms (or groups), L, in one dimension, which we denote as L = ....... An-1 B n-1 C n-1 D n-1 An Bn Cn Dn An+1 B n+1 C n+1 D n+1 ........ (3.5) Just like n-pentane in Eq. (3.2) we can fragment this lattice at (say) Level = 2 to write

L → ....... + Cn−1 Dn−1 An + Dn−1 An Bn + An Bn Cn + Cn Dn An+1 + Dn An+1 Bn+1 ... (3.6) .... − Dn−1 An − An Bn − BnCn − Cn Dn − Dn An+1 ...... ∞

→

∑ n=−∞

 Dn−1 An Bn + An BnCn + Cn Dn An+1  − D A − A B − B C − C D   n−1 n n n n n n n

The terms in the brackets in Eq. (3.6) constitute one unit cell (ABCD), and the energy of this lattice, per unit cell, is given by ∞

E=

∑ n=−∞

 E ( Dn−1 An Bn ) + E ( An BnCn ) + E ( Cn Dn An+1 )     −E ( Dn−1An ) − E ( An Bn ) − E ( BnCn ) − E ( Cn Dn ) 

(3.7)

Notice that the fragments in Eqs (3.6) and (3.7) contain groups in more than one unit cell, as groups are bonded across unit cell "boundaries". The extension of SMFA to a periodic lattice is achieved by taking a set of groups that contains all the groups in one "central" unit cell and all the groups that are within Level +1 bonds of them. This collection of groups is fragmented using the same algorithm used for molecules, and then all fragments that contain only groups outside the "central" unit cell are discarded. Similarly, any fragment that is just a lattice translation of another is discarded. This leaves the set of fragments that reproduce the infinite periodic lattice by translation, as in the second line of Eq. (3.6). The algorithm applies to any three dimensional crystal structure, or to a cleaved crystal which is truncated at a crystal plane. 74


14

Page 15 of 55

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 same periodic version of SMFA can be directly applied to a simulation of a fluid under periodic boundary conditions. The simulation "box" is the "central" unit cell. As an example, consider just one configuration of 270 water molecules from the simulation described above. Each water molecule is a group in SMFA, where, in this case, the groups are connected by hydrogen bonds. Fragmentation breaks these hydrogen bonds. For that configuration, at Level = 2, there are 1118 fragments that describe the 270 water molecules under periodic boundary conditions. Each fragment contains a small number of water molecules .The distribution of fragment sizes are shown in Figure 2.

Figure 2. The distribution of fragment sizes in 1118 fragments that represent the 270 water molecules at one (arbitrary) configuration in the simulation.


15


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 16 of 55

As the figure shows, fragments containing 2, 3, 4 or 5 monomers are common. However, a small number of fragments can contain 6 monomers. An example is sketched in Figure 3.

Figure 3. An example of a fragment of the water structure that contains six monomers. The individual water molecules are heavily connected in liquid water; many monomers are simultaneously hydrogen bonded to 4 other water molecules. This high degree of interconnectedness results in about four times as many fragments as there are monomers. In a more open configuration, as in a simple chain, the number of fragments is about twice the number of groups. Moreover, the computer time required to evaluate the composition of the fragments is much greater for a heavily connected configuration than for a more open structure. 3.2. Modifying SMFA fragments In the context of accurate quantum chemistry calculations for molecules, the evaluation of the molecular fragments using SMFA is fast, in the sense that the computer time required is tiny by comparison with the time for the ab initio calculations. However, as discussed in Section 1, every time step in a MD simulation


16

Page 17 of 55

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


must be executed in very little computer time if 106 - 108 steps or more are required. In this context, the complete SMFA procedure is very slow. Fortunately, the structure of a molecule-solvent system does not change much in a time of only 1 fsec. Hence, the compositions of the fragments do not change much in one MD time step. As we shall see, almost all of the fragments retain the same composition from one step to another. Hence, all we need to do is slightly modify the set of fragments at each time step. To modify the fragmentation, we use the fact that no fragment can contain groups that are separated by more than Level bonds. Recall that for SMFA a group was defined in section 3.1 to correspond to the usual definition of a chemical functional group in a covalent molecule; such groups are connected with other groups by covalent bonds, or by hydrogen bonds. In the case of liquid water, each H2O molecule is a group. Let us consider two groups in the system that either form a new bond or break their bond from one time step to the next. Denote these as groups 1 and 2. Each of these groups may be connected to other groups in the system. There will be a subset of connected groups that contains these two groups and all groups that lie up to Level bonded connections away from these two groups. Denote this set of groups (and connections) as S(1,2). For example, Figure 4 depicts a portion of a moleculesolvent system, with groups denoted by balls and bonds by straight lines. The dashed line between groups 1 and 2 denotes that a bond is made or broken between these two groups from one time step to the next.


17


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 18 of 55

Figure 4. A schematic representation of a subset of groups and bonds within some molecule-solvent system.

For Level = 2, the groups A, B, C, D, E, and F would be in S(1,2), while groups g, h, i, and j are not. The groups labelled g, h, i, and j cannot be in any fragment that contain groups 1 or 2. The set of fragments that contain groups 1 and 2 must come only from the subset S(1,2). In Figure 4, S(1,2) = ABC12DEF. Conversely, the fragments containing groups g, h, i, j, and the remainder of the whole system cannot be affected by the bonding between groups 1 and 2. Now suppose that groups 1 and 2 are bonded, then to denote this connection explicitly, we write S(1,2) as S(1,2;


18

Page 19 of 55

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


bonded). Alternatively, if no bond exists between groups 1 and 2, we denote S(1,2) as S(1,2; not bonded). Suppose now that, at some time, groups 1 and 2 are bonded and that the fragmentation, for Level = 2, is given by Eq. (3.3). Now we further suppose that at the next time step, the bond between groups 1 and 2 no longer exits (the distance between the groups is then longer than the allowed maximum). Then the change in the fragmentation caused by breaking this bond is ∆M:

∆M = S(1,2; not bonded) - S(1,2; bonded) .

(3.8)

Note that while ∆M contains no net groups, the two terms in ∆M differ in their connections. The combination of connected groups in ∆M can be fragmented using SMFA to obtain a change in fragments, ∆F:

∆M → ∆F

(3.9)

Note that ∆M contains only a small set of connected groups, independent of the size of the whole system. Since the number of groups is small, the computer time required in Eq. (3.9) is very short. Then the modified fragmentation is given by N frag

M → ∑ fn Fn + ∆F.

(3.10)

n=1

For the example of Figure 4, using Level = 2:

∆F = F2D + 2DE + 2D + 12 + A1 - 12DF - 2DE - A12 - DE

(3.11)

Generally, some fragments contained in ∆F will cancel with fragments of the same composition but opposite sign in the summation in the original sum of fragments in Eq. (3.10). For example, the terms -12DF and -A12 in Eq. (3.11) will cancel corresponding terms in the previous fragmentation, which occurred for Level = 2 because groups 1 and 2 were bonded.


19


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 20 of 55

Similarly, if groups 1 and 2 were originally not bonded, but became bonded at the next time step, then the earlier fragmentation would be modified as in Eq, (3.10), except that all the fragment signs in ∆F would be reversed. If several bonds are broken or formed from one time step to the next, then Eqs (3.8) to (3.10) are simply evaluated for the several pairs of atoms involved in the bond formation or breaking. Periodic boundary conditions For crystals or fluids under periodic boundary conditions, fragments contain groups from both the "central unit cell" or simulation box and from neighbouring cells. Each group is identified by a vector of 4 quantities (n,l1,l2,l3), where n denotes the number of the group within a unit cell, and l1, l2, l3 are lattice indices that enumerate the unit cells in the whole system. In a liquid simulation, the cubic unit cell normally contains hundreds of groups and has a side of tens of Angstrom in length. Fragments must contain at least one group in the central unit cell. Hence, for typical values of the parameter Level (2, 3, 4), l1, l2, l3 are restricted to values between -1 and 1. The groups in the fragmentation expression (3.3) are expressed in this notion. After Eq. (3.9) is evaluated, the fragments in ∆F are expressed in this lattice notion. Hence, when evaluating Eq. (3.10), some fragments in ∆F may cancel with fragments in the original summation because they are equivalent (with opposite sign) to within a lattice translation. 3.3. Liquid simulation We can now set out an algorithm or procedure for evaluating the forces (and energy) in a molecule-solvent system during a direct ab initio MD simulation. (1)

We assume that the fragmentation of the whole system has been carried out

before the MD simulation begins. The identities of the atoms that comprise each


20

Page 21 of 55

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


group are stored for later use. The composition of the fragments is stored in the (n,l1,l2,l3) group notation. We denote the set of fragments as {F(t)}. Moreover, we also store the current bonding of each group: For each group we record a list of groups [in (n,l1,l2,l3) notation] to which it is bonded. We denote this list of bonded connections as B(t). The set of elemental symbols and Cartesian coordinates of each fragment is stored, along with the total charge of the fragment and its spin multiplicity. (2)

Ab initio calculations are carried out for each fragment. The forces on each

atom are evaluated and combined with the forces from all calculations to give an approximation to the force on each atom in the whole system. The Cartesian coordinates of each atom are updated, according to the MD algorithm (say the velocity-Verlet algorithm). (3)

Some groups may move out of the box, when the coordinates are updated. The

usual MD algorithm is employed to translate such groups to within the opposite side of the box. Under normal conditions of temperature, pressure and size of the simulation box, the number of groups leaving the box is expected to be proportional to the surface area of the box, and hence to increase with increasing system size at less than a linear rate. (4)

For each such translation in (3), the lattice notation for each group (n,l1,l2,l3) is

updated on the assumption that fragment compositions and the lists of bonded groups are otherwise unaltered. For example, if group (n,0,0,0) leaves the box in the forward x axis direction, then the Cartesian coordinates of (n,0,0,0) are translated close to where (n,-1,0,0) was at the last time step. Hence, the list of groups to which (n,0,0,0) is bonded is modified to reduce the first lattice component of each neighbouring group by 1, as that is the list of bonded neighbours of (n,-1,0,0). Moreover, where


21


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 22 of 55

(n,l1,l2,l3) occurs in a list of neighbours of some other group, it is replaced by (n,l1+1,l2,l3). The fragment compositions are similarly modified: If (n,l1,l2,l3) occurs in a fragment, it is replaced by (n,l1+1,l2,l3). Corresponding modifications of the lists of bonded neighbours and fragment compositions are made if a group leaves the box backwards along the x axis, or along the y or z axes in either direction. Again, the computer time required for this is expected to increase less than linearly with increasing system size. (5)

As in many MD approaches, the simulation box is decomposed into a three

dimensional grid of "gridboxes".75 Here the grid box length is taken to be just longer than the maximum possible bond length so that each group can only be bonded to groups in its own gridbox or to groups in the immediately adjacent gridboxes. For the box of 270 water molecules of Section II, we use a 7 x 7 x 7 grid (343 gridboxes). Each group is assigned to a grid box, based on its location. This requires a computation time that is linear in the size of the system, denoted here as O(N), where N is the total number of groups. In simulations of very large systems, it may be an advantage to also assign each fragment to a grid box, based on some measure of the fragment's average position. The utility of this step will be discussed below. Again, since the number of fragments is only O(N), the computer time for this procedure is only O(N). (6)

For each atom in a group, the atom-atom distances are calculated to all atoms

in other groups that are within the same gridbox or nearest-neighbour gridboxes (only). For the whole system, the computational cost is only O(N). For the simulation of water, atom-atom distances are only calculated if the two atoms are hydrogen and oxygen. If the calculated distance is below the maximum allowed for a bond (including hydrogen bonds), then a bond is assigned, and the list of groups bonded to


22

Page 23 of 55

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


each group is stored for this configuration of the system. We denote this list as B(t+δt). (7)

Comparing the old list of bonding connections, B(t), with the new list,

B(t+δt), identifies bonds that have been broken or formed. Comparison is made between the two lists for each group, one group at a time, so the computation time is also only O(N). (8)

Create a difference of fragments, ∆M, as in Eq. (3.8) for each broken bond,

and the corresponding ∆M, with opposite sign, for each newly formed bond. The number of bonds broken or formed in one time step should be proportional to the size of the system. Hence, the computation time for creating all the ∆M terms should only be O(N). As we shall see below, even when the bonds are only hydrogen bonds, as in water, the total number of bonds formed and broken in one time step is very small compared to the total number of bonds. Covalent bonds are only formed or broken in chemical reactions, so that such events would be rare over time and very few in number at any time. (9)

Using SMFA, fragment the ∆M to ∆F, as in Eq. (3.9). Since the number of

sets of groups in ∆M for the whole system is only O(N) for a solvent (or independent of N for a covalent bond break/formation), the computer time for this process is also only O(N). (10)

Express each fragment in ∆F in the periodic (n,l1,l2,l3) group notation. The

number of fragments in ∆F is O(N), so this process requires only O(N) computation time. (11)

Combine ∆F with the set of fragments {F(t)} of step (1) above, as in Eq.

(3.10). Since both ∆F and {F(t)} contain O(N) fragments, the process of searching for cancelling fragments will require O(N2) computation time, although the "prefactor" is


23


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 24 of 55

small since the number of fragments in ∆F is small. However, if we have assigned the fragments in {F(t)} to the gridboxes, as mentioned in (5) above, and similarly assigned the fragments in ∆F , then the computation time should be only O(N), as each fragment in ∆F need only be compared with fragments in {F(t)} that are in the same gridbox or in nearby gridboxes. The result of this process is a new fragmentation, which we denote as {F(t+δt)}. (12)

We replace {F(t)} by {F(t+δt)} and replace the list of bonded groups B(t) by

B(t+δt), and return to step (2) to continue the trajectory.

4. Results We have implemented the procedure of section 3.3 using the sample of MD configurations of 270 water molecules under periodic boundary conditions, described in section 2. That is, we have not evaluated the forces on each fragment ab intio at each time step and updated the coordinates accordingly, but simply taken the model MD configurations at each time step. These configurations result from the MD force field implemented in the Gromacs program (as detailed in section 2) and recorded at 1 fsec intervals. In this first investigation, we have also not assigned the fragments in {F(t)} and ∆F to a grid, leaving a small sub-process which is formally O(N2) in computation time. The purpose of this model investigation is to highlight the scale of some important processes and determine if this application of SMFA is feasible. Some important scales can be identified. (i) Over 10,000 time steps, only 1163 water molecules left or entered the box from one time step to the next, or one molecule every 9 steps on average. This suggests that the translation of the group periodic coordinates, in the fragment and bonding lists, might only be updated approximately every 9 steps to minimise computation time.


24

Page 25 of 55

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


(ii) On average each configuration of the 270 water molecules has about 1100 bonds which involve at least one monomer in the central unit cell. At each time step, some bonds are broken and some new bonds are formed. Figure 5 presents the distribution of such bond changes at each time step over a sample of 1000 steps.

Figure 5. A histogram of the total number of bonding changes, bonds broken plus bonds formed, at each time step of 1 fsec during an interval of 1000 timesteps.

On average, 1.6 hydrogen bonds are broken at each step, and 1.6 new bonds are formed, giving 3.2 (± 1.8 standard deviation) bonding changes at each step. That is, only about 0.3% of bonds change at each step. This justifies the assumption that only a small modification of the fragmentation is required at each time step. As one might expect for a sample of configurations near equilibrium, the number of bonds breaking is very close to the number of bonds forming over time. (iii) The making and breaking of bonds results in a fluctuation in the total number of fragments over time. Figure 6 presents the variation of the number of fragments at each timestep, over a period of 1 psec. Figure 7 presents a histogram of the number of


25


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 26 of 55

fragments at each time step. The mean number of fragments in this time interval is 1023, with a standard deviation of 34.

Figure 6. The total number of fragments, at each timestep, is shown as a function of time over 1 psec.


26

Page 27 of 55

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


Figure 7. A histogram of the number of fragments at a single time step, accumulated over 1000 timesteps.

(iv) The procedure to update the fragmentation takes on average 0.12 sec on an Intel Xeon 2.6 GHz processor, with the current code. This time might well be reduced very substantially with more optimised code than has been achieved in this first implementation. Very little floating point arithmetic is actually involved, outside the calculation of atom-atom distances. Nevertheless, by reducing the time required for fragmentation from 101 - 103 secs to just 0.12 sec per timestep, the algorithm


27


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 28 of 55

presented herein ensures that SMFA can be feasibly applied to direct ab initio MD simulations. The major remaining computational task is then the ab initio calculations of the forces in the relatively small SMFA fragments. (v) A major advantage of the fragmentation approach is that the necessary ab initio calculations are independent and could be executed using different processors. In the example used herein, the number of independent fragments is about 1000. This does not imply that 1000 processors are required for optimum efficiency and reduction of the quantum chemistry calculation time. The optimum number of processors might be determined as follows. The time required for ab initio calculations at each time step is set by the time required to evaluate the forces for the largest fragment. The time required for the ab initio calculation of each fragment will depend on the parallelisation that can be achieved, for the particular ab initio method, on multiple cores within an individual processor and the parallelisation possible across multiple processors. This time will be a complicated function of the number of basis functions. However, we can assume that the fragment with the largest number of basis functions will require the longest calculation time. The remaining (smaller) fragments can be expected to require less processing time. It should be possible to estimate, a priori, the processing time required for each fragment ab initio calculation, taking all the factors above into account. Once the calculation time for each fragment has been estimated, there is an important "coarse-grained" parallelisation that follows from the characteristics of the SMFA fragment size distribution: the whole set of fragment calculations must be distributed over "processors" (in practice multiple cores and possibly multiple processors) to "balance the load" of the calculations and minimise the actual time required for the whole set of calculations. To illustrate this point, we have used a very


28

Page 29 of 55

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


simple model, assuming that the ab initio processing time for a fragment is proportional to the number of water molecules in the fragment, NH2O, raised to a power, say NH2O 4. A simple procedure to distribute the computation "load" might then be as follows. First, allocate the largest fragment to one processor. Then, taking the remaining fragments in order of increasing processing time, we allocate fragments to a second processor until the total processing time for all such fragments would exceed the time for the first fragment. At that point, we begin allocating fragments to a third processor, and so on. Any such simple algorithm should require negligible time to execute. At each time step, the number of processors required will vary. For the 1000 time steps in Figure 6, Figure 8 presents a histogram of the number of processors required at each step. In practice each "processor" may represent multiple cores in an individual processor unit, or even multiple processing units. The average number of processors in the figure is 98. As Fig. 8 shows, the number of processors required does vary considerably, but is much less that the total number of independent ab initio calculations, that is, the number of fragments.


29


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 30 of 55

Figure 8. A histogram of the number of processors required for ab initio calculations at a single time step, accumulated over 1000 time steps. 5. Discussion The SMFA method has been implemented in the context of ab initio molecular dynamics, under periodic boundary conditions. The test case used to demonstrate the approach was liquid water, simulated with 270 water molecules at a temperature of 300 K and a pressure of 1 atmosphere, under periodic boundary conditions. This test case was chosen, rather than a covalently bonded molecule, because the hydrogen bonds between water molecules are very much more labile than covalent bonds.


30

Page 31 of 55

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


Hence, the set of fragmentations, which represent the solvent in SMFA, change at virtually every 1 fsec time step in a simulation, whereas covalent bonding only changes in the rare event of a chemical reaction. Hence, the most difficult task for SMFA (or any fragmentation method) is to efficiently update, at each time step, the set of fragments that represent the solvent in a solvent-molecule system. It is interesting to note that the high connectivity of each water molecule is responsible for the long computer time required to apply the original SMFA procedure to the fragmentation of hundreds or thousands of water molecules. This same high connectivity is responsible for the fact that only a few connections (hydrogen bonds) are made or broken in a time as short as 1 fsec. The algorithm for updating the fragmentation then relies for its efficiency on this fact that few bonds are made or broken in one time step. As a consequence, an average of only 0.12 sec was required for the fragmentation up-dating scheme at each time-step. While this time is likely to be two orders of magnitude less than the ab initio computation time at each step, more efficient programing might reduce this time further. The test case employed herein shows that only incremental changes in the set of fragments occur at each time step, which ensures that the modification procedure of section 3.2 is an efficient approach. However, it is also clear from Figure 6 that over time the set of fragments changes significantly. Hence, in order to maintain an accurate estimate of the energy and forces, any fragmentation scheme must modify the composition of the fragments frequently, probably as often as every time step. A major advantage of fragmentation schemes is that all ab initio calculations are independent, so that the time required for ab initio calculations is independent of the size of the system, given enough processors. However, the number of processors required is much less than the total number of fragments, because the computational


31


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 32 of 55

time for the ab initio calculations is determined by the size of the largest fragment. A simple ansatz for how the computation time varies with fragment size was used to illustrate how the number of required processors might be determined. Note that if we had considered a much larger simulation box, the number of processors required would simply increase linearly with the number of water molecules. In previous applications of SMFA to molecules and water clusters, the interactions between water monomers or molecular groups that are well separated in distance, so-called "non-bonded" interactions, have been accounted for using both ab initio calculations and perturbation theory. Non-bonded interactions within a "cut-off" distance are evaluated using ab initio calculations. Interactions at longer distances are either treated using perturbation theory, or by evaluating the ab initio energy of the SMFA fragments in the presence of point charges that model the electrostatic and induction interactions with these long range parts of the whole system. Using such embedded charges ensures that SMFA with Level = 2 gives cluster binding energies within about 1%.73 In the context of direct dynamics, these point charges, associated with solvent molecules or molecular fragments at a distance, might be approximated using the equilibrium charge distribution of these molecules and fragments. This is similar to the approach adopted in traditional molecular dynamics. There would then be negligible computational time required to calculate the embedded charges for each fragment ab initio calculation at each time step. In the case of non-bonded interactions between solvent molecules within the cut-off distance, a simple approximation should suffice. Since SMFA fragments for only Level = 2 accurately account for the relatively strong, very close, cooperative, solvent-solvent interactions, it should be adequate to treat longer range effects using pairwise interactions. There would be negligible computational time required to identify changes in these pairwise


32

Page 33 of 55

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


interactions. These interactions are identified directly from the calculation of inter monomer distances in section 3.3 (5). Although this paper has considered a single "whole" system of water under periodic boundary conditions, it should be possible to treat a QM/MM or ONIOM

76

approach in the same way. That is, one might consider a very large molecule-solvent system (under periodic boundary conditions, as usual), but only treat a relatively small part of that system using direct ab initio molecular dynamics, while the remainder of the system employed standard molecular dynamics force fields. It is important to note that an efficient implementation of direct dynamics using SMFA would require that the SMFA procedure, the ab initio calculations code (say GAMESS) 77, and the molecular dynamics simulation codes be integrated in one package. This would avoid time consuming input and output of the atomic coordinates and velocities, the atomic coordinates and electronic structure details for each fragment ab initio calculation, and other data needed for the SMFA procedure. The amalgamation of quantum chemistry and quasiclassical trajectory codes has already been achieved for direct ab initio reaction dynamics in the gas phase. 78 It is also important to note that as the set of SMFA fragments and any set of longer range pairwise interactions changes with time, the net approximation to the forces on each atom in the system will change discontinuously. In earlier work, it was found that the error in the forces due to the fragmentation approximation can be very small, relative to the usual accuracy employed in ab initio calculations. For example, the error in the forces may be only of the order of magnitude that is deemed negligible in determining stationary geometries.

71

However, important questions must be

answered as to whether a fragmentation approximation to a PES provides the basis for a sufficiently accurate estimate of observable properties with classical trajectory


33


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 34 of 55

methods. These properties include "static" quantities, such as free energies, and "dynamic" properties, such as diffusion coefficients or reaction rates. The fragmentation approximation may possibly introduce systematic errors in the PES, and most likely introduces random errors (including discontinuities) in the PES. The magnitude of these errors (both systematic and random) can be reduced by using higher values of the parameter Level or by more accurately accounting for long-range interactions (including the use of embedded charges). Systematic evaluation of the convergence of property values with the convergence of the magnitude of errors in the forces will be necessary. These important questions must be left to later investigations. However, the algorithm presented herein now makes such investigations feasible.

Supporting Information: (i) The covalent radii used; (ii) the definition of single bonds; (iii) the definition of multiple bonds, and (iv) the Van der Waals radii used.

Acknowledgements The author gratefully acknowledges the assistance of his colleague Dr Megan O'Mara for supplying the water simulation data, and gratefully acknowledges helpful discussions with his colleague Dr Stephen Williams. The author also gratefully acknowledges a grant of computer time from the Australian National Computational Infrastructure National Facility.


34

Page 35 of 55

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)

Corchado, J. C.; Vhuang, Y.-Y.; Fast, P. L.; Hu, W.-P.; Liu, Y.-P.; Lynch, G.

C.; Nguyen, K. A.; Jackels, C. F.; Ramos, A. F.; Ellingson, B. A.; et al, POLYRATE 9.7: Computer Program for the Calculation of Chemical Reaction Rates for Polyatomics, 2007. (2)

Althorpe, S. C.; Clary, D. C. Quantum Scattering Calculations on Chemical

Reactions. Annu. Rev. Phys. Chem. 2003, 54, 493-529. (3)

Hu, W.; Schatz, G. C. Theories of Reactive Scattering. J. Chem Phys. 2006,

125, 132301. (4)

Zhang, W.; Zhou, Y.; Wu, G.; Lu, Y.; Pan, H.; Fu, B.; Shuai, Q.; Liu, L.; Liu,

S.; Zhang, L.; et al Depression of Reactivity by the Collision Energy in the Single Barrier H + CD4 → HD + CD3 Reaction. PNAS 2010, 107, 12782 -12785. (5)

Zhou, Y.; Zhang, D. H. Eight-Dimensional Quantum Reaction Rate

Calculations for the H+CH4 and H2+CH3 Reactions on Recent Potential Energy Surfaces. J. Chem. Phys. 2014, 141, 194307. (6)

Porter, R. N. Molecular Trajectory Calculations. Annu. Rev. Phys. Chem.

1974, 25, 317-355. (7)

Dynamics of Molecular Collisions, Miller, W. H., Ed.; Plenum Press: New

York, 1976. (8)

Heazlewood, B. R.; Jordan, M. J. T.; Kable, S. H.; Selby, T. M.; Osborn, D.

L.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Roaming is the Dominant Mechanism for Molecular Products in Acetaldehyde Photodissociation. Proc. Nat. Acad. Sci. 2008, 105, 12719-12724. (9)

Collins, M. A. Molecular Potential-Energy Surfaces for Chemical

Reaction Dynamics. Theor. Chem. Acc. 2002, 108, 313-324.


35


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

(10)

Page 36 of 55

Braams, B. J.; Bowman, J. M. Permutationally Invariant Potential

Energy Surfaces in High Dimensionality. Int. Rev. Phys. Chem. 2009, 28, 577-606. (11)

Xu, X.; Chen, J.; Zhang, D. H. Global Potential Energy Surface for the

H + CH4 ↔ H2 + CH3 Reaction using Neural Networks. Chin. J. Chem. Phys. 2014, 27, 373-379. (12)

Li, J.; Guo, H. An Accurate Full 15 Dimensional Permutationally

Invariant Potential Energy Surface for the OH + CH4 → H2O + CH3 Reaction. J. Chem Phys. 2015, 143, 221103. (13)

Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.;

Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. M.; Kollman, P. A. AMBER. J. Am. Chem. Soc. 1995, 117, 5179-5197. (14)

Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D.; Nilsson, L.;

Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al CHARMM: The Biomolecular Simulation Program,. J. Comp. Chem. 2009, 30, 15451615. (15)

Spoel, D. V. d.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.;

Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701-1718. (16)

Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions:

Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227-249. (17)

Warshel, A. Molecular Dynamics Simulations of Biological Reactions.

Acc. Chem. Res. 2002, 35, 385-395. (18)

Svensson, M.; Humbel, S.; Froese, R. D. J.; Toshiaki Matsubara;

Sieber, S.; Morokuma, K. ONIOM: A Multilayered Integrated MO + MM Method


36

Page 37 of 55

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


for Geometry Optimizations and Single Point Energy Predictions. A Test for Diels−Alder Reactions and Pt(P(t-Bu)3)2 + H2 Oxidative Addition. J. Phys. Chem. 1996, 100, 19357-19363. (19)

Vreven, T.; Morokuma, K.; Farkas, O.; Schlegel, H. B.; Frisch, H. B.

Geometry Optimization with QM/MM, ONIOM, and other Combined Methods. I. Microiterations and Constraints. J. Comput. Chem. 2003, 24, 760-769. (20)

Canfield, P.; Dahlbom, M. G.; Hush, N. S.; Riemers, J. R. Density-

Functional Geometry Optimization of the 150 000-Atom Photosystem-I Trimer. J. Chem. Phys. 2006, 124, 024301. (21)

Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko,

L. V.; Boatz, J. A. Accurate Methods for Large Molecular Systems. J. Phys. Chem. B 2009, 113, 9646-9663. (22)

Mullin, J. M.; Roskop, L. B.; Pruitt, S. R.; Collins, M. A.; Gordon, M.

S. Systematic Fragmentation Method and the Effective Fragment Potential: An Efficient Method for Capturing Molecular Energies. Journal of Physical Chemistry A 2009, 113, 10040-10049. (23)

Gao, J. Towards a Molecular Orbital Derived Empirical Potential for

Liquid Simulations. J. Phys. Chem. B 1997, 101, 657-663. (24)

Gao, J. A Molecular-Orbital Derived Polarization Potential for Liquid

Water. J. Chem. Phys. 1998, 109, 2346-2354. (25)

Xie, W.; Gao, J. Design of a New Generation Force Field: The X-POL

Potential. J. Chem. Theory Comput. 2007, 3, 1890-1900. (26)

Xie, W.; Orozco, M.; Truhlar, D. G.; Gao, J. X-POL Potential: an

Electronic Structure-Based Force Field for MD Simulation of a Solvated Protein in Water. J. Chem. Theory Comput. 2009, 5, 459-467.


37


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

(27)

Page 38 of 55

Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L.

Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632-672. (28)

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. (29)

Fedorov, D. G.; Kitaura, K. The Fragment Molecular Orbital Method:

Practical Applications to Large Molecular Systems; CRC Press: Boca Raton, 2009. (30)

Zhang, D. W.; Zhang, J. Z. H. Molecular Fractionation with Conjugate

Caps for Full Quantum Mechanical Calculation of Protein-Molecule Interaction Energy. Journal of Chemical Physics 2003, 119, 3599-3605. (31)

Zhang, D. W.; Zhang, J. Z. H. Full Ab Initio Computation of Protein-

Water Interaction Energies. Journal of Theoretical & Computational Chemistry 2004, 3, 43-49. (32)

Mei, Y.; Ji, C. G.; Zhang, J. Z. H. A New Quantum Method for

Electrostatic Solvation Energy of Protein. J. Chem. Phy. 2006, 125, 094906. (33)

Duan, L. L.; Mei, Y.; Zhang, Q. G.; Zhang, J. Z. H. Intra-Protein

Hydrogen Bonding is Dynamically Stabilized by Electronic Polarization. J. Chem. Phys. 2009, 130, 115102. (34)

Deev, V.; Collins, M. A. Approximate Ab Initio Energies by

Systematic Molecular Fragmentation. J. Chem. Phys. 2005, 122, 154102. (35)

Collins, M. A.; Deev, V. A. Accuracy and Efficiency of Electronic

Energies from Systematic Molecular Fragmentation. J. Chem. Phys. 2006, 125, 104104.


38

Page 39 of 55

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


(36)

Addicoat, M. A.; Collins, M. A. Accurate Treatment of Non-Bonded

Interactions within Systematic Molecular Fragmentation. J. Chem. Phys. 2009, 131, 104103. (37)

Bettens, R. P. A.; Lee, A. M. A New Algorithm for Molecular

Fragmentation in Quantum Chemical Calculations. J. .Phys. Chem. A 2006, 110, 8777-8785. (38)

Bettens, R. P. A.; Lee, A. M. On the Accurate Reproduction of Ab

Initio Interaction Energies between an Enzyme and Substrate. Chem. Phys. Lett. 2007, 449, 341-346. (39)

Le, H. A.; Lee, A. M.; Bettens, R. P. A. Accurately Reproducing Ab

Initio Electrostatic Potentials with Multipoles and Fragmentation. J. Phys. Chem. A 2009, 113, 10527-10533. (40)

Le, H. A.; Tan, H. J.; Ouyang, J. F.; Bettens, R. P. A. Combined

Fragmentation Method: A Simple Method for Fragmentation of Large Molecules. J. Chem. Theor. Comp. 2012, 8, 469-478. (41)

Ganesh, V.; Dongare, R. K.; Balanarayan, P.; Gadre, S. R. Molecular

Tailoring Approach for Geometry Optimization of Large Molecules: Energy Evaluation and Parallelization Strategies. J. Chem. Phys. 2006, 125, 104109. (42)

Sahu, N.; Yeole, S. D.; Gadre, S. R. Appraisal of Molecular Tailoring

Approach for Large Clusters. J. Chem. Phys. 2013, 138, 104101. (43)

Li, S. H.; Li, W.; Fang, T. An Efficient Fragment-Based Approach for

Predicting the Ground-State Energies and Structures of Large Molecules. J. Am. Chem. Soc. 2005, 127, 7215-7226.


39


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

(44)

Page 40 of 55

Li, W.; Li, S. H.; Jiang, Y. S. Generalized Energy-Based

Fragmentation Approach for Computing the Ground-State Energies and Properties of Large Molecules. J. Phys. Chem. A 2007, 111, 2193-2199. (45)

Hua, S. G.; Hua, W. J.; Li, S. H. An Efficient Implementation of the

Generalized Energy-Based Fragmentation Approach for General Large Molecules. J. Phys. Chem. A 2010, 114, 8126-8134. (46)

Jiang, N.; Ma, J.; Jiang, Y. Electrostatic Field-Adapted Molecular

Fractionation with Conjugated Caps for Energy Calculations of Charged Biomolecules. J. Chem. Phys. 2006, 124, 114112. (47)

Ramabhadran, R. O.; Raghavachari, K. Theoretical Thermochemistry

for Organic Molecules: Development of the Generalized Connectivity-Based Hierarchy. J. Chem. Theor. Comp. 2011, 7, 2094-2103. (48)

Mayhall, N. J.; Raghavachari, K. Molecules-in-Molecules: An

Extrapolated Fragment-Based Approach for Accurate Calculations on Large Molecules and Materials. J. Chem. Theor. Comp. 2011, 7, 1336-1343. (49)

Mayhall, N. J.; Raghavachari, K. Many-Overlapping-Body (MOB)

Expansion: A Generalized Many Body Expansion for Nondisjoint Monomers in Molecular Fragmentation Calculations of Covalent Molecules. J. Chem. Theor. Comp. 2012, 8, 2669-2675. (50)

Richard, R. M.; Herbert, J. M. A Generalized Many-Body Expansion

and a Unified View of Fragment-Based Methods in Electronic Structure Theory. J. Chem. Phys. 2012, 137, 064113. (51)

Collins, M. A.; Bettens, R. P. A. Energy-Based Molecular

Fragmentation Methods. Chem. Rev. 2015, 115, 5607-5642.


40

Page 41 of 55

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


(52)

Collins, M. A. Systematic Fragmentation of Large Molecules by

Annihilation. Phys. Chem. Chem. Phy. 2012, 14, 7744-7751. (53)

Collins, M. A.; Cvitkovic, M. W.; Bettens, R. P. A. The Combined

Fragmentation and Systematic Molecular Fragmentation Methods. Acc. Chem. Res. 2014, 47, 2776-2785. (54)

Collins, M. A. Molecular Potential Energy Surfaces Constructed from

Interpolation of Systematic Fragment Surfaces. J. Chem. Phys. 2007, 127, 024104. (55)

Frankcombe, T. J.; Collins, M. A. Growing Fragmented Potentials for

Gas-Surface Reactions: The Reaction between Hydrogen Atoms and HydrogenTerminated Silicon (111). J. Phys. Chem. C 2012, 116, 7793-7802. (56)

D’Arcy, J. H.; Jordan, M. J. T.; Frankcombe, T. J.; Collins, M. A. H2

Adsorption in a Porous Crystal: Accurate First-Principles Quantum Simulation. J. Phys. Chem. A 2015, 119, 12166-12181. (57)

Leforestier, C. Classical Trajectory Studies Using the Full Ab Initio

PES H- + CH4 → CH4 + H-. J. Chem. Phys. 1978, 68, 4406-4410. (58)

Chen, W.; Hase, W. L.; Schlegel, H. B. Ab Initio Classical Trajectory

Study of H2CO → H2 + CO Dissociation. Chem. Phys. Lett. 1994, 228, 436-442. (59)

Millam, J. M.; Bakken, V.; Chen, W.; Hase, W. L.; Schlegel, H. B. Ab

Initio Classical Trajectories on the Born–Oppenheimer Surface: Hessian-Based Integrators using Fifth-Order Polynomial and Rational Function Fits. J. Chem. Phys. 1999, 111, 3800-3805. (60)

Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and

Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471-2474. (61)

Garton, D. J.; Minton, T. K.; Troya, D.; Pascual, R.; Schatz, G. C.

Hyperthermal Reactions of O(3P) with Alkanes: Observations of Novel Reaction


41


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 42 of 55

Pathways in Crossed-Beams and Theoretical Studies. J. Phys. Chem. A 2003, 107, 4583-4587. (62)

Sun, L.; Hase, W. L. Born-Oppenheimer Direct Dynamics Classical

Trajectory Simulations. In Reviews in Computational Chemistry; Lipkowitz, K. B., Larter, R., Cundari, T. R., Eds.; Wiley: New York, 2003; Vol. 19; pp 79-146. (63)

Xie, J.; McClellan, M.; Sun, R.; Kohale, S. C.; Govind, N.; Hase, W.

L. Direct Dynamics Simulation of Dissociation of the [CH3--I--OH]- Ion-Molecule Complex. J. Phys. Chem. A 2015, 119, 817-825. (64)

Nishizawa, H.; Nishimura, Y.; Kobayashi, M.; Irle, S.; Nakai, H. Three

Pillars for Achieving Quantum Mechanical Molecular Dynamics Simulations of Huge Systems: Divide-and-Conquer, Density-functional Tight-Binding, and Massively Parallel Computation. J. Comp. Chem. 2016, 37, 1983-1992. (65)

Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F. v.; Hermans, J.

Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel: Dordrecht, 1981; pp 331-342. (66)

Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An Nlog(N)

Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. (67)

Allen, F. H.; Kennard, O.; Watson, D. G.; Brammer, L.; Orpen, A. G.

Tables of Bond Lengths determined by X-Ray and Neutron Diffraction. Part I . Bond Lengths in Organic Compounds. J. Chem. Soc. Perkin Trans. 1987, S1-S19. (68)

Nixon, J. F. The Coordination Chemistry of Compounds Containing

Phosphorus-Carbon Multiple Bonds. Chem. Rev. 1988, 88, 1327-1362. (69)

Netzloff, H. M.; Collins, M. A. Ab Initio Energies of Non-Conducting

Crystals by Systematic Fragmentation. J. Chem. Phys. 2007, 127, 134113.


42

Page 43 of 55

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


(70)

Collins, M. A. Ab Initio Lattice Dynamics of Nonconducting Crystals

by Systematic Fragmentation. J. Chem. Phys. 2011, 134, 164110. (71)

Collins, M. A. Molecular Forces, Geometries and Frequencies by

Systematic Molecular Fragmentation Including Embedded Charges. J. Chem Phys. 2014, 141, 094108. (72)

Reid, D. M.; Collins, M. A. Calculating Nuclear Magnetic Resonance

Shieldings Using Systematic Molecular Fragmentation by Annihilation. Phys. Chem. Chem. Phys. 2015, 17, 5314-5320. (73)

Pruitt, S. R.; Addicoat, M. A.; Collins, M. A.; Gordon, M. S. The

Fragment Molecular Orbital and Systematic Molecular Fragmentation Methods Applied to Water Clusters. Phys. Chem. Chem. Phys. 2012, 14, 7752-7764. (74)

Frankcombe, T. J.; Collins, M. A. Potential Energy Surfaces for Gas-

Surface Reactions. Phys. Chem. Chem. Phys. 2011, 13, 8379-8391. (75)

Quentrec, B.; Brot, C. New Method for Searching for Neighbours in

Molecular Dynamics Computations. J. Comput. Phys. 1975, 13, 430-432. (76)

Vreven, T.; Byun, K. S.; Komáromi, I.; Dapprich, S.; J. A.

Montgomery, J.; Morokuma, K.; Frisch, H. B. Combining Quantum Mechanics Methods with Molecular Mechanics Methods in ONIOM. J. Chem. Theory & Comput. 2006, 2, 815-826. (77)

Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure

Theory: GAMESS a Decade Later. In Theory and Applications of Computational Chemistry, the First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G., Eds.; Elsevier: Amsterdam, 2005, pp 1167-1189. (78)

Lourderaja, U.; Sun, R.; Kohale, S. C.; Barnes, G. L.; Jong, W. A. d.;

Windus, T. L.; Hase, W. L. The VENUS/NWChem Software Package. Tight


43


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 44 of 55

Coupling between Chemical Dynamics Simulations and Electronic Structure Theory. Comput. Phys. Commun. 2014, 185, 1074−1080.

Email Address: [email protected] Telephone number: 61 2 6125 3254


44

Page 45 of 55


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 ACS Paragon Plus Environment

45


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 46 of 55

Figure Captions Figure 1. A figurative representation of the positions of all 270 water molecules at one time step in the simulation.

Figure 2. The distribution of fragment sizes in 1118 fragments that represent the 270 water molecules at one (arbitrary) configuration in the simulation.

Figure 3. An example of a fragment of the water structure that contains six monomers.

Figure 4. A schematic representation of a subset of groups and bonds within some molecule-solvent system.

Figure 5. A histogram of the total number of bonding changes, bonds broken plus bonds formed, at each time step of 1 fsec during an interval of 1000 timesteps.

Figure 6. The total number of fragments, at each timestep, is shown as a function of time over 1 psec.

Figure 7. A histogram of the number of fragments at a single time step, accumulated over 1000 timesteps.

Figure 8. A histogram of the number of processors required for ab initio calculations at a single time step, accumulated over 1000 time steps.


46

Page 47 of 55


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 ACS Paragon Plus Environment

47


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

Figure 1. A figurative representation of the positions of all 270 water molecules at one time step in the simulation. 67x61mm (72 x 72 DPI)


Page 48 of 55

Page 49 of 55

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


82x156mm (72 x 72 DPI)



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

Figure 3. An example of a fragment of the water structure that contains six monomers. 82x60mm (281 x 281 DPI)


Page 50 of 55

Page 51 of 55

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


Figure 4. A schematic representation of a subset of groups and bonds within some molecule-solvent system. 82x147mm (193 x 193 DPI)



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

Figure 5. A histogram of the total number of bonding changes, bonds broken plus bonds formed, at each time step of 1 fsec during an interval of 1000 timesteps. 82x82mm (72 x 72 DPI)


Page 52 of 55

Page 53 of 55

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


Figure 6. The total number of fragments, at each timestep, is shown as a function of time over 1 psec. 82x156mm (72 x 72 DPI)



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

Figure 7. A histogram of the number of fragments at a single time step, accumulated over 1000 timesteps. 82x156mm (72 x 72 DPI)


Page 54 of 55

Page 55 of 55

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


Figure 8. A histogram of the number of processors required for ab initio calculations at a single time step, accumulated over 1000 time steps. 82x156mm (72 x 72 DPI)


Can Systematic Molecular Fragmentation Be ... - ACS Publications

Recommend Documents