Energy minimisation

The microiterative minimisation scheme in DL-FIND follows program logic similar to the original HDLCOpt implementation. The program flow of a microiterative minimisation with electrostatic embedding in DL-FIND is as follows:

1. Perform a full QM/MM calculation (macroiteration).
2. Fit point charges to the electron density in the QM region using the electrostatic potential (ESP). This is used to approximate the QM-MM interaction during the microiterations without requiring a QM calculation.
3. If not the first cycle, test convergence of the macroiteration (rejecting the step if the energy has increased). Terminate if convergence criteria are met.
4. Take a step in the inner region using the gradient from step 1.
5. Enter the microiterations:
   1. Calculate an ESP-approximated energy and gradient.
   2. Test convergence of the microiteration (if not the first). Continue to step 6 if converged. Reject the step if the approximated energy has increased and reduce the step size
   3. Take a step in the outer region using the ESP-approximated gradient
   4. Continue at step 5a
6. Continue with macroiterations at step 1.

Note that in the original HDLCOpt implementation the macroiterative step energy (i.e. the exact energy before the inner region step is taken) was used as a reference energy to test acceptance of the first microiterative step. This is acceptable for a minimisation because the energy should decrease in both the macroiterations and microiterations. However, the DL-FIND implementation is designed to work for methods where the energy may not decrease after a macroiterative step, so a reference ESP-approximated energy is calculated before any outer region steps are taken.

The main DL-FIND loop cycles once per energy evaluation. A flag (glob%imicroiter) is used to switch between macroiterations and microiterations.

At the macroiterative step, a full QM/MM calculation is carried out, as in a standard optimisation, by a call to the ChemShell QM/MM driver, giving an energy $E^0$ and gradient $\mathbf{g}^0$. In the same call the QM code is then used to calculate the electrostatic potential on a van der Waals surface around the QM region. A least-squares fit algorithm in ChemShell is used to fit point charges at the QM atom positions to reproduce the ESP. A second call to ChemShell is then made to recalculate the QM/MM energy and gradient using the point charge representation for the QM atoms, giving $E^{0}_{\text{ESP}}$ and $\mathbf{g}^0_$ESP.

Following Ref. [9], the corrected gradient $\mathbf{g}^1$ and energy $E^1$ used in the microiterations are:

$$\mathbf{g}^{1} = \frac{\partial E^{1}_{\text{ESP}}}{\partial \mathbf{x}_{\text{outer}}} + \mathbf{g}^{0}_{\text{corr}}$$ (1)

and

$$E^{1} = E^{1}_{\text{ESP}} + E^0 - E^0_{\text{ESP}} + \mathbf{g}^{0}_{corr} \cdotp ( \mathbf{x}^1_{\text{outer}} - \mathbf{x}^0_{\text{outer}} )$$ (2)

where

$$\mathbf{g}^{0}_{\text{corr}} = \left. \frac{\partial E^{0}}{\part... ...{\partial \mathbf{x}_{\text{outer}}} \right\vert _{\mathbf{x}^0_{\text{outer}}}$$ (3)

where superscript 1 indicates the current coordinates and 0 the respective coordinates from the last macroiterative calculation. The expressions for the corrected energies and gradients ensure a smooth connection between the microiterative and macroiterative cycles, as in the case of identical geometries the equations reduce to the exact answer. In DL-FIND the corrections are actually calculated for the full set of coordinates (not just the outer coordinates), but the energy expression collapses to the outer coordinates because the inner region is fixed, and the gradient correction for the inner coordinates is not used by the microiterative optimiser.

DL-FIND supports multiple instances of the low memory BFGS (L-BFGS) minimisation algorithm, and two instances are used during a microiterative minimisation. The first L-BFGS instance is used for macroiterative steps and covers all coordinates, although only inner region steps are permitted during macroiterations. The reason it covers all coordinates and not just the inner region is because the gradients on all atoms can then be taken into account when calculating the inner region step, and therefore the outer region does not have to be completely converged during the microiterative cycles for the optimisation to succeed. The second L-BFGS instance is initialised at the start of every microiterative loop and covers the outer region coordinates only. The corrected QM/MM energy and gradient is used to relax the environment while the inner region is kept fixed. The second L-BFGS data is kept in a separate set of arrays in the microiter module.

Unlike HDLCOpt, DL-FIND supports multiple coordinate systems and therefore special care must be taken during the microiterative cycles that the correct subset of coordinates are used. When the Cartesian coordinates used externally are transformed to the internal list of coordinates used by the optimiser, the new list is ordered into inner region coordinates followed by outer region coordinates. This makes it straightforward to operate over a single region. When Cartesians are used internally, this mapping is trivial. The delocalised coordinate (DLC) and total connection (TC) coordinate systems, however, may not be used, as they are delocalised over the entire system and therefore cannot be decomposed by region. Hybrid delocalised coordinates (HDLCs) and HDLC-TC may be used as they are only delocalised over discrete subsets of coordinates called residues, providing that no residues cross the inner and outer regions. The input is checked to ensure that this does not happen. In the HDLC case the list of internal coordinates is ordered by residue within each region.

A number of keywords were added to the DL-FIND input to specify that a microiterative optimisation should be carried out, the list of atoms or residues in the inner region, whether to carry out an ESP fit (a full QM calculation can also be used during the microiterations if requested), and the maximum number of microiterative cycles to carry out before taking another macroiterative step.

The microiterative minimisation scheme was tested for correctness on a QM glycine molecule solvated in water described by a classical forcefield. To check that the implementation gave similar performance improvements to the original HDLCOpt implementation, it was further tested on the water sphere system described in Ref. [9]. This consists of 4404 water molecules cut from equilibrated bulk water with a radius of 25 Å, as shown in Figure 1.

Figure 1: Water sphere used to test the DL-FIND implementation of microiterative QM/MM minimisation.

Three water molecules in the centre of the sphere were defined as the QM region and inner microiterative region. QM calculations were performed with the MNDO semi-empirical QM package using the AM1 Hamiltonian. Active regions of radius 5, 10, 15 and 20 Å were defined, within which the water molecules were described using the TIP3P force field. The atoms outside the active region were frozen during the optimisation.

Results for the optimisations are shown in Table 1.

Table 1: Comparison of number of optimisation steps required for convergence of a QM/MM water sphere optimisation using microiterative (QM=macro cycles and MM=micro cycles) and standard minimisation algorithms. DOF = degrees of freedom. $\Delta E$ is the difference in converged energy between the microiterative and standard runs.

Active		Microiterative opt.		Standard opt.	$\Delta E$ /
radius/Å	DOF	QM	MM	QM	kcal/mol
5	2493	114	2018	639	-2
10	6966	176	3230	2038	+17
15	14382	311	4682	1372	-3
20	26136	207	3949	1990	+10

Note that the DL-FIND tests were carried out in Cartesian coordinates (and therefore without bond constraints), as microiterative optimisation with HDLC coordinates was not implemented at the time. The number of degrees of freedom is therefore greater than in Ref. [9] and so the number of optimisation steps are not directly comparable. However, very similar trends are observed. As in the previous HDLCOpt tests, the optimisation takes place over a complex energy surface with many close lying minima and so different runs may take different paths and reach slightly different converged energies. For the same reason the number of cycles does not necessarily increase monotonically with the number of degrees of freedom. Nevertheless, the improvement in performance from standard to microiterative optimisation is clear. Over all sizes of active region, a reduction in the number of QM evaluations is seen of between 5 and 12 times. An illustration of the significant speed-up of convergence in terms of number of QM steps is shown in Figure 2.

Figure 2: Maximum gradient component (one of four convergence criteria) during microiterative and standard optimisation of a water sphere with an active region of 20 Å.

With a small, semi-empirical QM region this will not result in an overall reduction of computational time, but it demonstrates the potential for significant savings in typical production QM/MM calculations where the QM calculation dominates.

Microiterative optimisations with HDLC coordinates and HDLC constraints were subsequently tested using the glycine-water system. The converged HDLC minima were in good agreement with those found using Cartesian coordinates.