Molecular Dynamics

Molecular dynamics is a good method to perform thermodynamical averages, and to look at dynamical properties. It is also a good starting point for many other more advanced techniques.

During an MD one expects the density to change more or less continuously, thus it is possible to use the solutions of the previous steps to create the new initial guess for the density (and wavefunctions for OT). Indeed cp2k can use different EXTRAPOLATION techniques.

For MD a good extrapolation is PS with an EXTRAPOLATION_ORDER of 3. If you are doing something where from one step to the other there is little continuity, (geometry optimization, montecarlo), then something with little continuity like LINEAR_PS or just the previous density (USE_PREV_P) might be better. If you change the particle number then pick USE_GUESS that restarts from scratch with the restart method choosen in SCF_GUESS. The ASPC extrapolation is used in connection with langevin and a reduced SCF convergence (just on SCF step) and history restart to simulate big systems.

Trajectory

With the printkey TRAJECTORY you can control the output of the trajectory. The name of the file is by default projectname-pos-1.xyz and projectname is whatever you have written in PROJECT_NAME. In the prinkey you can also change the format from xyz to something else.

An interesting file to check is the projectname-1.ener file, in it you can find the following columns:

md_step, time[fs], e_kin [hartree], temp[K], e_pot [hartree], e_tot [hartree], cpu_time_per_step [s]

You should always check it and look at how the system equilibrates.

The .ener file and other interesting trajectory files are all controlled with from the PRINT section.

MD Convergence

If the MD has to be trusted then one should be sure that the trajectory can be trusted. Actually a simple, and very sensitive test that there are no big technical errors is to perform an NVE trajectory and look at the energy conservation. The energy conservation has normally two features, short time oscillations (that are larger when the system is not yet equilibrated) and a long range drift. If you express these in K, then you can compare them with the temperature that you are simulating at. Another (equivalent) possibility is to express them as fraction of the kinetic energy. The oscillation and drift (normally per ps, but it also depends on how many picoseconds you want to simulate, and if you want an NVE trajectory) should be small with respect to the kinetic energy (1% or less is a good value, but obviously it depends on the accuracy that you want to achieve, more might be acceptable, or less needed).

To improve the energy conservation one can either improve the forces with EPS_SCF and EPS_DEFAULT, improve \(\tilde n\) increasing the cutoff, and reduce the timestep. For the timestep a good value is normally around 1 fs, but if you are simulating high temperature or special system you might need to reduce it. Normally having good forces and larger timestep is advantageous. Obviously the timestep cannot be larger than the period of the highest oscillation of the system (typically H atoms, that is the reason why sometimes D atoms are used).

Conserving the total energy well is not enough if the system is very heterogeneous and the interesting part is small. In that case there is the risk that even a large error on it might pass unnoticed. If this is to be excepted (for example the atom with the sharpest gaussian basis set is in the interesting part) then checking that part of the system (oscillations, kinetic energy,…) is advisable.

Equilibration

For the result to be meaningful the system should be equilibrated, and not in a very unlikely state. If one is doing shooting or transition path sampling then this is less true, but still the system should not be in a very high-energy state.

So at the beginning you should build your system in some meaningful way, or from classical simulations. To have a better starting point you can optimize the energy (if your system is small), you anneal it, but it is not always needed. Then you can equilibrate it at a given temperature.

To equilibrate a system one can use a crude velocity rescale when it is too far away from the goal temperature (as established by TEMP_TOL). Doing an NVE simulation with TEMP_TOL is better way to equilibrate than using the NVT ensamble that uses the Nose-Hoover chain thermostats and might give spurios effects if you are far from equilibrium, as it tryes to conserve some extended quantity.the thermostat can store energy that he will give back at some later point. It should be taken into account that the temperature is not constant, but does oscillate in the NVE ensamble, these oscillations are connected with the specific heat and are inversely proportional with sqrt(N) where N is the system size.

A faster way to equilibrate the system is to use LANGEVIN as ENSEMBLE, and use a GAMMA of 0.001 [1/fs] (or larger if you want to equilibrate faster). Langevin introduces a viscous dissipative force and a random forces that are in equilibrioum at the given temperature. For equilibration purposes (and not to simulate removed degrees of freedom, like a solvent), the smaller gamma the longer it takes to equilibrate, and the closer the trajectory is to an NVE trajectory, the larger gamma the more the “environment thermostat” influences the system. Also With langevin if your system is in a bad initial configuration it is a good idea to set a TEMP_TOL, or first anneal a little the system.