Skip to content

Ab Initio Nanoreactor

The ab initio nanoreactor1 couples fast TeraChem dynamics with a virtual piston that periodically compresses and decompresses the system, dramatically accelerating the rate at which reactive events are sampled. The piston is implemented by alternating between two sets of spherical boundary conditions on fixed schedules.

The technique relies on the same spherical boundary conditions described in the AIMD overview. What makes the nanoreactor distinct is that the inner and outer radii (and force constants) are swapped back and forth with a regular period: a large, weakly confining sphere for several hundred timesteps, then a smaller, more aggressively confining sphere for a shorter window, then back to the large sphere, and so on. The collision rate during the compression phase is what drives reactivity.

Time-dependent boundary conditions

The schedule is set by two keywords:

Keyword Description
mdbc_t1 Number of timesteps for which the first set of BCs (md_r1, md_k1) is active per cycle
mdbc_t2 Number of timesteps for which the second set of BCs (md_r2, md_k2) is active per cycle

After mdbc_t1 steps under the first set, TeraChem switches to the second set for mdbc_t2 steps, then back to the first set for another mdbc_t1 steps, and so on for the rest of the trajectory.

By convention the first set is the larger, equilibration sphere and the second set is the smaller, compressing sphere — but the keywords don't care about the ordering; whichever sphere is "second" is the one that becomes active on the second leg of each cycle.

Two additional mdbc_* switches are commonly turned on for nanoreactor runs:

  • mdbc_hydrogen yes — by default the spherical confining potential is not applied to hydrogen atoms; nanoreactor runs typically want it applied to hydrogens as well.
  • mdbc_mass_scaled yes — scale the confining potential by atomic mass. Without this, compression can strip hydrogen atoms from heavy fragments. See the original paper for the rationale.

The four BC parameters md_r1, md_k1, md_r2, md_k2 and the two cycle times mdbc_t1, mdbc_t2 all need to be tuned to the system being studied.

Example

The example below sets up an unrestricted HF nanoreactor run at 1500 K, with a 6 Å radius weakly-confining sphere for 750 steps alternating with a 4 Å more strongly-confining sphere for 250 steps:

nanoreactor.in
coordinates       nanoreactor.xyz
basis             3-21g
method            uhf
charge            0
dispersion        no

# Thermostat: Langevin at 1500 K, starting from 1200 K
run               md
tinit             1200
thermostat        langevin
t0                1500
lnvtime           200

# SCF aids for radical chemistry
scf               diis+a
convthre          0.005
levelshift        yes
levelshiftvala    0.3
levelshiftvalb    0.1

timings           yes
nstep             30000
maxit             300

# Nanoreactor BCs: large sphere then small sphere, repeating
mdbc              spherical
md_r1             6.0
md_k1             3.0
md_r2             4.0
md_k2             5.0
mdbc_hydrogen     yes
mdbc_mass_scaled  yes
mdbc_t1           750
mdbc_t2           250
end

The level shift on the alpha and beta channels helps keep the UHF SCF from collapsing onto a closed-shell-like solution during transient radical encounters, which are routine inside the piston.

Analysis of nanoreactor trajectories — clustering reactive events and refining minimum-energy reaction paths — is handled by a separate companion package that is not (yet) distributed with TeraChem.

  1. L.-P. Wang, A. Titov, R. McGibbon, F. Liu, V. S. Pande, and T. J. Martinez, Nature Chem. 6, 1044 (2014). ↩