A Way To Speed Up Nonequilibrium Molecular ... - ACS Publications

Subscriber access provided by UNIV OSNABRUECK

Article

Elastic barrier dynamical freezing in free energy calculations: a way to speed up nonequilibrium molecular dynamics simulations by orders of magnitude Edoardo Giovannelli, Gianni Cardini, and Riccardo Chelli J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b01117 • Publication Date (Web): 15 Jan 2016 Downloaded from http://pubs.acs.org on January 19, 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.

Journal of Chemical Theory and Computation 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 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

Journal of Chemical Theory and Computation

Elastic barrier dynamical freezing in free energy calculations: a way to speed up nonequilibrium molecular dynamics simulations by orders of magnitude Edoardo Giovannelli, Gianni Cardini, and Riccardo Chelli∗ Dipartimento di Chimica, Universit` a di Firenze, Via della Lastruccia 3, I-50019 Sesto Fiorentino, Italy E-mail: [email protected]

Abstract An important issue concerning computer simulations addressed to free energy estimates via nonequilibrium work theorems, such as the Jarzynski equality [Phys. Rev. Lett. 1997, 78, 2690], is the computational effort required to achieve results with acceptable accuracy. In this respect, the dynamical freezing approach [Phys. Rev. E 2009, 80, 041124] has been shown to improve the efficiency of this kind of simulations, by blocking the dynamics of particles located outside an established mobility region. In this report, we show that dynamical freezing produces a systematic spurious decrease of the particle density inside the mobility region. As a consequence, the requirements to apply nonequilibrium work theorems are only approximately met. Starting from these considerations, we have developed a simulation scheme, called “elastic barrier dynamical freezing”, according to which a stiff potential-energy barrier is enforced at ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


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 boundaries of the mobility region, preventing the particles from leaving this region of space during the nonequilibrium trajectories. The method, tested on the calculation of the distance-dependent free energy of a dimer immersed into a Lennard-Jones fluid, provides an accuracy comparable to the conventional steered molecular dynamics, with a computational speed-up exceeding a few orders of magnitude.

1

Introduction

The free energy difference between thermodynamic states of a system represents a fundamental quantity in physics, chemistry and biology. At a coarse level, methods for free energy calculations can be classified according to their sampling conditions, which can be of equilibrium or nonequilibrium type. Concerning the second category of methodologies, two nonequilibrium work relationships, the Jarzynski Equality 1 (JE) and the Crooks Fluctuation Theorem 2,3 (CFT), are the basis for the implementation of effective algorithms. 4–12 Both theorems connect the free energy difference between two thermodynamic states to the work performed in externally driven trajectories that switch the system between such states. The external control is retained through a time-dependent parameter entering the Hamiltonian and correlated with a collective variable of the system, whose definition, albeit arbitrary, is dictated by the specific chemical-physic aspects of the process under study. The basic difference between JE and CFT is that the former exploits one set of driven trajectories leading from an established state to the other state, while in CFT an additional set of trajectories in the reverse direction of the process is also employed. JE and CFT are of practical use both in molecular dynamics and Monte Carlo simulations. 7,13–19 A typical example is represented by Steered Molecular Dynamics (SMD) simulations, 20 in which some mechanical quantity of the system, such as an interatomic distance 15,21 or a torsion angle, 13 is taken as driven collective variable. However, also more exotic nonequilibrium processes can be accounted for by using JE and CFT, in which thermodynamical variables (e.g., temperature or pressure) or potential-energy scaling factors undergo external control. 6,22 2


Page 2 of 35

Page 3 of 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


Owing to the flexibility and the proven efficacy of nonequilibrium sampling supported by JE and CFT post-analysis, a continuous research is devoted to reduce the computational cost of these methods. To this aim, several approaches have been developed, including biased path sampling, 23–25 optimal protocol strategies, 26 time-step increase, 27 premature interruption of dissipative paths 28,29 and adaptive SMD. 30,31 One of the major problems associated with these advanced nonequilibrium techniques is related to the computational effort needed to evaluate the interatomic forces, especially for large size molecular systems. In this respect, the dynamical freezing scheme 32 was proven to improve, even of several orders of magnitude, the efficiency of SMD simulations of large systems in computing the free energy as a function of collective variables. A method conceptually analogous to dynamical freezing, called configurational freezing, 33,34 was also devised for Monte Carlo simulations. The underlying idea of dynamical freezing consists in blocking the dynamical evolution of those particles involved marginally in the dissipation mechanisms of the driven trajectories. Particle freezing is realized through an isokinetic scaling of the atomic velocities accompanied by a simultaneous change of the masses. This procedure avoids time-demanding calculations of the forces between dynamically frozen particles. It is worth noting that, in general, the identification of the space domain where the driven paths take place, and hence dissipation occurs 34 (also called mobility region, hereafter MR), is not a simple task. In spite of this, it is easy to envisage potential applications in which defining the MR is straightforward; two important examples are the estimates of the binding free energy of a protein-ligand complex 14 and the adsorption free energy of a molecule-surface system. 35 The aim of this study is to analyze the systematic source of error occurring when nonequilibrium paths are produced by using the dynamical freezing scheme of Ref. 32, here referred to as Isokinetic Dynamical Freezing (IDF). In particular, we discuss the effects of the diffusion of dynamically active particles away from the MR. We find that such a diffusion, without a balancing flow of particles entering the same phase-space region (indeed, this inward particle flow is null because the particles outside the MR are frozen), prevents a correct sampling of

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

the system configurations during the nonequilibrium trajectories. Thus, an average density decrease inside the MR invalidates the nonequilibrium work theorems. As a matter of fact, JE and CFT are subject to the condition that the nonequilibrium trajectories are generated through a Markovian process, according to which an canonical phase-space sampling should be obtained upon fixing the control parameter. 36,37 In Sec. 2.2, we discuss the drop of this requirement in employing the IDF protocol. Based on the reasons of the IDF failure, in Sec. 2.3, we propose an alternative scheme, termed Elastic Barrier Dynamical Freezing (EBDF), which prevents the particles from leaving the MR. This is realized by applying a harmonic potential-energy barrier at the boundaries of the MR. In the approximation of a stiff potential-energy barrier, the mobile particles moving towards the MR boundaries, undergo elastic collisions, thus remaining confined within the MR and avoiding the density decrease. In Sec. 4, the EBDF scheme is tested on the calculation of the distance-dependent Potential of Mean Force 38–40 (PMF) of a dimer immersed into a Lennard-Jones fluid (system and simulation details are described in Sec. 3). The outcomes of the modified dynamical freezing method are compared to those achieved by IDF and standard SMD simulations, mainly focusing on the accuracy of the PMF profiles and on the computational cost. Conclusive remarks are given in Sec. 5.

2 2.1

Theory Nonequilibrium molecular dynamics simulations and free energy

Let us consider a system with Hamiltonian H(x), the vector x denoting coordinates q and momenta p of the atoms. According to a generic SMD simulation protocol, the system is driven out of equilibrium by means of a time-dependent control parameter λ(t) entering a harmonic potential U(q, t) = kharm [ζ(q) − λ(t)]2 added to the Hamiltonian. 41 The 4


Page 4 of 35

Page 5 of 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


parameter λ(t) is associated with a collective variable ζ(q) (such as interatomic distance, torsion angle, coordination number, etc.), depending generally on the atomic coordinates. An amount of trajectories lasting a time τ is produced, during which the control parameter evolves from the starting value λ(0) = λA to the final value λ(τ ) = λB following an established time schedule. Therefore, the end states A and B correspond to the Hamiltonians HA (x) = H(x) + U(q, 0) and HB (x) = H(x) + U(q, τ ), respectively. The initial configurations of the pulling trajectories are picked from the equilibrium distribution e−βHA (x) /QA , with β being the inverse temperature and QA the partition function of the state A. Of course, the direction of the process, from A to B, is arbitrary; analogous treatment could be done for the reverse direction. The PMF along the collective variable ζ(q), defined as R Φ(z) = −β −1 ln[ δ(ζ(q) − z)e−βH(x) dx], can be computed by using the work associated with the trajectories that switch the system from A to B and/or viceversa. For very large values

of kharm , the stiff spring approximation applies 20 and hence Φ(z) equals Φ(λ), namely the PMF expressed as a function of the control parameter. In this case, JE 1 and CFT-based relationships can be applied to compute Φ(z) ≡ Φ(λ), 4,42,43 as well as the free energy difference between A and B, 5 Φ(λB ) − Φ(λA ). If the stiff spring approximation does not hold, unidirectional 44 and bidirectional 4,45 reweighting techniques are adopted to estimate Φ(z) (the relevant relationships are summarized in Ref. 45).

2.2

Isokinetic dynamical freezing

In the context of nonequilibrium molecular dynamics simulations, the idea of dynamical freezing consists in distinguishing a MR, where the particles undergo free relaxation, from the remaining region of space, where the particles are kept in a fixed position under the assumption that their direct or indirect contribution to the dissipated work is less important. Specifically, in the IDF scheme, an isokinetic scaling of masses and velocities is applied at time 0 to the particles located outside the MR. The masses of the external particles, say m, are instantaneously increased to a very large value, M ≫ m, with a simultaneous velocity 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

scaling, v → v

p m/M , aimed at preserving the kinetic energy of the single particles. As

the total kinetic energy is unchanged, no additional work is done in the scaling process. The velocities of the external particles are thus reduced of a large amount and, due to the huge value of M , the accelerations become negligible in the equations of motion. This practically freezes the dynamics of these particles, so that the expensive calculation of the mutual forces can be avoided during the pulling trajectory. Mobile (unfrozen) particles, instead, evolve under the time-dependent Hamiltonian H(x, t) = H(x) + U(q, t). In order to limit the dissipation during nonequilibrium paths, the MR must contain the atoms or molecules mostly affected by the external control parameter, under the assumption that the system undergoes a local perturbation. 34 Although an intuitive identification of the perturbed portion of system could be relatively easy, effective definition of the MR may not be straightforward. A large size MR, for example, would allow for a better relaxation of the system and hence for a lower dissipation, ultimately leading to more accurate PMF estimates. However, an excessively large MR would increase unavoidably the computational cost of the pulling simulations. On the other side, a small MR would offer cheaper calculations, but at the prize of a large dissipation and hence inaccurate PMF estimates. Different implementations of the IDF can be designed. For example, the list of frozen particles can be fixed during the pulling trajectory or updated with a given frequency. Also the definition of MR can be time dependent, to account for the change of shape of the reactive site during the driven process. 32 Regardless the specific time evolution of the MR and the frequency of update in the list of frozen particles, a final isokinetic mass-velocity scaling is applied to the frozen particles to restore the original masses. In Ref. 32, it has been shown that the mass-velocity scaling protocol produces no additional contribution to the work W performed on the system during a pulling trajectory. Therefore, the total work Rτ is computed with the standard formula: W = 0 ∂t H(x, t) dt.

Concerning the time evolution of the system under IDF, we will show that, owing to freez-

ing, the outward particle flow through the MR boundaries is not balanced by the inward flow,

6


Page 6 of 35

Page 7 of 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


yielding an incorrect phase-space sampling and consequently a failure of the nonequilibrium work theorems. This is basically due to the fact that these theorems rely on the condition that the dynamics must preserve the canonical distribution upon fixing the control parameter λ, 36,37 a requirement which is not fulfilled in IDF. This underlying canonical distribution may however change according to the time dependence of the Hamiltonian H(x, t). The initial configuration of each SMD trajectory is randomly picked from an equilibrium simulation at fixed control parameter with no IDF scheme in place. Thus, the distribution of the initial configurations is canonical by construction. It is straightforward to show that such a distribution is not preserved if the dynamics is started by activating IDF, while keeping fixed the control parameter. To this purpose, we focus on the particle flow through the MR boundaries at the initial time. The inward and outward particle flows correspond, respectively, to the numbers of particles entering and leaving the MR in the unit of time, averaged on the configurational distribution. The canonical distribution would be conserved only if the average inward and outward flows are balanced. For example, observation of a net outward flow would reveal an average decrease of the density inside the MR. The outward particle flow at time 0 is defined as limt→0 [N in (0) − N in (t)], where N in (t) is the number of particles inside the MR at time t averaged over the system configurations resulting from evolving the initial configurations according to the Hamiltonian H(x, 0). The quantity N in (0) is clearly in accord to the canonical distribution. Since dynamically active (unfrozen) particles are inside the MR at time 0, their diffusion can only contribute to decrease, or to leave unchanged, the quantity N in (t) with respect to N in (0), giving rise to an average outward flow. Moreover, we assume the limit m/M → 0+ , since the ratio between the masses before and after the mass-velocity scaling can be chosen arbitrarily small. In these conditions, the velocities of the external particles are zero and therefore they do not contribute to the inward flow. Globally, an average outward flow (arising from the unfrozen particles) is observed at time 0, which alters the starting canonical distribution also for a system subject to stationary Hamiltonian equations of motion. The fact that an unfrozen particle, after

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

leaving the MR, may or may not undergo mass-velocity scaling according to the specific IDF scheme employed, does not matter. So far, we have examined the particle flow at time 0, thus avoiding to account for the retro-diffusion of unfrozen particles towards the MR. This can however happen in the specific IDF implementation in which mass-velocity scaling is applied only at the initial and final steps of the pulling trajectories, while leaving the unfrozen particles free to move through the sample during the trajectory. Nevertheless, we may argue that, blocking the dynamics of the external particles in the initial configuration, the particle density inside the MR at any time t will be systematically less than the equilibrium density into the MR itself. In fact, the unfrozen particles will distribute in some ratio between MR and the rest of the system through percolation phenomena. This is schematically represented in Fig. 1A. If the list of frozen particles is updated periodically, matter will accumulate near the MR boundaries. In fact, by using an updating list protocol together with an invariant MR, freezing would result irreversible. In a condensed-phase system, the increased density around the MR boundaries will create a potential-energy barrier which will stop the outward particle flow, finally reaching a stationary (non canonical) condition. This situation is sketched in Fig. 1B. From the above discussion, we conclude that the IDF is an approximated method relying on the assumption that the diffusion of unfrozen particles away from the MR is small enough to not affect significantly the underlying canonical distribution in the time of the pulling simulations. If this condition is satisfied, the use of nonequilibrium work theorems in IDF processes may be considered acceptable. This usually happens in SMD simulations with rather fast-rate pulling realizations and large MR.

2.3

Elastic barrier dynamical freezing

To tackle the problem described in Sec. 2.2, we outline the necessity of avoiding a net particle flow through the MR boundaries. Constraining the unfrozen particles to move within the MR 8


Page 8 of 35

Page 9 of 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


B

A

Figure 1: Panel A: Schematic representation of the distribution of frozen and mobile particles (blue and red filled circles, respectively) expected in stationary conditions (long-time limit) when the list of frozen particles is never updated during the simulation (IDF1-type simulations). Panel B : Schematic representation of the distribution of frozen and mobile particles expected in stationary conditions when the list of frozen particles is periodically updated during the simulation (IDF2-type simulations). In both panels, a bidimensional box is displayed for clarity reasons. The black circle represents the MR boundaries. Larger circles are used for mobile particles to highlight their distribution. prevents the density decrease, eventually preserving the canonical distribution as the stationary underlying distribution. This is the idea behind the configurational freezing method, 33,34 which is the Monte Carlo simulation analog of dynamical freezing. In configurational-freezing simulations, the configurations of the system, initially sampled according to a canonical distribution, evolve without leaving an established MR: no moves involving external (frozen) particles are generated and no moves bringing an internal (unfrozen) particle outside the MR are accepted. All other attempted moves are accepted according to the standard Metropolis criterion. 46 In principle, configurational freezing is able to generate a stationary canonical distribution on a single trajectory, thought this may only occur in an infinite time. Under the previous protocol, the system becomes practically non-ergodic, even if the pulling trajectories were realized very slowly. In spite of this, JE and CFT are valid because the dynamics is Markovian and the detailed balance is fulfilled. At the same time, the sampling occurring in the MR allows to enhance the system relaxation during the nonequilibrium process, improving the accuracy of PMF estimates. 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 35

To devise a molecular dynamics version of configurational freezing, we have to focus on the deep reason why configurational freezing works. This relies on the fact that configurational freezing is nothing but a limit case of the preferential sampling approach. 47 In preferential sampling, a high-mobility region is first identified. A parameter p defines how often we wish to move the out (external) particles relative to the in (internal) ones: p lies between 0 and 1, values close to 0 corresponding to much more frequent moves of the in particles. A Monte Carlo move consists of the following steps: (1) a particle is chosen at random; (2) if it is in, a trial move is made; (3) if it is out, a random number is generated uniformly between 0 and 1. If p is greater than the random number then a trial move is made. If not, then we return to step (1) without accumulating any averages. Trial moves are accepted with a probability dependent on p as follows: acc(iin → jin ) = min(1, pj /pi ) acc(iin → jout ) = min(1, a pj /pi )

(1)

acc(iout → jin ) = min(1, b pj /pi ) acc(iout → jout ) = min(1, pj /pi ), where pi and pj are the canonical probabilities of the microstates i and j, and a and b are a = p{1 + (p − 1)/[pN + (1 − p)Nin ]}−1

(2)

b = p−1 {1 − (p − 1)/[pN + (1 − p)Nin ]}−1 , where N is the total number of particles and Nin is the number of in particles in the microstate i. In the preferential sampling machinery, selecting a out particle is less probable than selecting a in particle (p < 0.5). This makes the probability of generating a iout → jin move smaller than a iin → jout move. Equilibrium between outward and inward flows is established if the probability of accepting an attempted iin → jout move is made smaller than that of accepting a iout → jin move. Equation 1 just sets the quantitative conditions to satisfy the above requirement, which corresponds to enforce the detailed balance. Configurational10


Page 11 of 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


freezing algorithm results from setting p = 0, namely, from selecting only in particles. Preservation of the detailed balance in dynamical freezing can thus be obtained by designing a way of equilibrating the inward and outward particle flows through the MR boundaries. To explain this balancing scheme we resort to heuristic arguments. We first suppose to enforce a moderate increase of the masses of the external particles, by applying an isokinetic mass-velocity scaling to them by some factor M > m to make their dynamics slower, but not frozen (by setting M not much larger than the natural masses m of the particles). This operation would reduce the diffusion motion of the external particles, at the same time keeping their dynamics Hamiltonian. An imbalance between inward and outward flows would occur because of the artificially slower diffusion of the external particles. To restore the flow equilibrium, we would need to decrease somehow the outward particle flow without affecting the Hamiltonian dynamics of the internal particles. This could be realized through the introduction of a potential-energy step function, which holds 0 inside the MR and U > 0 outside the MR. In the practice, however, in order to have a derivable potential-energy function everywhere in the space, a Fermi-like function could be used. The resulting energy barrier would decrease the outward particle flow. A proper tuning of the parameters M (affecting the inward flow) and U (affecting the outward flow) can lead to a correct flow balancing. The particles leaving the MR would undergo mass scaling by a factor M/m, while the particles entering the MR by m/M (mass restoration). The scaling could be continuous within the region where the Fermi-like potential is not constant. At the end of the trajectory, all the particles with mass M should undergo an isokinetic mass-velocity scaling to restore the original masses. To preserve the balancing between inward and outward flows, an increase of the M parameter would correspond to an increase of U and viceversa. Of course, by setting U = 0 everywhere in the space and M = m, one would generate a conventional (dynamical-freezing free) dynamics. Thanks to the overall Hamiltonian dynamics, as well as to the balancing of the particle flow through the MR surface obtained by means of a proper setting of M and U , one would design a dynamical analog of the preferential sampling Monte

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

Carlo method. Following this analogy, we can therefore switch to the EBDF scheme (from this sort of dynamical preferential sampling) as it has been done in switching from preferential sampling Monte Carlo to configurational freezing. 33 Avoiding sampling of the particles external to the MR in configurational freezing corresponds to stop the dynamics outside the MR, namely to enforce an isokinetic mass-velocity scaling at the beginning of the trajectories by setting M ≫ m. As a consequence, in order to preserve the particle flow, we need to introduce a potential-energy step (Fermi-like) function with U ≫ 0. The virtually infinite value of U “transforms” the Fermi-like potential into a stiff energy wall acting on the internal particles as an elastic barrier. At the beginning of the trajectories, a work equal to (N − Nin )U is performed on the system because of the introduction of the Fermi-like potential, while the isokinetic mass-velocity scaling does not produce additional work. Under these conditions, the number of frozen and mobile particles does not change during the dynamics, and hence the same work (N − Nin )U is subtracted at the end of the trajectory, when removal of the Fermi-like potential is applied simultaneously to the inverse mass-velocity scaling. Globally, no work arising from the EBDF scheme is yielded during a pulling trajectory. Finally, we note that the MR, so far assumed stationary along the path, can also be taken time dependent without affecting the detailed balance condition at each step of the process, and hence the validity of nonequilibrium work theorems. Indeed, the potentialenergy barrier, even though the MR surface changes during the path, will prevent the exiting of particles outside the momentary MR. In spite of this, the implementation of a MR with variable shape presents various technical difficulties, mainly arising from enforcing constanttemperature conditions. Thus, in order to preserve simplicity in the algorithm, we do prefer to consider a MR with fixed shape and size.

12


Page 12 of 35

Page 13 of 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


3

System and simulation details

The system, also employed in Ref. 32 to illustrate the IDF method, is formed by a bistable dimer immersed in a dense Lennard-Jones fluid. The dimer consists of two particles interacting through a double-well potential, Udim (r) = 3[(r − 1)2 − 0.1](r − 3)2 .

(3)

Udim (r) has two minima at the interparticle distances r ≃ 1 and r ≃ 3, corresponding to compact and extended configurations of the dimer, respectively. All quantities into Eq. 3 and in the rest of the article are in reduced units. The solute (dimer) and solvent particles have the same masses and evolve in a cubic box with standard periodic boundary conditions. The Lennard-Jones potentials for solvent-solvent and solute-solvent interactions are identical and vanish in the distance range 3.0-3.5 through a cubic switching function. The equations of motion are integrated with a time-step of 5 × 10−3 . A neighbor list, 48 updated every 92 time-steps (with cutoff radius 5.5), is used to compute interparticle forces. In order to illustrate better the computational advantages provided by EBDF in simulating large size systems, simulations have been performed with different box dimensions, fixing the particledensity at the value of 0.85. Specifically, samples of 6800, 10200, 13600, 17000 and 20400 particles have been analyzed. In SMD simulations, the harmonic potential U(r, t) = 5000[r − λ(t)]2 is used to switch the dimer between the compact (λA = 0.5) and extended (λB = 3.5) configurations along an established direction, thus avoiding the Jacobian correction to the PMF. 2000 forward realizations (λA → λB , corresponding to the dimer dissociation process) and 2000 backward realizations (λB → λA , corresponding to the dimer association process) are performed at constant pulling speed of 0.15 with constant-volume constant-energy (NVE) equations of motion, by starting from microstates sampled at equilibrium. These microstates are picked every 100 time-steps from two constant-volume constant-temperature (NVT) simulations in 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 35

which the λ parameter is fixed in turn to λA and λB (initial microstates for forward and backward realizations, respectively). The large force constant adopted into U(r, t) allows us to compute the PMF in the regime of stiff spring approximation. 20 As stated above, the particle density is 0.85, while the temperature is kept at the constant value of 0.8 by means of a Nosé-Hoover chain thermostat. 49 In both equilibrium and SMD simulations, the center of mass of the dimer is restrained around a fixed position through a harmonic potential with 104 force constant. The above setup is common to all SMD schemes under study. In the EBDF case, the MR for a given SMD trajectory is a sphere centered in a fixed point in space corresponding to that of the dimer center of mass in the initial microstate. According to this choice, the elastic energy barrier at the MR boundaries is enforced by P −2 adding the potential energy N i=1 Ubar (di ) to the Hamiltonian, where the sum is over the N − 2 solvent particles and Ubar (di ) is an energy term depending on the distance of the ith particle from the MR center: 1 Ubar (di ) = K(di − rcut )2 2

if

di ≥ rcut

Ubar (di ) = 0

if

di < rcut ,

(4)

where rcut = 6 is the radius of the sphere representing the MR. In the limit K → ∞, the total energy is conserved upon collision of a particle with the MR, and hence the potential Ubar (di ) acts as a perfect elastic-wall device, exactly as required to apply EBDF (see discussion in Sec. 2.3). In practice, to get a stable integrator for the equations of motion, K must take a large but finite value, so that the potential-energy barrier is only approximately elastic. Moreover, in order to use affordable time-steps, K should not be too large. As a matter of fact, there is not a precise recipe to determine the optimal value of K. To fulfill the elastic-barrier requirement in EBDF, preserving, at the same time, the integrability of the equations of motion with relatively large time-steps, the force associated to the MR barrier must be larger than the intermolecular forces, but not much larger than the forces typical

14


Page 15 of 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


of bond stretching motions. In this respect, multiple time-steps algorithms could also be implemented to avoid the use of small time-steps in the treatment of interparticle forces. In the present calculations, we have set K = 104 , which corresponds to the force constants used to perform the SMD trajectories and to restrain the position of the dimer center of mass. It is worth noting that, owing to the finite stiffness of the energy barrier, a non null outward flow of particles through the MR boundaries is expected at the initial time (the flow would be exactly zero for a perfectly stiff barrier). However, this outward flow would be balanced almost immediately by an inward flow arising from the unavoidable retro-diffusion of unfrozen particles induced by the energy barrier. Thus, an equilibrium would be reached between outward and inward flows, resulting in a slight decrease of the particle density inside the MR and a complementary density increase within a thin shell external to the MR. At the end of the SMD trajectory, the latter effect would yield a spurious energy exchange between system and barrier device, corresponding to a net dissipated work. However, in the assumption of a stiff energy barrier, we may neglect this additional dissipation. The results obtained from EBDF are compared to those of the IDF and standard schemes, by applying the same operative conditions and using the same sets of initial microstates. Concerning the IDF schemes, we consider three simulation setups, denoted as IDF1, IDF2 and IDF3. In the IDF1 scheme, mass-velocity scaling of the particles external to the MR is applied only at the first step of the SMD trajectories (the scaling enforced at the end of the trajectories is immaterial). The MR is a sphere, whose center and radius are defined as done for the EBDF case, while no update of the list of frozen particles is operated during the pulling trajectories. We remind that, for a SMD simulation performed in a quasi-static pulling regime, the IDF1 scheme would behave as represented in Fig. 1A. In the IDF2 scheme, the MR is defined as in IDF1, but the list of frozen particles is updated every 2.5 time units, according to their positions with respect to the stationary MR. Enforcing this scheme in quasi-static pulling regimes, one would roughly observe the

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 35

situation sketched in Fig. 1B. In the IDF3 scheme, the shape of the MR changes during the SMD trajectories, being the MR defined as the union of spheres of radius rcut = 6 centered on the two dimer particles. Moreover, the list of frozen particles is updated every 2.5 time units. Finally, standard SMD simulations, hereafter denoted as ST, are carried on by keeping all solvent particles free to relax.

4

Simulation tests

The purpose of the present tests is to carry out a comparative analysis between the EBDF protocol and the IDF1, IDF2, IDF3 and ST ones, addressing accuracy in PMF estimates and computational cost, above all. A qualitative view of the accuracy featuring the simulation schemes under study is given in Fig. 2, where the PMF as a function of the distance between the dimer particles is reported for samples with different size (6800, 13600 and 20400 particles). The calculation has been done using the bidirectional nonequilibrium work approach of Ref. 4. In a regime of stiff spring approximation, as in our case, the PMF is nF e−βWF (λ) Φ(λ) = −β ln nF + nB eβ[∆F −WF (λB )] nB eβ[WB (λA )−WB (λ)] + , nF + nB eβ[∆F +WB (λA )] back −1

+ forw

(5)

where WF (λ) is the work done on the system to switch the control parameter from λA to λ in the forward realizations (λA → λB direction) and WB (λ) is the work done on the system to switch the control parameter from λB to λ in the backward realizations (λB → λA direction). The angular brackets indicate path-ensemble averages in the forward and backward directions calculated over nF and nB realizations, respectively (in our case, 2000 realizations for each direction). Finally, ∆F = FB − FA is the free energy difference between the end states, calculated according to a Bennett-like 50 relationship. 5 In Fig. 2, the PMF curves Φ(λ) are

16


10

Φ(λ) vs Φref(λ)

8 6 4

IDF1 IDF2

2 0

IDF3 EBDF ST

N = 6800

10 8

Φ(λ) vs Φref(λ)

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


6 4

IDF1 IDF2

2 0

IDF3 EBDF ST

N = 13600

10

Φ(λ) vs Φref(λ)

Page 17 of 35

8 6 4

IDF1 IDF2

2 0

IDF3 EBDF ST

N = 20400

0.5

1.0

1.5

λ

2.0

2.5

3.0

3.5

Figure 2: PMF as a function of the dimer interparticle distance (corresponding to the control parameter λ), calculated from IDF1, IDF2, IDF3, EBDF and ST simulation schemes (solid lines). Each Φ(λ) profile is compared to the reference one, Φref (λ) (dashed lines), obtained from thermodynamic integration. For the sake of clarity, the PMF curves are shifted (labels in the panels indicate the adopted method). The three panels refer to calculations performed on samples with different size (see N , i.e., the number of particles, in the panels). fitted to a reference (accurate) PMF, Φref (λ), recovered from thermodynamic integration 38 exploiting a series of NVT molecular dynamics simulations lasting 200 time units each, and distributed with a λ resolution of 0.05 for numerical integration purposes. Fit is done by shifting Φ(λ) by an offset obtained with a least-squares procedure aimed at minimizing the root mean square deviation 42

σ = |λB − λA |

−1

Z

λB

2

[Φ(λ) − Φref (λ)] dλ λA

17


12

.

(6)


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 35

Regardless of the sample size, clear indications on the relative performances of the various schemes cannot be easily inferred from the figure. From our point of view, the similarity among the Φ(λ) profiles and, especially, the small differences between dynamical-freezing and ST outcomes are extremely positive aspects. Actually, a detailed comparative inspection of the IDF and EBDF curves, perhaps reveals EBDF as the most performing dynamicalfreezing scheme (for the N = 6800 sample, EBDF outperforms even ST). The scenario does not change significantly using samples of 10200 and 17000 particles, as it can be appreciated in Fig. S1 of the Supporting Information. The above qualitative considerations are supported by the deviation of Φ(λ) from Φref (λ) expressed in terms of root mean square deviation σ (Eq. 6). The σ parameter related to IDF1, IDF2, IDF3, EBDF and ST schemes is reported in Fig. 3 as a function of the number of particles. In addition to the results achieved by the bidirectional method (Fig. 3A), we also report those obtained through the JE applied to forward and backward realizations (Figs. 3B and 3C, respectively): ΦF (λ) = −β lnhe−βWF (λ) iforw ,

(7)

ΦB (λ) = −β lnhe−βWB (λ) iback .

(8)

Here, the angular brackets have the same meaning as in Eq. 5. We point out that, owing to the different definitions of WF (λ) and WB (λ) (Eq. 5), ΦF (λ) and ΦB (λ) are shifted. In any case, the offset of ΦF (λ) and ΦB (λ) applied to compute σ is determined in order to satisfy the equalities ΦF (λA ) = Φref (λA ) and ΦB (λB ) = Φref (λB ). Comparisons between ΦF (λ), ΦB (λ) and Φref (λ) are displayed in Fig. S2 of the Supporting Information. While for the bidirectional method of Eq. 5 almost all σ values fall within the interval (0.05-0.25), the JE (Eqs. 7 and 8) gives deviations that may exceed 0.8, and, in any case, systematically above 0.2. Considering the biased sampling problems affecting the finite-sample exponential averages, 51–53 the worse accuracy of the JE is not surprising. No systematic dependence

18


A

ST EBDF IDF1 IDF2 IDF3

0.8

σ

0.6 0.4 0.2 0 0.8 0.6

σ

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.4 0.2

B

0 0.8 0.6

σ

Page 19 of 35

0.4 0.2 0

C 6800

10200

13600

17000

20400

N

Figure 3: Root mean square deviation σ (Eq. 6) as a function of the number of particles in the sample, N , computed from IDF1, IDF2, IDF3, EBDF and ST simulations (see legend). Panel A: Φ(λ) in Eq. 6 is from the bidirectional method (Eq. 5). Panel B : Φ(λ) ≡ ΦF (λ) in Eq. 6, with ΦF (λ) from Eq. 7. Panel C : Φ(λ) ≡ ΦB (λ) in Eq. 6, with ΦB (λ) from Eq. 8. Lines linking the circles are guides for eyes. of σ from the sample size can be inferred from Fig. 3, and, most importantly, the noisy trend of the data does not point to significant differences among the various schemes. The only systematic feature can, perhaps, be noted in the general loss of accuracy of the IDF3 scheme when combined with the bidirectional method and the JE in the forward direction. Although a clear explanation to such a behavior cannot be supplied easily, one should not exclude the simple accidental origin. Overall, we may conclude that all dynamical-freezing schemes, and in particular the EBDF one, behave comparably well with respect to the standard approach. From the physical standpoint, this basically means that the chosen MR 19



well embeds the region of space where the driven dynamics of the dimer occurs, avoiding excessive dissipation via effective system relaxation. After verifying that EBDF and the other IDF schemes do not alter significantly the PMF, a comparative analysis between dynamical-freezing and ST protocols is required to quantify the possible advantages of the former approaches in terms of computational cost. Such an information is illustrated in Fig. 4, where we report the wall-clock time necessary to perform simulation trajectories lasting 25 time units by using the IDF1, IDF2, EBDF and ST methods. Results for the IDF3 scheme are not reported because the intent is to evaluate the 2500 ST EBDF IDF1 IDF2

2000

Wall-clock time (s)

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

1500

150

1000

135 120

500

105 6800

0

6800

10200

13600

10200 13600 17000 20400

17000

20400

N

Figure 4: Wall-clock time (in s units) needed to perform IDF1, IDF2, EBDF and ST simulations (see legend) lasting 25 time units as a function of the number of particles in the sample, N . The dimer particles are restrained in the compact configuration (λ = λA ). A zoomed view of the EBDF, IDF1 and IDF2 data is displayed in the inset. Lines are guides for eyes. The simulations have been performed with a personal computer powered by a AMD Phenom 9650 Quad-Core processor. computational cost of dynamical-freezing schemes in conditions of comparable MR size. In this respect, we remind that IDF3 is the only protocol in which the MR changes in shape and size according to the externally guided evolution of the dimer particles. In any case, in these tests, no particle-particle steering is applied, so that the dimer particles are restrained about the compact dimer configuration, i.e., λ = λA . In such simulation conditions, therefore, IDF3 is equivalent to IDF2. For the simulations supplied with dynamical freezing (EBDF, 20


Page 20 of 35

Page 21 of 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


IDF1 and IDF2), the MR setup is described in Sec. 3. Moreover, the same initial system configuration has been used for all simulations. The power of dynamical freezing mainly lies in the fact that CPU time is pretty independent on the system size. 32 Therefore, in order to highlight this important aspect, in Fig. 4 we report the wall-clock times as a function of the number of particles in the sample. To better reveal few details of the dynamical-freezing results, a zoomed view of the EBDF, IDF1 and IDF2 data is also reported. Indeed, the speedup obtained by dynamical freezing is evident. To accomplish a simulated trajectory using ST, a computational time ranging from 500 to 2500 s, in dependence of N , is necessary, whereas dynamical-freezing methods require a much more modest effort (100-160 s). In ST simulations, the wall-clock time as a function of N is approximately linear, which is consistent with the use of a neighbor-list algorithm. 48 A linear increasing trend of the wallclock time with N is also observed when using dynamical freezing, even if the slopes are more than 30 times smaller with respect to ST (37, 72 and 83 times for IDF2, IDF1 and EBDF, respectively). Indeed, while in ST simulations the linear trend is yielded by the increase of the number of particle-particle interactions, in dynamical freezing the effect is substantially ascribed to the periodical update of the neighbor list, which takes place every 92 time-steps. Among the dynamical-freezing schemes, the greater slope noted in IDF2 is due to a different implementation of the neighbor-list algorithm, ultimately correlated to the need of updating periodically the list of frozen particles. The previous results reveal a close similarity between EBDF and IDF as the computational cost is concerned. Also the accuracies in PMF estimation are comparable, and both approaches well reproduce the reference PMF profile. Nevertheless, in Sec. 2.2 it has been argued that the altered particle distribution at the MR boundaries upon applying IDF compromises the validity of nonequilibrium work theorems. The theorems can, however, be used if the spurious density decrease arising from IDF does not affect significantly the nonequilibrium work distribution, which typically happen when pulling simulations are short and the MR is large. According to these considerations (see also discussion in Sec. 2.2), the

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

comparable behaviors of EBDF and IDF can be ascribed to an appropriate choice of the MR and/or to the fast-rate pulling protocol employed. As we will see below (discussion of Fig. 7), the MR size is certainly one of the two conditions which turns out to be grounded. Quantitative insights on the outward particle flow through the MR boundaries can be gained computing ∆Nin = Nin (0) − Nin (τ ), which corresponds to the difference between the numbers of particles inside the MR 54 at the beginning and at the end of the trajectory (Nin (0) and Nin (τ ), respectively). Since in dynamical-freezing simulations the number of particles inside the MR can only decrease, ∆Nin can take values in the range [0, Nin (0)]. In order to evaluate the diffusion effects regardless of the time-dependent perturbation of the system, ∆Nin is computed using the MR setup described in Sec. 3, but restraining the dimer in its compact configuration, namely fixing λ to λA . Moreover, ∆Nin is computed by averaging over 100 trajectories lasting τ = 25 time units, each starting from a canonically sampled configuration. Fig. 5A displays ∆Nin as a function of the number of particles in the sample, calculated from IDF1, IDF2 and EBDF simulations. The outcomes of IDF3 simulations are not reported because, as explained when discussing Fig. 4, IDF3 and IDF2 schemes coincide in stationary conditions. The trend of ∆Nin is nearly independent on the system size, because the diffusion phenomenon is limited to the neighborhood of the MR surface, and hence unaffected by the overall size of the simulation box. Particle-density decrease inside the MR discloses an average outward particle flow for all dynamical-freezing schemes. The effect is more pronounced for IDF2 (∆Nin ≃ 60), while it is minimum for EBDF (∆Nin ≃ 5). The larger decrease noted by using IDF2 with respect to IDF1 (∆Nin ≃ 25) is ascribed to the periodic update of the list of frozen particles, which is done every 2.5 time units. When the list of frozen particles is updated, the mobile particles external to the MR are isokinetically frozen. Thus, the freezing mechanism prevents back diffusion of a variable number of particles each time it is applied, transforming the MR boundaries into a sort of “black hole” for mobile particles. Immediately after the freezing event, the inward particle flow becomes null, while

22


Page 22 of 35

Page 23 of 35

12 60

3

10 8 6

4

K = 10

30

4

20

5

K = 10

10 0

B

∆Nin

40

K = 10

A

EBDF IDF1 IDF2

50

∆Nin

6800

10200 13600 17000 20400 6800

2

0 10200 13600 17000 20400

N

N

Figure 5: Panel A. Decrease of the number of particles inside the MR, ∆Nin , as a function of the number of particles in the sample, N , for IDF1, IDF2 and EBDF simulations (see legend) lasting 25 time units, with the dimer particles being restrained in the compact configuration. In EBDF simulations, the force constant entering the potential-energy barrier is K = 104 (see Eq. 4). Panel B. Decrease of the number of particles inside the MR, ∆Nin , as a function of N for EBDF simulations performed using three different force-constant values, K = 103 , 104 , 105 (see labels in the panel). Lines are guides for eyes. a residual outward flow persists. The resulting effect is a rapid decrease of density inside the MR, a process which can be well recognized in Fig. 6, where ∆Nin as a function of time is reported for IDF1, IDF2 and EBDF simulations. 55 In fact, in the IDF2 simulation, 70

EBDF IDF1 IDF2

60 50

∆Nin

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


40 30 20 10 0

0

5

10

15

20

25

Simulation Time

Figure 6: Decrease of the number of particles inside the MR, ∆Nin , as a function of time, for IDF1, IDF2 and EBDF simulations (see legend). In EBDF simulations, the force constant entering the potential-energy barrier is K = 104 (Eq. 4). The curves are averaged on 1000 simulations in stationary conditions using the setup specified in the text.

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 35

large changes in slope occur at regular time intervals in correspondence of frozen-particle list updates, highlighting the annihilation of the inward particle flow. A plateau is reached at long times owing to the formation of a high-density particle layer, external to the MR boundaries, which prevents further outward diffusion (this will be more evident in the following discussion of Fig. 7). Such an effect is also observed in IDF1 simulations. In this case, however, the high-density particle layer is due to a dynamical equilibrium between the particles internal to the MR and the particles percolated outside the MR. In the EBDF simulations, equilibrium is reached almost instantaneously due to the stiff potential-energy barrier. A further interesting feature of Figs. 5A and 6 is the residual, albeit small, value of ∆Nin observed in EBDF simulations. This is due to a not perfect stiffness of the potentialenergy barrier, which is regulated by the force constant K into Eq. 4. In fact, as discussed in Sec. 2.3, an ideal potential-energy barrier would be featured by infinite force constant K. Accordingly, an enhancement of K should be accompanied by a decrease of ∆Nin . This is indeed confirmed in Fig. 5B, where ∆Nin as a function of N is reported for EBDF simulations carried out using three different values of K, namely 103 , 104 and 105 (with usual simulation setup). By increasing K, ∆Nin decreases from ∼11 to less than 2: in conclusion, we reasonably do expect that limK→∞ ∆Nin = 0. The quantity ∆Nin , while allowing to detect a global density change, does not provide information on the spatial localization of the phenomenon. To this regard, a more exhaustive picture is achievable through the distribution function g(s), quantifying the average particledensity on a spherical shell centered on the MR center, with thickness ∆s = 0.05 and radius s, relative to the bulk particle-density,

g(s) =

V n(s) . 4πs2 ∆s N

(9)

In Eq. 9, n(s) is the average number of particles in the spherical shell of volume 4πs2 ∆s, N is the total number of particles and V is the volume of the simulation box. The g(s)

24


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


functions for IDF1, IDF2, EBDF and ST are computed on the final configurations of 1000 trajectories, 25 time units long, performed on a sample of 10200 particles in the stationary conditions specified above. 56 The g(s) functions are shown in Fig. 7. All profiles display a relevant structuration, as ST EBDF IDF1 IDF2

2.8 2.4

2.0 1.6 1.2 0.8

2.0

0.4

g(s)

Page 25 of 35

1.6

0.0

5.6 5.8

6

6.2 6.4

1.2 0.8 0.4 0.0

0

1

2

3

s

4

5

6

7

Figure 7: Distribution functions g(s) of the distance s of the particles from a reference point (Eq. 9), computed from final configurations of ST, EBDF, IDF1 and IDF2 simulations (see legend). In the case of dynamical-freezing simulations, the reference point is the MR center, while, for the ST simulation it corresponds to the dimer center of mass. In the EBDF simulation, the force constant entering the potential-energy barrier is K = 104 (Eq. 4). We remind that MR boundaries correspond to s ≡ rcut = 6. A zoomed view around the MR boundary zone (vertical dashed line) is reported in the inset. revealed by the fifth solvation shell of the dimer visible in the ST simulation, which represents the reference. In the current context, the most important feature is the large deviation from the ST reference curve yielded by IDF1 and IDF2 schemes, especially near the distance featuring the MR boundaries (s ≃ 6). In all cases, the negative and positive deviations from the ST curve, occurring around s ≃ 6, are evidence of the net outward particle flow featuring the dynamical-freezing mechanism. According to the change of density observed in Figs. 5 and 6, the deviation is more pronounced for IDF2, which spread throughout the second peak of g(s). In the IDF1 case, only the broad peaks beyond 3.5 (fourth and fifth solvation shells) display significant deviation from the reference curve. We notice, however, that formation of a high-density layer in IDF1 simulations does not adhere to a statistical 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

representation of the equilibrium state. Rather, the expected stationary situation would be the one sketched in Fig. 1A, where a fraction of mobile particles is distributed homogeneously either inside or outside the MR. The limit distribution of Fig. 1A would yield a step in the g(s) function: specifically, g(s) < 1 for s < 6 and g(s) > 1 for s > 6. Clearly, the time needed to reach a stationary distribution would depend on the sample density outside the MR, being longer for higher densities. 57 At variance with IDF1 and IDF2, the EBDF curve fits very well the ST one, apart from a small bump at s ≃ 6. We notice, however, that this slight density deviation is due to the approximated stiffness of the potential-energy barrier. In fact, the bump practically disappears when using a larger force constant K (see Fig. S3 in the Supporting Information).

5

Concluding remarks

Dynamical freezing has been revealed to enhance the efficiency of free energy calculations based on nonequilibrium path sampling within the context of Jarzynski 1 and Crooks 2 work theorems. In this article, a critical analysis of the Isokinetic Dynamical Freezing (IDF) simulation approach 32 is provided, addressing weakness and advantages of the method. Under dynamical freezing, a consistent computational gain in nonequilibrium steered molecular dynamics simulations is obtained blocking the time evolution of particles located outside an established region of space. The choice of this Mobility Region (MR) is somehow arbitrary. If, on one side, the MR must be large enough to prevent excessive dissipation during the driven simulations, on the other side, it must not be too large in order to ensure a significant simulation speed-up. Indeed, the expensive calculation of interatomic forces among frozen particles becomes irrelevant for the system evolution, and hence it can be safely avoided. This leads to a relevant saving in CPU time with respect to the standard approach, even of several orders of magnitude. Nevertheless, we have outlined an important problematic aspect of the IDF scheme,

26


Page 26 of 35

Page 27 of 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


which is essentially related to a spurious decrease of particle density inside the MR during the driven process. This systematic drop in density eventually invalidates the approaches adopted to compute the potential of mean force with respect to collective variables of the system, i.e. the Jarzynski 1 and Crooks 2 work theorems. Inspired by the Monte Carlo version of dynamical freezing, called configurational freezing, 33,34 an alternative strategy to enforce dynamical freezing to particles outside the MR is proposed. We have shown that the condition necessary to guarantee validity of nonequilibrium work theorems can be set through the introduction of a potential-energy barrier, which defines the MR boundaries preventing a net outward particle flow. This method, called Elastic Barrier Dynamical Freezing (EBDF), has been tested on a pulling process applied to a two-particle dimer immersed into a Lennard-Jones fluid. Insertion of a rigid barrier for defining the MR has been approximated by a stiff harmonic potential. The potential of mean force as a function of the dimer particle-particle distance has been computed with EBDF and three different IDF schemes and the results have been compared to those of the standard steered molecular dynamics. In spite of the considerable speed-up obtained by using EBDF and IDF, especially for large size systems, computational accuracies have been found comparable to that of the standard approach. The good performances of the IDF schemes can be ascribed to a proper choice of the MR, which is evidently large enough to prevent a significant deviation of the work distribution upon density depletion. In fact, while the nonequilibrium work theorems can be safely used in the EBDF framework, the IDF technique can only be applied under the assumption that the change of particle density does not alter significantly the work distribution. An important issue to be considered is the possibility of extending the EBDF scheme to constant-pressure constant-temperature (NPT) thermodynamical ensembles, to allow for Gibbs free energy calculations. This improvement is not obvious, because EBDF is intrinsically a constant-volume based simulation approach. In this respect, we may however anticipate that it is possible to switch from NPT to NVT or NVE equations of motion at

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

the beginning of the pulling trajectories. Such a switch would be based on an isokinetic scaling of the inertial factors and momenta associated with thermostat and barostat devices (switch to NVE) or with the barostat device alone (switch to NVT). We notice that, based on the same principles, a NVT to NVE switch has been realized in the present calculations. The formal treatment of this extension of the method and the related numerical tests will be object of further investigations. Furthermore, we outline the almost full compatibility of EBDF with early path-sampling techniques developed in the framework of nonequilibrium simulation schemes. 27–29,58 In this respect, the relative simplicity of the method may easily allow implementation into modern simulation programs devised for performing steered molecular dynamics simulations. 59 The method here presented can be extended to molecular systems straightforwardly. The idea of treating a molecular system as formed by two subsystems, namely a core to which we are interested to and a surrounding environment considered in much less detail, is encountered often in computational chemistry (in dynamical freezing, the “detail” we refer to is of dynamical nature; indeed, dynamical freezing is based on neglecting the dynamics of molecules not involved in dissipation mechanisms). The situation resembles, for instance, that encountered in studies based on quantum-mechanics molecular-mechanics techniques 60 or in ab initio investigations of molecules where the solvent is accounted for by a polarizable continuum model, 61,62 or by reaction field approaches. 63–65 We can however envisage important cases for which dynamical freezing may not be a suitable route to free energy calculations. For example, computing folding free energies is clearly a challenging task, because, in this kind of processes, the dissipation can be distributed quite homogeneously through the sample. This implies that, in order to obtain acceptable accuracies, the mobility region should correspond to a significant amount of simulation volume, at the price of a modest gain in computer speed. In other cases involving proteins, such as protein-ligand binding free energy calculations, 14 dissipation is instead expected to be localized around the binding pocket, and hence, limiting the mobility region to the residues of the

28


Page 28 of 35

Page 29 of 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


binding domain and to the surrounding area, may provide important simulation speedup. We deserve a final remark to systems featured by long range interactions, such as electrostatic interactions. In this respect, regardless of the adopted approach, be it the standard Ewald 66 or a particle-mesh based 67 Ewald method, the most critical aspect is certainly the treatment of the reciprocal lattice contribution to forces and energy. Although this may present a considerable programming difficulty, accounting for dynamical freezing is expected to lead to computational gain, as the sum corresponding to the charge-based structure factor can be restricted to the mobile atoms alone, while taking into CPU memory the contribution arising from dynamically frozen atoms. These and other issues aimed at improving the methodology will be pursued in forthcoming studies.

Associated content Fig. S1. PMF as a function of the dimer interparticle distance, calculated from IDF1, IDF2, IDF3, EBDF and ST simulation schemes with 10200 and 17000 particles. Fig. S2. PMF as a function of the dimer interparticle distance, calculated from IDF1, IDF2, IDF3, EBDF and ST simulation schemes with 6800, 10200, 13600, 17000 and 20400 particles by using the Jarzynski equality in the forward and backward directions. Fig. S3. Distribution functions of the distance of the particles from a reference point, computed from final configurations of ST and EBDF simulations of 10200 particles. This information is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgement This work was supported by the Italian Ministero dell’Istruzione, dell’Università e della Ricerca.

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

References (1) Jarzynski, C. Phys. Rev. Lett. 1997, 78, 2690–2693. (2) Crooks, G. E. J. Stat. Phys. 1998, 90, 1481–1487. (3) Crooks, G. E. Phys. Rev. E 2000, 61, 2361–2366. (4) Minh, D. D. L.; Adib, A. B. Phys. Rev. Lett. 2008, 100, 180602. (5) Shirts, M. R.; Bair, E.; Hooker, G.; Pande, V. S. Phys. Rev. Lett. 2003, 91, 140601. (6) Giovannelli, E.; Cardini, G.; Chelli, R. Phys. Rev. E 2015, 92, 043310. (7) Trang, N. D.; Paolo, C.; Varani, G.; bussi, G. J. Chem. Theory Comput. 2013, 9, 1720–1730. (8) Giovannelli, E.; Cardini, G.; Gellini, C.; Pietraperzia, G.; Chelli, R. J. Chem. Theory Comput. 2015, 11, 3561–3571. (9) Zerbetto, M.; Piserchia, A.; Frezzato, D. J. Comput. Chem. 2014, 35, 1865–1881. (10) Gong, Z. P.; Quan, H. T. Phys. Rev. E 2015, 92, 012131. (11) Xiao, G.; Gong, J. Phys. Rev. E 2014, 90, 052132. (12) Gundermann, J.; Siegert, S.; Kantz, H. Phys. Rev. E 2014, 89, 032112. (13) Piserchia, A.; Zerbetto, M.; Frezzato, D. Phys. Chem. Chem. Phys. 2015, 17, 8038– 8052. (14) Nicolini, P.; Frezzato, D.; Gellini, C.; Bizzarri, M.; Chelli, R. J. Comput. Chem. 2013, 34, 1561–1576. (15) Siders, P. D. Mol. Simul. 2015, DOI: 10.1080/08927022.2015.1083101.

30


Page 30 of 35

Page 31 of 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


(16) Forti, F.; Boechi, L.; Estrin, D. A.; Marti, M. A. J. Comput. Chem. 2011, 32, 2219– 2231. (17) Bastug, T.; Chen, P. C.; Patra, S. M.; Kuyucak, S. J. Chem. Phys. 2008, 128, 155104. (18) Ngo, V.; Stefanovski, D.; Haas, S.; Farley, R. A. PLoS One 2014, 9, e86079. (19) Ozer, G.; Quirk, S.; Hernandez, R. J. Chem. Theory Comput. 2012, 8, 4837–4844. (20) Park, S.; Schulten, K. J. Chem. Phys. 2004, 120, 5946–5961. (21) Procacci, P.; Marsili, S.; Barducci, A.; Signorini, G. F.; Chelli, R. J. Chem. Phys. 2006, 125, 164101. (22) Cossins, B. P.; Foucher, S.; Edge, C. M.; Essex, J. W. J. Phys. Chem. B 2009, 113, 5508–5519. (23) Ytreberg, F. M.; Zuckermann, D. M. J. Chem. Phys. 2004, 120, 10876–10879. (24) Sun, S. X. J. Chem. Phys. 2003, 118, 5769–5775. (25) Geissler, P. L.; Dellago, C. J. Chem. Phys. B 2004, 108, 6667. (26) Schmiedl, T.; Seifert, U. Phys. Rev. Lett. 2007, 98, 108301. (27) Lechner, W.; Oberhofer, H.; Dellago, C.; Geissler, P. L. J. Chem. Phys. 2006, 124, 044113. (28) Chelli, R.; Gellini, C.; Pietraperzia, G.; Giovannelli, E.; Cardini, G. J. Chem. Phys. 2013, 138, 214109. (29) Giovannelli, E.; Gellini, C.; Pietraperzia, G.; Cardini, G.; Chelli, R. J. Chem. Phys. 2014, 140, 064104. (30) Ozer, G.; Valeev, E. F.; Quirk, S.; Hernandez, R. J Chem. Theory Comput. 2010, 6, 3026–3038. 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

(31) Ozer, G.; Keyes, T.; Quirk, S.; Hernandez, R. J. Chem. Phys. 2014, 141, 064101. (32) Nicolini, P.; Chelli, R. Phys. Rev. E 2009, 80, 041124. (33) Nicolini, P.; Frezzato, D.; Chelli, R. J. Chem. Theory Comput. 2011, 7, 582–593. (34) Chelli, R. J. Chem. Theory Comput. 2012, 8, 4040–4052. (35) Zhang, Z.; Wu, T.; Wang, Q.; Pan, H.; Tang, R. J. Chem. Phys. 2014, 140, 034706. (36) Jarzynski, C. Phys. Rev. E 1997, 56, 5018–5035. (37) Schöll-Paschinger, E.; Dellago, C. J. Chem. Phys. 2006, 125, 054105. (38) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300–313. (39) Chipot, C., Pohorille, A., Eds. Free Energy Calculations: Theory and Applications in Chemistry and Biology; Springer, Berlin, 2007; Vol. 86. (40) McQuarrie, D. A. Statistical Mechanics; HarperCollinsPublishers, New York, USA, 1976. (41) The use of a second-order polynomial for U(q, t) is very common, even if any arbitrary functional form can be used to bound ζ(q) to the control parameter λ(t). (42) Chelli, R.; Procacci, P. Phys. Chem. Chem. Phys. 2009, 11, 1152–1158. (43) Chelli, R.; Marsili, S.; Procacci, P. Phys. Rev. E 2008, 031104. (44) Hummer, G.; Szabo, A. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 3658–3661. (45) Nicolini, P.; Procacci, P.; Chelli, R. J. Phys. Chem. B 2010, 114, 9546–9554. (46) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. N.; Teller, E. J. Chem. Phys. 1953, 21, 1087–1092. (47) Owicki, J. C.; Scheraga, H. A. Chem. Phys. Lett. 1977, 47, 600–602. 32


Page 32 of 35

Page 33 of 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


(48) Frenkel, D.; Smit, B. Understanding Molecular Simulations: From Algorithms to Applications; Academic Press: San Diego, California, 2002. (49) Tuckerman, M. E.; Liu, Y.; Ciccotti, G.; Martyna, G. J. J. Chem. Phys. 2001, 115, 1678–1702. (50) Bennett, C. H. J. Comput. Phys. 1976, 22, 245–268. (51) Oberhofer, H.; Dellago, C.; Geissler, P. L. J. Phys. Chem. B 2005, 109, 6902–6915. (52) Shirts, M. R.; Pande, V. S. J. Chem. Phys. 2005, 122, 144107. (53) Gore, J.; Ritort, F.; Bustamante, C. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 12564– 12569. (54) More exactly, these are the particles having a distance from the center of the MR smaller than rcut = 6; see Eq. 4. (55) The simulations have been performed on samples of 10200 particles with the setup specified above. (56) In the ST simulation, the reference point for g(s) is the center of mass of the dimer particles, whose distance is restrained around the value of λA (compact dimer configuration). (57) In condensed-phase systems, equilibrium could not be reached in times typical of a molecular dynamics simulations. (58) Zerbetto, M.; Frezzato, D. Phys. Chem. Chem. Phys. 2015, 17, 1966–1979. (59) Tribello, A. G.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. Comput. Phys. Commun. 2014, 185, 604–613. (60) Marx, D.; Hutter, J. Ab initio molecular dynamics; Cambridge University Press, New York, 2009; pp 267–285. 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

(61) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999–3093. (62) Cammi, R. Molecular response functions for the polarizable continuum model ; Springer, 2013. (63) Friedman, H. L. Mol. Phys. 1975, 29, 1533–1543. (64) Brancato, G.; Rega, N.; Barone, V. J. Chem. Phys. 2006, 124, 214505. (65) Mancini, G.; Brancato, G.; Barone, V. J. Chem. Theory Comput. 2014, 10, 1150–1163. (66) Ewald, P. Ann. Phys. 1921, 64, 253–287. (67) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J Chem. Phys. 1995, 103, 8577–8593.

34


Page 34 of 35

Page 35 of 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


TOC graphic

35


A Way To Speed Up Nonequilibrium Molecular ... - ACS Publications

Recommend Documents