Constrained Isothermal-Isobaric Molecular Dynamics with Full Atomic Virial † Giovanni Ciccotti Centre Europe ´ en de Calcul Atomique et Mole ´ culaire, 46 Alle ´ e d’Italie, 69007, Lyon, France, and INFM, Department of Physics, UniVersita ´ La Sapienza, P. za A. Moro, 5, 00185, Rome, Italy Glenn J. Martyna Department of Chemistry, Indiana UniVersity, Bloomington, Indiana 47405-4001 Simone Melchionna Department of Chemistry, UniVersity of Cambridge, Lensfield Road, Cambridge CB2 1EW, U.K. Mark E. Tuckerman* Centre Europe ´ en de Calcul Atomique et Mole ´ culaire, 46 Alle ´ e d’Italie, 69007, Lyon, France, and Department of Chemistry and Courant Institute of Mathematical Sciences, New York UniVersity, New York, New York 10003 ReceiVed: February 15, 2001; In Final Form: May 7, 2001 We consider the isobaric-isothermal molecular dynamics method in a system subject to a set of holonomic constraints with a full atomic description of the virial. By applying the non-Hamiltonian statistical mechanical theory recently developed by Tuckerman et al. [Tuckerman, M. E.; Mundy, C. J.; Martyna, G. J. Europhys. Lett. 1999, 45, 149], the Kneller-Mu ¨lders equations [Kneller, G. R.; Mu ¨lders, T. Phys. ReV.E 1996, 54, 6825] are analyzed, and it is determined that they sample the desired ensemble only under certain circumstances. In general, a bias arising from the conservation of total momentum needs to be corrected. Thus, a new set of equations of motion is presented and the phase space generated by these equations is shown to be correct under all circumstances. 1. Introduction The macroscopic properties of a condensed matter system can be described in different sets of variables that characterize the thermodynamic state. For example, a one component system subject to the action of a piston and in contact with a heat bath is described thermodynamically by the pressure, temperature and number of particles. The system, then, samples states in the isothermal-isobaric, or NpT, ensemble. The isothermal-isobaric ensemble best reflects the conditions under which condensed phase experiments are often performed. Thus, it is important to be able to mimic these conditions in molecular dynamics (MD) simulations of finite-size systems. The advantages of employing the NpT ensemble can be summarized in the following three points. First, contact with experiments can be made directly. Second, the density and energy of the system adjust automatically so that the internal pressure and temperature match the corresponding external values. These conditions are very difficult to control in other ensembles. Third, some internal metastabilities, due to the finite size of the simulation cell, are removed due to the presence of volume fluctuations. This fact leads to a speed-up of calculations involving global spatial rearrangements and phase transitions. In past years, therefore, a number of methods have been proposed in the MD literature to simulate systems in the NpT ensemble. 1,2,5-7 For MD simulations in the canonical and isothermal-isobaric ensembles, use is often made of equations of motions not directly derived from a Hamiltonian. Examples include the so- called Nose ´ -Hoover (NH), 8,2 Nose ´ -Hoover chain (NHC) 3 and generalized Gaussian moment 4 equations of motion, and the equations of motion of refs 1, 2, and 5-7. In all of these examples, the phase space is extended to include a new set of dynamical variables in addition to the usual (r,p) variables, which act to control the temperature (“thermostat” variables) and/or pressure (“barostat” variables) fluctuations in the system. The thermostat variables act on the system momenta, while the barostat variables act on the system positions and momenta. The use of non-Hamiltonian dynamics poses the problem of describing the statistical distribution sampled by the dynamical variables in the full phase space. The statistical behavior of non- Hamiltonian systems cannot be derived simply using the standard text-book treatment employed in the Hamiltonian case. In fact, as Tuckerman, Mundy, and Martyna have recently shown, to define an invariant measure in phase space, one needs to introduce a Riemannian manifold with an associated metric. 9 In other words, the phase space of a non-Hamiltonian dynamical system is viewed as a manifold with a local curvature, and phase space has structure with global curvature. In the Hamiltonian case, the phase space is clearly flat. In MD simulations, it is often necessary to impose a set of holonomic constraints on a system. This may be done for the purpose of freezing a subset of internal degrees of freedom, such as bond lengths or bend angles, in order to eliminate from the dynamics, 10 high-frequency modes which can limit the integration time step. Moreover, constraints are also employed † Part of the special issue “Bruce Berne Festschrift”. 6710 J. Phys. Chem. B 2001, 105, 6710-6715 10.1021/jp010601s CCC: $20.00 © 2001 American Chemical Society Published on Web 06/23/2001