10.7.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
$v_{x,1}{v}_{y,1}{v}_{z,1}$
$v_{x,2}{v}_{y,2}{v}_{z,2}$
$v_{x,N}{v}_{y,N}{v}_{z,N}$
$end

Here $v_{x,i}$, $v_{y,i}$, and $v_{z,i}$ are the $x$, $y$, and $z$ Cartesian velocities of the $i$th 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 10.7.5).

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

$mass
$m_1{m}_2$
$m_3$ ! mass of 3rd atom
$m_4$
$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 $x_1$, $y_1$, $z_1$, $x_2$, $y_2$, $z_2$, …, $x_N$, $y_N$, $z_N$ 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 $\mathrm{\Delta }𝐅=𝐅(t+\delta 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. 370.

• •

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

• •

  ElecGradNorm: Records the norm and largest magnitude element of the electronic gradient matrix $\mathrm{𝐅𝐏}-\mathrm{𝐏𝐅}$ 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:

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,³⁷¹ based on a suggestion by Pulay and Fogarasi,⁷⁶⁶ is available for this purpose.

The Fock matrix elements $F_{\mu \nu }$ 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$:

The $N+1$ extrapolation coefficients $c_n$ 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. 371, 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:

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.³⁰ While these equations are linear, their dimensionality necessitates an iterative solution,^{749, 486} which, in practice, looks much like the SCF equations. Extrapolation is again useful here,⁸⁷⁸ and the syntax for $Z$-vector (response) extrapolation is similar to Fock extrapolation.

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 $\mu$, and the value of $\mu$ controls both the accuracy and the efficiency of the method. In the limit of small $\mu$ the nuclei and the orbitals propagate adiabatically, and ELMD mimics true classical dynamics. Larger values of $\mu$ 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,³⁷⁰ 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 ${\tau }_{\mathrm{nuc}}$ is set by the period of the fastest vibrational frequency, and electronic oscillations, whose time scale ${\tau }_{\mathrm{elec}}$ may be estimated according to³⁷⁰

A conservative estimate, suggested in Ref. 370, is that essentially exact classical dynamics is attained when ${\tau }_{\mathrm{nuc}}>$ 10 ${\tau }_{\mathrm{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 ${\tau }_{\mathrm{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.^{371, 766}

The following $rem variables are required for ELMD jobs: