3.12 Molecular dynamics

For molecular dynamics runs within the microcanonical ensemble (“$NVE$”) or for deterministic thermostats, the wave function can be extrapolated in order to reduce the computational effort to reach self-consistency. This section specifies what quantity is actually extrapolated and how.

Following the approach of Kühne et al. [181], we extrapolate the contra-covariant density matrix $\mathrm{𝐏𝐒}$, the product of the ordinary (purely contravariant) one-particle density matrix $𝐏$ and the overlap matrix $𝐒$. The contra-covariant density matrix is then used to project new wave function coefficients from the last iteration ${𝐜}_{\mathrm{new}}={(\mathrm{𝐏𝐒})}_{\mathrm{extra}}{𝐜}_{\mathrm{𝐨𝐥𝐝}}$. As a last step, the one-particle coefficients are orthonormalized.

The extrapolation scheme we use is sketched in Fig. 3.1 for the example of an ordinary function $f(t)$ in a real variable $t$. For a p$no$ scheme (chosen with wf_extrapolation polynomial $n$ $o$) $n$ old iterations $f(t_0-k\mathrm{\Delta }t)$, $k=1,\mathrm{\dots },n$ are stored. An ansatz function with $n$ degrees of freedom is then used to fit all of these previous iterations and evaluated at $t_0$ and the value used to intialize the SCF procedure.

The choice of the fitting function should aim at two targets: accurate extrapolation and time-reversal symmetry. The parameter $o$ specifies the order of the extrapolation. The higher $o$, the better the extrapolation gets with decreasing time stemps $\mathrm{\Delta }t$. Time-reversal symmetry is useful to avoid energetic drifts for not-so-tight SCF accuracy settings. By adding odd terms (odd with respect to $t_0+t\leftrightarrow t_0-t$), time-reversal symmetry is enhanced [175].

Therefore, in contrast to [255, 136], the $(n-o-1)$ remaining degrees of freedom are not resolved within a least-squares fit but instead to add fitting functions ${(t-t_0)}^{2k+1}$ to enhance time-reversal symmetry and thus energy conversion. Please note, however, that time-reversal symmetry itself only enhances energy conservation and not necessarily dynamical properties of the trajectories.

If you are interested in energy conservation, “p31” is in general a good choice to start with. If a good initial guess for the SCF procedure is desired, you should give “p32” a try. Please note that deactivation of extrapolation (wf_extrapolation none) is actually the same as the “p10” extrapolation (see Fig. 3.1) and simply uses the one-particle coefficients of the last iteration are to initialize the SCF procedure.

In geometry.in, the line containing the velocity must follow the line containing the atom that the velocity refers to. This can be used, e.g., to initialize a molecular dynamics run, or for analysis purposes later. Note that, for molecular dynamics, the FHI-aims standard output prints this information in the proper format, as part of a particular geometry associated with a molecular dynamics trajectory.

A negative number signals that the ending criterion is not checked, in fact, the default setting is N=-1.

This keyword is for a rough initialization only, and is overridden by any successful calls of the MD restarting procedure through MD_restart. The default initialization for all velocities is zero.

This option is useful for non-periodic systems, allowing to weed out some residual numerical noise in the forces. However, seeming rotations of the unit cell can easily appear in periodic systems for completely normal, physical reasons. Since this led to some confusion, this option is now actively disabled (code stops) for periodic systems. If you have a good reason to use this option in periodic systems, it can be re-introduced by hand by setting the keyword force_MD_clean_rotations to .true..

Possible values for option are .true., .false., time_restart, or filename. The default is .false..

The data for the molecular dynamics integrator is always written to a file aims_MD_restart.dat after each time step (regardless of whether or not keyword MD_restart is used). If MD_restart is set, this keyword controls the reading of such data from a previous run (if that exists), to either

• •

  continue a previous run (.true.), or

• •

  to use previous position / velocity information but reset all timings (time_restart), or

• •

  to begin a new run from scratch (.false.)

The default file name is used unless the user specifies a specific file filename from where the input is to be read.

Important: If the starting structure and velocities are taken from the restart file, any exact positions noted in geometry.in are ignored. However, the species identification, possible initial charges, and other per-atom information in geometry.in remain valid and must be present as always. If the option .true. is set, the molecular dynamics clock is also read from file, in all other restarts the clock is reset to zero.

Possible values for option are .true. and .false.. The default is .true.. Use this flag to switch to the ASCII format in the MD restart file.

This keyword is the key control for all molecular dynamics. The runtimes time are specified in ps. There are five different possible ensembles, each of which require different options - as described in a separate subsection below.

When a schedule of different temperatures, thermostats etc is required within the same run (e.g., initialization followed by NVE, change of temperature, etc.), use the MD_schedule keyword instead of MD_run.

Must be followed by an arbitrary number of lines MD_segment, where different temperatures, thermostats, etc. may be specified for each segment of the run.

Keyword MD_segment must only appear in consecutive lines after an MD_schedule keyword. Instead of a single set of MD thermostats, temperatures etc. throughout the simulation, this keyword allows to set specific values for only a segment of the full run.

The option polynomial 3 1 strongly reduces a possible energy drift in NVE runs even for moderate force accuracy settings and additionally lowers the number of SCF cycles per time step.

Warning: This feature is experimental. It most probably will not enhance convergence of metallic systems and is not parallelized. The niklasson06 schemes are not suitable for long runs of nontrivial physical systems because a regular reinitialization would be necessary.

The same warning as for wf_extrapolation applies.

When used with a molecular dynamics schedule of time segments with different thermostats, the following subkeywords should appear behind an MD_segment keyword instead of MD_run.

This option is a useful addition to the structural relaxation, which can be done in principle with a molecular dynamics run as well. In almost all cases, however, the BFGS algorithm as called with the relax_geometry keyword is preferable.

Andersen’s [9] simple definition of a thermostat that randomly resets the velocity of individual atoms to a Maxwell-Boltzmann distributed value with a given frequency. This is an example of a rather harsh thermostat.

Note that the Berendsen thermostat does NOT resemble any physical ensemble, but that it is a useful standard to initialize a molecular dynamics simulation. The relaxation time tau must be chosen by the user, there is no default. Remember that the special case $\tau =\mathrm{\Delta }t$ reproduces the correct temperature exactly at every time step and can be used for a very rough initialization.

This is the Bussi-Donadio-Parrinello thermostat, as described in Ref. [49] It a) properly samples the canonical distribution and b) preserves time correlations. The relaxation time tau (in ps) must be chosen by the user, there is no default. The performance of the thermostat is claimed in Ref. [49] to be practically independent of the value of $\tau$, for condensed phases. For clusters, though, a proper tuning of tau might be necessary. When this ensemble is chosen, the conserved pseudo-Hamiltonian (as described in Ref. [49]) is printed in the output.

This is a flexible thermostat based on the Generalized Langevin Equation, that adds extra degrees of freedom to the equations of motion and has a frequency dependent memory kernel, so that one can adjust the performance of the thermostat for different degrees of freedom in the system. It requires also matrices as inputs, for which we refer the reader to the keywords MD_gle_A and MD_gle_C. Please read references [55, 58, 56] and references therein before using this. It can be used to sample the canonical ensemble or to break detailed balance and simulate other conditions, including approximate quantum effects. It has a conserved quantity in the same spirit of the BDP thermostat.

Probably the most popular thermostat in the literature.

Due to numerical issues, this thermostat has been disabled for the time being.

Units should be 1/ps, and the matrix can be downloaded from http://epfl-cosmo.github.io/gle4md/. If you wish to model different dynamics (quantum effects, or thermalization of PIMD), one can also specify a C matrix with the keyword MD_gle_C.

This matrix is optional for the usage of the GLE thermostats. If sampling the canonical ensemble it is not needed. Otherwise, it is. Units should be K, and the matrix can be downloaded from http://epfl-cosmo.github.io/gle4md/.

The best way to perform path integral molecular dynamics and other more advanced dynamics techniques in FHI-aims is through the i-PI python wrapper [57]. This code is available free of charge. Information about the code can be found in http://ipi-code.org/ and we provide a quick tutorial on how one can make it work with FHI-aims in the tutorials. Through this interface, one can perform all types of dynamics (classical or quantum) with a wide range of thermostats, barostats, and different path integral acceleration techniques. i-PI uses internet sockets for the communication with the client codes, making it also easy to join different codes in the same simulation (e.g. a thermodynamic integration), as long as they can all communicate with i-PI. Note that NPT and NST (constant stress) simulations are currently supported only for functionals where the analytical stress tensor is available. If you wish to perform an NPT simulation, please add also compute_analytical_stress .true. to your control.in file.

The basic keywords that can appear in the control.in file of FHI-aims are listed below. Examples of how to run FHI-aims bound to i-PI are available in the folder aimsfiles/examples/ipi with the corresponding i-PI and FHI-aims input files, as well as a short explanation on how to run the programs. Examples of FHI-aims being used with i-PI can also be found in the i-PI distribution. Please refer also to the i-PI manual and contact the developer responsible for the FHI-aims interface (Mariana Rossi) if you have any doubts about the usage of FHI-aims with this code.

Using TCP/IP sockets (Internet sockets) is the most straightforward and consistent way of using the i-PI wrapper with FHI-aims, particularly on HPC systems where UNIX sockets are unavailable. However, IP/TCP ports are assigned dynamically on a system, and therefore there can be no guarantee that a previously assigned port in an input file will be free when the calculation starts. To avoid this problem we recommend using the utility script get_free_port.py in the utilities directory to set the port to right before initializing the calculation either locally or inside the submission script on an HPC system. This script will automatically check if a requested port and change it if it is not free or find a free port, and update both the relevant information in the i-PI and FHI-aims input files. To use the script simply, run

  > get_free_port.py -x inputs.template.xml \
                     -ox inputs.xml \
                     -c control.in \
                     -oc {system_name_descriptor}/control.in

or for cases when multiple calculators are used

  > get_free_port.py -x inputs.template.xml \
                     -ox inputs.xml \
                     -c control.in \
                     -oc {system_1_name_descriptor}/control.in \
                         {system_2_name_descriptor}/control.in

In cases where you simply want to overwrite the template input files do not set the -ox or -oc variables. Additionally if you are running a system with a variable IP address, you can update the hostaddress with the following command:

  > get_free_port.py --host HOST_ADDRESS \
                     -x inputs.template.xml \
                     -ox inputs.xml \
                     -c control.in \
                     -oc {system_name_descriptor}/control.in

where ----host default will set the host to be the default for the system you are on. For a complete description of the functionality of the script run

  > get_free_port.py --help