1、Chapter 3 Molecular Dynamics Simulation3.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.N
2、ewtons 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
4、approach 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 usedEquilibrate the system After equilibration, perform
5、 the actual 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
6、 measure an 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
7、over all degrees of freedom, Nf, we have: 21( )( )22NiiBfimv tk T tN21( )( )NiiiBfmv tT tk NThe 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%.1/fNAverage over many fluctuations In which c
8、ase should we worry about quantum effects?When we consider the the translational 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 co
9、urse, you can read the references on the coming slide. References for Molecular 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, Protei
10、ns. A Theoretical Perspective of Dynamics, Structure and Thermodynamics. Advances in chemical Physics. Volume LXXI. JM Haile, 1992, Molecular Dynamics Simulation. Elementary Methods.4.2 Molecular Dynamics: A ProgramThe best introduction to MD is to consider a simple program. We keep it as simple as
11、possible to illustrate some important features of MD. It is constructed as: read in the parameters that 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
12、 particles; Integrate Newtons equations of motion; measure the quantities in the present time;After completion of the central loop, compute and print the averages of measured quantities, and stop. Central loop, the core of the simulationA simple MD programProgram MD simple MD programcall init initia
13、lizationt=0 set timedo while (t.lt.tmax) MD loop call force(f,en) determine the forces call integrate(f,en) integrate equations of motion t=t+delt update time call sample sample averagesenddostopSubroutine init, force, and integrate will be described later. Subroutine sample is used to calculate ave
14、rages of quantities of interest like pressure and temperature, etc.4.2.1 Initiallization To start the simulation, we should 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 res
15、ult in an appreciable overlap; choose positions compatible with the structure that we are aiming to simulate.As the equilibrium 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 so
16、lid state of a model system, it is logical to prepare the system in the crystal structure of interest. If we are interest 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
17、Initial positions randomly; use final configuration of a liquid at a higher temperature or lower density.How to avoid it ?Initial 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 fo
18、r a new run and adjust the temperature and density to the desired values.Well-equilibrated The observed values should not depend on the initial condition. If does, there are two possibilities: The system really behaves non-ergodically.In glassy materials or low-temperature, disordered alloys. Sampli
19、ng configuration space is inadequate-not yet reach equilibrium.Lattice positions(FCC)C CELL = LENGTH / REAL( NC ) CELL2 = 0.5 * CELLC * SUBLATTICE 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
20、(3) = CELL2C * SUBLATTICE D * RX0(4) = CELL2 RY0(4) = 0.0 RZ0(4) = CELL2ADBCCELLLENGTHC * CONSTRUCT THE LATTICE FROM THE UNIT CELL * M = 0 DO 100 IZ = 1, NC DO 100 IY = 1, NC DO 100 IX = 1, NC DO 110 IREF = 1, 4 RX0(IREF+M) = RX0(IREF) + CELL * REAL ( IX - 1 ) RY0(IREF+M) = RY0(IREF) + CELL * REAL (
21、 IY - 1 ) RZ0(IREF+M) = RZ0(IREF) + CELL * REAL ( IZ - 1 )110 CONTINUE M = M + 4100 CONTINUEMove one cell in one directionInitial velocities There are two ways to initial velocities: one is to use random number distributed uniformly in the 0.5,0.5, and the other is done by randomly selecting from a
22、Maxwell-Boltzmann distribution. In the first way, the initial velocity distribution is Maxwellian neither in shape nor even in width. Shift all the velocitiestotal momentum is zeroScale all the velocitiesWith a factor (T/T(t)1/2T(t)=TDesired temperatureOn equilibration it will become a MB distributi
23、onInitialization of a MD ProgramSubroutine 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) velocity center of mass SUMV2=SUMV2+VX(I)*2+VY(I)*2+VZ(I)*2ENDDOSUMV=SUMV/FLOAT(N
24、PART)SUMV2=SUMV2/FLOAT(NPART)FS=SQRT(3.*T/SUMV2) Scale factor of the velocitiestemperatureInitialization of a MD Program-continueDO I=1, NPART VX(I)=(VX(I)-SUMV)*FS shift and scale RXM=RX(I)-VX(I)*DT Position previous time stepEnddoReturnendWe do not really use the velocities in program to solve equ
25、ations 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 this sentence to generate previous positions.4.2.2 The Force Calculation
26、 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 distance. If we use no tricks, the time needed for the evaluation of the forc
27、es 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 coordinates, Hamiltons equations become/iiimrp( )iiV r irpf If a given
28、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 forcexyz irijk( )( )( )ijijijxururrrxrx fThe Force Calculation-continue222rxyz222rxxxyz( )ijijxurxrr fFor a Lennard-Jones system
29、 (in reduced units),21264811( )(0.5)xxfrrrrNote that here x stands for xi-xj, and r is the distance between i and j.Program for the calculation of the forceSubroutine 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,npar
30、t loop over all pairs RXIJ=RX(I)-RX(J) RXIJ=RXIJ-BOXL*ANINT(RXIJ/BOXL) RIJSQ=RXIJ*2+RYIJ*2+RZIJ*2Periodic boundary conditionProgram 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 for
31、ce F(J)=F(J)-FF*RXIJ EN=EN+4*R6I*(R6I-1)-ECUT update energy ENDIF ENDDOENDORETURN A simple Lennard-Jones force routineEND can be found in f17.for4.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 in
32、troduce the so-called Verlet algorithm, which is 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,324( )()( )( )( )()23!tttttttttOtm frrvr Similarly,324( )()( )( )( )()23!tttttttttOtm frrvr Integrat
33、ing the Equations of Motion-Verlet algorithm Summing the two equations in the above slide, we obtain: 24( )()()2 ( )()tttttttOtmfrrr2( )()2 ( )()tttttttmfrrr The estimate of the new position contains an error that is of order t4. We does not use the velocity to compute the new position. However, it can be derived2()()( )()2tttttOttrrvVerlet algorithm-calculation procedurervattttttttttttt Update timeCalculate v and a at present from present positionCalculate position at next time step from present and last positions, and present force.
