T rot = L1wr / grot'B
3.6 Verification of the simulations
69
We can never guarantee that aur eode is free of errors, nor that the results meets with the theoretical restrietions. There are however several ways to gain increased eonfi-dence in the validity of the output. We have been coneemed with elimination of eode errors (of coursel), verification of liquid equilibrium, and investigation of the ergodic behaviour of the system. In what follows, we will describe the methods and various tests we have applied to assess the reliability of the results. It has not been possible to verify eaeh and every aspect of the simulation, but we have seleeted what we believe are the most important and the most interesting, balanced against the complexity in performanee.
3.6. 1 Code errors
We have compiled the eode with three different compilers, the f77 compiler at DEC Alpha 3000/400, the fe compiler at Convex 220 and tlle Convex Application Com-piler. The last compiler performs interprocedural analysis that tracks the flow of data and control between procedures. It is thus capable of, for instance, detecting errors that arise between different procedures, or inconsistencies in common blocks.
A system run under the same conditions at the Alpha and Convex computers would not generate identical trajectories as the organization of operations is different. Henee the accumulation of round-off errors is also different. We have compared the same eode at twodifferent computers and at two different levels of optimization at the same computer. Any differences in mean values will provide us with an insight to round-off errors. Table 3.2 give results for mixture temperatureT,mixture pressureP, and mixture potential energyEpfor a equimolar mixture of water and ethanol.
Table 3.2
NVT-results for a water-ethanol mixture of 128 molecules of each, run for a total of 50000 steps at steplength 1fs and with properties averaged for the last 35000 steps.
Sta.ndard deviationøfor each property is reported in parenthesis. All simulations are double precision. Ethanol potential slightly different from the OPLS-model applied in Chapter 6.
Computer Compiler <T> [K] <p> [MPa] <U> [J/g]
Convex 200 No optimization 292.7 (0.4) 81.5 (1.8) -1287.0 (0.5) DEC Alpha No optimization 292.7 (0.4) 81.9 (1.9) -1281.6 (0.5) Optimization 293.8 (0.4) 79.2 (1.9) -1291.5 (0.6)
70 Chapter3Simulation details
The deviations relative the non-optimized DEC Alpha simulation are largest for the pressure by nearly 3%, while the deviation in potential energy is less than 1%. We notice however that the deviation in potential energy is well outside ± 20'.Haile [28]
present an analysis of a l-dimensional Lennard-Jones system of100 atoms simulated for 10000and 20000 step, and our findings are of the same order as his. As we have simulated twice as lang and also use a more complicated potential, we find our results to be satisfying, and conclude that the deviations are of numerical nature.
Different mean values show that the system follow different trajectories for different computers or compilers. From inspection of the time development of potential energy we find (not presented) that within 1500 timestep the differences has grown larger than the fourth digit for different computers. The same difference does not show up until15000 timesteps for different compilers for the same computer.
From this brief investigation we also notice that the differences are larger with differ-ent optimization levels than with differdiffer-ent computers. Note that the non-optimized and the optimized simulations with DEC Alpha stations yield different temperatures, which may in part be responsible for the differences in pressure and configurational energy.
All Dur production runs are performedwith DEC Alpha stations at the highest optimi-zation level.
We have rewritten central parts of the codeto detect errors, and we have run our pro-gram for simple, atomic ane component systems. The mixture simulation has been tested out for various compositions of water in water, and alcohol in alcohol. The re-sults are the same as with ane component.
As total energy of the extended system, Equation (2.29), and the Hamiltonian of the isolated system are strictly eonstant in an ideal simulation, monitoring these constants of the motion is a powerful check of eode errors. To have a constant total energy, fluctuations in all contributions to total energy must cancel. Errors tend to show up quickly as a divergenee. Same minor departures from ideality must be tolerated though, as the finite difference equations only are correct to a definite order in h. We are satisfied with an aecuracy in the constants of motion to the fourth digit. The role of truncation can be evaluated by decreasing the timestep, then the accuracy should improve.
For each simulation we inspect the constant of the motion and estimate the drift as difference between lowest and highest va1ues, normally appearing at the beginning and end of simulation, respectively.
Chapter3 Simulation details 71 3.6.2 Equilibrium liquid state and stability, general requirements
To have thermodynamical equilibrium, the system should be mechanically and ther-mal stable. That is, all thermodynamic properties of the system should tluctuate around stable mean values, and the 1st law of thermodynamics should be satisfied. It is also important that the system is in the intended state (liquid).
For all simulations we have therefore controlied the performance of the following quantities:
• The translational order parameterpgiven by [15]
N
pek)=
kL
cos (k· ri) (3.37)i=l
where N is the number of molecules,riis the position of the ith molecule, andk is a reciprocal vector of the initial lattice. We start from a solid lattiee where phas the value unitY, but in the liquid state we expectrto tluctuate around zero with an amplitude oflIN.
• Energy conservation, the values of the constants of the motion defined in Equation (2.9) for the NVE simulations and in Equation (2.29) for the NVT simulations should be not onIy stable but also showing no tluctuations (see Subsection 3.6.1).
• Stability of time development of pressure and temperature,
• The velocity distribution initially taken from the Maxwell distribution should not depart from this. The x-component of linear velocity is sampled every 150th step, and the distribution is compared to the theoretical at the end of simulation.
(3.38) where miis mass of molecule i,vixis velocity of molecule i in the x-direction and Tis calculated temperature.
• The kinetie and rotational energies should be equally partitioned among all coordi-nate directions, and we save the square of rotational and translational velocities every 150th step,
• The total linear momentum is conserved in molecular dynamics simulations, the angular momentum is not. A sum over all components of linear and angular mo-menta is performed each 150th step, and saved to disk.
72
These and other tests are thoroughly discussed in [28].
Chapter3Simulation detaiis
3.6.3 Particular requirements for the extended system simulations
For the thermostatted simulations it is particularly interesting to see whether the simu-lations compare to the canonical ensemble. As discussed previously in Section 2.5, this is not guaranteed in advance. Also the extended system must satisfy the 1st law of thermodynamics, and there should be thermal equilibrium between the system and the reservoirs. We have thus looked at some quantities specific to the thermostatted simulation where the theoretical outcome is known.
• The reservoir parameters Stransand Srot are initialized to 1 at the beginning of the thermostatted simulation. Their average should not depart from this value, though they naturally will fluctuate, The parameters are saved to disk every 100th step, and coarse grain averages and standard deviations are calculated within the pro-gram.
• To have thermodynamical equilibrium between the system and the reservoirs, the rates of energy transferlltransand "rotmust fluctuate around zero. If not, there will be a net energy flux across the system boundaries. We therefore monitortlle time-development ofl1transand Tlror Not only should their mean values be close to zero separately, but the sum of their averages should also be closer to zero than the dividual ones. Coarse grain averages and standard deviations are reported, and in-stantaneous values are saved to disk every 100th step
• From the equipartition principle, the kinetie energy of a thermostat should be pro-portional to lhkBT.The potential energies of the reservoirs should ideally fluctuate around an average of zero.
• The running averages of the higher moments of translational, rotational, reservoir translational alld reservoir rotational energies. The first moment is simply the mean. If taken about the mean, the second is the variance, the third moment repre-sents skewness, and the fourth is kurtosis which measures the degree of peaked-ness [89]. Also the fifth moment is included. In a canonical simulation they are expected to approach the values given below [47]:
Chapter3 Simulation details
(~) = gj'ikBT
(KJ>
= ~(kBT)2(Kl>
= gj(kBT) 3(Kj> = ~
[3;i+6](kBT)4(Kl>
= gj[5gj+12](kBT) 573
(3.39)
j specify the particular energy as listed above. gj is the corresponding number of degrees of freedom, see Subsection 2.5.2. It takes the value unitY for the reservoir energies. The simulational m'th moment is calculated from expansion of
(3.40)
• Finally we analyse the distribution of simulation box energies (internal energies).
In a canonical simulation the internal energy should be consistent with a Gaussian distribution.
(3.41)
Ergodicity in connection with the Nose-Hoover thermostat was diseussed in Section 2.5. To have reliable averages, the simulation of both the system and its sorroundings must sample a representative part of phase space. But even if the mean values are sta-ble and give exellent values, the fluctuations around the mean can be highly non-canonica1. The outcome of the two final tests above signals the canonical behaviour of the system.