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 t...

Повний опис

Збережено в:

Бібліографічні деталі
Дата:	2008
Автори:	Kim, B., Chong, S.-H., Ishizuka, R., Hirata, F.
Формат:	Стаття
Мова:	English
Опубліковано:	Інститут фізики конденсованих систем НАН України 2008
Назва видання:	Condensed Matter Physics
Онлайн доступ:	http://dspace.nbuv.gov.ua/handle/123456789/119006
Теги:	Додати тег Немає тегів, Будьте першим, хто поставить тег для цього запису!
Назва журналу:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Цитувати:	An attempt toward the generalized Langevin dynamics simulation / B. Kim, S.-H. Chong, R. Ishizuka, F. Hirata // Condensed Matter Physics. — 2008. — Т. 11, № 1(53). — С. 179-190. — Бібліогр.: 12 назв. — англ.

Репозитарії

Digital Library of Periodicals of National Academy of Sciences of Ukraine

id	irk-123456789-119006
record_format	dspace
spelling	irk-123456789-1190062017-06-04T03:03:09Z An attempt toward the generalized Langevin dynamics simulation Kim, B. Chong, S.-H. Ishizuka, R. Hirata, F. 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. Представлено спробу узагальнити метод моделювання ланжевен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льне середовище розчинника, що розроблялась Ванґом та Уленбеком. 2008 Article An attempt toward the generalized Langevin dynamics simulation / B. Kim, S.-H. Chong, R. Ishizuka, F. Hirata // Condensed Matter Physics. — 2008. — Т. 11, № 1(53). — С. 179-190. — Бібліогр.: 12 назв. — англ. 1607-324X PACS: 60 DOI:10.5488/CMP.11.1.179 http://dspace.nbuv.gov.ua/handle/123456789/119006 en Condensed Matter Physics Інститут фізики конденсованих систем НАН України
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection	DSpace DC
language	English
description	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.
format	Article
author	Kim, B. Chong, S.-H. Ishizuka, R. Hirata, F.
spellingShingle	Kim, B. Chong, S.-H. Ishizuka, R. Hirata, F. An attempt toward the generalized Langevin dynamics simulation Condensed Matter Physics
author_facet	Kim, B. Chong, S.-H. Ishizuka, R. Hirata, F.
author_sort	Kim, B.
title	An attempt toward the generalized Langevin dynamics simulation
title_short	An attempt toward the generalized Langevin dynamics simulation
title_full	An attempt toward the generalized Langevin dynamics simulation
title_fullStr	An attempt toward the generalized Langevin dynamics simulation
title_full_unstemmed	An attempt toward the generalized Langevin dynamics simulation
title_sort	attempt toward the generalized langevin dynamics simulation
publisher	Інститут фізики конденсованих систем НАН України
publishDate	2008
url	http://dspace.nbuv.gov.ua/handle/123456789/119006
citation_txt	An attempt toward the generalized Langevin dynamics simulation / B. Kim, S.-H. Chong, R. Ishizuka, F. Hirata // Condensed Matter Physics. — 2008. — Т. 11, № 1(53). — С. 179-190. — Бібліогр.: 12 назв. — англ.
series	Condensed Matter Physics
work_keys_str_mv	AT kimb anattempttowardthegeneralizedlangevindynamicssimulation AT chongsh anattempttowardthegeneralizedlangevindynamicssimulation AT ishizukar anattempttowardthegeneralizedlangevindynamicssimulation AT hirataf anattempttowardthegeneralizedlangevindynamicssimulation AT kimb attempttowardthegeneralizedlangevindynamicssimulation AT chongsh attempttowardthegeneralizedlangevindynamicssimulation AT ishizukar attempttowardthegeneralizedlangevindynamicssimulation AT hirataf attempttowardthegeneralizedlangevindynamicssimulation
first_indexed	2025-07-08T15:04:41Z
last_indexed	2025-07-08T15:04:41Z
_version_	1837091606185902080
fulltext	Condensed Matter Physics 2008, Vol. 11, No 1(53), pp. 179–190 An attempt toward the generalized Langevin dynamics simulation B.Kim1,2, S.-H.Chong2, R.Ishizuka2, F.Hirata2 1 Department of Physics, Changwon National University, Changwon 641–773, Korea 2 Department of Theoretical Studies, Institute for Molecular Science, Okazaki 444–8585, Japan Received December 10, 2007 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. Key words: generalized Langevin theory, MD simulation, RISM, 3D-RISM PACS: 60 1. 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 thermody- namics 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 solute- solvent 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 calcu- lating 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. c© B.Kim, S.-H.Chong, R.Ishizuka, F.Hirata 179 B.Kim et al. Considerable efforts have been exerted to overcome such a shortcoming inherent to the molec- ular 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 interaction- site 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. 2. Projection operator method: summary Since the projection operator method is well-known [6–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 180 An attempt toward the generalized Langevin dynamics simulation operator iL whose expression will be given hereinafter: dA(t) dt ≡ iLA(t) ≡ {A,H}PB(t), (2.1) 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 A(t) = eiLtA(0) ≡ eiLtA. (2.2) Now the projection operator P is defined as P ( · · · ) ≡ ( A, · · · )( A,A )−1 A. (2.3) The inner product (a,b) denotes an average of the canonical distribution exp ( −H/kBT ) : (a,b) ≡ 〈 a∗b 〉 = Z−1 ∫ dΓa∗b exp ( −H(Γ)/kBT ) , (2.4) 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 P2 = P. After projecting A-component out of the microscopic degrees of freedom, the exact time evo- lution equation for A(t) is given by dA(t) dt = iΩ · A(t) − ∫ t 0 dsK(t − s) · A(s) + f(t). (2.5) Here the frequency matrix iΩ, the memory matrix K(t), and the fluctuating force vector f(t) are given by iΩ = ( A, Ȧ ) · ( A,A )−1 , K(t) = ( f , f(t) ) · ( A,A )−1 , f(t) = exp ( tQiL ) QȦ, (2.6) 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) dC(t) dt = iΩ · C(t) − ∫ t 0 dsK(t − s) · C(s). (2.7) 3. 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 H ≡ H0 + H1 + H2, H0 = N ∑ i=1 [ p2 i 2m + ∑ j>i U0(\|ri − rj \|) ] , H1 = Nu ∑ i=1 [ P2 i 2M + ∑ j>i U1(\|Ri − Rj \|) ] , H2 = N ∑ i=1 Nu ∑ j=1 Uint(\|ri − Rj \|), (3.1) 181 B.Kim et al. where m and M respectively denote the mass of solute particle and solvent atom. The Hamiltonian of the solvent is denoted by H0 where ri and pi are respectively the position and momentum of the ith atom of the liquid, and U0(rij) (rij ≡ \|ri−rj \|) is the pair potential energy of the solvent atoms. H1 is the Hamiltonian of the Nu solute particles, and Ri and Pi are the position and momentum of the ith solute particle, and Uint(\|ri −Rj \|) is the interaction potential energy between the solute and the solvent. The associated Liouville operator iL is given by iL ≡ iL0 + iL1 + iL2, iL0 ≡ N ∑ i=1 [ 1 m pi · ∂ ∂ri − ∑ j>i ∂U0(rij) ∂ri · ∂ ∂pi ] , iL1 ≡ Nu ∑ i=1 [Pi M · ∂ ∂Ri + F (u) i · ∂ ∂Pi ] , iL2 ≡ − N ∑ i=1 Nu ∑ j=1 ∂Uint(\|ri − Rj \|) ∂ri · ∂ ∂pi + Nu ∑ i=1 F (v) i · ∂ ∂Pi , (3.2) where rij ≡ \|ri − rj \| and Rij ≡ \|Ri −Rj \|. In (3.2), F (u) i and F (v) i are respectively the force exerted on the ith solute particle by the other solute particles and by the solvent liquid, and are given by F (u) i = − ∑ j 6=i ∂U1(Rij) ∂Ri , F (v) i = − N ∑ j=1 ∂Uint(\|Ri − rj \|) ∂Ri . (3.3) Our dynamic variable A(t) is chosen to be A(t) =     R(t) P(t) δρk(t) JL k (t)     . (3.4) For convenience, we regard R(t) and P(t) 3Nu 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: δρ(r, t) ≡ ∑ i δ ( r − ri(t) ) − ρ0, δρk(t) = ∫ dreik·rδρ(r, t) = ∑ i eik·ri(t) − (2π)3ρ0δ(k). (3.5) Likewise, JL k (t) is the Fourier component of the longitudinal part of the current of the solvent liquid: ρ̇k(t) = 1 m ∑ i ik · pi(t)e ik·ri(t) ≡ ikJL k (t), JL k (t) = 1 m ∑ i k̂ · pi(t)e ik·ri(t) , (3.6) 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 182 An attempt toward the generalized Langevin dynamics simulation the average over the canonical distribution exp ( − βH(Γ) ) with β ≡ 1/(kBT ): ( a,b ) ≡ 〈 a∗b 〉 = 1 Z ∫ dΓa∗(Γ)b(Γ) e−βH(Γ) , (3.7) where Z is the partition function Z ≡ ∫ dΓ exp ( − βH(Γ) ) . 3.1. The correlation matrix ( A, A ) The correlation matrix C = ( A,A ) is given by ( A,A ) =     (R,R) (P,R) (δρk,R) (JL k ,R) (R,P) (P,P) (δρk,P) (JL k ,P) (R, δρk) (P, δρk) (δρk, δρk) (JL k , δρk) (R, JL k ) (P, JL k ) (δρk, JL k , ) (JL k , JL k )     . (3.8) We first identify the vanishing elements. The following elements vanish: (R,P) = 0, (R, JL k ) = 0, (P,R) = 0, (P, δρk) = 0, (P, JL k ) = 0, (δρk,P) = 0, (δρk, JL k ) = 0, (JL k ,R) = 0, (JL k ,P) = 0, (JL k , ρk) = 0. (3.9) They vanish since the momentum integrations ∫ dpi pi exp ( −β ∑ i p 2 i /2m ) = 0 or ∫ dPP exp ( − βP2/2M ) = 0. We look at the nonvanishing elements. The momentum and current correlation are easy to compute. The momentum correlation is given by (P,P) = 1 ZP ∫ dPPPe−βP 2/2M = 1 ZP ∫ dPP ( − MkBT ) ∂ ∂P e−βP 2/2M = MkBT1,(3.10) where ZP ≡ ∫ dPP exp ( −βP2/2M ) and 1 is the (3Nu×3Nu) unit matrix. The equation (3.10) is nothing but the equipartition theorem. Likewise the correlation of the longitudinal current involves the momentum integration 1 Zp ∑ α,γ k̂αk̂γ ∫ dpN pα i pγ j exp ( − β ∑ l p2 l /2m ) = mkBTδij ∑ α,γ k̂αk̂γδαγ = mkBTδij , (3.11) where ∫ dpN ≡ ∫ dp1dp2 · · · dpN . Using this result, one can easily obtain the longitudinal current correlation as (JL k , JL k ) = 〈 JL −k JL k 〉 = NkBT m . (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, (3.13) where L is a (3Nu × 3Nu) matrix. We then consider (R, δρk). First note that (R, δρk) = 〈 Rρk 〉 − (2π)3ρ0δ(k) 〈 R 〉 = 〈 Rρk 〉 since 〈 R 〉 = 0. Therefore Bk ≡ (R, δρk) = 〈 Rρk 〉 = 1 Z c ∫ dRNudrN R ( ∑ i eik·ri ) e−βH, (3.14) 183 B.Kim et al. where Zc ≡ ∫ dRNudrNe−βH. It is shown in the Appendix that the correlations Bk and BT k vanish when the system is homogeneous: Bk = BT k = 0. (3.15) Finally, we have the static structure factor of solvent: Sk ≡ 1 N (δρk, δρk) = 1 N 〈 δρ−kδρk 〉 . (3.16) Summing up the above results, we have the following block-diagonal matrix for (A,A). ( A,A ) =     L O 0 0 O MkBT1 0 0 0T 0T NSk 0 0T 0T 0 NkBT m     , (3.17) where O denotes the (3Nu × 3Nu) zero matrix, 0 is the (3Nu × 1) matrix, and the superscript T is the transpose matrix. 3.2. Inverse of (A, A) The inverse of the correlation matrix (A,A) is then given by ( A,A )−1 =     L−1 O 0 0 O 1 MkBT 1 0 0 0T 0T 1 NSk 0 0T 0T 0 m NkBT     . (3.18) 3.3. The frequency matrix iΩ Here we compute the frequency matrix iΩ which is defined as iΩαβ ≡ ∑ γ (Aγ , Ȧα) [ (A,A)−1 ] γβ . (3.19) We first look at the elements of the matrix (Aγ , Ȧα): iΩ =     (R, Ṙ) (P, Ṙ) (δρk, Ṙ) (JL k , Ṙ) (R, Ṗ) (P, Ṗ) (δρk, Ṗ) (JL k , Ṗ) (R, ρ̇k) (P, ρ̇k) (δρk, ρ̇k) (JL k , ρ̇k) (R, J̇L k ) (P, J̇L k ) (δρk, J̇L k ) (JL k , J̇L k )     · (A,A)−1. (3.20) First we obtain some elements of Ȧ using the Liouville operator (3.2). Ṙi = iLRi = Pi M , Ṗi = iLPi = Fi ≡ F (u) i + F (v) i , ρ̇k = iLρk = ikJL k , (3.21) where Fi 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 Ṗ or J̇L k , it is more convenient to use the integration by parts. It is useful to remember that whereas Ṙ and ρ̇k involve single momentum (P or pi), Ṗ and J̇L k involve zero or two momentums (two pi). Using this fact, we can easily identify the vanishing elements: (R, Ṙ) = 0, (R, ρ̇k) = 0, (R, J̇L k ) = 0, (P, Ṗ) = 0, (P, ρ̇k) = 0, (P, J̇L k ) = 0, (δρk, Ṙ) = 0, (δρk, Ṗ) = 0, (δρk, ρ̇k) = 0, (JL k , Ṙ) = 0, (JL k , Ṗ) = 0, (JL k , J̇L k ) = 0. (3.22) 184 An attempt toward the generalized Langevin dynamics simulation The nonvanishing elements are (R, Ṗ) = − 1 M (P,P) = −kBT1, (P, Ṙ) = 1 M (P,P) = kBT1, (δρk, J̇L k ) = ikNkBT m , (JL k , ρ̇k) = ikNkBT m . (3.23) Taking all these into account, we obtain iΩ =     O kBT1 0 0 −kBT1 O 0 0 0T 0T 0 ikNkBT m 0T 0T ikNkBT m 0     · (A,A)−1. (3.24) Using the inverse correlation matrix (3.18), we compute iΩ as iΩ =     O 1 M 1 0 0 −kBTL−1 O 0 0 0T 0T 0 ik 0T 0T ikkBT mSk 0     . (3.25) 3.4. The reversible part From (2.5), the reversible part of the Langevin equation is given by iΩ ·A(t). Using (3.25), we obtain iΩ · A(t) =     P(t)/M −kBTL−1 · R(t) ikJL k (t) ikkBT mSk δρk(t)     . (3.26) 3.5. The fluctuating force The fluctuating force at t = 0 from (2.6) is given by f = QȦ = (1 − P)Ȧ = Ȧ − iΩ · A, (3.27) where we used PȦ = (A, Ȧ) · (A,A)−1A = iΩ · A. We first obtain Ȧ =     Ṙ Ṗ ρ̇k J̇L k     =     P/M F ikJL k J̇L k     (3.28) and iΩ · A =     P/M −kBTL−1 · R ikJL k ikkBT mSk δρk     (3.29) which is obtained by setting t = 0 in (3.26). Using the above two results, we obtain for the fluctuating force: f =     0 W 0 Rk     , f(t) = eitQL     0 W 0 Rk     , (3.30) 185 B.Kim et al. where W ≡ F + kBTL−1 · R, Rk ≡ J̇L k − ikkBT mSk δρk . (3.31) 3.6. The memory matrix The memory function matrix K(t) is calculated as Kαβ(t) ≡ ∑ γ ( fγ , fα(t) )[( A,A )−1] γβ =     O O 0 0 O 1 MkBT (W,W(t)) 0 m NkBT (Rk,W(t)) 0T 0T 0 0 0T 1 MkBT (W,Rk(t)) 0 m NkBT (Rk,Rk(t))     , (3.32) where W(t) ≡ exp ( itQL ) W and Rk(t) ≡ exp ( itQL ) Rk. 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 Nu remains finite, the term (m/NkBT ) ( Rk,W(t) ) will vanish since the ensemble average ( Rk,W(t) ) will remain finite. The other term m NkBT ( Rk,Rk(t) ) will not vanish since the ensemble average ( Rk,Rk(t) ) is proportional to N . Therefore, only the latter term survives in the thermodynamic limit. In Appendix B, we show that (W,Rk(t)) = 0. (3.33) Therefore, the final expression of the memory matrix is given by Kαβ(t) =     O O 0 0 O 1 MkBT (W,W(t)) 0 0 0T 0T 0 0 0T 0T 0 m NkBT (Rk,Rk(t))     . (3.34) 3.7. 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) dR(t) dt = P(t) M , dP(t) dt = −kBTL−1 · R(t) − 1 MkBT ∫ t 0 ds 〈 W ei(t−s)QLW 〉 · P(s) + eitQLW, dρk(t) dt = ikJL k (t), dJL k (t) dt = ikkBT mSk δρk(t) − m NkBT ∫ t 0 ds 〈 R−k ei(t−s)QLRk 〉 JL k (s) + eitQLRk, W ≡ F + kBTL−1 · R, Rk ≡ J̇L k − ikkBT mSk δρk . (3.35) 4. 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 186 An attempt toward the generalized Langevin dynamics simulation 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 kBTL−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 solute- solvent 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 〈 Wi Wj(t) 〉 has off-diagonal components. 5. 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 inho- mogeneity 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 phenomeno- logical 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 187 B.Kim et al. 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. Acknowledgements BK acknowledges financial support by the Institute for Molecular Science. BK also would like to express his gratitude to Prof. Kyozi Kawasaki for his support, guidance, and many inspiring discussions on various themes of nonequilibrium statistical physics. A. Calculation of Bk We show that the correlation Bk ≡ 〈 Rρk 〉 vanishes in the homogeneous system. We consider the ith component of Bk which is explicitly given by 〈 Riρk 〉 = 1 ZC ∫ dRNu ∫ drN Ri ∑ l eik·rle−βH . (A.1) If we shift all the integration variables as rl → rl + Ri and Rj → Rj + Ri except Ri variable (i.e., i 6= j) in the integrand of (A.1), we can separate the Ri-integration as follows: 〈 Riρk 〉 = ∫ dRi Ri eik·Ri ∫ dRi · ∫ dRNu−1 ∫ drN ∑ l e ik·rle−βH′ ∫ dRNu−1 ∫ drNe−βH′ , (A.2) where the Hamiltonian H′ is obtained from H with the shift of the variables, and does not involve Ri. We consider the first integration factor in (A.2). For simplicity we suppress the index i. Its X-component is given by ∫ dRX eik·R ∫ dR = 1 L3 ( ∫ L/2 −L/2 Xeik1XdX )( ∫ L/2 −L/2 eik2Y dY )( ∫ L/2 −L/2 eik3ZdZ ) = 1 L3 · 2 k1 [ 1 k1 sin (k1L 2 ) − L 2 cos (k1L 2 ) ] · 2 k2 sin (k2L 2 ) · 2 k3 sin (k3L 2 ) .(A.3) Note that this integral vanishes in the thermodynamic limit L → ∞ due to the oscillating term eik·R. The second integration factor in (A.2) remains finite. Therefore, we conclude that 〈 Xρk 〉 = 0 in the thermodynamic limit. Since this will hold for the other components, we conclude that Bk = 〈 Rρk 〉 = 0. (A.4) B. Calculation of 〈 We itQL Rk 〉 Here we show that the memory matrix 〈 W eitQLRk 〉 vanishes. We consider its ith component 〈 Wi eitQLRk 〉 : 〈 Wi eitQLRk 〉 = 〈 Wi ( 1 + tQiL + t2 2! QiLQiL + · · · ) Rk 〉 . (B.1) For simplicity of notation, we write Wi and Rk from (3.31) as Wi = Fi + CijRj , Rk = J̇L k − c δρk , (B.2) where the summation is implied for repeated indices, and the matrix C ≡ kBTL−1 and c ≡ ikkBT/mSk. 188 An attempt toward the generalized Langevin dynamics simulation We first consider the first term of (B.1) 〈 WiRk 〉 = 〈 ( Fi + CijRj )( J̇L k − c δρk ) 〉 = 〈 FiJ̇ L k 〉 − c 〈 Fiδρk 〉 + Cij 〈 Rj J̇ L k 〉 − cCij 〈 Rjδρk 〉 . (B.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. 〈 Fiδρk 〉 = 1 Z C ∫ dRNudrN ∑ l eik·rlFie −βH = 1 Z C ∫ dRNudrN ∑ l eik·rl ( − ∂H ∂Ri ) e−βH = 1 Z C ∫ dRNudrN ∑ l eik·rlβ−1 ∂ ∂Ri e−βH = 0, (B.4) where the last equality results upon doing the Ri-integration by parts. Obviously, we will have the same result for the first term, 〈 FiJ̇ L k 〉 = 0, since J̇L k does not contain Ri. Therefore, we have shown that 〈 WiRk 〉 = 0. (B.5) Next we consider 〈 WiQiLRk 〉 = 〈 WiQṘk 〉 . We need to compute the term PṘk which is given by the from (2.3) PṘk = 1 kBT Clm 〈 RlṘk 〉 Rm + 1 MkBT 〈 PlṘk 〉 Pl + 1 NSk 〈 δρ−kṘk 〉 δρk + m NkBT 〈 JL −kṘk 〉 JL k = 1 NSk 〈 δρ−kṘk 〉 δρk + m NkBT 〈 JL −k Ṙk 〉 JL k , (B.6) where the last equality holds since 〈 RlṘk 〉 = − 〈 PlRk 〉 = 0 and 〈 PlṘk 〉 = 0. Using (B.6), we have QṘk = (1 − P)Ṙk = Ṙk − 1 NSk 〈 δρ−kṘk 〉 δρk − m NkBT 〈 JL −kṘk 〉 JL k . (B.7) It is important to note that QṘk only involves the solvent coordinates and momenta. Then, since 〈 Wiδρk 〉 = 0 and 〈 WiJ L k 〉 = 0, we obtain 〈 WiQiLRk 〉 = 〈 WiQṘk 〉 = 〈 WiṘk 〉 = 〈( Fi + CijRj ) Ṙk 〉 = 〈 FiṘk 〉 + Cij 〈 RjṘk 〉 = 0, (B.8) 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 Rk will never generate the solute- variable components, and hence 〈 Wi(QiL)nRk 〉 = 0. Therefore, we obtain the final result 〈 WetQiLRk 〉 = 0. (B.9) References 1. Alder B.J., Wainwright T.E., J. Chem. Phys., 1957, 27, 1208. 2. Ermak D.L., McCammon J.A., J. Chem. Phys., 1978, 69, 1352. 3. Lamm G., Szabo A., J. Chem. Phys., 1986, 85, 7334. 4. Hirata F. (ed.) Molecular Theory of Solvation. Kluwer Academic Publishers, 2003. 5. Ishizuka R., Chong S.-H., Hirata F., J. Chem. Phys. in press. 6. Balucani U., Zoppi M. Dynamics of the Liquid State. Oxford University Press, 1994. 7. Mazenko G.F. Nonequilibrium Statistical Mechanics. Wiley-VCH Verlag, 2006. 8. Reichman D.R., Charbonneau P. Mode Coupling Theory. J. Stat. Mech P05013 2005. 189 B.Kim et al. 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 dynamics – here 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м1,2, С.-Х.Чонґ2, Р.Iшiдзука2, Ф.Хiрата2 1 Фiзичний факультет, Нацiональний унiверситет Чанґвона, Чанґвон, Корея 2 Вiддiл теоретичних дослiджень, Iнститут молекулярної науки, Оказакi, Японiя Отримано 10 грудня 2007 р. Представлено спробу узагальнити метод моделювання ланжевен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я Ланжевена, моделювання методом МД, RISM, 3D-RISM PACS: 60 190

An attempt toward the generalized Langevin dynamics simulation

Репозитарії

Схожі ресурси