Unfortunately, we are faced with this choice. In fact, if we are interested in dynamic properties the Monte Carlo we will not suit 1). If we are not interested in dynamic properties, then we face a very uncertain choice.
In principle, the dynamics is a generic method, but there are several problems. First, the choice of the thermostat and barostat. Depending on it we can lose thermodynamic or dynamic properties, or both. Secondly, the choice of integration step. The step of dynamic algorithms is limited (a typical limit of 5 fs). Near this border there are large fluctuations of energy, it may be a systematic drift, distortion of other properties.
The Monte Carlo method guarantees the canonical distribution and consequently accurate thermodynamic properties. The properties are not distorted when we make long steps. However, with increasing step decreases the probability of its acceptance. As a result, the effectiveness of the methods is approximately the same.
As a result both methods give similar results in the quality and in the speed of calculations. If we want to be assured of quality, the more reliable is the Monte Carlo method. Molecular dynamics in some cases may be a little faster and it is indispensable in the calculation of dynamic properties.
Probably as the most versatile can be considered Nosé–Hoover thermostat. It is suitable for calculations of both thermodynamic and dynamic properties. However, it is sometimes advantageous to use other thermostats.
Berendsen thermostat is very stable. It allows you to quickly reach the desired temperature. Thermostat keeps the moment, but not the inertia moment. Thanks to these properties it is convenient to use it in the initial stages of work, especially in the construction of models. However, it does not provide a canonical distribution, so it should not be used for precision work.
Andersen thermostat provides a canonical distribution. Therefore, it is useful in studying the thermodynamic properties. Its main drawback is that it introduces additional frictional force, distorting the dynamic properties. Another disadvantage is that it does not keep an immobility of the model mass center.
Lowe thermostat similar to the Andersen thermostat, but it keeps the momentum and the inertial moment of the model. It also has an effect on the dynamic properties, but to a lesser extent. Its particular drawback is that it acts uniformly only on the homogeneous model. If the model is inhomogeneous, or worse has a non-spherical shape, the action of the thermostat will be different in different directions.
A Hybrid Monte Carlo can be considered as an ideal “thermostat” if we forget about the dynamic properties.
Dependence of the water diffusion coefficient from the applied thermostat2):
| Thermostat | d, 10−5 * cm2/s |
|---|---|
| Experiment | 2.13-2.66 |
| Nose-Hoover, τ = 0.3 ps | 2.3 |
| Nose-Hoover, τ = 0.03 ps | 2.3 |
| Nose-Hoover, τ = 0.01 ps | 2.4 |
| Andrsen, τ = 0.1 ps | 2.3 |
| Andrsen, τ = 1.0 ps | 0.8 |
| Lowe, τ = 1.0 ps | 2.2 |
| Berendsen, τ = 0.4 ps | 2.1 |
There are many algorithms to generate molecular dynamics trajectories, but the most widely used Verlet algorithm3). It has nice features:
This combination of properties makes it out of competition among the one-step algorithms.
Limiting the step size is determined by the most frequent vibrations in the existing model. In order to the algorithm working it is necessary to calculate not less than 4-5 points during the period of most rapid vibrations, the better 10 or more. The most rapid vibrations in molecular dynamics are usually vibrations of valence OH bond. Their period is about 10 fs. Accordingly the boundary of the stability of the Verlet algorithm is 2 fs. It’s not possible to make calculations on the border, the model will fall apart. Therefore, from a practical viewpoint the limit value is 1 fs. If we want to be sure of the accuracy of the calculation it is better to work with a step of 0.5 fs. This is a typical value used in the practice.
For highly accurate work the fourth-order algorithm VEFRL can be recommended. In step 0.5 fs, it will be 1000 times more accurate (energy conservation) than Verlet algorithm, and with a step 0.2 fs VEFRL will work with the utmost precision attainable on the processor. But we must bear in mind that VEFRL requires 4 calculation of the force on each step and thus runs 4 times slower. Because of this, for routine simulations it is unsuitable.
For routine simulations, we need not the exact algorithm, but one that allows us to work with the big steps. These are multi-step algorithms, such as RESPA. RESPA allows the routine simulation with 2 fs steps.
Step of any algorithm can be increased by 2-3 times if we attribute a large mass of the hydrogen atoms (4-10 a.e.). Because of this the dynamic properties of the model will be distorted, but the thermodynamic properties remain unchanged. Accordingly RESPA allows us to work with a 5 fs step.
For simple liquids equilibration time is 200-1000 ps. This time may be recommended for most simulations.
Folding time of the fastest proteins in implicit water is 2-10 ns, and in the explicit one is 4 μs.
For use dynamic as an optimization the 1-100 K temperature and ~100 ps time can be recommend.