From Classical to Quantum and Back: A ... - ACS Publications

Article pubs.acs.org/JCTC

From Classical to Quantum and Back: A Hamiltonian Scheme for Adaptive Multiresolution Classical/Path-Integral Simulations Karsten Kreis,†,‡ Mark E. Tuckerman,*,§,∥,⊥ Davide Donadio,†,# Kurt Kremer,† and Raffaello Potestio*,† †

Max Planck Institute for Polymer Research, Ackermannweg 10, 55128 Mainz, Germany Graduate School Materials Science in Mainz, Johannes Gutenberg University Mainz, Staudinger Weg 9, 55128 Mainz, Germany § Department of Chemistry, New York University (NYU), New York, New York 10003, United States ∥ Courant Institute of Mathematical Sciences, NYU, New York, New York 10012, United States ⊥ NYU−East China Normal University Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China # Department of Chemistry, University of California at Davis, One Shields Avenue, Davis, California 95616, United States ‡

ABSTRACT: Quantum delocalization of atomic nuclei affects the physical properties of many hydrogen-rich liquids and biological systems even at room temperature. In computer simulations, quantum nuclei can be modeled via the pathintegral formulation of quantum statistical mechanics, which implies a substantial increase in computational overhead. By restricting the quantum description to a small spatial region, this cost can be significantly reduced. Herein, we derive a bottom-up, rigorous, Hamiltonian-based scheme that allows molecules to change from quantum to classical and vice versa on the fly as they diffuse through the system, both reducing overhead and making quantum grand-canonical simulations possible. The method is validated via simulations of low-temperature parahydrogen. Our adaptive resolution approach paves the way to efficient quantum simulations of biomolecules, membranes, and interfaces.

1. INTRODUCTION Nuclear quantum delocalization plays a crucial role in lowtemperature systems, for example, helium or hydrogen,1−5 which can undergo a superfluid transition, and it affects in nontrivial ways a large variety of processes at more standard thermodynamic conditions. This is the case for proton transfer in biomolecules and membranes and in DNA oxidation,6−12 the thermodynamics of ice,13 the structure of water adlayers on catalysts,14,15 the structure and dynamics of bulk water at room temperature,16−19 and aqueous proton and hydroxide transport.20−25 To account for these effects in computer simulations, Feynman’s path-integral (PI) formulation of quantum statistical mechanics2,26,27 can be employed, which enables an accurate description of quantum delocalization via Monte Carlo (MC) or molecular dynamics (MD) simulations,2,27 albeit at an increased computational cost. A strategy to overcome this limitation is to restrict the PI description of the atoms to a (small) spatial region where quantum effects are important and to model the remaining atoms as classical particles interacting via an appropriately chosen effective potential. Molecules diffusing across the boundary separating these two regions change their representation “on the fly” from classical to quantum and vice versa. This approach, which is obviously viable only for sufficiently short De Broglie wavelengths, is © 2016 American Chemical Society

beneficial particularly when the region that needs to be modeled quantum mechanically is small compared to the size of the full system. Examples of such systems include liquid−solid or liquid−liquid interfaces,28−30 aqueous solutions,20−25 and proteins.31−33 The simplified model in the classical region also allows the number of molecules in the system to change during the simulation,34,35 thereby generating a grand-canonical distribution. More generally, an approach in which a classical and a PI model are concurrently used would allow a substantial computational gain and enable the simulation of larger systems for significantly longer sampling times. A first step in this direction was taken in the framework of the adaptive resolution simulation (AdResS) scheme36−38 by merging quantum and classical effective forces.39−42 These works and subsequent advancements43,44 demonstrated the possibility of investigating the properties of a system of light particles by explicitly considering their quantum nature locally without affecting the overall thermodynamic balance between quantum and classical regions. The AdResS scheme, however, is based on an interpolation of forces that does not admit a global Hamiltonian structure;45 the PI formulation, on the contrary, relies on a configurational potential energy. Hence, a Received: March 4, 2016 Published: May 23, 2016 3030

DOI: 10.1021/acs.jctc.6b00242 J. Chem. Theory Comput. 2016, 12, 3030−3039

Article

Journal of Chemical Theory and Computation theoretically rigorous formulation of a hybrid quantum− classical system should be based on a Hamiltonian function defined everywhere in the simulation domain. In this paper, we provide a quantum−classical coupling protocol via a global Hamiltonian, allowing us to treat a system of particles quantum mechanically only in a restricted region of space and classically everywhere else. Additionally, the particles can freely diffuse across the simulation domain and switch the nature of their interactions according to their position in space. The validation of the proposed approach is carried out by means of MC simulations. We have organized the paper as follows: In section 2, we review the classical Hamiltonian version of AdResS, denoted by H-AdResS, and show how to extend it to a quantum system with variable mass. The quantization of such a system leads to a condition on the minimal size of the hybrid region. In section 3, we provide details of our validation studies, including system setup, free energy compensation terms, an MC algorithm for sampling the PI H-AdResS distribution, and specific simulations that demonstrate the performance of the method. In section 4, we discuss the speedup of the new approach over a full PI calculation. In section 5, we place the PI H-AdResS scheme in the context of, and compare it to, other methods designed to reduce the cost of PI calculations. Finally, in section 6, we present a brief summary and conclude the article.

Figure 1. Illustration of the simulation setup for quantum−classical simulations. Red and blue correspond to a larger and a smaller radius of gyration, respectively. The green line shows the interpolation function λ, smoothly changing from 0 in the classical region to 1 in the quantum region. The smooth transition from extended to collapsed molecules demonstrates the transition from quantum mechanical to classical behavior. The particles freely move between the regions and change their descriptions accordingly.

2. METHODOLOGY Here, we make use of the H-AdResS method,46−48 which was developed to perform adaptive resolution MD/MC simulations on the basis of a global Hamiltonian and is thus the most suitable framework to correctly model each subsystem with appropriate interactions. A typical H-AdResS system is partitioned into two regions connected via a hybrid buffer region. The resolution of a particle depends on the value of a representative coordinate R and is parametrized by a continuous function λ(R) smoothly switching from 0 to 1 in the hybrid region (see Figure 1). The total potential energy of each molecule is obtained by interpolating the two resolutions. The classical H-AdResS Hamiltonian H of a system of pointlike particles reads

harmonic springs of frequency ωP = P /β ℏ, where β = 1/ kBT. The exact quantum behavior is recovered only for P → ∞; however, sufficiently accurate results are obtained with finite values of P, which typically range from 16 to 48 beads for standard PI simulations.11,19,49−53 Owing to the H-AdResS setup, the interaction between the ring polymers changes adaptively in space according to the potential energy interpolation as shown in eq 1. The quantum behavior, dictated by the strength of the springs connecting the beads of each ring, is nevertheless the same everywhere. A strategy for switching between quantum and classical descriptions is to modify the masses of the atoms as larger masses correspond to stiffer springs with constant mω2P; a large mass causes the ring polymers to collapse, and the particles approach their classical limit. We thus define m → μ(λ) = λm + (1 − λ)M, where μ(λ) smoothly switches from a mass μ(0) = M to a mass μ(1) = m ≪ M and λ = λ(Rα) with Rα being the particle positions. When μ is equal to the physical mass m, the particles are light and the quantum zero-point motion becomes important. The mass M should be chosen large enough so that the particles become essentially classical. Assigning the masses in this way, which causes them to become position dependent, introduces the problem of quantizing a system with coordinate-dependent masses, and we turn to this problem now (as a reference for the regular approach for constant-mass particles, see Tuckerman27). Consider first the Hamiltonian operator for a particle of mass μ(x) in one dimension subject to a potential V (x). The approach easily generalizes to N Boltzmann particles in any number of dimensions. When the mass depends on x, we have the problem of representing the Hamiltonian as a Hermitian operator, which we can obtain by writing it in the following form: 1 /̂ = p ̂μ−1(x)̂ p ̂ + V (x)̂ (2) 2

N

H=2+

∑ [λαVα1 + (1 − λα)V α0 − ΔH(R α)] α=1

(1)

where 2 is the kinetic energy, α indexes the N particles, and λα = λ(Rα). Single-particle potentials VRes α (with Res = 0, 1) are the sum of all intermolecular potentials acting on particle α, properly normalized so that double counting is avoided.46,47 The term ΔH, referred to as the free energy compensation (FEC),46,47 is an external field acting in the hybrid transition region to eliminate the density imbalance that naturally occurs in such dual-resolution systems. In fact, different models of the same physical system exhibit a free energy difference that needs to be neutralized to enforce identical thermodynamical and/or structural properties (e.g., density) everywhere in the simulation domain. The FEC, whose calculation is described in section 3.2, aptly compensates for these free energy imbalances. The incorporation of the H-AdResS Hamiltonian in the PI formalism is straightforward. The effective potential energy obtained from the PI quantization of the Hamiltonian as shown in eq 1, assuming Boltzmann statistics, describes a set of N ring polymers, each containing P points or “beads” connected by 3031


Article

Journal of Chemical Theory and Computation where μ−1(x̂) is the inverse mass operator and p̂ is the

⎛ β ⎞ ⟨xk|exp⎜ − p ̂μ−1(x)̂ p ̂⎟|xk + 1⟩ ⎝ 2P ⎠

momentum operator. Using this Hamiltonian, we seek to formulate the partition function Q = Tr[exp{−β /̂ }] as a PI.

∞

=

Introducing the usual set of P − 1 resolutions of the identity

⎡ dμ−1 β⎛ exp⎢ − ⎜p2 μ−1(xk + 1) + iℏp ⎢ 2P ⎝ dx ⎣

operator and the Trotter factorization of the Boltzmann operator, we can write the trace in the coordinate basis as P

Q = lim

P →∞

⎛

⎞

∫ dx1⋯dxP∏ ⟨xk|exp⎝− Pβ V (x)̂ ⎠ ⎜

⎡ β ℏ2 dμ−1 ⎢(xk − xk + 1) − 2P dx ⎣

xP + 1= x1

P

= lim

P →∞

∫

⎛ β ⎞ dx1⋯dxP∏ exp⎜ − V (xk)⎟⟨xk| ⎝ ⎠ P k=1

⎛ β ⎞ exp⎜ − p ̂μ−1(x)̂ p ̂⎟|xk + 1⟩ ⎝ 2P ⎠

To derive the matrix elements in eq 3, we introduce the momentum identity resolution as

⎛ dμ−1 ⎞ ⎟ ⎜ ⎝ dx ⎠ x

⎛ β ⎞ ⟨xk|exp⎜ − p ̂μ−1(x)̂ p ̂⎟|xk + 1⟩ ⎝ 2P ⎠ ⎛

∞

=

⎜

⎟

⎛ dμ−1 ⎞ ⎟ ⎜ ⎝ dx ⎠ x

(4)

Given that the limit P → ∞ is ultimately taken, we can work with an infinitesimal version of the exponential operators by

⎛ dμ ⎞ ⎜ ⎟ ⎝ dx ⎠ x

⎛ β ⎞ ⎛ ⎞ β −1 p ̂μ (x)̂ p ̂⎟|x⟩ ⟨p|exp⎜ − p ̂μ−1(x)̂ p ̂⎟|x⟩ ≈ ⟨p|⎜1 − ⎝ 2P ⎠ ⎝ ⎠ 2P

Now, we introduce the commutator [μ−1(x̂), p̂] and write

(9)

⎛ dμ ⎞ 1 ⎜ ⎟ μ2 (xk + 1) ⎝ dx ⎠xk +1

= k+1

≪ k+1

⟨Δx⟩ ≡

μ (x)̂ p ̂ = p ̂μ (x)̂ + [μ (x), ̂ p ̂]

(10)

2Δxk , k + 1 Λ2μ(xk + 1)

μ(xk + 1) (11)

β ℏ2 /(Pμ(x)) . Because

1 ∑ Δxl2,l+ 1 P l=1

= Λμ (P − 1)/P ≈ Λμ (12)

for a free ring of constant mass μ and typical values of P, we can approximate Δxk,k+1 ≈ Λμ(xk+1) and write

(6)

dμ(x) 2μ(x) ≪ dx Λμ(x)

Substituting eq 6 into eq 5 yields ⎛ iℏβ dμ−1 ⎞ β 2 −1 p ̂ μ (x)̂ − p̂ ⟨p|⎜1 − ⎟| x ⟩ 2P 2P dx ̂ ⎠ ⎝ ⎛ β p 2 −1 iℏβ dμ−1 ⎞ ⎟ p μ (x ) − = ⟨p|x⟩⎜1 − 2P 2P dx ⎠ ⎝ ⎡ β⎛ dμ−1 ⎞⎤ ≈ ⟨p|x⟩exp⎢ − ⎜p2 μ−1(x) + iℏp ⎟⎥ ⎢⎣ 2P ⎝ dx ⎠⎥⎦

β ℏ2

k+1

P

−1

dμ−1 = p ̂μ−1(x)̂ + iℏ dx ̂

2P Δxk , k + 1

using the definition Λμ(x) ≡

(5)

−1

≪

the condition becomes

expanding the exponential to first order. Thus, we obtain

−1

(8)

where we defined Δxk,k+1 = |xk − xk+1|. Because

⎞

∫−∞ dp⟨xk|p⟩⟨p|exp⎝− 2βP p̂μ−1(x)̂ p̂⎠|xk+1⟩

⎤2 ⎫ ⎪ ⎥⎬ ⎦⎪ xk + 1 ⎭

for the matrix elements in eq 3, where the last equality has been obtained by introducing the matrix elements ⟨x|p⟩ = exp(ipx /ℏ)/ 2π ℏ and performing the momentum integration by completing the square. From eq 8, we see that the inverse mass derivative term can be neglected if the following condition holds:

(3)

xP + 1= x1

⎞⎤ ⎟⎥ ⎠⎥ xk + 1 ⎦

⎛ μ(xk + 1)P ⎞1/2 ⎧ ⎪ βμ(xk + 1)P exp⎨− =⎜ 2 ⎟ ⎪ 2(β ℏ)2 ⎝ 2πβ ℏ ⎠ ⎩

⎟

k=1

⎛ β ⎞ exp⎜ − p ̂μ−1(x)̂ p ̂⎟|xk + 1⟩ ⎝ 2P ⎠

∫−∞ dp⟨xk|p⟩⟨p|xk+1⟩

(13)

for an arbitrary position x. The inequality in eq 13 places a lower bound on the width of the hybrid region (where λ(x) ≠ const): the interpolation must be sufficiently smooth so that the mass derivative can be safely neglected. The derivation and the criterion generalize to an arbitrary number N of interacting particles in three dimensions. In fact, typical potentials do not significantly change the radius of gyration of the ring polymers; additionally, it is necessary to consider only the components of the bead−bead distances parallel to the direction in which the mass changes. The size of the hybrid region compatible with eq 13 is still generally small

(7)

where the operators are now replaced by the corresponding eigenvalues. Substituting eq 7 into eq 4 gives 3032


Article

Journal of Chemical Theory and Computation (1−2 nm in our simulations for parahydrogen at 20 K, even smaller at room temperature), thus having no significant impact on the computational overhead. Introducing the H-AdResS potential energy, we obtain the following partition function for N interacting Boltzmann particles in three dimensions: 3/2 ⎡ P N ⎛ mP ⎞ Q = lim ⎢ ∏ ∏ ⎜ ⎟ 2 P →∞⎢ ⎣ k = 1 α = 1 ⎝ 2πβ ℏ ⎠

∫

⎤ μ drα , k⎥e−βVP ⎥⎦

The mass m is set to the molecular hydrogen mass mH2 = 2.001 au. In the classical (CL) region, the increased mass is chosen as M = 100 mH2. In the quantum region, we employ the Silvera− Goldman potential56,57 with a cutoff at 0.9 nm for the intermolecular interaction potential V1, but in the classical region, we make use of a shifted, purely repulsive Weeks− Chandler−Andersen (WCA) potential:58 ⎧ ⎡⎛ 12 ⎤ ⎛ σ ⎞6 σ ⎞ 1⎥ ⎪ ⎢ ⎪ ϵ − + 4 : r ≤ Rc ⎜ ⎟ ⎜ ⎟ V 0(r ) = ⎨ ⎢⎣⎝ r − r0 ⎠ 4 ⎥⎦ ⎝ r − r0 ⎠ ⎪ ⎪0 : r > Rc ⎩

(14)

with 2 ⎧ μα , k ⎪ μα , k ωP 3 ⎨ log =∑∑⎪ |rα , k − rα , k + 1|2 − 2 2β m k=1 α=1 ⎩ P

VPμ

+

N

⎫ ⎪ 1 [λα , k Vα1 , k + (1 − λα , k )V α0, k − ΔH(rα , k)]⎬ ⎪ P ⎭

(17)

where r = |rα,k − rβ,k|, (α ≠ β) denotes the distance between beads of the same imaginary time slice k in different molecules α and β. Furthermore, we choose ϵ = 1.0 kJ/mol, σ = 0.14 nm, and r0 = 0.15 nm. The cutoff is given by Rc = 21/6σ + r0. The two potentials are graphically presented in Figure 2. The WCA

(15)

and μα,l = μ(rα,l), where the representative coordinates Rα,k in the functions λα,k are the bead positions rα,k. In eq 15, the position-dependent prefactor has been explicitly introduced in the potential VμP as a logarithmic function of the bead masses so that it can be treated as a conventional energy term and fully removed from the Hamiltonian by means of the FEC function ΔH in eq 1.54 The light mass m has been used as the reference mass scale. The ring polymers described by the energy function VμP (eq 15) are expanded in the region where the mass is small and collapse to nearly classical point-like particles in the largemass region.

3. VALIDATION To validate the proposed quantum−classical coupling scheme, adaptive PI MC simulations of liquid parahydrogen at 20 K with P = 16 are performed. Liquid hydrogen between 14 and 25 K exhibits pronounced quantum behavior;4,5,55 therefore, it constitutes a probing test case. A Trotter number of P = 16 is likely not large enough for full convergence of physical observables in this system under the aforementioned conditions; however, our focus is less on the exact properties of parahydrogen and more on the validation of the proposed PI H-AdResS approach. 3.1. System Setup. We consider a system composed of 1458 hydrogen molecules in a slab of dimensions of 11.0 nm × 2.5 nm × 2.5 nm (molecular density 28.4 cm3/mol) with periodic boundary conditions in all directions. The width of the low-mass quantum (QM) region is set to dQM = 3.0 nm, and the thickness of each hybrid (HY) transition region is dHY = 2.0 nm. To assign to a bead its position-dependent resolution λ, its distance from the boundary between the quantum and the hybrid region is computed, that is, δxα,k = |xα,k| − dQM/2, where xα,k denotes the X coordinate of the bead in a coordinate system with its origin at the center of the simulation box. This quantity is then employed in the resolution function λ(δxα,k). This function is given generically by ⎧1 :x≤0 ⎪ ⎪ 2⎛ π x ⎞ ⎟ : 0 < x < dHY λ(x) = ⎨ cos ⎜ ⎝ 2 dHY ⎠ ⎪ ⎪0 : x ≥ dHY ⎩

Figure 2. Nonbonded intermolecular interaction potentials used in the adaptive quantum−classical simulations. The red curve is the Silvera− Goldman potential, which is employed in the low-mass quantum region. The blue curve shows the shifted WCA potential, which is used in the high-mass classical region.

potential in the classical region is not interpreted as a classical model for low-temperature parahydrogen. Rather, it is parametrized to reproduce approximately only the hard-core radius of the reference quantum particles and to serve as a crude model of particles occupying a reservoir with which the quantum region can exchange matter. We purposely avoid fitting the classical potential to the structure of the reference to demonstrate the generality of the protocol; indeed, the potential could be chosen arbitrarily, thereby allowing various options and applications. One could, for example, focus on reproducing specific properties in the classical region and/or employ a computationally advantageous model such as an ideal gas (noninteracting particles).59 The mass of the particles in the classical region is 100 times larger than that of the molecules in the quantum region, and the chosen set of parameters also satisfies eq 13. Finally, we stress that in the classical region the WCA interaction between ring polymers is computed using only the centers of mass of the rings, thus gaining an effective reduction of the computational

(16) 3033


Article

Journal of Chemical Theory and Computation cost. This simplification is allowed by the essentially point-like structure of the rings in the classical regions, as can be seen from the radius of gyration profile (Figure 3). The number of

Figure 4. Normalized density profiles obtained when a correction is applied in the hybrid region based on KTI only (blue) and when applying the iterative protocol described in the text (red).

FEC with the radius of gyration profiles (see Figure 3). Indeed, the density profile is mostly distorted in the area where the radius of gyration of the ring polymers changes more steeply, that is, where the local environment of molecules changes sharply. Therefore, to remove also the remaining fluctuations in the obtained density profile after applying the KTI-based FEC, an iterative approach similar to the one devised by Fritsch et al. is employed.38 The scheme reads ΔH i + 1(x) = ΔH i(x) − Δhi(x)

(18)

with ΔH (x) = ΔHKTI(x), the initial correction term obtained by KTI converted from a function of resolution into a function of position. Furthermore, 0

Figure 3. Top: Radius of gyration rg and the normalized density ρ profiles. The blue line corresponds to the radius of gyration of molecules in a fully quantum reference simulation. Bottom: Ring polymers’ center of mass (CoM) and bead−bead RDFs of the quantum−classical simulation calculated in the quantum region and of a fully quantum reference simulation. The black curve is the RDF of a fully classical (P = 1) system of particles interacting via the Silvera− Goldman potential.

⎧ 1 i ⎪− ln{ρ (x0)} : x ≤ x0 β ⎪ ⎪ 1 i i Δh (x) = ⎨− ln{ρ (x)} : x0 < x < x1 β ⎪ ⎪ 1 i ⎪− β ln{ρ (x1)} : x ≥ x1 ⎩

computations per pair of molecules is reduced from P = 16 to P = 1. A snapshot of simulations is shown in Figure 1, which shows the gradual change in size of the ring polymers, indicating a transition between classical and quantum mechanical behavior. 3.2. Free Energy Compensation. To modulate the thermodynamic imbalance between the classical high-mass and the quantum low-mass subsystems, an FEC is applied.46,47 First, we compute the compensation term ΔHKTI(λ) by Kirkwood thermodynamic integration (KTI)60 of a smaller system of 360 molecules in a box with dimensions of 2.570 nm × 2.570 nm × 2.570 nm. However, using only the correction term derived by KTI is insufficient to reach a flat density profile (see Figure 4). The interpolation between subsystems with very different masses and intermolecular interaction potentials leads to large thermodynamic imbalances and strong forces pushing molecules from one region to the other. The KTI-based FEC alone is incapable of sufficiently correcting these imbalances between the classical and the quantum domains because it provides only a mean field estimate of the necessary compensation. This becomes especially clear when comparing the density profile obtained when applying only the KTI-based

(19)

where x0 =

dQM 2

− 0.5 nm

(20)

dQM

+ dHY + 0.5 nm (21) 2 and x is measured as the absolute distance from the center of the simulation box. Here, ρi(x) denotes the density profile obtained when applying the FEC ΔHi(x) of iteration i [therefore, ρ0(x) = ρKTI(x), where ρKTI(x) corresponds to the initial density profile obtained from simulations in which only the KTI-based FEC term is applied]. The term Δhi(x) acts in an extended hybrid region, defined as the regular hybrid region plus additional 0.5 nm within both the quantum and the classical regions. As we have two hybrid transition regions in the presented simulation setup above, averages over both regions are used to derive one consistent and symmetric correction, which is applied in the same way in both regions. The protocol converges by construction when a flat density profile is achieved. However, to speed up the convergence, on x1 =

3034


Article

Journal of Chemical Theory and Computation

The different values for the σi’s of all simulations are presented in Table 1. In the dual-resolution simulation, they are

top of the iterative density-based correction, we also add a linear ramp term i (x ) Δhramp

Table 1. Widths of the Gaussian Distributions Employed To Draw the Random Displacements and Rotations for the KTI, for the Reference Simulations, and for the Adaptive Quantum−Classical Simulations

⎧0 : x ≤ x0 ⎪ = ⎨−(ai /β) × (x − x0)/(x1 − x0) : x0 < x < x1 ⎪ − (a / β ) : x ≥ x1 ⎩ i (22)

for some of the later iterations when the overall density difference between the two regions has become small. Note that this additional term is not strictly required but can help to speed up the iterative procedure. The values of the coefficients ai that have been employed here to favor convergence are reported in the following section. The protocol is executed in an iterative fashion until a sufficiently flat density profile is obtained (see Figure 4). Note that the converged profile is also significantly smoother and more uniform than the ones obtained with a related approach by Agarwal and Delle Site.42 3.3. Monte Carlo Sampling. To sample the system’s configuration space, we employ a standard Metropolis MC algorithm.27 For the KTI of the small system, we run 16 simulations with 105 sweeps each. The λ parameter increases linearly every sweep by 10−5. The results are averaged after the simulations. Employing the KTI FEC term thus obtained, we then run 14 iterations of simulations applying the protocol set out above to refine the density profile. The parameters ai for the ramp term in the correction protocol are chosen as a7 = 0.06, a8 = 0.04, a9 = 0.05, a10 = 0.05, a11 = 0.02, a12 = 0.03, and a14 = 0.02. For all other iterations, no ramp term is added, that is, ai = 0. Each iteration consists of 32 parallel simulations, each of these running 1.5 × 104 equilibration sweeps and another 7.5 × 103 sweeps during which the density profile is measured. Also here, after each iteration, the results are averaged. Having reached a sufficiently smooth density profile, we then utilize the FEC from the KTI and the iterative protocol to perform the main production simulations. For these, we perform 32 simulations in parallel, each running 1.5 × 104 sweeps after another 1.5 × 104 equilibration sweeps. Afterwards, the results [i.e., the radial distribution functions (RDFs), the density profiles, the radius of gyration profiles, and the particle number histograms] are once again averaged over all simulations. Each sweep is constituted by N × P attempted MC moves on randomly chosen molecules, with N being the total number of molecules and P being the Trotter number (N = 1458 and P = 16 in the production-run simulations). Three different kinds of moves are randomly performed: Whole-molecule displacements: The chosen molecule is displaced as a whole by moving its center of mass. The direction is chosen randomly from a uniform spherical distribution, and the distance is drawn from a Gaussian distribution with zero mean and width σCoM. Molecule rotations: The chosen molecule is rotated as a whole around a randomly oriented axis passing through its center of mass. The angle is chosen randomly from a Gaussian distribution with zero mean and width σrot. Individual Trotter-bead moves: An individual bead of the molecule is randomly chosen and displaced. The direction is chosen randomly from a uniform spherical distribution, and the distance is drawn from a Gaussian distribution with zero mean and width σbead.

simulation

σCoM (nm)

σrot (rad)

σbead (nm)

KTI adaptive simulation quantum reference classical reference

0.1 0.1 0.1 0.1

0.5 0.5 0.5

0.03 0.03 0.07

chosen such that they result in adequate acceptance ratios for the moves both in the classical high-mass and in the quantum low-mass regions. When picking a molecule for an MC move, the probabilities for performing whole-molecule displacements or -molecule rotations are 1 each, whereas the probability for 13

Trotter-bead moves is 11 . This choice leads to a convenient 13 balance between whole-molecule motions and Trotter-bead fluctuations. 3.4. Reference Simulations. To evaluate the results of the adaptive quantum−classical simulations, we perform fullquantum as well as full-classical reference simulations of liquid parahydrogen for comparison. The same setup is used for the adaptive simulations, and likewise the temperature is set to T = 20 K and the Silvera−Goldman potential is employed. For the full-quantum simulations, we choose P = 16 as in the adaptive simulations, whereas the classical simulations are performed with P = 1. In both cases, 32 simulations are run in parallel, each one for 1.5 × 104 equilibration sweeps and another 1.5 × 104 production sweeps during which properties are computed. Afterwards, the results are averaged. The values used for the σi’s in the reference simulations are presented in Table 1. In the classical simulations, all moves are as described above with the obvious exception of bead and rotating moves, which do not exist for individual classical particles. The calculations of the particle-number fluctuations in the quantum simulations are performed in exactly the same manner and using the same area as those in the adaptive simulations. The RDFs of the classical and quantum reference simulations are calculated in the complete simulation domains. 3.5. Results. The radius of gyration of the ring polymers in the quantum region, shown in Figure 3, perfectly reproduces that of a corresponding full-quantum simulation, whereas it drops in the classical region by ≈90%, indicating nearly classical behavior. At the same time, the FEC enforces a uniform density profile everywhere in the system. A quantitative measure of the fluid structure is provided by the ring’s center of mass and bead−bead (intermolecular bead pairs with the same Trotter index) RDFs. In spite of the remarkable differences between the quantum fluid and the classical model, the RDFs measured in the quantum region perfectly match those from the fullquantum reference simulation. Furthermore, the particle-number fluctuations in the inner quantum region (the subdomain of the quantum region where particles are 1.0 nm far from the hybrid region, that is, a distance larger than the 0.9 nm cutoff) are calculated and compared to the same quantity in the full-quantum simulations (see Figure 5). The fluctuations in the adaptive simulation 3035


Article


total volume ratio, for example, employing a small spherical quantum region in a large, three-dimensional box, the computational gain could be close to Smax. Such a situation could, for example, occur in studies of proton-transfer reactions in membranes, where only a small part of the membrane needs to be described quantum mechanically, but the remainder of the system could be treated via a simpler classical model. Similarly, in simulations of enzymes, it is likely sufficient to treat a large region of the molecule classically in a computationally efficient manner, whereas quantum nuclear effects can be taken into account at the active site where they play a significant role.33 The method also enables efficient sampling of a quantum grand-canonical ensemble by embedding a quantum subsystem in a particle reservoir, which could be significantly larger than the quantum domain to allow for realistic particle fluctuations in the latter. At the same time, the structure and dynamics of the reservoir, being unimportant, could be described by a much coarser and numerically more efficient classical model. In all such cases, an adaptive quantum−classical coupling scheme, as proposed in this work, would provide a significant advantage. To demonstrate the improved computational efficiency, we perform four sets of simulations for different box sizes, with each set consisting of a full-quantum and an adaptive simulation in which the quantum region as well as both hybrid regions have a width of 2.0 nm each. The total box sizes as well as the corresponding molecule numbers for the simulations are presented in Table 2. The temperature and density are the

Figure 5. Particle-number probability distribution of the inner part of the quantum region in the adaptive simulation (blue crosses) and of the same volume in a reference full-quantum simulation (red circles). The black curve is a Gaussian fit to the latter (μ = 132.4, σ = 4.383).

match the reference closely. This means that the quantum subdomain exchanges particles with its surroundings in the same manner as it would if embedded in a full-quantum system. Therefore, by coupling to a sufficiently large particle reservoir, the scheme can be employed to reliably and rigorously simulate the quantum grand-canonical ensemble with only a little additional computational cost due to the numerically efficient classical particle reservoir.

4. SPEEDUP OVER FULL-QUANTUM SIMULATIONS As mentioned earlier, in the proposed quantum−classical coupling scheme, interactions in the classical region do not need to be calculated P times (with P being the Trotter number). Rather, because of the collapse of the polymer rings, they are computed only once between the centers of mass of the (quasi-point-like) rings. Additionally, a numerically simpler potential with a shorter cutoff can be used in the classical region. In general, the overall speedup of the simulations strongly depends on the details of the implementation of the algorithm and is therefore platform dependent. For example, if there is a high overhead in the code, then the overall computational gain by more efficient potential energy calculations will be small. If the program spends most of its time with these calculations, then a significant improvement is possible. The expected speedup S for the energy calculations of shortrange intermolecular potentials, compared to full-quantum simulations, can be estimated as PN S∼ P⟨NQM⟩ + cCLpot⟨NCL⟩ (23)

Table 2. Number of Molecules and Box Geometries for the Different Sets of Simulations for the Calculation of the Computational Gain of Adaptive Quantum−Classical over Full-Quantum Simulations number of molecules

Lx (nm)

Ly (nm)

Lz (nm)

1325 3974 6623 9272

10.0 30.0 50.0 70.0

2.5 2.5 2.5 2.5

2.5 2.5 2.5 2.5

same as before. In all cases but the first, the classical regions are significantly larger than the quantum ones. For each setup, we perform eight simulations and average the results. Furthermore, the simulations are run for 500 sweeps, and in the case of the adaptive simulations, the previously derived FEC is applied. The set of MC moves is chosen for both the quantum and the adaptive simulations to be the one previously used for the adaptive simulations. To obtain platform-independent results, we measure only the time our code spends with potential energy calculations. These times are plotted in Figure 6. Additionally, the corresponding speedups, defined as Tquantum/Tadaptive with Tadaptive (Tquantum) being the time spent for the energy calculations in the adaptive (quantum) simulations, are presented. Furthermore, we plot Vtotal/VQM, the ratio of the total box volume and the volume of the quantum plus the hybrid regions. In the latter part, all quantum mechanical energy calculations associated with the ring polymers are performed explicitly. This quantity is essentially the speedup derived above because Vtotal/VQM≈N/ ⟨NQM⟩. It can be seen that the adaptive simulations are significantly faster than their corresponding full-quantum counterparts. For the largest box, the energy calculations in the adaptive

where N is the total number of particles in the system, ⟨NQM⟩ is the average number of particles in the quantum and hybrid regions, ⟨NCL⟩ is the average number of particles in the classical and hybrid regions, and cCLpot is a constant depending on the molecular interaction potential in the classical region. If the same potentials are used in the quantum and classical subsystems, then cCLpot = 1, whereas a computationally cheaper potential in the classical region would result in a smaller cCLpot with the extreme being a noninteracting ideal gas, in which case cCLpot = 0. Hence, the speedup would approach Smax = N/ ⟨NQM⟩, its upper bound, in a good implementation of the scheme with a simple interaction in the classical region. Therefore, in simulation setups featuring a small quantum-to3036


Article


long-range interactions in PI simulations can be safely evaluated only on the ring-polymer centroids, while still yielding highly accurate results.62 Because, in our methodology, the centroid is well defined in both the classical and the quantum regions, the computation of long-range forces could then be seamlessly carried out at the cost of a classical simulation on the whole system. More recently, it was shown that the ring-polymer contraction scheme can also be extended to ab initio density functional theory simulations.63,64 We note that reactive systems often require the use of electronic structure methods in the reactive region, and these can be combined with PI approaches, as was done by Wang et al.33 Combining PI HAdResS with such an approach would require an ab initio version of the scheme, which has not yet been developed but is the subject of future work. We therefore defer the discussion of quantum mechanics/molecular mechanics (QM/MM) approaches within H-AdResS for such a future publication. A closely related method to the ring-polymer contraction scheme is the mixed time slicing approach,65 in which different degrees of freedom, that is, particles, are modeled with a different number of Trotter beads to incorporate nuclear quantum effects only for the relevant light atoms. This technique could also be combined with our work resulting, for example, in a scheme where different Trotter bead numbers are used for different atoms in the quantum region, whereas in the classical domain all PIs would collapse and behave completely classically. Another method for performing efficient PI simulations was introduced by Ceriotti et al., who showed that carefully parametrized colored noise thermostats can speed up the convergence in PI simulations by allowing the use of significantly fewer required beads, thus reducing the overall computational cost.66−69 This scheme is compatible with our approach as well and could be utilized to reduce the computational complexity further within the quantum subregion. Finally, we mention higher order Trotter factorizations,70 either directly sampled via MC71 or combined with a posteriori reweighting of averages 72,73 or a cumulant expansion74 within PI MD. When a fourth-order factorization scheme is employed, P can be generally reduced by a factor of 4 while maintaining the same level of accuracy. It is important to note, however, that in dual-resolution simulation approaches this computational speedup would add to the advantage of allowing short-range or even no interactions in the lowresolution (the classical or, more generally, the coarse-grained) region, which, by itself, leads to substantial savings. A feature that differentiates our method from those mentioned here is that it is local, in that it gains its overall speedup by restricting the expensive quantum treatment to a limited region in space. The other approaches, on the other hand, are global in the sense that they reduce the computational cost associated with the PI computations themselves in a full PI simulation setup. Therefore, the proposed method can be regarded as complementary to these techniques because any of the schemes that we discussed in this section can be used to treat the PI region in the PI H-AdResS approach. This would reduce the overall computational complexity even further, which could prove especially useful for large simulations where only a very limited subregion requires a PI treatment. Finally, although we use an MC approach to validate our scheme, it is worth pointing out that the approach can also be applied in the framework of regular PI MD as well as centroid and ring-polymer MD75−77 techniques to calculate approximate

Figure 6. Simulation times required for the potential energy calculations in the quantum and adaptive simulations of slab systems having four different box lengths (Lx). Inset: Speedup of the adaptive quantum−classical simulations together with the upper bound on the speedup. The linear increase of the speedup follows the increase in size of the full box (when the size of the quantum region is kept fixed). The lines are a guide to the eye. If all box sides were to be varied in a hypothetical setup in which the quantum region were spherical, one would expect the curve to exhibit a cubic increase.

quantum−classical simulations are faster by a factor of ≈9 than those in the full-quantum simulations. Furthermore, it can be seen that the time required for the adaptive simulations stays nearly constant for the different box sizes. The reason for this is that the computational cost of the interactions between classical molecules is negligible compared to the time required for the computation of the potential energies in the quantum region. Therefore, as shown in the graph, the speedup in our simulations is indeed close to its upper bound.

5. RELATION TO OTHER PATH-INTEGRAL METHODS In this section, we describe the present methodology in the context of, and in relation to, other established techniques that also reduce the computational effort of traditional PI simulations. As a first observation, we recall from section 1 that the first attempts to couple a PI and a classical model of a liquid were performed in the framework of the force-based, nonHamiltonian AdResS method.39−42 This scheme, which is introduced af ter the system has already been quantized, provides a pragmatic coupling of the two descriptions. The methodology presented here represents an advance in the classical−PI hybrid methodology as it stems from a unique Hamiltonian function that governs the entire system and couples the different models before quantization is performed. This enables a more fundamental description of the setup and a more consistent treatment of the various parts of the system. Furthermore, the scheme presented can be seamlessly employed in tandem with other methods specifically developed to reduce the computational cost of PI simulations, as will be highlighted below. Consider, first, the ring-polymer contraction scheme by Markland and Manolopoulos.61,62 This approach leads to a significant decrease in numerical complexity by contracting lowfrequency nonbonded and long-range interactions to only a few imaginary time slices or even a single time slice. For example, 3037


Article


(17) Paesani, F.; Yoo, S.; Bakker, H. J.; Xantheas, S. S. J. Phys. Chem. Lett. 2010, 1, 2316. (18) Ceriotti, M.; Cuny, J.; Parrinello, M.; Manolopoulos, D. E. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 15591. (19) Fritsch, S.; Potestio, R.; Donadio, D.; Kremer, K. J. Chem. Theory Comput. 2014, 10, 816. (20) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601. (21) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361. (22) Tuckerman, M. E.; Marx, D.; Parrinello, M. Nature 2002, 417, 925. (23) Wu, Y. J.; Chen, H. N.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 467. (24) Berkelbach, T. C.; Lee, H.-S.; Tuckerman, M. E. Phys. Rev. Lett. 2009, 103, 238302. (25) Marx, D.; Chandra, A.; Tuckerman, M. E. Chem. Rev. 2010, 110, 2174. (26) Feynman, R.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (27) Tuckerman, M. E. Statistical Mechanics: Theory and Molecular Simulation; Oxford University Press: New York, 2010. (28) Gelfand, M. P.; Fisher, M. E. Phys. A (Amsterdam, Neth.) 1990, 166, 1. (29) Aguado, A.; Scott, W.; Madden, P. A. J. Chem. Phys. 2001, 115, 8612. (30) Schmitz, F.; Virnau, P.; Binder, K. Phys. Rev. E 2014, 90, 012128. (31) Hwang, J.-K.; Warshel, A. J. Am. Chem. Soc. 1996, 118, 11745. (32) Engel, H.; Doron, D.; Kohen, A.; Major, D. T. J. Chem. Theory Comput. 2012, 8, 1223. (33) Wang, L.; Fried, S. D.; Boxer, S. G.; Markland, T. E. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 18454. (34) Mukherji, D.; Kremer, K. Macromolecules 2013, 46, 9158. (35) Wang, H.; Hartmann, C.; Schütte, C.; Delle Site, L. Phys. Rev. X 2013, 3, 011018. (36) Praprotnik, M.; Delle Site, L.; Kremer, K. J. Chem. Phys. 2005, 123, 224106. (37) Praprotnik, M.; Delle Site, L.; Kremer, K. Annu. Rev. Phys. Chem. 2008, 59, 545. (38) Fritsch, S.; Poblete, S.; Junghans, C.; Ciccotti, G.; Delle Site, L.; Kremer, K. Phys. Rev. Lett. 2012, 108, 170602. (39) Poma, A. B.; Delle Site, L. Phys. Rev. Lett. 2010, 104, 250201. (40) Poma, A. B.; Delle Site, L. Phys. Chem. Chem. Phys. 2011, 13, 10510. (41) Potestio, R.; Delle Site, L. J. Chem. Phys. 2012, 136, 054101. (42) Agarwal, A.; Delle Site, L. J. Chem. Phys. 2015, 143, 094102. (43) Agarwal, A.; Zhu, J.; Hartmann, C.; Wang, H.; Delle Site, L. Molecular dynamics in a grand ensemble: Bergmann-Lebowitz model and adaptive resolution simulation. New J. Phys. 2015, 17 (8), 083042. (44) Agarwal, A.; Delle Site, L. Grand-Canonical Adaptive Resolution Centroid Molecular Dynamics: Implementation and application. Comput. Phys. Commun. DOI: 10.1016/j.cpc.2016.05.001, in press. (45) Delle Site, L. Phys. Rev. E 2007, 76, 047701. (46) Potestio, R.; Fritsch, S.; Español, P.; Delgado-Buscalioni, R.; Kremer, K.; Everaers, R.; Donadio, D. Phys. Rev. Lett. 2013, 110, 108301. (47) Potestio, R.; Español, P.; Delgado-Buscalioni, R.; Everaers, R.; Kremer, K.; Donadio, D. Phys. Rev. Lett. 2013, 111, 060601. (48) Kreis, K.; Donadio, D.; Kremer, K.; Potestio, R. EPL 2014, 108, 30007. (49) Tuckerman, M. E.; Marx, D. Phys. Rev. Lett. 2001, 86, 4946. (50) Habershon, S.; Markland, T. E.; Manolopoulos, D. E. J. Chem. Phys. 2009, 131, 024501. (51) Wang, L.; Ceriotti, M.; Markland, T. E. J. Chem. Phys. 2014, 141, 104502. (52) Wilkins, D. M.; Manolopoulos, D. E.; Dang, L. X. J. Chem. Phys. 2015, 142, 064509. (53) Miller, T. F., III; Manolopoulos, D. E. J. Chem. Phys. 2005, 122, 184503. (54) Potestio, R. Eur. Phys. J. B 2014, 87, 245.

quantum dynamical properties and time-correlation functions. This is the focus of an ongoing study.

6. CONCLUSIONS We have derived a protocol for concurrent multiscale PI simulations on the basis of a rigorous, bottom-up Hamiltonian formulation. In view of the reduced computational complexity in the classical subdomain, our method enables a computationally more efficient sampling of configurations compared to a full-quantum simulation. This, in turn, allows an extension of the accessible time and length scales. Additionally, the possibility to spatially switch off and on the quantum treatment of the system makes this scheme a powerful tool to investigate the role played by nuclear delocalization in soft matter.78−80 Furthermore, the technique is compatible with other approaches, alleviating the computational complexity of PI simulations. Future applications of the proposed approach are diverse and include, for example, simulations of the quantum grand-canonical ensemble as well as adaptive quantum− classical simulations of interface systems and biologically relevant systems such as membranes and proteins.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (M.E.T.) *E-mail: [email protected] (R.P.) Funding

K. Kreis is a recipient of a fellowship funded through the Excellence Initiative (DFG/GSC 266). R.P., D.D., and K. Kremer acknowledge funding from SFB-TRR146 of the German Research Foundation (DFG). M.E.T. acknowledges support from CHE-1301314 of the National Science Foundation. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Azuah, R. T.; Stirling, W. G.; Glyde, H. R.; Boninsegni, M.; Sokol, P. E.; Bennington, S. M. Phys. Rev. B 1997, 56, 14620. (2) Ceperley, D. M. Rev. Mod. Phys. 1995, 67, 279. (3) Nozières, P.; Pines, D. The Theory of Quantum Liquids, Advanced Book Classics; Pines, D., Ed.; Perseus: Cambridge, MA, 1999. (4) Scharf, D.; Martyna, G.; Klein, M. Fiz. Nizk. Temp. 1993, 19, 516. (5) Lindenau, T.; Ristig, M. L.; Gernoth, K. A.; Dawidowski, J.; Bermejo, F. J. Int. J. Mod. Phys. B 2006, 20, 5035. (6) Smirnov, A. Y.; Mourokh, L. G.; Nori, F. J. Phys.: Condens. Matter 2011, 23, 234101. (7) Jeuken, L. J. C.; Bushby, R. J.; Evans, S. D. Electrochem. Commun. 2007, 9, 610. (8) Haines, T. H. Prog. Lipid Res. 2001, 40, 299. (9) Löwdin, P.-O. Rev. Mod. Phys. 1963, 35, 724. (10) Rein, R.; Harris, F. E. Science 1964, 146, 649. (11) Perez, A.; Tuckerman, M. E.; Hjalmarson, H. P.; von Lilienfeld, O. A. J. Am. Chem. Soc. 2010, 132, 11510. (12) Jacquemin, D.; Zúñiga, J.; Requena, A.; Céron-Carrasco, J. P. Acc. Chem. Res. 2014, 47, 2467. (13) Pamuk, B.; Soler, J. M.; Ramı ́rez, R.; Herrero, C. P.; Stephens, P. W.; Allen, P. B.; Fernández-Serra, M.-V. Phys. Rev. Lett. 2012, 108, 193003. (14) Li, X.-Z.; Probert, M. I. J.; Alavi, A.; Michaelides, A. Phys. Rev. Lett. 2010, 104, 066102. (15) Nagata, Y.; Pool, R. E.; Backus, E. H. G.; Bonn, M. Phys. Rev. Lett. 2012, 109, 226101. (16) Morrone, J. A.; Car, R. Phys. Rev. Lett. 2008, 101, 017801. 3038


Article

Journal of Chemical Theory and Computation (55) Leachman, J. W.; Jacobsen, R. T.; Penoncello, S. G.; Lemmon, E. W. J. Phys. Chem. Ref. Data 2009, 38, 721. (56) Silvera, I. F.; Goldman, V. V. J. Chem. Phys. 1978, 69, 4209. (57) Silvera, I. F. Rev. Mod. Phys. 1980, 52, 393. (58) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237. (59) Kreis, K.; Fogarty, A. C.; Kremer, K.; Potestio, R. Eur. Phys. J.: Spec. Top. 2015, 224, 2289. (60) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300. (61) Markland, T. E.; Manolopoulos, D. E. J. Chem. Phys. 2008, 129, 024105. (62) Markland, T. E.; Manolopoulos, D. E. Chem. Phys. Lett. 2008, 464, 256. (63) Marsalek, O.; Markland, T. E. J. Chem. Phys. 2016, 144, 054112. (64) John, C.; Spura, T.; Habershon, S.; Kühne, T. D. Phys. Rev. E 2016, 93, 043305. (65) Steele, R. P.; Zwickl, J.; Shushkov, P.; Tully, J. C. J. Chem. Phys. 2011, 134, 074112. (66) Ceriotti, M.; Bussi, G.; Parrinello, M. Phys. Rev. Lett. 2009, 102, 020601. (67) Ceriotti, M.; Parrinello, M.; Markland, T. E.; Manolopoulos, D. E. J. Chem. Phys. 2010, 133, 124104. (68) Ceriotti, M.; Manolopoulos, D. E.; Parrinello, M. J. Chem. Phys. 2011, 134, 084104. (69) Ceriotti, M.; Manolopoulos, D. E. Phys. Rev. Lett. 2012, 109, 100604. (70) Takahashi, M.; Imada, M. J. Phys. Soc. Jpn. 1984, 53, 3765. (71) Suzuki, K.; Tachikawa, M.; Shiga, M. J. Chem. Phys. 2010, 132, 144108. (72) Jang, S.; Jang, S.; Voth, G. A. J. Chem. Phys. 2001, 115, 7832. (73) Perez, A.; Tuckerman, M. E. J. Chem. Phys. 2011, 135, 064104. (74) Poltavsky, I.; Tkatchenko, A. Chem. Sci. 2016, 7, 1368. (75) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 100, 5093. (76) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 100, 5106. (77) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2004, 121, 3368. (78) Lambeth, B. J.; Junghans, C.; Kremer, K.; Clementi, C.; Delle Site, L. J. Chem. Phys. 2010, 133, 221101. (79) Fritsch, S.; Junghans, C.; Kremer, K. J. Chem. Theory Comput. 2012, 8, 398. (80) Fogarty, A. C.; Potestio, R.; Kremer, K. J. Chem. Phys. 2015, 142, 195101.

3039


From Classical to Quantum and Back: A ... - ACS Publications

Recommend Documents