Computational Physics (Lecture 20 )

1. Computational Physics(Lecture 20) PHY4370

2. Constant pressure, temperature, and bond length • when we want to compare the simulation results with the experimental measurements: • environmental parameters, • such as temperature or pressure, are the result of thermal or mechanical contact of the system with the environment or the result of the equilibrium of the system. • For a macroscopic system, we deal with average quantities by means of statistical mechanics, • so it is highly desirable for the simulation environment to be as close as possible to a specific ensemble.

3. The example ofatomic cluster systems is closely related to the microcanonical ensemble, • which has the total energy of the system conserved. • It is important that we also be able to find ways to deal with constant-temperature and/or constant-pressure conditions.

4. There have been many efforts to model realistic systems with simulation boxes by introducing some specific procedures. • However, all these procedures do not introduce any new concepts in physics. • They are merely numerical techniques to make the simulation boxes as close as possible to the physical systems under study.

5. Constant pressure: the Andersen scheme • The scheme for dealing with a constant-pressure environment was devised by Andersen (1980) with the introduction of an environmental variable, the instantaneous volume of the system, in the effective Lagrangian. • When the Lagrange equation is applied to the Lagrangian, the equations of motion for the coordinates of the particles and the volume result. • The constant pressure is then a result of the zero average of the second-order time derivative of the volume.

6. where the last two terms are added to deal with the constant pressure from the environment.

7. The parameter M can be viewed here as an effective inertia associated with the expansion and contraction of the volume Ω, • and P_0 is the external pressure, which introduces a potential energy P_0 Ωto the system under the assumption • that the system is in contact with the constant-pressure environment.

8. The coordinate of each particle r_iis rescaled with the dimension of the simulation box, L = Ω ^1/3, because the distance between any two particles changes with L and the coordinates independent of the volume are then given by xi = ri /L, which are not directly related to the changing volume.

9. This effective Lagrangianis not the result of any new physical principles or concepts but is merely a method of modeling the effect of the environment in realistic systems. • Now if we apply the Lagrange equation to the above Lagrangian, we obtain the equations of motion for the particles and the volume Ω,

11. which are usually evaluated immediately after the volume and the coordinates are updated.

12. In practice, we first need to set up the initial positions and velocities of the particles and the initial volume and its time derivative. • The initial volume is determined from the given particle number and density, • and its initial time derivative is usually set to be zero. • The initial coordinates of the particles are usually arranged on a densely packed lattice, • for example, a face-centered cubic lattice, and the initial velocities are usually drawn from the Maxwell distribution. • One should test the program with different Ms in order to find the value of M that minimizes the fluctuation.

13. A generalization of the Andersen constant-pressure scheme was introduced by Parrinello and Rahman (1980; 1981) to allow the shape of the simulation box to change as well. • This generalization is important in the study of the structural phase transition. • With the shape of the simulation box allowed to vary, the particles can easily move to the lattice points of the structure with the lowest free energy.

14. The idea of Parrinello and Rahman

15. Constant temperature: the Nosėscheme • The constant-pressure scheme discussed above is usually performed with an ad hoc constant-temperature constraint, which is done by rescaling the velocities during the simulation to ensure the relation between the total kinetic energy and the desired temperature in the canonical ensemble.

16. This rescaling can be shown to be equivalent to a force constraint up to first order in the time step τ . • The constraint method for the constant-temperature simulation is achieved by introducing an artificial force −η p_i, • which is similar to a frictional force if η is greater than zero or to a heating process if η is less than zero.

17. The equations of motion are modified under this force to

18. The most popular constant-temperature scheme is that of Nose (1984a; 1984b), • who introduced a fictitious dynamical variable to take the constant-temperature environment into account. • The idea is very similar to that of Andersen for the constant-pressure case. • In fact, one can put both fictitious variables together to have simulations for constant pressure and constant temperature together. • Herewe will briefly discuss the Nose scheme.

19. where s and v_sare the coordinate and velocity of an introduced fictitious variable that rescales the time and the kinetic energy in order to have the constraint of the canonical ensemble satisfied. • The rescaling is achieved by replacing the time element dtwith dt/s and holding the coordinates unchanged, that is, x_i= r_i. • Thevelocity is rescaled with time: • We can then obtain the equation of motion for the coordinate xi and the variable s by applying the Lagrange equation. • Hoover(1985; 1999) showed that the Nose Lagrangian leads to a set of equations very similar to the result of the constraint force scheme discussed at the beginning of this subsection.

20. The Nose equations of motion are given in the Hoover version by with s0 as the initial value of s at t = t0.

21. We can discretize the above equation set easily with either the Verlet algorithm or one of the Gear schemes. • Note that the behavior of the parameter s is no longer directly related to the simulation; • It is merely a parameter Nose introduced to accomplish the microscopic processes happening in the constant-temperature environment.

22. We can also combine the Andersen constant-pressure scheme with the Nose constant-temperature scheme in a single effective Lagrangian

23. Constant bond length • Another issue we have to deal with in practice is that for large molecular systems, such as biopolymers, the bond length of a pair of nearest neighbors does not change very much even though the angle between a pair of nearest bonds does. • If we want to obtain accurate simulation results, we have to choose a time step much smaller than the period of the vibration of each pair of atoms. • This costs a lot of computing time and might exclude the applicability of the simulation to more complicated systems, such as biopolymers.

24. A procedure commonly known as the SHAKE algorithm (Ryckaert, Ciccotti,andBerendsen, 1977; van Gunsteren and Berendsen, 1977) was introduced to deal with the constraint on the distance between a pair of particles in the system. • The idea of this procedure is to adjust each pair of particles iteratively to have in each time step.

25. Here dijis the distance constraint between the ith and j thparticles and Δis the tolerance in the simulation. • The adjustment of the position of each particle is performed after each time step of the molecular dynamics simulation. • Assume that we are working on a specific pair of particles and for the l thconstraint and that we would like to have

26. where rij= r j − riis the new position vector difference after a molecular time step starting from r(0)ij and the adjustments for the l − 1 constraints have been completed. • Here δrij= δrj− δriis the total amount of adjustment needed for both particles.

27. which is reasonable, because gi j is a small number during the simulation.

28. More importantly, by the end of the iteration, all the constraints will be satisfied as well; • all gi j go to zero at the convergence. • Equation (8.100) is used to estimate each gi j for each constraint in each iteration. • After one has the estimate of gij for each constraint, the positions of the relevant particles are all adjusted. • Theadjustmentshave to be performed several times until the convergence is reached. • For more details on the algorithm, see Ryckaertet al. (1977). • This procedure has been used in the simulation of chain-like systems as well as of proteins and nucleic acids.

29. Structure and dynamics of real materials • A numerical simulation of a specific material starts with a determination of the interaction potential in the system. • In most cases, the interaction potential is formulated in a parameterized form, which is usually determined separately from the available experimental data, first principles calculations, and condition of the system under study. • The accuracy of the interaction potential determines the validity of the simulation results. • Accurate model potentials have been developed for many realistic materials, for example, the Au (100) surface (Ercolessi, Tosatti, • and Parrinello, 1986) and Si3N4 ceramics (Vashishta et al., 1995).

30. We then need to set up a simulation box under the periodic boundary condition. • Because most experiments are performed in a constant-pressure environment, we typically use the constant-pressure scheme developed by Andersen (1980) or its generalization (Parrinello and Rahman, 1980; 1981). • The size of the simulation box has to be decided together with the available computing resources and the accuracy required for the quantities to be evaluated.

31. The initial positions of the particles are usually assigned at the lattice points of a closely packed structure, • for example, a face-centered cubic structure. • The initial velocities of the particles are drawn from the Maxwell distribution for a given temperature. • The temperature can be changed by rescaling the velocities. • This is extremely useful in the study of phase transitions with varying temperature, • such as the transition between different lattice structures, glass transition under quenching, or liquid–solid transition when the system is cooled down slowly.

32. The advantage of simulation over actual experiment also shows up when we want to observe some behavior that is not achievable experimentally due to the limitations of the technique or equipment. • For example, the glass transition in the Lennard–Jones system is observed in molecular dynamics simulations but not in the experiments for liquid Ar, because the necessary quenching rate is so high that it is impossible to achieve it experimentally.

34. where g(r) is the pair distribution and S(k) is the angular average of S(k).