模拟物理导论模拟物理导论凝聚态物质的数值模拟方法(V)马红孺http:/.cn2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型Molecular systems:In most cases the interaction part can be approximated by pair interactions:One famous example is the Lennard-Jones potential2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型A very important quantity in statistical mechanics is the pair correlation function g( r,r0 ), defined aswhereIt may also be written as2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型For a homogeneous system the pair correlation function depends only on the distance between r and r0 . In this case we denote it as g(r) .The g( r,r0 ) is proportional to the probability that given a particle at point r and find another particle at point r0 . At large distance g(r) tends to 1, we may define the total correlation functionThe Fourier transform of the above function gives the static structure function (or structure factor)2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型The structure function is defined as the correlation function of Fourier component of density fluctuations The density is defined as: ,and the density fluctuation is: ,and its Fourier component is: 2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型当体积趋于无限时, 红颜色的部分可以略去.2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型The structure factor can be measured directly by scattering experiments and can also be calculated by simulations.Many physical quantities can be expressed in terms of the pair correlation functions, for example the energy in NVT ensemble isThe pressure is2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型The compressibilityThis expression can be derived from the fluctuations of particle numbersSince so2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型On the other hand, it can be proved thatWe have the final result.The time correlation function is the correlations of two physical quantities at different times,For systems at equilibrium the time correlation function is a function of the time difference only and can be written as2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型The velocity auto correlation function of the ith particle isThis can be derived from the definition (we will back to this point) Which is related to the diffusion constant of the particle. which holds for large t.2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型In general, transport coefficient is defined in terms of the response of a system to a perturbation.where is the transport coefficient, and A is a physical variable appearing in the perturbation Hamiltonian. There is also an Einstein relation associated with this kind of expressionwhich holds for large t , (t , where is the relaxation time of ).2003-10-21上海交通大学理论物理研究所马红孺分子模型分子模型The shear viscosity is given byorHereThe negative of P is often called stress tensor.2003-10-21上海交通大学理论物理研究所马红孺Monte Carlo 模拟模拟Monte Carlo simulation of Particle Systems粒子系统的Monte Carlo 模拟和自旋系统原则上是一样的。Metropolis 算法为：1，随机或顺序选取一个粒子，其位置矢量为 ，对此粒子做移动2，计算前后的能量差，决定是否接受移动。3，在达到平衡后，收集数据，计算物理量。2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟Molecular dynamics simulationsMD method is essentially the integration of the equation of motion of the classical many-particle system in a period of time. The trajectories of the system in the phase space are thus obtained and averages of the trajectories give various physical properties. Since we work on real dynamics in MD simulations we can also study the dynamic properties of the system such as relaxation to equilibrium, transport etc.2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟Consider a rectangular volume of L1 L2 L3, with Nclassical particles put in. The particles are interact with each other. In principle, the interaction include pair interactions, three body interactions as well as many body interactions. For simplicity we will consider here only pair interactions. In this case each particle feel a force by all other particles and we further assume the force is depend only on distances from other particles and for each pair the force directed along the line join the pair of particles. So the force on the ith particle iswhere is an unit vector along rj-ri.2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟Periodic boundary condition(PBC)where L are vectors along the edges of the rectangular system volume and the sum over n is with all integers n . Usually this sum is the most time consuming part in a simulation.2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟General procedure of MD (NVE ensemble)1. Initialize;2. Start simulation and let the system reach equilibrium;3. Continue simulation and store results.2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟Initialize: 1, Specify the number of particles and interaction;2, Setup the simulation box;3, Specify the total energy of the system;4, Assign position and momenta of each particle. a, In many cases we assign particles in a FCC lattice,If we use cubic unit cell and cube BOX then the number of particles per unit cell is 4, and the total number of particles are a 4M3, M=1,2,3, . That is we may simulation systems with total number of particles N=108, 256, 500, 864, . b, The velocities of particles are draw from a Maxwell distribution with the specified temperature. This is accomplished by drawing the three components of the velocity from the Gaussian distribution.2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟The distribution of the x-component of velocity isDraw numbers from a Gaussian:Consider:Thenwhere v2=vx2+vy2 and2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟So the distribution of vx and vy may be obtained from v and . The distribution of v:The distribution of is uniform in the interval 0,2 . 2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟Generate random numbers for a given distributionFor a given distribution P(y) we consider how to get a random number y draw from P(y) from a random number x draw from uniform 0,1, i.e., we are going to find a function f(x), from which for a series of random numbers x distributed uniformly in the interval 0,1, y=f(x) will distributed according to P(y).2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟thenSinceExponential distribution2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟The distribution of v:2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟Draw random numbers uniformly distributed in the interval 0,2. Another method of draw random numbers in the Gaussian distribution is through the following empirical methods. Consider the distribution2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟According to the central limit theorem, if we draw uniform random numbers ri in interval 0,1, and define a variablewhen n! 1 the distribution of is the Gaussian distributionIf we take n=12, we get2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟After the generation of the velocity of each particle, we may shift the velocity so that the total momentum is zero.The standard Verlet algorithm is the first successful method in history and still wide used today in different forms. It isTo start the integration we need r(h), given by2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟Variations of this method areandBoth of these variations are mathematically equivalent to the original one but more stable under finite precision arithmetic.2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟The temperature of the system is given by the equal partition theorem, that is the average of kinetic energy of each degree of freedom is half kBT, The N-1 is due to the conservation of the total momentum which reduce the degree of freedom by 3.To reach the desired temperature we may scale the velocity at every few steps of integration2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟After the system reach to equilibrium the integration continue in the same method as above without scaling of velocity. The data are stored or accumulated for the calculating physical properties. The static properties of physical quantity A is given by time average2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟here A is the value of A at th time step. Usually the data stored in each step include:1, the kinetic energy2, the potential energy 3, the virial 2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟We also needs data to calculate the pair correlation function, this is done by divide the interval 0,r into sub intervals ir, (i+1)r, at each stage of updating, add the number of pairs with separation in the interval ir, (i+1)r, to an array n(i) and find the average value after simulation, the pair correlation function given by2003-10-21上海交通大学理论物理研究所马红孺分子动力学模拟分子动力学模拟练习:1, Write programs for the two methods to generate Guassian random numbers.2, Compare the two methods for efficiency and quality.3, Generate random numbers with exponential distribution by means of the transformation method described before and check the quality. 2003-10-21上海交通大学理论物理研究所马红孺