1、2.3.4 Potential truncation and its corrections The chemical potential may also be related to g(r) 13200ln()4( ) ( , )Bk Tdr v r g rdrParameter coupling the two atomThermal de Broglie wavelenth In the MD or MC, in order to save time, we introduce a spherical cutoff. The cutoff distance should be suff

2、icient large to ensure that it is a small perturbation, but it must be no greater than BOXL/2 for the consistency with minimum image convention. Long-range corrections Therefore simulation calculates only the part of integration in energy, pressure and chemical potential equations (from 0 to rC). Th

3、ey should be corrected. Assuming the radial distribution function as unit in the distance r rC, the energy, pressure and chemical potential can be corrected by2full2( )CCrEENr v r dr23full(2/3)( )/CCrPPr dv rdr dr2full4( )CCrr v r drSubscript C refers to the quantities calculated from simulation.2.3

4、.5 Time correlation functions and transport coefficients Correlation between two different quantities A and B are measured in the usual statistical sense, via the correlation coefficient cAB/( ) ( )ABcA BAB ensAAA2222( )ensensensAAAA The absolute value of cAB lies between 0 and 1, with values close

5、to 1 indicating a high degree of correlation. By considering A and B to be evaluated at two different times, the resulting quantity is a function of the time difference and called time correlation function cAB(t).Time correlation functions and transport coefficients-continue For identical phase func

6、tions, cAA(t) is called an autocorrelation function and its time integral (from t =0 to t =) is a correlation time tA. These functions are of great interest in computer simulation because1) They give a clear picture of the dynamics in a fluid2) tA is related to macroscopic transport coefficients3) F

7、ourier transform is related to experimental spectra. The non-normalized correlation function is defined as( )( )(0)( ( )( (0)ABensensCtA tBAtBTime correlation functions and transport coefficients-continue So that ( )( )/( ) ( )ABABctCtAB2( )( )/( )( )/(0)AAAAAAAActCtACtC Transport coefficients are d

8、efined in terms of the response of a system to a perturbation.0( ) (0)dt A t ATransport coefficientA variable appearing in the perturbation term in the Hamiltonian. Einstein relation22( ( )(0)tA tADiffusion coefficient We can use time correlation function or Einstein relation to calculate transport

9、coefficients in computer simulations. Or go back to first principles and conducting a suitable non-equilibrium simulation. For using equilibrium MD, we give here the equations in the microcanonical ensemble, for a fluid composed of N identical molecules. The diffusion coefficient is given by Green-K

10、ubo formula or Einstein relationship01( )(0)3iiDdttvvThe center-of-mass velocity of a single molecules212( )(0)3iitDtrrMolecular position at time t.3directionsVelocity autocorrelation function and mean squared displacement 01234-0.20.00.81.0 =0.863 g/cm3 =1.396 g/cm3Cvv(t)t / psVelocity aut

11、ocorrelation functions for liquid argonVariation in mean squared displacement during a MD simulation of argonShear viscosity0( )(0)BVdt Ft Fk T Green-Kubo Formula:1/iiiiiiiFp pmr fV,xy yz zx1/iiiijijiij iFp pmrfVAn off-diagonal element of the pressure tensor Einstein relationship22( )(0)BVtLtLk T1ii

12、iLr pVBulk viscosity Green-Kubo formula04( )(0)3VBVdtFtFk T( )( )( )FtFtFFtPIn calculation, we use 11( )Tr ( )33F tF tFand( )( )F tF tP Einstein relationship242 ()( )(0)3VBVtLtLPtk TThermal conductivity20( )(0)TBVdt jt jk TA component of the energy current, i.e., the time derivative of .1()iiiirVMak

13、e no contribution if ri=0, as in the case in a normal one-component MD simulation. Green-Kubo Formula Einstein relationship222( )(0)TBVttk THow to get i ?How to get the energy per molecule? Assuming pairwise potentials, the potential energy of two molecules is taken to be divided equally between the

14、m:212/2( )iiiijj ipmv rMomenta of molecule i. The quantities F, F and are multi-particle properties, properties of the system as a whole, and so additional averaging over the N particles is not possible. Consequently viscosity and thermal conductivity are subject to much greater statistical imprecis

15、ion than D.2.3.6 Quantum corrections Most of this course will deal with the computer simulation of system within the classical approximation. Even within the limitations of a classical simulation, it is still possible to estimate quantum corrections of thermodynamic functions. This is achieved by ex

16、panding the partition function in powers of Plancks constant = h/2.22311exp( ) 1( )!24NNVTNiQdVVNmirrrrTruncated at forth-order and the first- and third-orders are zero because the sum of forces is zero.Quantum corrections-continue The expansion in the previous slide accounts for the leading quantum

17、-mechanical diffraction effects; other effects, such as exchange, are small for most of interest. It leads to the correction to Helmholtz free energy 2221trans24/iANfmMean square force on one atom in the simulation. Better statistics are obtained by averaging over all N atoms.22222tr20( )2( )( )( )(

18、 )4812NNd v rdv rAd g rv rr g rdrdrrdrr For hard spheres, it amounts to replacing the hard sphere diameter by +/81/2.Quantum corrections-molecular systems The quantum correction equation in previous slide is also the translational correction for a system of N molecules, where it is understood that m

19、 stands for molecular mass and is the mean square force acting on the molecular center of mass. For a molecular system, additional corrections must be applied to take account of rotational motion. For linear molecules:2221rot24/6iANINIMean square torque Moment of inertia2.3.7 Monitoring the equilibr

20、ation Any molecular simulations need equilibration after initialization and before production. The purpose of the equilibration phase is to enable the system to evolve from the starting configuration to reach equilibrium. How to monitor the equilibration ? Using mean squared displacement; Using radi

21、al distribution function Using order parameter A solid latticeMonitoring the equilibration by the mean squared displacement and RDF Mean squared displacement (MSD) For a fluid, with no underlying regular structure, MSD gradually increases with time. For a solid, however, the MSD oscillates about a m

22、ean value. RDF Particularly useful for detecting the presence of two phases. It is characterised by a larger than expected first peak and by the fact that g(r) does not decay towards a value of 1 at long distances.Order parameter One way to measure translational oeder in a system initially in a FCC

23、lattice was suggested by Verlet, whose order parameter is ()/3xyz141cosNixixNaLength of one edge of the unit cell Initially, all of the coordinates are multiples of a/2 and so =1. As the simulation proceeds the should gradually decrease to a value of zero, indicating that the atoms are distributed r

24、andomly.Order parameterVariation in Verlet order parameter during the equilibration phase of a MD simulation.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-b

25、ody system.Means that the nuclear motion of the constituent particles obeys the laws of classical mechanics.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 ex

26、periments, we proceed: Preparing a sample of the material studied; Connecting the sample to a measuring instrument; 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. I

27、n a MD simulation, we follow exactly the same approach.e.g., a thermometer, manometer, or viscometer, etc.MD 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.Inter

28、action energy potential, pair potential is frequently usedEquilibrate the system After equilibration, perform 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,

29、the system undergoes an irreversible change during the experiment, or we do not measure what we think. How to 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 over all degrees of freedom, Nf, we have: 21( )( )22NiiBfimv tk T tN21( )( )NiiiBfmv tT tk NT