An attempt toward the generalized Langevin dynamics simulation

An attempt to generalize the Langevin dynamics simulation method is presented based on the generalized Langevin theory of liquids, in which the dynamics of both solute and solvent is treated by the generalized Langevin equations, but the integration of the equation of motion of solute is made in the manner similar to the ordinary molecular dynamics simulation with discretized time steps along a trajectory. A preliminary result is derived based on an assumption of the uniform solvent density. The result is regarded to be a microscopic generalization of the phenomenological Langevin theory for the harmonic oscillator immersed in a continuum solvent developed by Wang and Uhlenbeck.


Introduction
It has been five decades since the molecular dynamics simulation scored its first step in the study of liquids and solutions [1].Accelerated by the increasing power of computer, the method is now enjoying the status of a standard tool to explore the molecular aspects of physical, chemical, and biological processes in liquids and solutions.However, the method is facing a high barrier which may not be overcome by the improvement of computing power alone.One of these problems is the phenomena related to the thermodynamic limit essentially such as the phase transition and separation.These phenomena are induced by a subtle balance of the two quantities of thermodynamics in the phases, entropy and enthalpy, both of which do make sense in the thermodynamic limit, namely in the infinite limit of V and N where V is the volume of container and N is the number of molecules in the container.The other problem is a large-scale fluctuation and dynamics of a solute in solution.A typical example of such phenomena is the protein folding in solutions.The phenomenon is not only governed by both the local force field inside protein and the solutesolvent interactions, but also is connected with the large-scale density fluctuations of solvent, or the collective density excitations.The problem is concerned with the long time limit, or the low frequency limit of the dynamic behavior as well as with the thermodynamic limit.So, it will be extremely difficult to carry out such a simulation using the direct molecular dynamics method.There is another problem of biological interest in which the molecular simulation will be hindered in the sense of "molecular recognition" or ligand binding process.The molecular simulation can be applied to the process in some limiting cases, but it will never be predictive.Why is it so?It is because the molecular simulation is designed essentially for sampling limited area in the phase space.The molecular recognition, which is a rare event, requires the sampling of the area in the phase space which the molecular simulation cannot access.At best, the method is capable of calculating the free energy change from bulk to the host cavity using the thermodynamic perturbation technique.However, the method requires in advance the experimental information regarding the position and the number of ligands inside the host cavity as well as the number of water molecules to be replaced.Unfortunately, the information is the one which should be predicted by the theory.

B.Kim et al.
Considerable efforts have been exerted to overcome such a shortcoming inherent to the molecular simulation method.A typical one is the Langevin dynamics simulation, in which protein dynamics including the conformational fluctuation is described by the Langevin equation, or by stochastic differential equation.The frictional force on the atomic motion of protein is represented by the phenomenological ones such as the Oseen tensor [2,3].
The method is implemented in the calculation of dynamics of solution in a variety of phenomena including ion channels.The applicability of the method is severely limited, however, due to the oversimplified model of liquid and solutions.The over-simplification is manifested in two aspects, i.e., structure and dynamics.The solution is represented by the structureless continuum.Therefore, the method cannot describe the interaction induced by solvent, such as hydrogen-bond.The friction is just local in time, which means that no memories along the trajectory are included.It makes the short-time dynamics of solute unrealistic.
A legitimate way of improving the method is to generalize the Langevin equation with respect to the two points stated above: liquid structure and dynamics.For the past decade, we have proposed theories to describe the dynamics of molecular liquids, which combines the reference interactionsite method (RISM) with the generalized Langevin equation (GLE): the liquid structure by RISM and the dynamics by GLE [4].The theory has described the collective density excitation in water, and its response to the ion dynamics.However, the theory is not good enough to describe the dynamics of a complex system like a protein in aqueous solutions.The reason is obvious.It is because the treatment of such a system requires the structure and dynamics of nonuniform or inhomogeneous liquids.Just imagine a liquid in contact with a protein and that in bulk solution.The conventional treatment of GLE already involves the higher order correlations in its formalism, which are conventionally treated by employing the superposition approximation, because it is concerned with the time evolution of the higher order moment of the density fluctuation.On the other hand, in case of the nonuniform liquid, it concerns the lowest order moment of the density fluctuation, or the average density being not constant but dependent on position.For example, the average density of liquid in the cavity of protein is never the same as that in the bulk solution.So, it is an essential requirement to describe the liquid dynamics subject to the field from protein in order to be able to handle the higher order moments of the density fluctuations.
Concerning the problems stated above, we have recently made a large step toward the final solution.That is the generalization of the RISM theory to the inhomogeneous density regime, which we refer to as "Inhomogeneous RISM (IRISM)" [5].The IRISM has described the site-site correlation functions of a pair of water molecules subject to the field of an ion, which is nothing but the three-body correlation functions, "ion-water-water".The new theory has a potential capability of describing the higher order moment of the density fluctuations if it is incorporated into the generalized Langevin treatment of liquids.That is our motivation to formulate the "generalized Langevin dynamics simulation" in which the dynamics of both solute and solvent is treated on equal footing by the generalized Langevin equations, but the integration of the equation of motion of solute is made in the manner similar to the ordinary molecular dynamics simulation with discretized time steps along a trajectory.
In the present article, we present the idea of the generalized Langevin dynamics simulation by taking the simplest model in which a finite number of solute molecules is immersed in the infinite number of solvent molecules.This paper should be regarded as a progress report, because we have not yet reached the final goal.We first develop a general formalism to treat a solute-solvent system entirely based on the projection operator method.We then consider the limit of uniform solvent density in order to make contact with the classical Langevin treatment.

Projection operator method: summary
Since the projection operator method is well-known [6][7][8], we here only summarize the general results of the method.It gives the time evolution equation of a dynamic variable A(t) which is a function of microscopic variables.Its microscopic time evolution is governed by the Liouville operator iL whose expression will be given hereinafter: where {a, b} PB is the Poisson bracket, and H is the Hamiltonian of the system.The formal solution of (2.1) is given by Now the projection operator P is defined as 3) The inner product (a, b) denotes an average of the canonical distribution exp − H/k B T : where Γ denotes all microscopic degrees of freedom in the system.The operator projects out only the "component" of A from the object (• • •).Then obviously PA = A holds.It also has the idempotent property P 2 = P.
After projecting A-component out of the microscopic degrees of freedom, the exact time evolution equation for A(t) is given by Here the frequency matrix iΩ, the memory matrix K(t), and the fluctuating force vector f (t) are given by where Q ≡ 1 − P. One can easily show that the fluctuating force f (t) does not have A-component, i.e., (A, f (t)) = 0. Using this feature and the linearity of the equation, we immediately obtain the following dynamic equation for the auto-correlation function of A(t), C(t) (2.7)

Generalized Langevin equations
We consider a system of finite number of spherical solute particles immersed in a simple liquid solvent.The Hamiltonian of our system is given by where m and M respectively denote the mass of solute particle and solvent atom.The Hamiltonian of the solvent is denoted by H 0 where r i and p i are respectively the position and momentum of the ith atom of the liquid, and U 0 (r ij ) (r ij ≡ |r i − r j |) is the pair potential energy of the solvent atoms.H 1 is the Hamiltonian of the N u solute particles, and R i and P i are the position and momentum of the ith solute particle, and U int (|r i − R j |) is the interaction potential energy between the solute and the solvent.The associated Liouville operator iL is given by where are respectively the force exerted on the ith solute particle by the other solute particles and by the solvent liquid, and are given by Our dynamic variable A(t) is chosen to be For convenience, we regard R(t) and P(t) 3N u component-column vectors.Here δρ k (t) is the Fourier component of the density fluctuation δρ(r, t) ≡ ρ(r, t) − ρ 0 (ρ 0 is the average number density) of the solvent liquid: Likewise, J L k (t) is the Fourier component of the longitudinal part of the current of the solvent liquid: where k = |k| and k = k/k.Balucani and Zoppi [6] have worked out the special case where only the momentum of a solute particle was chosen as a dynamic variable.We now proceed to obtain the specific expressions of the equations (2.5) or (2.6).The first one has to compute the correlation matrix A, A and its inverse A, A −1 .The inner product denotes the average over the canonical distribution exp − βH(Γ) with β ≡ 1/(k B T ): where Z is the partition function Z ≡ dΓ exp − βH(Γ) .

The correlation matrix A, A
The correlation matrix C = A, A is given by We first identify the vanishing elements.The following elements vanish: (R, P) = 0, (R, J L k ) = 0, (P, R) = 0, (P, δρ k ) = 0, (P, They vanish since the momentum integrations dp i p i exp − β i p 2 i /2m = 0 or dPP exp − βP 2 /2M = 0. We look at the nonvanishing elements.The momentum and current correlation are easy to compute.The momentum correlation is given by where Z P ≡ dPP exp − βP 2 /2M and 1 is the (3N u × 3N u ) unit matrix.The equation (3.10) is nothing but the equipartition theorem.Likewise the correlation of the longitudinal current involves the momentum integration where dp N ≡ dp 1 dp 2 • • • dp N .Using this result, one can easily obtain the longitudinal current correlation as (3.12) The remaining elements (R, R), (R, δρ k ), (δρ k , R), and (δρ k , δρ k ) involve the spatial coordinates only.We consider them in order.In the present work we will not specify particular form for correlation of initial position of solute particles since here we are interested in laying out a general structure of the dynamics: R, R ≡ L, where L is a (3N u × 3N u ) matrix.
We then consider (R, δρ k ).First note that (R, where Z c ≡ dR Nu dr N e −βH .It is shown in the Appendix that the correlations B k and B T k vanish when the system is homogeneous: Finally, we have the static structure factor of solvent: Summing up the above results, we have the following block-diagonal matrix for (A, A).
where O denotes the (3N u × 3N u ) zero matrix, 0 is the (3N u × 1) matrix, and the superscript T is the transpose matrix.

Inverse of (A, A)
The inverse of the correlation matrix (A, A) is then given by (3.18)

The frequency matrix iΩ
Here we compute the frequency matrix iΩ which is defined as We first look at the elements of the matrix (A γ , Ȧα ): (3.20) First we obtain some elements of Ȧ using the Liouville operator (3.2).

Ṙi = iLR
where F i is the total force exerting on the ith solute particle by the solvent as well as by other solute particles.Actually when we compute the elements involving Ṗ or JL k , it is more convenient to use the integration by parts.It is useful to remember that whereas Ṙ and ρk involve single momentum (P or p i ), Ṗ and JL k involve zero or two momentums (two p i ).Using this fact, we can easily identify the vanishing elements: (R, Ṙ) = 0, (R, ρk ) = 0, (R, JL k ) = 0, (P, Ṗ) = 0, (P, ρk ) = 0, (P, JL k ) = 0, (δρ k , Ṙ) = 0, (δρ k , Ṗ) = 0, (δρ k , ρk ) = 0, The nonvanishing elements are Taking all these into account, we obtain (3.24) Using the inverse correlation matrix (3.18), we compute iΩ as (3.25)

The reversible part
From (2.5), the reversible part of the Langevin equation is given by iΩ • A(t).Using (3.25), we obtain (3.26)

The fluctuating force
The fluctuating force at t = 0 from (2.6) is given by where we used P Ȧ = (A, Ȧ) which is obtained by setting t = 0 in (3.26).Using the above two results, we obtain for the fluctuating force: (3.31)

The memory matrix
The memory function matrix K(t) is calculated as In (3.32), the two terms exhibit an explicit N -dependence.In the thermodynamic limit in which N is taken to be infinite while the number of solute particles N u remains finite, the term (m/N k B T ) R k , W(t) will vanish since the ensemble average R k , W(t) will remain finite.The other term m N kBT R k , R k (t) will not vanish since the ensemble average R k , R k (t) is proportional to N .Therefore, only the latter term survives in the thermodynamic limit.
In Appendix B, we show that Therefore, the final expression of the memory matrix is given by (3.34)

The explicit form of the exact dynamic equations
With the explicit results of the previous sections, we here write down an explicit form for the time evolution equation (2.5)

Discussions
Our main results in the present work are summarized in (3.35).The first two equations are concerned with the dynamics of solute particles, whereas the last two are describing the solvent dynamics.The equations were derived under the assumption that the average liquid (solvent) density is uniform everywhere.Due to this assumption, the static effect of solute field on the solvent density is disregarded entirely.Nevertheless, some static as well as dynamic coupling between solute and solvent is retained.
The equations for solute have a typical expression of the original Langevin equation, but each term in the right hand side has microscopic description in contrast to the original one [9].The first term is related to the variance-covariance matrix of mean square displacement of solute, which signifies the conformational fluctuation of the molecule.The factor k B T L −1 is related to a frequency matrix of the fluctuation, diagonalization of which gives rise to an "effective normal mode" of the fluctuation.Here, the word "effective" means that the conformation and the normal mode is influenced by solvent.The second term is the damping or drag term.In the original Langevin equation, this term is local in time and is described by the phenomenological expressions such as Oseen tensor or Rotne-Prager tensor, which models the so-called "hydrodynamic interactions" [3].In contrast, our drag term has a memory, the kernel of which is expressed by a time correlation function of the fluctuating force.The fluctuating force consists of the two terms; one is an explicit force acting on a solute particle from solvent as well as from other particles in solute, and the other one is due to the conformational fluctuation of solute.The third term stands for the fluctuation force acting on the solute, which is orthogonal to the dynamic variables at time zero.Since the equations are meant to be integrated directly in the real space, it will become a challenge how to integrate the last two terms [10].
The last two equations in (3.35) describe the time evolution of the solvent density and current fluctuations in the form of the generalized Langevin equation, which is the origin of the drag force acting on the solute particle.Due to our assumption of the uniform solvent density, the equations are essentially independent of the solute influence.The equations can be solved routinely in terms of the memory equation for the dynamic structure factor, or the van Hove space-time correlation function, from which one should be able to evaluate the drag force acting on the solute particles, or the "hydrodynamic force".
The equation derived here is a generalization of the phenomenological Langevin equation for a system of harmonic oscillator derived by Wang and Uhlenbeck [9] and employed by Lamm and Szabo [3] for their "Langevin mode analysis" of macromolecules.The generalization is made in two ways.Firstly, not only the intramolecular interaction of solute but also solvent-solvent and solutesolvent interactions are explicitly taken into account.Secondly, the solvent dynamics is explicitly incorporated in the solute dynamics through the friction kernel which is non-local in time.The equations also provide a molecular basis for the phenomenological Rouse-Zimm model [11], because the equation (3.35) can describe the hydrodynamic interaction between solute particles mediated by the solvent since the memory matrix W i W j (t) has off-diagonal components.

Concluding remarks
In the present study, we have presented a result concerning the way of building a new theory for the generalized Langevin dynamics simulation.In this preliminary study, we neglected the inhomogeneity of the solvent density induced by solute.Based on this assumption, we have derived the equations which can be a molecular basis for the conventional treatments using the phenomenological Langevin equations.The results can be readily extended to a system in which the solvent consists of polyatomic molecules.The recipes for such an extension have been already provided in our earlier publications [12].
However, the theory presented here will not be sufficient in describing the realistic dynamics of molecules with biological interests, say, ions and water dynamics in an ion channel.It is because such molecules in confined space are under strong effect of the solute field, and the densities of the molecules are entirely different from those in bulk.It is obvious that our assumption of the uniform solvent density breaks down in that situation.The extension of the present theory to such cases requires an account for the density profile and the density pair correlation functions of solvent around the solute.An extension of the RISM theory to such inhomogeneous liquids has been already presented in our recent work [5].A study for the purpose of incorporating the inhomogeneity of liquid density into the generalized Langevin dynamics simulation is in progress in our group.
We first consider the first term of (B.1) 3) We have already shown that the last two terms vanish (Appendix A and (3.22)).We now show that the first two terms vanish as well.
where the last equality results upon doing the R i -integration by parts.Obviously, we will have the same result for the first term, F i JL k = 0, since JL k does not contain R i .Therefore, we have shown that Next we consider W i QiLR k = W i Q Ṙk .We need to compute the term P Ṙk which is given by the from (2.3) where the last equality holds since R l Ṙk = − P l R k = 0 and P l Ṙk = 0. Using (B.6), we have It is important to note that Q Ṙk only involves the solvent coordinates and momenta.Then, since W i δρ k = 0 and W i J L k = 0, we obtain ) where in the last line the first term vanishes from the argument shown in (B.4), and we already showed that the second term vanishes.
It is now clear that the repeated applications of QiL on R k will never generate the solutevariable components, and hence W i (QiL) n R k = 0. Therefore, we obtain the final result (B.9) 9. Wang M.C., Uhlenbeck G.E., Rev. Mod.Phys., 1945, 17, 323.10.At this point it seems appropriate to mention that there has been a call from a prominent statistical physicist some years ago to look into the projected dynamics from a computational point of view in his lecture in S. T. Cho Memorial Workshop held in Seoul, Korea in 1997.We here quote below the relevant paragraph from K. Kawasaki's lecture appeared in Progress in Statistical Physics, edited by Sung W. et al, World Scientific, 1998." • • • I expect that in future a method will be found to deal with this kind of projected dynamicshere projected dynamics more precisely means the one governed by a Liouville operator like (1 − P)iL where the relevant part is projected out-directly by computer, and in fact I would like to encourage those familiar with both computer and formal manipulation to look into this problem.The reason is that, generally speaking, projected dynamics has a simpler character despite its appearance than the original unprojected dynamics because complex dynamical behavior that often shows up at long times are projected out.I would like to take this opportunity to mention that there will be a chance to open up a new approach to extend the use of computer enormously if one can find a way to handle projector by computer.Computer can be used to obtain a quantity involving projected dynamics such as "memory function" in some general sense which usually has a short decay time and needs less computer time.This result can be fed into an identity which relates a quantity we need which may exhibit complicated long time behavior with memory function which can be studied efficiently by computer."11.Doi M., Edwards S.F.The Theory of Polymer Dynamics.Oxford University Press, 1986.12. Chong S.-H., Hirata F., Phys.Rev. E, 1998, 57, 1691; 1998, 58, 6188; 1998, 58, 7296.Спроба в напрямку комп'ютерного моделювання узагальненої ланжевенiвської динамiки Представлено спробу узагальнити метод моделювання ланжевенiвської динамiки на основi узагальненої ланжевенiвської теорiї рiдин, в яких динамiка i домiшки i розчинника розглядається узагальненими рiвняннями Ланжевена, але iнтегрування рiвнянь руху домiшки робиться подiбним чином до звичайного методу молекулярної динамiки з дискретизованим часовим кроком вздовж траєкторiї.Отримано попереднiй результат, що грунтується на припущеннi однорiдної густини розчинника.Результат трактується як мiкроскопiчне узагальнення феноменологiчної ланжевенiвської теорiї для гармонiчного осцилятора, помiщеного у суцiльне середовище розчинника, що розроблялась Ванґом та Уленбеком.