Starting from the perfect crystal at equilibrium lattice constant a0, we can assign initial velocities to the atoms and perform MD simulations. In the simplest simulation, no thermostat is introduced to regulate the temperature, and no barostat is introduced to regulate the stress. The simulation then corresponds to the NVE ensemble, where the number of particles N, the cell volume V (as well as shape), and total energy E are conserved. This simulation is usu­ally performed as a benchmark to ensure that the numerical integrator is implemented correctly and that the time step is small enough.

The instantaneous temperature Tinst is defined in terms of the instantaneous kinetic energy K through the relation K = (3N/2)kBTinst, where kB is Boltzmann’s constant. Therefore, the velocity can be initialized by assigning random numbers to each component of every atom and scaling them so that Tinst matches the desired temperature. In practice, Tinst is usually set to twice the desired temperature for MD simulations of solids, because approximately half of the kinetic energy flows to the potential energy as the solids reach thermal equilibrium. We also need to subtract appropriate constants from the x, y z components of the initial velocities to make sure the center-of-mass linear momentum of the entire cell is zero. When the solid contains surfaces and is free to rotate (e. g., a nanoparticle or a nanowire), care must be taken to ensure that the center-of-mass angular momentum is also zero.

Figure 7(a) plots the instantaneous temperature as a function of time, for an MD simulation starting with a perfect crystal and Tinst = 600 K, using the Velocity Verlet integrator13 with a time step of A t = 1 fs. After 1 ps, the temperature of the simulation cell is equilibrated around 300 K. Due to the finite time step At, the total energy E, which should be a conserved quantity in Hamiltonian dynamics, fluctuates during the MD simulation. In this simula­tion, the total energy fluctuation is <2 x 10-eV per atom, after equilibrium has been reached (t> 1 ps). There is also zero long-term drift of the total energy. This is an advantage of symplectic integrators11, and also indicates that the time step is small enough.

The stress of the simulation cell can be computed by averaging the Virial stress for time between 1 and 10 ps. A hydrostatic pressure P = — (sxx + + ozz)/3 =

1.33 ± 0.01GPa is obtained. The compressive stress develops because the crystal is constrained at the zero-temperature lattice constant. A convenient way to find the equilibrium lattice constant at finite tem­perature is to introduce a barostat to adjust the vol­ume of the simulation cell. It is also convenient to introduce a thermostat to regulate the temperature of the simulation cell. When both the barostat and thermostat are applied, the simulation corresponds to the NPT ensemble.

The Nose-Hoover thermostat11,33,34 is widely used for MD simulations in NVT and NPT ensem­bles. However, care must be taken when applying it to perfect crystals at medium-to-low temperatures, in which the interaction between solid atoms is close to harmonic. In this case, the Nose-Hoover thermo­stat has difficulty in correctly sampling the equi­librium distribution in phase space, as indicated by periodic oscillation of the instantaneous temperature.

2	2
03 0
0 20 40 60 t (Ps)
0
80	100
(b)

Figure 7 (a) Instantaneous temperature Tinst and Virial pressure p as functions of time in an NVE simulation with initial temperature at 600 K. (b) Tinst and P in a series of NVT at T = 300 K, where the simulation cell length L is adjusted according to the averaged value of P.

The Nose-Hoover chain35 method has been devel­oped to address this problem.

The Parrinello-Rahman19 method is a widely used barostat for MD simulations. However, periodic oscillations in box size are usually observed during equilibration of solids. This oscillation can take a very long time to die out, requiring an unreasonably long time to reach equilibrium (after which meaningful data can be collected). A viscous damping term is usually added to the box degree of freedom to accel­erate the speed of equilibration. Here, we avoid the problem by performing a series of NVT simulations, each one lasting for 1 ps using the Nose-Hoover chain method with Velocity Verlet integrator and At = 1 fs. Before starting each new simulation, the simulation box is subjected to an additional hydro­static elastic strain of e = (P)/B0, where (P) is the average Virial pressure of the previous simulation, where B = 2000GPa is an empirical parameter.

The instantaneous temperature and Virial pressure during 100 of these NVT simulations are plotted in Figure 7(b). The instantaneous temperature fluctuates near the desired temperature (300 K) nearly from the beginning of the simulation. The Virial pressure is well relaxed to zero at t = 20 ps. The average box size from 50 to 100 ps is L = 16.5625 A, which is larger than the initial value of 16.5290 A. This means that the normal strain caused by thermal expansion at 300 K is exx = 0.00203. Hence, the coefficient of thermal expan­sion is estimated to be a = exx/T = 6.8 x 10_6K_1.35