9.8 Ab Initio Molecular Dynamics

9.8.1 Overview and Basic Job Control

Initial Cartesian coordinates and velocities must be specified for the nuclei. Coordinates are specified in the $molecule section as usual, while velocities can be specified using a $velocity section with the form:

$velocity
vx,1vy,1vz,1
vx,2vy,2vz,2
vx,Nvy,Nvz,N
$end

Here vx,i, vy,i, and vz,i are the x, y, and z Cartesian velocities of the ith nucleus, specified in atomic units (bohrs per atomic unit of time, where 1 a.u. of time is approximately 0.0242 fs). The $velocity section thus has the same form as the $molecule section, but without atomic symbols and without the line specifying charge and multiplicity. The atoms must be ordered in the same manner in both the $velocity and $molecule sections.

As an alternative to a $velocity section, initial nuclear velocities can be sampled from certain distributions (e.g., Maxwell-Boltzmann), using the AIMD_INIT_VELOC variable described below. AIMD_INIT_VELOC can also be set to QUASICLASSICAL, which triggers the use of quasi-classical trajectory molecular dynamics (see Section 9.8.5).

The nuclear mass can be initialized by a $mass section with the form:

$mass
m1m2
m3 ! mass of 3rd atom
m4
$end

The total number in the $mass section must be equal to number of atoms. Unit is the atomic unit, e.g., mass of hydrogen-1 atom is 1.00783. If the $mass section is not initialized, the default mass will be used.

Although the Q-Chem output file dutifully records the progress of any ab initio molecular dynamics job, the most useful information is printed not to the main output file but rather to a directory called “AIMD” that is a subdirectory of the job’s scratch directory. (All ab initio molecular dynamics jobs should therefore use the –save option described in Section 2.7.) The AIMD directory consists of a set of files that record, in ASCII format, one line of information at each time step. Each file contains a few comment lines (indicated by “#”) that describe its contents and which we summarize in the list below.

  • Cost: Records the number of SCF cycles, the total CPU time, and the total memory use at each dynamics step.

  • EComponents: Records various components of the total energy (all in hartree).

  • Energy: Records the total energy and fluctuations therein.

  • MulMoments: If multipole moments are requested, they are printed here.

  • NucCarts: Records the nuclear Cartesian coordinates x1, y1, z1, x2, y2, z2, …, xN, yN, zN at each time step, in either bohrs or Ångstroms.

  • NucForces: Records the Cartesian forces on the nuclei at each time step (same order as the coordinates, but given in atomic units).

  • NucVeloc: Records the Cartesian velocities of the nuclei at each time step (same order as the coordinates, but given in atomic units).

  • TandV: Records the kinetic and potential energy, as well as fluctuations in each.

  • View.xyz: Cartesian-formatted version of NucCarts for viewing trajectories in an external visualization program.

For ELMD jobs, there are other output files as well:

  • ChangeInF: Records the matrix norm and largest magnitude element of Δ𝐅=𝐅(t+δt)-𝐅(t) in the basis of Cholesky-orthogonalized AOs. The files ChangeInP, ChangeInL, and ChangeInZ provide analogous information for the density matrix 𝐏 and the Cholesky orthogonalization matrices 𝐋 and 𝐙 defined in Ref. 373.

  • DeltaNorm: Records the norm and largest magnitude element of the curvy-steps rotation angle matrix Δ defined in Ref. 373. Matrix elements of Δ are the dynamical variables representing the electronic degrees of freedom. The output file DeltaDotNorm provides the same information for the electronic velocity matrix dΔ/dt.

  • ElecGradNorm: Records the norm and largest magnitude element of the electronic gradient matrix 𝐅𝐏-𝐏𝐅 in the Cholesky basis.

  • dTfict: Records the instantaneous time derivative of the fictitious kinetic energy at each time step, in atomic units.

Ab initio molecular dynamics jobs are requested by specifying JOBTYPE = AIMD. Initial velocities must be specified either using a $velocity section or via the AIMD_INIT_VELOC keyword described below. In addition, the following $rem variables must be specified for any ab initio molecular dynamics job:

AIMD_METHOD
       Selects an ab initio molecular dynamics algorithm.
TYPE:
       STRING
DEFAULT:
       BOMD
OPTIONS:
       BOMD Born-Oppenheimer molecular dynamics. CURVY Curvy-steps Extended Lagrangian molecular dynamics.
RECOMMENDATION:
       BOMD yields exact classical molecular dynamics, provided that the energy is tolerably conserved. ELMD is an approximation to exact classical dynamics whose validity should be tested for the properties of interest.

TIME_STEP
       Specifies the molecular dynamics time step, in atomic units (1 a.u.  = 0.0242 fs).
TYPE:
       INTEGER
DEFAULT:
       None.
OPTIONS:
       User-specified.
RECOMMENDATION:
       Smaller time steps lead to better energy conservation; too large a time step may cause the job to fail entirely. Make the time step as large as possible, consistent with tolerable energy conservation.

AIMD_TIME_STEP_CONVERSION
       Modifies the molecular dynamics time step to increase granularity.
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       n The molecular dynamics time step is TIME_STEP/n a.u.
RECOMMENDATION:
       None

AIMD_STEPS
       Specifies the requested number of molecular dynamics steps.
TYPE:
       INTEGER
DEFAULT:
       None.
OPTIONS:
       User-specified.
RECOMMENDATION:
       None.

Ab initio molecular dynamics calculations can be quite expensive, and thus Q-Chem includes several algorithms designed to accelerate such calculations. At the self-consistent field (Hartree-Fock and DFT) level, BOMD calculations can be greatly accelerated by using information from previous time steps to construct a good initial guess for the new molecular orbitals or Fock matrix, thus hastening SCF convergence. A Fock matrix extrapolation procedure,374 based on a suggestion by Pulay and Fogarasi,794 is available for this purpose.

The Fock matrix elements Fμν in the atomic orbital basis are oscillatory functions of the time t, and Q-Chem’s extrapolation procedure fits these oscillations to a power series in t:

Fμν(t)=n=0Ncntn (9.16)

The N+1 extrapolation coefficients cn are determined by a fit to a set of M Fock matrices retained from previous time steps. Fock matrix extrapolation can significantly reduce the number of SCF iterations required at each time step, but for low-order extrapolations, or if SCF_CONVERGENCE is set too small, a systematic drift in the total energy may be observed. Benchmark calculations testing the limits of energy conservation can be found in Ref. 374, and demonstrate that numerically exact classical dynamics (without energy drift) can be obtained at significantly reduced cost.

Fock matrix extrapolation is requested by specifying values for N and M, as in the form of the following two $rem variables:

FOCK_EXTRAP_ORDER
       Specifies the polynomial order N for Fock matrix extrapolation.
TYPE:
       INTEGER
DEFAULT:
       0 Do not perform Fock matrix extrapolation.
OPTIONS:
       N Extrapolate using an Nth-order polynomial (N>0).
RECOMMENDATION:
       None

FOCK_EXTRAP_POINTS
       Specifies the number M of old Fock matrices that are retained for use in extrapolation.
TYPE:
       INTEGER
DEFAULT:
       0 Do not perform Fock matrix extrapolation.
OPTIONS:
       M Save M Fock matrices for use in extrapolation (M>N)
RECOMMENDATION:
       Higher-order extrapolations with more saved Fock matrices are faster and conserve energy better than low-order extrapolations, up to a point. In many cases, the scheme (N = 6, M = 12), in conjunction with SCF_CONVERGENCE = 6, is found to provide about a 50% savings in computational cost while still conserving energy.

When nuclear forces are computed using underlying electronic structure methods with non-optimized orbitals (such as MP2), a set of response equations must be solved.14 While these equations are linear, their dimensionality necessitates an iterative solution,777, 493 which, in practice, looks much like the SCF equations. Extrapolation is again useful here,914 and the syntax for Z-vector (response) extrapolation is similar to Fock extrapolation.

Z_EXTRAP_ORDER
       Specifies the polynomial order N for Z-vector extrapolation.
TYPE:
       INTEGER
DEFAULT:
       0 Do not perform Z-vector extrapolation.
OPTIONS:
       N Extrapolate using an Nth-order polynomial (N>0).
RECOMMENDATION:
       None

Z_EXTRAP_POINTS
       Specifies the number M of old Z-vectors that are retained for use in extrapolation.
TYPE:
       INTEGER
DEFAULT:
       0 Do not perform response equation extrapolation.
OPTIONS:
       M Save M previous Z-vectors for use in extrapolation (M>N)
RECOMMENDATION:
       Using the default Z-vector convergence settings, a (M,N)=(4,2) extrapolation was shown to provide the greatest speedup. At this setting, a 2–3-fold reduction in iterations was demonstrated.

Assuming decent conservation, a BOMD calculation represents exact classical dynamics on the Born-Oppenheimer potential energy surface. In contrast, so-called extended Lagrangian molecular dynamics (ELMD) methods make an approximation to exact classical dynamics in order to expedite the calculations. ELMD methods—of which the most famous is Car–Parrinello molecular dynamics—introduce a fictitious dynamics for the electronic (orbital) degrees of freedom, which are then propagated alongside the nuclear degrees of freedom, rather than optimized at each time step as they are in a BOMD calculation. The fictitious electronic dynamics is controlled by a fictitious mass parameter μ, and the value of μ controls both the accuracy and the efficiency of the method. In the limit of small μ the nuclei and the orbitals propagate adiabatically, and ELMD mimics true classical dynamics. Larger values of μ slow down the electronic dynamics, allowing for larger time steps (and more computationally efficient dynamics), at the expense of an ever-greater approximation to true classical dynamics.

Q-Chem’s ELMD algorithm is based upon propagating the density matrix, expressed in a basis of atom-centered Gaussian orbitals, along shortest-distance paths (geodesics) of the manifold of allowed density matrices 𝐏. Idempotency of 𝐏 is maintained at every time step, by construction, and thus our algorithm requires neither density matrix purification, nor iterative solution for Lagrange multipliers (to enforce orthogonality of the molecular orbitals). We call this procedure “curvy steps” ELMD,373 and in a sense it is a time-dependent implementation of the GDM algorithm (Section 4.5) for converging SCF single-point calculations.

The extent to which ELMD constitutes a significant approximation to BOMD continues to be debated. When assessing the accuracy of ELMD, the primary consideration is whether there exists a separation of time scales between nuclear oscillations, whose time scale τnuc is set by the period of the fastest vibrational frequency, and electronic oscillations, whose time scale τelec may be estimated according to373

τelec(μεLUMO-εHOMO)1/2 (9.17)

A conservative estimate, suggested in Ref. 373, is that essentially exact classical dynamics is attained when τnuc> 10 τelec. In practice, we recommend careful benchmarking to insure that ELMD faithfully reproduces the BOMD observables of interest.

Due to the existence of a fast time scale τelec, ELMD requires smaller time steps than BOMD. When BOMD is combined with Fock matrix extrapolation to accelerate convergence, it is no longer clear that ELMD methods are substantially more efficient, at least in Gaussian basis sets.374, 794

The following $rem variables are required for ELMD jobs:

AIMD_FICT_MASS
       Specifies the value of the fictitious electronic mass μ, in atomic units, where μ has dimensions of (energy)×(time)2.
TYPE:
       INTEGER
DEFAULT:
       None
OPTIONS:
       User-specified
RECOMMENDATION:
       Values in the range of 50–200 a.u. have been employed in test calculations; consult Ref. 373 for examples and discussion.