Release66:Qmmm NEB Calculations

From NWChem

Template:Release66:QMMM Optimization Links Experimental implementation of Nudged Elastic Band (NEB) method is available for reaction pathway calculations with QM/MM. The actual pathway/beads construction involves (by default) only the region containing QM and link atoms (referred to as qmlink). The rest of the system plays a passive role and is quenched/optimized each time a gradient on a bead is calculated.

The initial guess for NEB pathway can be generated using geometries of the starting and ending point provided by the .rst files. These are set in the input using the following directive

set qmmm:neb_path_limits xxx_start.rst xxx_end.rst

where xxx_start.rst xxx_end.rst refers to starting and ending point of the NEB pathway. Both rst files have to be present at the top level directory. It should be noted that only coordinates of qmlink region will be used from these two files. The initial coordinates for the rest of the system come from reference rst file provided in the MD block

md
  system xxx_ref
  ...
end

Typically this reference restart file ( xxx_ref.rst ) would be a copy of a restart file for starting or ending point.

The number of beads in the NEB pathway, initial optimization step size, and number of optimization steps are set using the following directives

set neb:nbeads 10
set neb:stepsize 10
set neb:steps 20

The calculation starts by constructing initial guess for the pathway (consisting of a sequence of numbered rst files) by combining linearly interpolated coordinates of the qmlink regions from starting and ending rst files and classical coordinates from the reference file. Next phase involves calculation of the gradients on qmlink region atoms for each of the beads. This involves two steps. First classical region around the qmlink region is relaxed following standard QM/MM optimization protocol. Aside the fact that optimization region cannot be qmlink, all other optimization directives apply and should be set in the QM/MM block following standard convention, e.g.

qmmm
 region  solvent
 maxiter 1000
 ncycles 1
 density espfit
end

or

qmmm
 region  mm_solute solvent
 maxiter 300          1000
 ncycles 3
 density espfit
end

In both examples presented above we utilized espfit option for density to speed up calculations. Note that optimization region cannot be qmlink!

After the optimization has been performed the gradient on qmlink region is calculated. The procedure is repeated for all the beads. After that the bead coordinates will be advanced following NEB protocol and the entire cycle will be repeated again.

In addition to interpolated initial guess, one can also specify custom initial path represented by numbered sequence of restart files stored in the perm directory. This behavior will be triggered automatically in the absence of qmmm:neb_path_limits directive. The default naming of the custom initial path is of the form <system>XXX.rst, where <system> is the prefix of reference restart file as defined in MD block and XXX is the 3-digit integer counter with zero blanks (001,002, ..., 010, 011, ..). If needed the prefix for the custom initial path can be adjusted using

set qmmm:neb_path <string prefix>

The progress of NEB calculation can be monitored by

grep gnorm <output file>

Experience shows that the value of gnorm less or around O(10^-2) indicates converged pathway. The current pathway in the XYZ format can be found in the output file (look for XYZ FILE string) and viewed as animation in some of the molecular viewers (e.g. JMOL)

Here is an example of QM/MM NEB calculations, which involves generation of the initial guess for the pathway based on restart files for starting (cloh_rs.rst) and ending (cloh_ps.rst) point. The reference restart file (cloh_neb.rst) is simply a copy of starting (cloh_rs.rst) restart file. As always all the restart and topology (cloh.top) files have to reside in the top level directory. It is very important to explicitly set permanent and scratch directory to be distinct from top level calculation directory.

This particular system consists of Cl- and OH radical species in aqueous solution. The reaction pathway will only involve coordinates of Cl- and OH, which constitutes the QM region. A total of 10 beads will be used, using 5 NEB refinement steps with initial step size of 1.0. The solvent part of the system will be always optimized with respect to current bead configuration, involving in this case 1000 steepest descent iterations (see region, maxiter directives). During solvent optimization, QM region will be represented by point ESP charges (see espfit directive). The entire example directory including output file can be downloaded here.

memory total 1500 Mb
start cloh
permanent_dir ./perm
scratch_dir   ./data
charge -1
basis "ao basis"
 * library "6-31+G*"
end
dft
mult 2
XC b3lyp
iterations 5000
end
md
 system cloh_neb
 cutoff 1.0
end
qmmm
bqzone 10.0
region  solvent
maxiter 1000
ncycles 1
density espfit
end
set qmmm:neb_path_limits cloh_rs.rst cloh_ps.rst
set neb:nbeads 10
set neb:stepsize 1.0
set neb:steps 5
task qmmm dft neb