SIAM J. SCI. COMPUT. c  2006 Society for Industrial and Applied Mathematics Vol. 27, No. 6, pp. 2121–2139 NUMERICAL STOCHASTIC INTEGRATION FOR QUASI-SYMPLECTIC FLOWS ∗ R. MANNELLA † Abstract. Algorithms for the numerical integration of Langevin equations obeying detailed bal- ance are introduced. The algorithms are derived fulfilling the requirements that they should become symplectic in the deterministic frictionless limit, and should reproduce the equilibrium distributions to some higher order in the integration time step when the equilibrium distribution exists. Exten- sions to the case when the system volume or pressure are kept constant are discussed. Comparisons with other integration schemes are carried out both for static and dynamical quantities. Key words. numerical integration, stochastic differential equations, symplectic integration AMS subject classifications. 65C30, 65Z05 DOI. 10.1137/040620965 1. Introduction. Stochastic processes are well known to be at the heart of many physical systems. 1 Several approaches have therefore been developed to understand the dynamics realized by specific models: among others, one of techniques which is widely used is Monte Carlo integration of the equations of motion. The literature on the numerical integration of stochastic differential equations (SDEs) is huge; we will limit ourselves to mentioning a couple of classical citations widely used in the physics community [2, 3, 4]. Additional comments and references can be found in [5]. The numerical algorithms presented in these works are very general ones and can be applied to basically any stochastic flow; however, they might not be the optimal ones for cases when additional information about the system under study is available. An important class for which specialized algorithms can be derived is given by the equations of motion ˙ x = v, ˙ v = −γv + F (x)+ ξ (t), (1.1) where ξ (t) is a random Gaussian noise, with zero average and standard deviation 〈ξ (t)ξ (s)〉 =2γTδ(t − s). (1.2) In the following, we will also use V (x), defined as F (x) ≡−V ′ (x). Note that although for simplicity we will first deal with only one Brownian walker, we will then extend the algorithms to the case of many walkers, where x and v are vectors, and F i (x), the force acting on the ith walker, is a function of all other walkers and/or an external field. The above equation is commonly found in the liquid state literature (for numerical schemes appropriate in the integration of the Brownian dynamics of a liquid, see among others [6, 7, 8, 9, 10, 11, 12, 13]), and some algorithms have been proposed, ∗ Received by the editors December 16, 2004; accepted for publication (in revised form) July 7, 2005; published electronically February 21, 2006. http://www.siam.org/journals/sisc/27-6/62096.html † Dipartimento di Fisica and INFM, UdR Pisa, Universit`a di Pisa, Largo Bruno Pontecorvo 3, I-56127 Pisa, Italy (mannella@df.unipi.it). 1 For a classical introduction, see [1]. 2121