Molecular dynamics with Langevin equation using local harmonics and Chandrasekhar's convolution

A numerical method for studying molecular systems subject to a random force field leading to a Gaussian velocity distribution and described by the Langevin equation is presented. Two basic elements constitute the formulation: local harmonic modes and Chandrasekhar's formula for the distribution function for a convolution involving a random function. First, by linearizing the governing Langevin equations locally and employing an orthogonal change of coordinates, an explicit solution for the displacement and velocity is constructed. Second, Chandrasekhar's formula is employed in deriving the probability distribution function of the displacements and the velocities coming from the random forces. The local mode analysis is essential for the use of the Chandrasekhar's formula, since we need the formal solution as a convolution of the random forces and the local Green's function. For an illustration of the method in a significant case representative of real problems, we study a one dimensional idealization of a long chain molecule possessing internal energy barriers and subjected to an applied tension. The results are compared with the predictions of a conventional approximate method where a finite number of random realizations are generated in each time step. This truncation constitutes an approximation to obtain the desired Gaussian probability distribution function for the velocities which is reached in the limit of an infinity of random realizations. The calculations show that the conventional approximations may be acceptable only for short times, small temperatures, and average values over very long times. In particular, these approximations fail to give accurate results for transient phenomena, show slow convergence with the increase in the number of random realizations, and predict large values for the variance even in the steady regime. The new proposed method on the other hand, (i) incorporates the mathematically and conceptually correct limit for the distribution function, (ii) is quite stable with respect to increases in the value of the time increment as well as in terms of fluctuations characterized by the variance, (iii) leads to considerable savings in computer time over the approximate method, and (iv) has the proper description during the transient regime, which is usually the most interesting phase of dynamical processes.