JOURNAL OF COMPUTATIONAL. PHYSICS 23, 321-341 (1977) Numerical integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes JEAN-PAUL RYCKAERT*, GIOVANNI CICCOTTI+, AND HERMAN J. C. BERENDSEN* Centre Europten de Calcul Atomique et Molhxlaire (CECAM), Bhtiment 506, Umbersit de Paris XI, 91405 Orsay, France Received July 19, 1976 A numerical algorithm integrating the 3N Cartesian equations of motion of a system of N points subject to holonomic constraints is formulated. The relations of constraint remain perfectly fulfilled at each step of the trajectory despite the approximate character of numeric- al integration. The method is applied to a molecular dynamics simulation of a liquid of 64 n- butane molecules and compared to a simulation using generalized coordinates. The method should be useful for molecular dynamics calculations on large molecules with internal degrees of freedom. 1. INTR~D~JCTI~N The method of molecular dynamics (MD), which has been widely used in the past for studying simple liquids and solids, has more recently been applied to molecular systems with internal degrees of freedom such as N, [l], H,O [2] and even C,H,, [3]. In applying the MD method three problems arise: (a) the choice of a suitable mechani- cal model, (b) the derivation of the equations of motion of the system and (c) the choice of an efficient algorithm for the numerical integration of these equations. In polyatomic molecules, the fast internal vibrations are usually decoupled from rotational and translational motions and can therefore be frozen by introducing a certain number of rigid bonds and angles in the skeleton of the molecule. For example, Nz becomes a rod, H,O a rigid triangle and C4H,, a nonrigid solid with one internal rotation [l, 2, 31. The classical way to treat such systems is in terms of generalized coordinates (Lagrange-Hamilton formalism), but as the number of internal degrees of freedom increases it rapidly becomes harder to write down explicitly the appropriate equations of motion. * Fact& des Sciences, Universite Libre de Bruxelles, Brussels, Belgium. + Gruppo Nazionale di Struttura della Materia, Consiglio Nazionale delle Ricerche, and Istituto di Fisica “G. Marconi,” Universita di Roma, Roma, Italy. *Laboratory of Physical Chemistry, The University of Groningen, Zemikelaan, NG8002 Groningen, the Netherlands. 327 Copyright Q 1977 by Academic Press. Inc. AlI rights of reproduction in any form reserved. ISSN 0021-9991328 RYCKAERT. CICCOTTI AND BERENDSEN Some time ago, Orban and Ryckaert [4] suggested the use of Cartesian equations of motion in order to describe the dynamical behavior of n-alkane chains, built up of n CH2 or CH, point units connected by n - 1 rigid bonds and n - 2 rigid angles between adjacent bonds. These equations are in fact the Lagrange equations of motion of the first kind [5], in which the forces of constraint appear explicitly; the dependence of these forces on the positions and velocities of the centers of force is obtained from the relations of constraint (see Section 2 or [4]). On integrating these equations numerically, Orban and Ryckaert obtained quite promising preliminary results. However, one difficulty remained in their approach. The constraints are satisfied exactly at some initial time, but not at later times, because in the numerical integration of a large set of coupled differential equations the computed trajectory deviates more and more from the true one as time proceeds. This is a consequence of the approximate character of the algorithm which is used, and in the present case means that the numerical values of the constrained bond lengths and angles gradually depart from their original values. It follows that after a sufficiently long time the character of the system is strongly modified. In principle, the time step can be reduced sufficiently to obtain an acceptable discrepancy in the constraints after the total time of integra- tion, but this would be very inefficient. In order to avoid these difficulties in the nume- rical computations we have developed a method which is still based on Cartesian coordinates, but leads now to a trajectory in which all constraints are fulfilled exactly at each step of the integration. This is obtained without any loss of precision. Consider a system of N interacting points subject to I holonomic constraints %(W) = 0 (k = I,..., I). (1) The force acting on each point can be divided into two contributions: the force Fi due to the potential energy and the force of constraint Gt due to all constraints ok involving the ith particle. Gi can be written [5] Gi = - ik X,(r) Via, , 1 where the {h,(r)} are a set of I Lagrangian multipliers depending only on time. Let us suppose that the integration algorithm used takes into account the time derivatives of forces up to order n - 2 (n is often equal to 2), and that it is equivalent to a Taylor expansion of the coordinates r(t + At) up to order (At)“, where At is the time step. In order to be consistent with the integration algorithm, the derivatives of G, and hence of (A} also have to be known up to order n - 2. To avoid the errors implied in the conservation of the constraints, which are accumulative, we employ a method by which the 1 relations of constraint are exactly fulfilled at time t + At. Thus instead of solving for {hrP2’ (t)}, a set of parameters (n} are obtained. In Section 2 we show that the substitution of (yk} for {h~-2) (t)} means that the calculation of coordinates at successive steps is stiI1 exact up to (dt)m, and thus can be applied in any molecular dynamics algorithm. In addition, the constraints are now automatically fulfilled.INTEGRATION OF A SYSTEM WITH CONSTRAINTS 329 In Section 3, we apply the procedure to the well-known algorithm used by Verlet [6]. This algorithm is particularly suitable for use with our method because the forces of constraint have not to be evaluated explicitly at time t. Their effect is obtained instead from the calculation of the set of 1 parameters {n}. The 1 equations giving them are quadratic in {y); they can be solved by an iterative