Ab Initio Steered Molecular Dynamics (AISMD)

Specific atoms can be subjected to constant external forces to generate force-modified potential energy surfaces (FMPESs) and to run molecular dynamics on them.¹ TeraChem supports two distinct kinds of steering, selected via the steering keyword:

Fixed steering — atoms are pulled toward lab-frame fixed points. The modified PES depends on the orientation of the molecule. This is the appropriate choice when running dynamics in the presence of pulling forces.
Adaptive steering — atoms are pulled toward (or pushed away from) other atoms in the molecular frame. The modified PES is invariant to rigid rotations of the molecule. This is the appropriate choice for geometry optimization and minimum-energy-path calculations under pulling forces.

Fixed steering

In fixed steering, each steered atom is attached to a Cartesian "anchor point" in the lab frame and a constant force pulls the atom toward that point. The total potential energy is

\[V_{total}(r) = V_{ab\text{-}initio}(r) + \sum_i^{N_{attach}} F_0^i \left( \lVert r_i^{fix} - r_i \rVert - \lVert r_i^{fix} - r_i^{0} \rVert \right)\]

where \(r_i^{0}\) is the initial position of the \(i\)th steered atom, \(r_i^{fix}\) is its fixed anchor point, and \(F_0^i\) is the steering force applied to it. Because the anchor points are specified in the lab frame, the energy depends on the orientation of the molecule (the anchors do not rotate with it).

Steering file format

Fixed steering is configured through a separate text file (default name: steering) located in the same directory as the input file. The first line gives the number of steered atoms \(N_{attach}\); subsequent lines have the form

<atom_index>  <fixed_point_x>  <fixed_point_y>  <fixed_point_z>  <force>

Coordinates of the anchor points are in bohr and the force is in atomic units of force (1 nN \(\approx\) 0.012138 a.u.). Atom indices are 0-based (the first atom in the coordinate file is atom 0). The maximum number of steered atoms is 16.

A two-atom example pulling atom 3 toward \((-50, 0, 0)\) bohr and atom 5 toward \((+50, 0, 0)\) bohr with force 0.04 a.u. each:

steering

2
3 -50 0 0 0.04
5  50 0 0 0.04

Activating fixed steering

In the input file:

Value of `steering`	Effect
`yes`	Read from the default file `steering` in the input directory
(any other string)	Read from that file path
`no`	Disabled

steered_md.in

coordinates       tDCC4nN.xyz
basis             6-31g*
method            ub3lyp
charge            0
dftd              d3
purify            no

run               minimize
new_minimizer     yes
min_coordinates   hdlc
convthre          1e-6
precision         mixed
nstep             18

steering          steering.xyz    # custom steering-file name
end

Adaptive steering

In adaptive steering, pulling forces are defined along the line between two specific atoms in the molecule, so the forces rotate with the molecule. The total potential energy is

\[V_{total}(r) = V_{ab\text{-}initio}(r) + \sum_i^{N_{pairs}} F_0^i \, \lVert r_i^{left} - r_i^{right} \rVert\]

where \(N_{pairs}\) is the number of atom pairs subject to pulling forces and \(r_i^{left}\), \(r_i^{right}\) are the positions of the two atoms in the \(i\)th pair. A positive force expands the distance between the pair (pulling them apart); a negative force compresses it (pushing them together).

Adaptive steering is not recommended for molecular dynamics; it is most appropriate for geometry optimization and minimum-energy-path calculations under external forces.

Adaptive steering keywords

Adaptive steering is configured entirely from input-file keywords (no separate file):

Keyword	Description
`steering adaptive`	Activate adaptive steering
`steeratom1`	Comma-separated list of "left" atom indices, one per pulling pair
`steeratom2`	Comma-separated list of "right" atom indices, one per pulling pair
`steerforce`	Comma-separated list of pulling forces (atomic units), one per pair

Atom indices for adaptive steering are 1-based (the first atom is atom 1). There must be no whitespace within the comma-separated lists. The three lists must all have the same length; the maximum number of pulling pairs is 5.

Example: two-pair benzene pulling

aismd_benzene.in

basis        6-31g
method       b3lyp
coordinates  c6h6.xyz
run          minimize

steering     adaptive
steeratom1   1,4
steeratom2   6,5
steerforce   0.012,-0.012
end

This input applies a pulling force of 1 nN (0.012 a.u.) between atoms 1 and 6 (tending to expand their separation) and a pushing force of 1 nN between atoms 4 and 5 (tending to compress their separation).

Example: single pair, QM/MM with Amber

aismd_qmmm.in

prmtop          molecule.prmtop
coordinates     initial.rst7
qmindices       qm_indices.txt

method          b3lyp
basis           3-21g
charge          0
purify          no

run             minimize
nstep           15
min_coordinates cartesian

steering        adaptive
steeratom1      1
steeratom2      3
steerforce      +0.03641340     # 3.0 nN x 0.012138 a.u./nN
end

M. T. Ong, J. Leiding, H. Tao, A. M. Virshup, and T. J. Martinez, J. Amer. Chem. Soc. 131, 6377 (2009). ↩