分子动力学论文--算法和参数设置1

分子动力学论文算法和参数设置1

欢迎来到我的博客！坚持比努力重要。

**1.molecular dynamics simulation is a technique for computing the equilibrium and transport properties of a classical many-body system.
分子动力是一个在计算 一个经典多体转换性能的系统 技术
eg:1.分子动力学模拟是一种计算经典多体系统平衡和输运性质的技术。
In this context,the word classical means that the nuclear motion of the constituent particles obeys the laws of classical mechanics
Chapter 4Molecular Dynamics Simulations
Molecular Dynamics simulation is a technique for computing the equilibrium and transport properties of a classical many-body system.
In this context, the word classical means that the nuclear motion of the constituent particles obevs the laws of classical mechanics.
This is an excellent approximation for a wide range of materials.
Only when we consider the translational or rotational motion of light atoms or molecules (He, H2, D,) or vibrational motion with a frequency v such that hv > kpT should we worry about quantum effects.Of course, our discussion of this vast subject is necessarily incomplete.Other aspects of the Molecular Dynamics techniques can be found in (19, 39-41].4.1 Molecular Dynamics: The IdeaMolecular Dynamics simulations are in many respects very similar to real experiments. When we perform a real experiment, we proceed as follows We prepare a sample of the material that we wish to study. We connect this sample to a measuring instrument (e.g., a thermometer, manometer, or viscosimeter), and we measure the property of interest during a certain time nterval. If our measurements are subject to statistical noise (as most measurements are), then the longer we average, the more accurate our measurement becomes. In a Molecular Dynamics simulation, we follow exactly the same approach. First, we prepare a sample: we select a model system consisting of N particles and we solve Newton’s equations of motion for this system until the properties of the system no longer change with time (we
第四章分子动力学模拟
分子动力学模拟是一种计算经典多体系统平衡和输运性质的技术。在这种情况下，“经典”一词意味着构成粒子的核运动符合经典力学定律。这是对各种材料的极好近似。只有当我们考虑光原子或分子(He, H2, D)的平动或旋转运动或频率为v的振动运动时，才会考虑hv >我们应该担心量子效应吗?当然，我们对这个大问题的讨论是不完整的。分子动力学技术的其他方面可以在(19,39 -41).4.1分子动力学:思想分子动力学模拟在许多方面与真实实验非常相似。当我们进行一个真实的实验时，我们如下进行我们准备了一个我们想要研究的材料的样品。我们将样品连接到测量仪器(例如温度计、压力计或粘度计)上，然后在特定时间内测量感兴趣的性质。如果我们的测量受到统计噪声的影响(就像大多数测量一样)，那么我们取平均值的时间越长，我们的测量就越精确。在分子动力学模拟中，我们采用完全相同的方法。首先，我们准备一个样本:我们选择一个由N个粒子组成的模型系统，我们求解这个系统的牛顿运动方程，直到系统的性质不再随时间变化(我们)
equilibrate the system). After equilibrafon, we perform the actual measurement. In fact, some of the most common mistakes that can be made when performing a computer experiment are verv similar to the mistakes that can be made in real experiments (e.g., the sample is not prepared correctlv, the measurement is too short, the system undergoes an irreversible change during the experiment, or we do not measure what we think).To measure an observable quantity in a Molecular Dynamics simulation, we must first of all be able to express this observable as a function of the positions and momenta of the particles in the svstem. For instance, a convenient definition of the temperature in a (classical) many-body system makes use of the equipartition of energv over all degrees of freedom that enter quadratically in the Hamiltonian of the system. In particular for the average kinetic energy per degree of freedom, we have
平衡系统)。平衡后，我们进行实际测量。事实上,一些最常见的错误,可以在执行计算机实验verv类似错误,可以在实际实验中(例如,样品没有准备correctlv,测量太短,在实验系统经历了一个不可逆转的变化,或者我们不衡量我们认为)。为了在分子动力学模拟中测量一个可观测量，我们首先必须能够将这个可观测量表示为svstem中粒子位置和动量的函数。例如，在一个(经典的)多体系统中，温度的一个方便的定义是利用在系统的哈密顿量中以二次方式进入的所有自由度上的能量均分。特别是对于每个自由度的平均动能，我们有
1/2mv(2)=1/2k(角标B)T
In a simulation, we use this equation as an operational definition of the temperature. In practice, we would measure the total kinetic energy of the system and divide this by the number of degrees of freedom Nt (=3N-3 for a system of N particles with fixed total momentum’). As the total kinetic energy of a system fluctuates, so does the instantaneous temperature:
在模拟中，我们使用这个方程作为温度的操作定义。在实践中，我们要测量系统的总动能，然后除以自由度Nt(对于总动量固定的N个粒子的系统=3N-3)。当系统的总动能波动时，瞬时温度也随之波动:
T(t)=(N,i=1求和)m（角标i）v上角标2下角标i/k[下B]N[下f]
The relative fluctuations in the temperature will be of order 1/N. As Nf is typiçally on the order of 102-103, the statistical fluctuations in the temperature are on the order of 5-10%. To get an accurate estimate of the temperature, one should average over many fluctuations.
温度的相对波动将是1/N阶。由于Nf通常在102-103的数量级，因此温度的统计波动在5-10%的数量级。为了准确估计温度，人们应该对许多波动进行平均。
4.2 Molecular Dynamics: A Program
The best introduction to Molecular Dvnamics simulations is to consider a simple program.
The program we consider is kept as simple as possible to illustrate a number of important features of Molecular Dvynamics simulations.
The program is constructed as follows:1. We read in the parameters that specify the conditions of the run (e.g., initial temperature, number of particles, density, time step).
4.2分子动力学:程序分子
对分子动力学模拟最好的介绍是考虑一个简单的程序。我们所考虑的程序尽量保持简单，以说明分子动力学模拟的一些重要特征。
程序构造如下:1. 我们读入指定运行条件的参数(例如。初始温度，粒子数，密度，时间步长)。
Algorithm 3 (A Simple Molecular Dynamics Program)
算法3(一个简单的分子动力学程序)
program md 简单的医学程序
call int 初始化
t=0
do while（t.lt , tmax） MD循环
call force（f , en）
call integrate(f , en) 运动积分方程
t=t+delt
call sample
enddo
stop
end
Comment to this algorithm:

1. Subroutines init, force, integrate, and sample will be described in Algorithms 4, 5, and 6, respectively. Subroutine sample is used to calculate averages like pressure or temperature
2. We initialize the system (i.e., we select initial positions and velocities).
3. We compute the forces on all particles.
4. We integrate Newton’s equations of motion. This step and the previous one make up the core of the simulation. Thev are repeated until we have computed the time evolution of the svstem for the desired lengthof time
5. After completion of the central loop, we compute and print the averages of measured quantities, and stop.
   Algorithm 3 is a short pseudo-algorithrm that carries out a Molecular Dynamics simulation for a simple atomic system. We discuss the different operations in the program in more detail.
   对该算法的评论:
   1。子例程init、force、integration和sample将分别在算法4、5和6中进行描述。子程序样例用于计算压力、温度等参数。
   2。我们初始化系统(例如，我们选择初始位置和速度)。
   3.计算所有粒子上的力。
   4.我们对牛顿的运动方程进行积分。这一步和previpus one构成了仿真的核心。它们不断重复，直到我们计算出所需时间长度下系统的时间演化。
   5.中央回路完成后，计算并打印测得量的平均值，然后停止。
   算法3是一个简短的伪算法，对一个简单的原子系统进行分子动力学模拟。我们将更详细地讨论程序中的不同操作。
   4.2.1 Initialization
   To start the simulation, we should assign initial positions and velocities to all particles in the system. The particle positions should be chosen compatible with the structure that we are aiming to simulate. In anv event, the particles should not be positioned at positions that result in an appreciable overlap of the atomic or molecular cores. Often this is achieved by initially placing
   4.2.1初始化准备
   为了开始模拟，我们应该给系统中的所有粒子分配初始位置和速度。粒子位置的选择应该与我们要模拟的结构相一致。在anv事件中，粒子不应该被放置在导致原子或分子核心明显重叠的位置上。这通常是通过最初的放置来实现的
   Algorithm 4 (Initialization of a Molecular Dynamics Program)
   算法4(分子动力学程序的初始化)
   subroutine init MD程序的初始化
   sumv=0
   sumv2=0
   do i=1，npart
   x(i)=lattice.pos(i) 把粒子放在晶格上
   v(i)=(ranf()-0.5) 给随机速度
   sumv=sumv+v(i) 速度质心
   sumv2=sumv2+v(i)**2 动能
   enddo
   sumv=sumv/npartn 质量的速度中心
   sumv2=sumv2/npart 均方速度
   fs=sqrt(3*temp/sumv2) 速度的比例因子
   do i=1,npart 设定所需动能和设定
   v(i)=(v(i)-sumv)*fs 速度转向零
   xm(i)=x(i)-v(i)*dt 位置前时间步长
   enddo
   return
   end
   Comments to this algorithm:
   对该算法的评论:
6. Function lattice-pos gives the coordinates of lattice position i and ranf () gives a uniformly distributed random number. We do not use a Maxwell-Boltzmann distribution for the velocities; on equilibration it will become a Maxwell-Boltzmann distribution.
7. In computing the number of degrees of freedom, we assume a three-dimensional system (in fact, wve approximate N, by 3N).
8. 函数grid -pos给出了格点位置i的坐标，ranf()给出了均匀分布的随机数。我们不用马克斯韦尔-波尔兹曼分布来表示速度;在平衡状态下，它将变成马克斯韦尔-波尔兹曼分布。
9. 在计算自由度的数目时，我们假设一个三维系统(实际上，我们约取N×3N)。
   the particles on a cubic lattice, as described in section 3.2.2 in the context of Monte Carlo simulations.In the present case (Algorithm 4), we have chosen to start our run from a simple cubic lattice. Assume that the values of the density and initial temperature are chosen such that the simple cubic lattice is mechanically unstable and melts rapidly. First, we put each particle on its lattice site and then we attribute to each velocity component of every particle a value that is drawn from a uniform distribution in the interval -0.5.0.51. This initial velocity distribution is Maxwellian neither in shape nor even in width. Subsequentlv, we shift all velocities, such that the total momentum is zero and we scale the resulting velocities to adjust the mean kinetic energy to the desired value. We know that, in thermal equilibrium, the following relation
   should hold:
   在本例中(算法4)，我们选择从一个简单的立方格开始运行。假设选择了密度和初始温度的值，使简单立方晶格在机械上不稳定并迅速融化。首先，我们把每个粒子放到它的格点上，然后我们给每个粒子的每个速度分量赋予一个值，这个值是从区间内的均匀分布中得到的-0.5.0.51。这种初速度分布在形状上甚至在宽度上都不是麦克斯韦式的。随后，我们改变所有的速度，这样总动量为零，我们调整得到的速度来调整平均动能到期望的值。我们知道，在热平衡中，有以下关系应:
   （v[上2下aerf]=k[下B]T/m
   where va is the a component of the velocity of a given particle. We can use his relation to define an instantaneous temperature at time t T(t):
   va是一个给定粒子的速度的一个分量。我们可以用his关系式定义t时刻的瞬时温度t (t):
   k[下B]T(t)=(N,i>1求和)mv上2下ai/Nf
   Clearly, we can adjust the instantaneous temperature T(t) to match the desired temperature T by scaling all velocities with a factor (T/T(t)1/2. Thisinitial setting of the temperature is not particularly critical, as the temperature will change anyway during equilibration.As will appear later, we do not really use the velocities themselves in our algorithm to solve Newton’s equations of motion. Rather, we use the positions of all particles at the present (x) and previous (xm) time steps, combined with our knowledge of the force (f) acting on the particles, to predict the positions at the next time step. When we start the simulation, we must bootstrap this procedure by generating approximate previous positions. Without much consideration for any law of mechanics but the conservation of linear momentum, we approximate x for a particle in a direction by xm(i) = x(i) -v(i)dt. Of course, we could make a better estimate of the true previous position of each particle. But as we are only bootstrapping the simulation, we do not worry about such subtleties.
   显然，我们可以调整瞬时温度，使其与期望的温度匹配，方法是用一个因子(T/T)1/2来缩放所有速度。这个温度的初始设定不是特别关键，因为在平衡过程中温度无论如何都会改变。正如后面将会出现的，我们在算法中并没有使用速度本身来求解牛顿运动方程。相反，我们使用所有粒子在当前(x)和先前(xm)时间步长的位置，结合我们对作用在粒子上的力(f)的知识，来预测下一个时间步长的位置。当我们开始模拟时，我们必须通过生成近似的先前位置来引导这个过程。除了线性动量守恒，不考虑任何力学定律，我们用xm(i) = x(i) -v(i)dt来近似一个质点的x。当然，我们可以更好地估计每个粒子之前的真实位置。但是，由于我们只是在引导模拟，所以我们并不担心这些细微之处。
   4.2.2受力计算
   What comes next is the most time-consuming part of almost all Molecular Dynamics simulations: the calculation of the force acting on every particle.If we consider a model system with pairwise additive interactions (as we do in the present case), we have to consider the contribution to the force on particle i due to all its neighbors. If we consider only the interaction between a particle and the nearest image of another particle, this implies that, for a system of N particles, we must evaluate N x (N-1)/2 pair distances.This implies that, if we use no tricks, the time needed for the evaluation of the forces scales as N2. There exist efficient techniques to speed up the evaluation of both short-range and long-range forces in such a way that the computing time scales as N, rather than N2. In Appendix F, we describe some of the more common techniques to speed up the simulations. Although the examples in this Appendix apply to Monte Carlo simulations, the same techniques can also be used in a Molecular Dynamics simulation. However, in the present, simple example (see Algorithm 5) we will not attempt to make
   接下来是几乎所有分子动力学模拟中最耗时的部分:计算作用在每个粒子上的力。如果我们考虑一个具有成对加性相互作用的模型系统(就像我们在本例中所做的那样)，我们必须考虑所有相邻粒子i对力的贡献。如果我们只考虑一个粒子和另一个粒子最近的像之间的相互作用，这意味着，对于一个有N个粒子的系统，我们必须计算N x (N-1)/2对距离。这意味着，如果我们不使用技巧，评估力所需要的时间为N2。现有的有效技术可以加速近程和远程力的评估，计算时间尺度为N，而不是N2。在附录F中，我们描述了加速模拟的一些更常见的技术。虽然本附录中的示例适用于蒙特卡洛模拟，但同样的技术也可以用于分子动力学模拟。然而，在目前，简单的例子(见算法5)，我们将不尝试作出
   Algorithm 5 (Calculation of the Forces)
   算法5(力的计算)
   subroutine force(f ,en) 确定的力
   en=0 和能量，
   do i=1,npart
   f(i)=0
   enddo
   do i=1,npart-1
   do j=i+1,npart 循 环遍历所有对
   xr=x(i)-x(j)
   xr=xr-boxnint(xr/box) 周期性边界cohditions
   r2=xr2
   if(r2.lt.rc2) then 最截止
   r2i=1/r2
   r6i=r2i3
   ff=48r2rr6r(r6i-1)-ecut Lennard-Jopes涌现
   f(i)=f(i)+ffxr 更新力
   f(j)=f(j)-ffxr
   en=en+4r6i(r6i-1)-ecut 更新能量
   endif
   enddo
   enddo
   return
   end
   Comments to thisalgorithm：
   评论thisalgorithm
10. For efficiencu reasons the factors 4 and 48 are usually taken out of the force loop and taken into account at the end of the calculation for the energy.
11. 由于效率的原因，系数4和48通常被从力环中取出，并在能量计算结束时加以考虑。
12. The term ecut is the value of the potential at r = rc; for the Lennard-Jones potential, we have
13. ecut是r = rc处的势的值;对于Lennard-Jones的潜力，我们有
    ecut=4(1/r[上12下c]-1/r[上6下c])
    the program particularly efficient and we shall, in fact, consider all possible pairs of particles explicitly.We first compute the current distance in the x, y, and zdirections betweer each pair of particles i and i. These distances are indicated by xr. As in the Monte Carlo case, we use periodic boundary conditions (see section 3.2.2). In the present example, we use a cutoff at a distance rc in the explicit calculation of intermolecular interactions, where rc is chosen to be less than half the diameter of the periodic box. In that case we can always limit the evaluatior
    这个程序特别有效，事实上，我们将明确地考虑所有可能的粒子对。我们首先计算每一对粒子i和i在x、y和z方向上的当前距离，这些距离用xr表示。在蒙特卡罗的情况下，我们使用周期边界条件(见3.2.2节)。在本例中，我们在分子间相互作用的显式计算中使用距离rc的截断，其中rc被选择为小于周期箱直径的一半。在这种情况下，我们总是可以限制求值程序
    of intermolecular interactions between i and i to the interaction between i and the nearest periodic image of j.In the present case, the diameter of the periodic box is denoted by box. If we use simple cubic periodic boundary conditions, the distance in any direction between i and the nearest image of i should alwavs be less (in absolute value) than box/2. A compact way to compute the distance between i and the nearest periodic image of i uses the nearest integer function (nint (x) in FORTRAN). The nint function simply rounds a real number to the nearest integer.2 Starting with the x-distance (say) between i and any periodio image of i, xr, we compute the x-distance between i and the nearest image of i as xr=xr-boxnint (xr/box). Having thus computed all Cartesian components of rii, the vector distance between i and the nearest image of j, we compute r (denoted by r2 in the program). Next we test if , is less than r2, the square of the cutoff radius. If not, we immediately skip to the next value of i. It perhaps is worth emphasizing that we do not compute Iril itself, because this would be both unnecessary and expensive (as it would involve the evaluation of a square root).If a given pair of particles is close enough to interact, we must compute the force between these particles, and the contribution to the potential energy. Suppose that we wish to compute the x-component of the force
    i与i之间的分子间相互作用与i与最近的j周期像之间的相互作用。在本例中，周期盒的直径用box表示。如果我们使用简单的三次周期边界条件，i和i的最近图像在任何方向上的距离都应该比box/2小(绝对值)。一种紧凑的计算i和i的最近周期图像之间的距离的方法使用最近的整数函数(nint (x)在FORTRAN)。nint函数只是将一个实数四舍五入为最接近的整数。2从i和i的任意周期像xr的x距离开始，我们计算i和i的最近像的x距离为xr=xr-boxnint (xr/box)。这样计算了rii的所有笛卡尔分量，即i与j的最近像之间的向量距离，我们计算r(在程序中用r2表示)。接下来我们测试，是否，小于r2，截断半径的平方。如果没有，我们立即跳到i的下一个值。或许值得强调的是，我们不计算Iril本身，因为这既不必要又昂贵(因为它将涉及求平方根)。如果一对给定的粒子距离足够近，可以相互作用，我们必须计算这些粒子之间的力，以及它们对势能的贡献。假设我们要计算力的x分量**