《分子动力学模拟基础汇总课件.ppt》由会员分享，可在线阅读，更多相关《分子动力学模拟基础汇总课件.ppt（25页珍藏版）》请在三一办公上搜索。
1、Chapter 3 Molecular Dynamics Simulation,3.1 Molecular Dynamics:The Idea,What is molecular dynamics?,It is a technique to compute the equilibrium and transport properties of a classical many-body system.,Means that the nuclear motion of the constituent particles obeys the laws of classical mechanics.
2、,Newtons law,Lagrangian equation,and Langevin equation.,This is an excellent approximation for a wide range of materials.,MD simulation is similar to real experiments,When we perform a real experiments,we proceed:,Preparing a sample of the material studied;,Connecting the sample to a measuring instr
3、ument;,Measuring the property of interest during a certain time;,If the measurements are subject to statistical noise,then the longer we average,the more accurate our measurement becomes.,In a MD simulation,we follow exactly the same approach.,e.g.,a thermometer,manometer,or viscometer,etc.,MD appro
4、ach,First,prepare a sample:select a model system consisting of N particles;,Second,solve Newtons equation of motion until the properties of the system no longer change with time.,Interaction energy potential,pair potential is frequently used,Equilibrate the system,After equilibration,perform the act
5、ual measurement.,Some of the most common mistakes in MD are similar to the mistakes that may be made in real experiments.,The sample is not prepared correctly,the time is too short,the system undergoes an irreversible change during the experiment,or we do not measure what we think.,How to measure an
6、 observable quantity?,To measure an observable quantity in a MD 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 system.,Let us take temperature as an example.Making use of the equipartition of energy over all degre
7、es of freedom,Nf,we have:,The relative fluctuations in the temperature will be of order.As Nf is typically of the order of 102-103,the statistical fluctuations are of the order of 5-10%.,Average over many fluctuations,In which case should we worry about quantum effects?,When we consider the the tran
8、slational or rotational motion of light atoms or molecules,or vibrational motion with a frequency such that h kBT.,He,H2,D2,etc.,Of course,our course of this vast subject is incomplete.If you need the knowledge beyond the course,you can read the references on the coming slide.,References for Molecul
9、ar Dynamics,MP Allen&DJ Tildesley,1987,Computer Simulation of Liquids.HC Berendsen&WF van Gunsteren,1984,Molecular Dynamics Simulations:Techniques&Approached,NATO ASI Series C123.CL Brooks et al.,1988,Proteins.A Theoretical Perspective of Dynamics,Structure and Thermodynamics.Advances in chemical Ph
10、ysics.Volume LXXI.JM Haile,1992,Molecular Dynamics Simulation.Elementary Methods.,4.2 Molecular Dynamics:A Program,The best introduction to MD is to consider a simple program.We keep it as simple as possible to illustrate some important features of MD.It is constructed as:,read in the parameters tha
11、t specify the conditions of the run(e.g.,initial temperature,number of particles,density,time step,etc.);Initialize the system(select initial ri and vi);Compute the energy and forces on all particles;Integrate Newtons equations of motion;measure the quantities in the present time;After completion of
12、 the central loop,compute and print the averages of measured quantities,and stop.,Central loop,the core of the simulation,A simple MD program,Program MD simple MD programcall init initializationt=0 set timedo while(t.lt.tmax)MD loop call force(f,en)determine the forces call integrate(f,en)integrate
13、equations of motion t=t+delt update time call sample sample averagesenddostop,Subroutine init,force,and integrate will be described later.Subroutine sample is used to calculate averages of quantities of interest like pressure and temperature,etc.,4.2.1 Initiallization,To start the simulation,we shou
14、ld assign initial positions and velocities to all particles in the system.,Often this is achieved by initially placing the particles on a cubic lattice.,Avoiding positions that result in an appreciable overlap;choose positions compatible with the structure that we are aiming to simulate.,As the equi
15、librium properties of the system do not depend on the choice of initial conditions,all reasonable initial conditions are in principle acceptable.,Initial positions,To simulate a solid state of a model system,it is logical to prepare the system in the crystal structure of interest.,If we are interest
16、 in the fluid phase,we simply prepare the system in any convenient crystal structure.,If the density close to the freezing,selection of crystal structure is unwise:meta-stable,Initial positions randomly;use final configuration of a liquid at a higher temperature or lower density.,How to avoid it?,In
17、itial positions-continue,In any cases,it is usually preferable to use the final configuration of an earlier simulation at a nearby state point as the starting configuration for a new run and adjust the temperature and density to the desired values.,Well-equilibrated,The observed values should not de
18、pend on the initial condition.If does,there are two possibilities:,The system really behaves non-ergodically.,In glassy materials or low-temperature,disordered alloys.,Sampling configuration space is inadequate-not yet reach equilibrium.,Lattice positions(FCC),C CELL=LENGTH/REAL(NC)CELL2=0.5*CELLC*S
19、UBLATTICE A*RX0(1)=0.0 RY0(1)=0.0 RZ0(1)=0.0C*SUBLATTICE B*RX0(2)=CELL2 RY0(2)=CELL2 RZ0(2)=0.0C*SUBLATTICE C*RX0(3)=0.0 RY0(3)=CELL2 RZ0(3)=CELL2C*SUBLATTICE D*RX0(4)=CELL2 RY0(4)=0.0 RZ0(4)=CELL2,A,D,B,C,CELL,LENGTH,C*CONSTRUCT THE LATTICE FROM THE UNIT CELL*M=0 DO 100 IZ=1,NC DO 100 IY=1,NC DO 10
20、0 IX=1,NC DO 110 IREF=1,4 RX0(IREF+M)=RX0(IREF)+CELL*REAL(IX-1)RY0(IREF+M)=RY0(IREF)+CELL*REAL(IY-1)RZ0(IREF+M)=RZ0(IREF)+CELL*REAL(IZ-1)110 CONTINUE M=M+4100 CONTINUE,Move one cell in one direction,Initial velocities,There are two ways to initial velocities:one is to use random number distributed u
21、niformly in the 0.5,0.5,and the other is done by randomly selecting from a Maxwell-Boltzmann distribution.,In the first way,the initial velocity distribution is Maxwellian neither in shape nor even in width.,Shift all the velocities,total momentum is zero,Scale all the velocitiesWith a factor(T/T(t)
22、1/2,T(t)=T,Desired temperature,On equilibration it will become a MB distribution,Initialization of a MD Program,Subroutine init initialization of MD programSUMVX=0 SUMV2=0DO I=1,NPART RX(I)=LATTICE-POS(I)Place the particles on a lattice VX(I)=(RANF()-0.5)Give random velocities SUMVX=SUMVX+VX(I)veloc
23、ity center of mass SUMV2=SUMV2+VX(I)*2+VY(I)*2+VZ(I)*2ENDDOSUMV=SUMV/FLOAT(NPART)SUMV2=SUMV2/FLOAT(NPART)FS=SQRT(3.*T/SUMV2)Scale factor of the velocities,temperature,Initialization of a MD Program-continue,DO I=1,NPART VX(I)=(VX(I)-SUMV)*FS shift and scale RXM=RX(I)-VX(I)*DT Position previous time
24、stepEnddoReturnend,We do not really use the velocities in program to solve equations of motion but positions at present RX(I)and previous RXM(I)time steps,combined with our knowledge of the force F acting on the particles,to predict the positions at the next time step.When start a simulation,we use
25、this sentence to generate previous positions.,4.2.2 The Force Calculation,The calculation of the force acting on every particle is the most time-consuming part of almost all MD simulations.,If we consider pairwise additive interaction,for a system of N particles,we must evaluate N(N-1)/2 pair distan
26、ce.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.,The Force Calculation-continue,For Cartesian c
27、oordinates,Hamiltons equations become,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.For the x-component of the force,The Force Calculation-continue,For a Lennard-Jones system(in reduced units),N
28、ote that here x stands for xi-xj,and r is the distance between i and j.,Program for the calculation of the force,Subroutine force(f,en)determine the force&energyEn=0 set energy to zeroDo i=1,npartf(i)=0 set forces to zeroEnddoDo I=1,npart-1 do j=I+1,npart loop over all pairs RXIJ=RX(I)-RX(J)RXIJ=RXI
29、J-BOXL*ANINT(RXIJ/BOXL)RIJSQ=RXIJ*2+RYIJ*2+RZIJ*2,Periodic boundary condition,Program for the calculation of the force-continue,IF(RIJSQ.LT.RCUTSQ)THEN test cut-off R2I=1./RIJSQ R6I=R2I*3 FF=48.*R2I*R6I*(R6I-0.5)LJ potential F(I)=F(I)+FF*RXIJ update force F(J)=F(J)-FF*RXIJ EN=EN+4*R6I*(R6I-1)-ECUT u
30、pdate energy ENDIF ENDDOENDORETURN A simple Lennard-Jones force routineEND can be found in f17.for,4.2.3 Integrating the Equations of Motion,Now that we have compute all forces between the particles,we can integrate Newtons equation of motion.I will introduce the so-called Verlet algorithm,which is
31、not only one of the simplest,but also usually the best.,To derive it,we start with a Taylor expansion of the coordinate of a particle,around time t,Similarly,Integrating the Equations of Motion-Verlet algorithm,Summing the two equations in the above slide,we obtain:,The estimate of the new position
32、contains an error that is of order t4.We does not use the velocity to compute the new position.However,it can be derived,Verlet algorithm-calculation procedure,Update time,Calculate v and a at present from present position,Calculate position at next time step from present and last positions,and present force.,
