====== Method choice ======

===== Dynamics or Hybrid Monte Carlo? =====

Unfortunately, we are faced with this choice.
In fact, if we are interested in dynamic properties the Monte Carlo we will not suit ((there are variants of Monte Carlo for dynamic properties too)). 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__.

===== Thermostat =====

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 thermostat((1 ns of molecular dynamics, with 1 fs steps of rRESPA four-level integrator, 216 molecules of the [[Water properties|flexible SPC]] water in periodic boundary conditions at the 1.0042 g/cm<sup>3</sup> density)):
^Thermostat^d, 10<sup>−5</sup> * cm<sup>2</sup>/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|





===== Step size in molecular dynamics =====

There are many algorithms to generate molecular dynamics trajectories, but the most widely used Verlet algorithm((he's a //leap frog// algorithm)).
It has nice features:
  - reversible
  - has no energy drift 
  - stable for large steps
  - requires only one calculation of the forces on the step
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.


===== Simulation time =====
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.