AIMD

TeraChem can carry out ab initio molecular dynamics calculations of various types. Any method which provides the analytic gradient can be used to provide potential energy surfaces and forces. Born-Oppenheimer molecular dynamics is used, where the forces are determined from a self-consistent calculation, i.e. the wavefunction coefficients are fully determined before computing the interatomic forces.

Initial conditions

A trajectory needs starting positions and velocities. TeraChem offers three ways to supply them:

Random velocities at a chosen temperature (the default). Set the temperature with tinit; if tinit is unset, t0 is used instead. Initial velocities are sampled from a Maxwell-Boltzmann distribution and the net linear and angular momenta are removed.
Velocities from a file by setting velocities <filename>. The file uses the same XYZ-style format that the trajectory output uses. This is the appropriate choice for restarting from a previous trajectory or for using precomputed conditions.
Wigner or Husimi sampling around a stationary point through run initcond, which diagonalizes the Hessian at the input geometry and samples positions and velocities from the harmonic distribution. See Initial Condition Generation for details, including the transition-state-sampling mode.

Choosing an ensemble

By default TeraChem runs in the microcanonical (NVE) ensemble — Newton's equations are solved without any feedback on the kinetic energy. NVE is the appropriate choice when computing dynamical properties: in the words of the legacy manual, "some would argue this is the only choice." For sampling, equilibration, or free-energy calculations, NVT can be more convenient. The ensemble is selected with thermostat:

NVE (no thermostat keyword)

If thermostat is omitted, Hamilton's equations are integrated by velocity Verlet and the total energy should be conserved up to integrator error. Time steps of 0.5 fs are typical for molecules built from light atoms; heavier or constrained systems can tolerate longer steps.

`thermostat rescale` — velocity rescaling

Hamilton's equations are integrated as in NVE, but every rescalefreq steps the velocities are rescaled to match the target temperature:

\[\alpha = \sqrt{\frac{K_{\text{desired}}}{K_{\text{actual}}}} \;, \qquad v_i^{\text{rescaled}} = \alpha\, v_i\]

with $K_{\text{desired}} = \tfrac{1}{2} N_{\text{DOF}} k_B T$ and $N_{\text{DOF}} = 3 N_{\text{atom}} - N_{\text{constraints}}$. Velocity rescaling does not exactly reproduce the canonical distribution and has been shown to introduce pathological behavior over long simulations,¹ so it is not recommended for production NVT work. Setting rescalefreq to a value much larger than the trajectory length is one way to run an effectively-NVE simulation.

`thermostat nhc` — Nosé-Hoover chain

Deterministic NVT integrator. The chain oscillator period is controlled by nhctime (fs).

`thermostat langevin` — Langevin dynamics

Stochastic NVT integrator with explicit friction and noise terms. The damping time is set by lnvtime (fs).

`thermostat bp` — Bussi-Parrinello stochastic velocity rescaling

The recommended thermostat for general NVT simulations. Like Langevin in that it samples the canonical distribution, but with weaker perturbation of the underlying dynamics. The relaxation time is set by lnvtime (fs).

The t0 keyword sets the target temperature in K, and can be a comma-separated list paired with nstep to define a temperature ramp (see the parameter table below).

Spherical boundary conditions

To prevent "evaporation" from a finite cluster or to confine the system to a chosen density, TeraChem can apply a spherical confining potential during MD. The potential is the sum of two harmonic terms,

\[U_{\text{constr}}(r) = k_1 \left( (r - R_{\text{center}}) - R_1 \right)^2 + k_2 \left( (r - R_{\text{center}}) - R_2 \right)^2\]

where $R_{\text{center}}$ is the center of mass of the system on the first MD frame and then held fixed. By default $k_1 = 10.0$ kcal/(mol Å²) and $k_2 = 0$.

The simplest way to apply spherical BCs is:

mdbc        spherical
md_density  1.0

Here md_density (g/mL) is used to compute $R_1$ automatically from the number of atoms in the system. To set the radius explicitly use md_r1 instead of md_density. Confining forces are not applied to hydrogen atoms by default; override this with mdbc_hydrogen yes.

The full list of mdbc_* keywords is in the Parameters for AIMD jobs table below. For time-dependent boundary conditions used to accelerate reactive sampling, see Nanoreactor.

Files produced by an AIMD calculation

Molden files - There are a number of these, with different orbitals
Singlet.x.molden: The natural orbitals corresponding to the xth singlet state
state_averaged_natural_orbital.molden:

Parameters for AIMD jobs

A minimal AIMD input requires run md plus a target temperature (t0), a time step (timestep), and a number of steps (nstep). By default initial velocities are randomly sampled at tinit (which defaults to t0), so a typical NVE or NVT run can be set up with just a handful of keywords.

Integrator and step control

Keyword	Type	Default	Description
`timestep`	float	1.0 fs	Velocity-Verlet integration time step
`nstep`	int	1	Number of MD steps. Can also be given as a comma-separated list paired with `t0` to define a temperature schedule (see below).
`t0`	float	300.0	Target temperature in K. Can be a comma-separated list of temperatures to ramp through; the corresponding `nstep` list must have the same length.
`tinit`	float	value of `t0`	Temperature used to sample initial random velocities
`velocities`	string	`random`	Either `random` (sample from a Maxwell-Boltzmann distribution at `tinit`) or the name of a file containing explicit velocities

Thermostat

Keyword	Type	Default	Description
`thermostat`	string	none (NVE)	Thermostat: `rescale`, `nhc` (Nose-Hoover chain), `langevin`, or `bp` (Bussi-Parrinello Langevin). Omit for NVE dynamics.
`rescalefreq`	int	as set internally	For `thermostat rescale`: rescale velocities every N steps
`nhctime`	float	as set internally	For `thermostat nhc`: oscillator period of the NHC, in fs
`lnvtime`	float	as set internally	For `thermostat langevin` and `thermostat bp`: damping time, in fs

Restarts and output

Keyword	Type	Default	Description
`restartmd`	filename	not set	Restart MD from this checkpoint file
`restartmdfreq`	int	100	Write a restart checkpoint every N steps
`md_binaryoutput`	bool	no	Write the trajectory in binary format instead of text
`md_outputfreq`	int	1	Write trajectory output every N steps

Boundary conditions and external forces

Keyword	Type	Default	Description
`mdbc`	string	none	Boundary condition: `spherical` (the most common; see above), `surface`, `disk`, or `hemisphere`
`md_r1`, `md_r2`	float	—	Radii (in Å) of the inner and outer confining spheres
`md_k1`, `md_k2`	float	$k_1$ = 10.0, $k_2$ = 0	Force constants of the inner and outer confining potentials in kcal/(mol Å²)
`md_density`	float	—	Density (g/mL) used to compute `md_r1` automatically; overridden by an explicit `md_r1`
`mdbc_hydrogen`	bool	no	Apply the confining potential to hydrogen atoms
`mdbc_mass_scaled`	bool	no	Scale the confining potential by atomic mass (used for nanoreactor compression to prevent hydrogen stripping)
`mdbc_t1`, `mdbc_t2`	int	—	Time-dependent BC schedule: alternate between the first and second sets of `md_r/md_k` parameters for these many steps each (used for nanoreactor runs)
`mdbc_no_centering`	bool	no	Disable re-centering of the confining sphere on the system's initial center of mass
`md_surface`	filename	—	XYZ file defining the atoms that make up the surface for `mdbc surface`
`md_sel1`, `md_sel2`	string	—	Atom selections that the boundary potential acts on
`md_diskheight`	float	—	Height (Å) of the confining disk for `mdbc disk`
`external_force`	filename	not set	File defining external forces applied during the trajectory

Fixing atoms during MD

Selected atoms can be held fixed during a TeraChem MD trajectory by listing their indices in a file called fixed_atoms placed in the same directory as the input file. The format is straightforward: the first line gives the number of fixed atoms, and each subsequent line gives the index of one fixed atom (one index per line). Atom indices are 0-based (the first atom in the coordinate file is atom 0).

fixed_atoms

The force on a fixed atom is zeroed out at every step with no modification to the forces on the other atoms. This is equivalent to setting the masses of the fixed atoms to infinity. Fixing atoms only takes effect during MD — it has no effect on geometry optimizations.

In QM/MM calculations, the fixed_atoms file applies to QM atoms only. To fix MM atoms, place their indices in a parallel file called fixed_mmatoms using the same format.

Steered molecular dynamics (AISMD)

Constant external forces can be applied to selected atoms to produce force-modified potential energy surfaces, either in the lab frame (for dynamics) or in the molecular frame (for optimizations under load). See Steered MD (AISMD) for the full description.

Nanoreactor

Reactive sampling can be dramatically accelerated by alternating between two sets of spherical boundary conditions on a fixed schedule, implementing a "virtual piston" that periodically compresses the system. This is the ab initio nanoreactor; see Nanoreactor for the full description and an example input.

Overriding atomic masses

The default atomic masses in TeraChem are those of the most abundant naturally occurring isotope. Frequencies and MD trajectories both depend on these masses, and either can be modified by including a $masses block in the input file:

$masses
C 13.0
15 H 2.0
$end

Two line formats are accepted:

<element> <mass> — sets all atoms of the named element to the given mass in amu.
<atom_index> <element> <mass> — sets the mass of one specific atom (the index is 1-based) and verifies that the element symbol matches the input geometry. If the element does not match, TeraChem aborts; this guards against off-by-one indexing errors.

Multiple lines are processed in order, so a later C 12.0 will overwrite an earlier C 13.0. The $masses block must end with $end.

S. C. Harvey, R. K.-Z. Tan, and T. E. Cheatham III, J. Comp. Chem. 19, 726 (1998). ↩