Molecular Dynamics

Molecular dynamics in Octopus is performed as part of a time-dependent calculation. In this mode the electronic system is propagated in time, and the ionic positions can also be updated. This page summarizes the molecular-dynamics functionality available in Octopus, and gives the main related variables.

The most common workflow is to first compute a ground state and then run a time-dependent calculation with moving ions. The ions are treated as classical particles with positions $\mathbf R_I(t)$, masses $M_I$, velocities $\mathbf v_I(t)$, and forces $\mathbf F_I(t)$. In its simplest form the ionic equation of motion is

$$ M_I \frac{d^2 \mathbf R_I}{dt^2} = \mathbf F_I(t), $$

where the force is computed from the electronic state and the ion-ion interaction. In a Born-Oppenheimer calculation the force is evaluated from the instantaneous electronic ground state, while in Ehrenfest dynamics it is evaluated from the time-dependent electronic state.

For a worked example of a Born-Oppenheimer molecular-dynamics run and the corresponding trajectory analysis, see the tutorial Born-Oppenheimer Molecular Dynamics.

Overview of supported molecular-dynamics functionality

The central switch for molecular dynamics is TDDynamics, used together with MoveIons. The following table summarizes the main types of ionic motion supported by Octopus.

Starting a molecular-dynamics run

A molecular-dynamics run is a time-dependent run, so the calculation mode must be set to td. The number of time steps and the time step itself are controlled in the usual way for time-dependent calculations. The variable TDTimeStep sets the propagation time step. It has to be chosen explicitly, and should be small enough that the time evolution remains stable and that the frequencies of interest can be resolved. The length of the propagation can be specified either as a number of steps (TDMaxSteps) or as a total propagation time (TDPropagationTime). The variables TDMaxSteps and TDPropagationTime are alternatives and should not be used at the same time.

To allow the ions to move, set

MoveIons = yes

The default of MoveIons depends on whether ionic velocities have been specified explicitly or implicitly: if initial velocities are provided, Octopus assumes that the ions should move; otherwise, the default is not to move the ions. In documentation and production inputs, it is usually clearer to set MoveIons explicitly.

Choice of dynamics

The variable TDDynamics selects the type of dynamics during the time propagation.

Ehrenfest dynamics

Ehrenfest dynamics (TDDynamics = ehrenfest) is the default value of TDDynamics. In this mode, the electronic orbitals are propagated in real time while the ions are propagated classically under the forces associated with the time-dependent electronic state. This mode is appropriate when the dynamics of the electrons and ions should be followed together, for example in real-time simulations where the electrons are driven away from the ground state.

When using Ehrenfest dynamics, the electronic time step is normally the limiting numerical parameter. A useful check is to propagate without an external perturbation and verify that the total energy remains stable over the intended propagation time.

Born-Oppenheimer molecular dynamics

In Born-Oppenheimer molecular dynamics (TDDynamics = bo), the ions are propagated while the electronic system is kept in its ground state at each ionic step. At each molecular-dynamics step Octopus performs a self-consistent calculation, then computes the forces used to propagate the ions. Schematically, the force on ion $I$ is $$ \mathbf F_I^{\mathrm{BO}}({\mathbf R}) = - \frac{\partial E_{\mathrm{KS}}^{\mathrm{gs}}({\mathbf R})}{\partial \mathbf R_I}. $$

The ionic positions are propagated with the velocity-Verlet algorithm. Written for a time step $\Delta t$, the update is

$$ \mathbf R_I(t+\Delta t) = \mathbf R_I(t) + \mathbf v_I(t)\Delta t + \frac{\mathbf F_I(t)}{2M_I}\Delta t^2, $$

$$ \mathbf v_I(t+\Delta t) = \mathbf v_I(t) + \frac{\mathbf F_I(t)+\mathbf F_I(t+\Delta t)}{2M_I}\Delta t. $$

Born-Oppenheimer molecular dynamics is therefore usually more expensive per ionic step than an Ehrenfest propagation with the same time step, but it is the natural choice for adiabatic ionic dynamics, thermal motion on the electronic ground-state surface, and vibrational analysis from a trajectory.

The time step should be chosen according to the ionic motion one wants to resolve, but the electronic ground-state calculation at each step must also be sufficiently converged for the forces to be reliable. For a practical example, see the Born-Oppenheimer molecular-dynamics tutorial linked above. The Born-Oppenheimer tutorial gives a complete worked example for a water molecule.

Accelerated ionic time scale

For advanced real-time calculations, TDIonicTimeScale can be used to change the relative time scale of ionic and electronic motion.

This variable plays a role analogous to the fictitious-mass parameter in Car-Parrinello-type ideas. Increasing it linearly accelerates the ionic time step relative to the electronic time step, but also increases the deviation from the Born-Oppenheimer surface. A value different from 1 means that the electrons no longer follow the physical electronic time scale.

This option should therefore be used as an acceleration or algorithmic device, not as a physical rescaling of time.

Initial velocities

Initial velocities determine whether the ions start from rest, from user-defined velocities, from a file, from a random thermal distribution, or from a phonon-mode displacement. If no velocities are provided, Octopus initializes the velocities to zero.

The initial velocities determine the initial ionic kinetic energy, $$ K_{\mathrm{ion}} = \frac{1}{2}\sum_I M_I |\mathbf v_I|^2, $$

and, for an unconstrained system with $N_{\mathrm{dof}}$ degrees of freedom, the corresponding instantaneous kinetic temperature is

$$ T_{\mathrm{ion}} = \frac{2K_{\mathrm{ion}}}{N_{\mathrm{dof}} k_{\mathrm B}}. $$

Octopus can initialize velocities in several ways.

Explicit velocity block

The most direct way to set initial ionic velocities is the Velocities block, for instance:

%Velocities
  "C" | -1.7 | 0.0 | 0.0
  "O" |  1.7 | 0.0 | 0.0
%

Each line gives the species label and the three Cartesian velocity components for one atom. The order of the entries must match the order of the atoms in the coordinate specification. If the order is inconsistent, the velocities will be assigned to the wrong atoms.

This can for instance be use to initiate an ion dynamics for a given vibration mode/photon mode.

Reading velocities from files

Velocities can also be read from files, using a format consistent with the corresponding coordinate-file variables, see XYZVelocities, XSFVelocities, and PDBVelocities.

As with the Velocities block, the ordering of atoms in the velocity file must match the ordering of atoms in the coordinate specification.

Random thermal velocities

The variable RandomVelocityTemp assigns random velocities at the initial time to the atoms, following a Boltzmann distribution at the temperature given by the variable, in Kelvin. The default value is 0 K. The random-number seed can be modified with the GSL_RNG_SEED environment variable.

For a Maxwell-Boltzmann distribution, each Cartesian velocity component has a variance that scales as

$$ \langle v_{I\alpha}^2 \rangle = \frac{k_{\mathrm B}T}{M_I}, $$

so lighter atoms receive larger typical velocities than heavier atoms at the same temperature.

Initial displacements from phonon modes

Octopus can also initialize positions and velocities from a selected phonon mode. The variable InitialDisplacementMode selects a non-zero mode index. When it is set, Octopus reads the phonon modes from PhononModesFile and displaces the ionic positions and velocities according to that mode.

The amplitudes of the initial displacement are controlled by

• InitialDisplacementAmplitudePos for the position displacement;
• InitialDisplacementAmplitudeVel for the velocity displacement.

The file specified by PhononModesFile is a plain-text file containing the frequencies and normal-mode eigenvectors. The atom order in this file must be the same as in the input geometry. The variable PhononModesZeroThreshold controls which low-frequency modes are treated as zero-frequency modes, such as translations or rotations, and ignored when generating initial displacements.

Thermostats and temperature control

The ionic thermostat is selected with Thermostat. The available values are:

The target temperature is not entered directly in Thermostat. Instead, TemperatureFunction gives the name of a function in the TDFunctions block. Temperature values are given in Kelvin.

For a constant 300 K thermostat target, one can write

Thermostat = nose_hoover
TemperatureFunction = "temperature"
ThermostatMass = 1.0

%TDFunctions
  "temperature" | tdf_cw | 300.0
%

The TDFunctions block supports several time-dependent function shapes, including constant (tdf_cw), Gaussian (tdf_gaussian), cosinoidal (tdf_cosinoidal), trapezoidal (tdf_trapezoidal), file-based (tdf_from_file), and expression-based (tdf_from_expr) functions. This makes it possible to use either a constant target temperature or a temperature protocol, for example a ramp or externally tabulated function.

Velocity scaling

With Thermostat = velocity_scaling, the velocities are rescaled to control the ionic temperature. In the simplest picture, if $T_{\mathrm{target}}$ is the target temperature and $T_{\mathrm{ion}}$ is the instantaneous kinetic temperature, a scaling factor

$$ \lambda = \sqrt{\frac{T_{\mathrm{target}}}{T_{\mathrm{ion}}}} $$

changes the velocities according to

$$ \mathbf v_I \leftarrow \lambda \mathbf v_I. $$

This method is useful for simple temperature control and for preparing trajectories, but it does not represent a full deterministic heat bath.

Nose-Hoover thermostat

With Thermostat = nose_hoover, Octopus uses a Nose-Hoover thermostat. The variable ThermostatMass sets the fictitious mass associated with the thermostat variable. Schematically, a Nose-Hoover thermostat augments the ionic equations of motion with a friction-like variable $\xi$,

$$ \dot{\mathbf v}_I = \frac{\mathbf F_I}{M_I} - \xi \mathbf v_I, $$

$$ Q\dot{\xi} = 2K_{\mathrm{ion}} - N_{\mathrm{dof}} k_{\mathrm B} T_{\mathrm{target}}, $$

where $Q$ is the thermostat mass and $T_{\mathrm{target}}$ is the value provided through TemperatureFunction and TDFunctions. Larger thermostat masses produce slower temperature feedback, while smaller masses produce faster coupling to the target temperature.

Prescribed ionic motion

Octopus also supports cases where the ionic motion is prescribed rather than determined by the forces.

Constant-velocity ions

The variable IonsConstantVelocity is experimental. If it is set to yes, ions move with the constant velocity specified by the initial conditions and are not affected by forces:

$$ \mathbf R_I(t) = \mathbf R_I(0) + \mathbf v_I(0)t. $$

This is useful for controlled model trajectories or for testing the electronic response to moving ionic potentials.

Time-dependent displacements

The experimental block IonsTimeDependentDisplacements specifies time-dependent displacements from equilibrium positions,

$$ \mathbf r_I(t) = \mathbf r_{I,0} + \Delta\mathbf r_I(t). $$

The block gives one line per atom to be displaced:

%IonsTimeDependentDisplacements
  atom_index | "dx(t)" | "dy(t)" | "dz(t)"
%

The functions dx(t), dy(t), and dz(t) must match function names defined in TDFunctions. When this block is set, the ions are not affected by forces. A simple oscillatory displacement can be defined using tdf_from_expr:

%TDFunctions
  "zero" | tdf_cw | 0.0
  "zosc" | tdf_from_expr | "0.05*sin(0.01*t)"
%

%IonsTimeDependentDisplacements
  1 | "zero" | "zero" | "zosc"
%

Periodic systems: cell dynamics and pressure

For periodic systems, Octopus can also propagate the simulation cell. The variable CellDynamics controls whether cell relaxation is performed during a time-propagation run. This is based on the Parrinello-Rahman equation of motion for the cell and is only available for periodic systems.

If the lattice vectors are grouped into a cell matrix $\mathbf h(t)$, ionic positions can be written in terms of reduced coordinates $\mathbf s_I$ as

$$ \mathbf R_I(t) = \mathbf h(t)\mathbf s_I(t). $$

Cell dynamics introduces equations of motion for $\mathbf h(t)$, so that the cell can respond to internal stresses. Schematic Parrinello-Rahman dynamics has the form

$$ W\ddot{\mathbf h} \sim V\left(\boldsymbol{\sigma}_{\mathrm{int}} - \boldsymbol{\sigma}_{\mathrm{ext}}\right), $$

where $W$ is a fictitious cell mass, $V$ is the cell volume, and $\boldsymbol{\sigma}$ denotes stress. For hydrostatic pressure, the external stress is proportional to the identity matrix,

$$ \boldsymbol{\sigma}_{\mathrm{ext}} = -P_{\mathrm{ext}}\mathbf I. $$

The external hydrostatic pressure is set with HydrostaticPressure. This variable can be used for geometry optimization and molecular dynamics in periodic systems.

Ground-state recalculation during time evolution

Some time-dependent outputs compare the propagated states with ground-state Kohn-Sham states. If the ions move, the ionic potential changes during the propagation, and one may want the reference ground-state orbitals to follow the instantaneous ionic geometry.

This is controlled by the variables RecalculateGSDuringEvolution, and RecalculateGSInterval .

When RecalculateGSDuringEvolution is enabled, the ground state is not recalculated at every step. Instead, RecalculateGSInterval sets the number of time-propagation steps between two ground-state recalculations. A smaller interval follows the instantaneous Hamiltonian more closely, while a larger interval reduces the computational cost of outputs that need the instantaneous ground-state basis.

Self-consistency during electronic propagation

For time-dependent Kohn-Sham propagation, the Hamiltonian depends on the density and the propagation step is, in principle, nonlinear. The variables TDStepsWithSelfConsistency and TDSCFThreshold control whether and how self-consistency is imposed during electronic time propagation.

TDStepsWithSelfConsistency gives the number of propagation steps for which self-consistency is imposed. The value all_steps imposes self-consistency for all steps, while 0 means that no propagation step is forced to be self-consistent. The convergence threshold for this self-consistency is TDSCFThreshold.

For many molecular-dynamics calculations, especially Born-Oppenheimer molecular dynamics, the convergence of the ground-state calculation at each ionic configuration is often the dominant accuracy parameter. For Ehrenfest dynamics or other coupled electron-ion calculations, the electronic time-step and self-consistency settings should be checked by monitoring conserved quantities such as the total energy in a suitable test propagation.

Analysis of trajectories

The main trajectory-based analyses are done with utilities. To create an XYZ trajectory that can be visualized with external programs, use oct-xyz-anim. To compute vibrational information from a molecular-dynamics trajectory, use oct-vibrational_spectrum.

The vibrational spectrum utility is based on the velocity-autocorrelation function. A typical definition is

$$ C_{vv}(t) = \left\langle \sum_I \mathbf v_I(t_0)\cdot\mathbf v_I(t_0+t) \right\rangle_{t_0}, $$

and the spectrum is obtained from a Fourier transform,

$$ I(\omega) \propto \int_{-\infty}^{\infty} C_{vv}(t)e^{-i\omega t},dt. $$

The variable VibrationalSpectrumTimeStepFactor can be used by oct-vibrational_spectrum to select how many molecular-dynamics time steps are used in the analysis. The variable AnimationSampling controls how many frames are used when creating movies with oct-xyz-anim.

Caveats

Several molecular-dynamics features are marked experimental in the variable documentation, including Born-Oppenheimer dynamics, constant-velocity ions, and time-dependent ionic displacements. Users should therefore validate energy conservation, temperature control, and force convergence for their systems.

The time step should always be convergence-tested. For thermostatted runs, the target temperature function, thermostat choice, and ThermostatMass should be tested to ensure that the dynamics are appropriate for the intended physical interpretation.