Shell model optimisation

In a standard electrostatic embedding QM/MM calculation, the classical MM atoms polarise the QM region by entering the QM Hamiltonian as point charges, but the QM atoms do not in turn polarise the MM region. A shell model force field can be used to enable polarisation of the MM region. Each MM atom consists of a core and a shell, which are connected by a spring force. Some of the MM charge is assigned to the shell, which can move under the influence of the electrostatic potential of the QM region. Because the polarisation of the QM region is in turn influenced by the positions of the shells, each QM/MM energy and gradient evaluation consists of a series of QM and MM calculations which are iterated to self-consistency.

In the context of microiterative optimisation, no special changes need to be made to the algorithms in DL-FIND, as the shells are treated as invisible from the point of view of the optimiser and only the final converged QM/MM energy and gradient is required. However, changes to ChemShell were made to support efficient microiterative shell model optimisation with ESP fitting. Whenever a full QM/MM energy and gradient evaluation is carried out, the polarisation of QM and MM regions is iterated to self-consistency as normal. This is the case for macroiterations, and also microiterations if ESP fitting is not used (although the QM region is fixed if it is in the inner region, the polarisation of the QM region by the MM region will not be fixed). However, if the QM region is approximated by ESP fitted charges, the polarisation of the QM region is also fixed for the duration of the microiterations. The ESP charges are then fed directly to GULP so that the shells can be relaxed in a single GULP calculation (taking into account how the forces on the shells change as they move, rather than feeding in the forces themselves, which would need to be iterated over).

Shell model microiterative minimisation was tested for correctness on the MgO surface defect example used to test the HDLCOpt implementation [9]. The ESP fitting procedure worked well and the optimisation converged to the same minimum as HDLCOpt.

In order to evaluate the performance of all the shell model microiterative optimisation algorithms on HECToR, a test system involving a reaction is required. The reaction we have used is the dissociation of hydrogen over Li-doped MgO, as shown in Figure 4.

Figure 4: QM/MM model of hydrogen dissociation over Li-doped MgO. Hydrogen shown in white, oxygen in red, magnesium in cyan, and point charges in purple. The lithium atom is below the surface.

The test system consists of a cluster of 6349 atoms plus 87 point charges around the edge of the cluster to mimic the effect of the bulk material's electrostatics. The QM region consisted of 33 atoms in the centre of the cluster, including the hydrogen molecule, the lithium dopant, 5 oxygen atoms, and 25 magnesium atoms, of which 4 were treated at a full QM level and 21 were modelled using pseudopotentials to form a boundary around the QM region, which ensures that electron density does not leak into the surroundings. There were 834 active atoms optimised in each run. For the microiterative calculations, the inner region was defined to be the same as the QM region (including boundary atoms).

The QM calculations were performed using GAMESS-UK using the B97-2 functional with a Def2 TZVP basis (274 spherical harmonic basis functions). The calculations were unrestricted with a multiplicity of 2. The MM calculations were performed in GULP using a shell model potential.

The ESP fitting procedure was found to fail for this system, with unphysical charges being assigned to many of the QM atoms. This is probably due to the indistinct nature of the QM-MM boundary, in which some MM oxygen atoms are within the region enclosed by the QM pseudopotential magnesium atoms. A number of remedies were attempted, including fitting only to the full QM atoms while constraining the charges on the boundary QM atoms, but they did not resolve the problem. Another approach would be to develop pseudopotentials for the oxygen atoms so the boundary can be well-defined, but this would be new research and beyond the scope of the current dCSE project. We therefore performed the tests without ESP fitting, which still gives a good indication of the savings that can be expected in terms of macroiterative cycles when the microiterative algorithms are applied. The tests also confirm that all the algorithms can be used with shell model polarised forcefields.

The tests consisted of minimisation of the physisorbed and dissociated hydrogen reaction endpoints, P-RFO and dimer optimisation of the transition state for the dissociation, and NEB optimisation of the full reaction path. The results are given in Table 2.

Table 2: Comparison of energies and number of optimisation steps required for convergence of the hydrogen/Li-doped MgO system for microiterative and standard minimisation algorithms. NEB energies are those of the highest energy image (no climbing image).

Opt.	E(std)/h	E(mic)/h	Standard opt.	Microiterative opt.
type			Cycles	Macro	Micro
H$_2$ min	-1842.669	-1842.669	27	3	22
HH min	-1842.689	-1842.689	50	4	50
P-RFO	-	-1842.662	-	22	95
dimer	-1842.662	-1842.662	47	35	309
NEB	-1842.664	-1842.666	44	32	31

Unlike the case of the water sphere, the minimisation runs here are much more constrained structurally and so it is not surprising that both the standard and microiterative optimisations result in highly consistent energies. In both the molecular hydrogen and dissociated cases, the number of macroiterative cycles is dramatically lower (roughly ten fold less) than the number of cycles in the standard optimisation. This is likely a function of the good starting geometries and regular nature of the system, and this level of improvement is unlikely to be seen on more complex energy surfaces. Nevertheless it is a very positive result for the microiterative algorithm.

For transition state optimisation, the improvements are more modest. It was not possible to run a full P-RFO optimisation due to the computational cost of calculating the Hessian that would be required. However, the good agreement in energies between microiterative P-RFO and both standard and microiterative dimer optimisations shows that the algorithms are working well. The number of macroiterative cycles for the dimer method is lower than the number of standard cycles, albeit not by as large a number as for the minimisations. In this case P-RFO performed better than the dimer method both in terms of macroiterative and microiterative cycles, but it should be emphasised that this may not be true in general, as there are many factors that influence the number of cycles required to converge. Furthermore, even the microiterative version of P-RFO still requires substantial computational expense to calculate a Hessian.

The excellent agreement in energies between the standard and microiterative dimer methods suggests that the microiterative degrees of freedom are spectators to the reaction, which is consistent with the results of the NEB optimisation of the reaction path. 10 NEB images were used (with no freezing of images or a climbing image), and a relatively loose convergence criterion was employed. This, together with task-farming parallelism (10 workgroups), was necessary to make the NEB calculation tractable on HECToR. The convergence criterion likely explains the small difference in optimised energies between the highest energy standard and microiterative images, rather than relaxation of the environment which would have shown up in the transition state optimisations. The reduction in cycles is similar to that of the dimer method, although much fewer microiterations are required in the case of NEB, which indicated that the initial interpolated reaction path is of very good quality. This is not surprising for a very constrained system with many spectator degrees of freedom.

Overall, the results show that all the microiterative algorithms are functioning correctly for shell model optimisation of materials systems on HECToR, and that microiterative relaxation of the environment does reduce the number of macroiterations required in the test system.