Implementation

To obtain the QM/MM potential energy every molecular system was partitioned into a reactive part (the core^4.3) treated quantum mechanically with the PM3 [32] or the AM1 [31] semiempirical Hamiltonians and a nonreactive part (the environment) treated by means of molecular mechanics with the AMBER force field[51].

The QM atoms are influenced by the partial charges of the MM atoms, and, in addition, bonding and van der Waals interactions between the two regions are included consistently. All the non-bonded interaction parameters used for $H_{QM/MM}$ (see section 1.2.3.1) are taken from AMBER force field. We used hydrogen link atoms to cap exposed valence sites due to bonds which cross the QM/MM boundary. The AMBER 5.0 (Roar-cp module)[235] program was used to carry out the QM/MM calculations. This potential energy method is the same as the one used in the Mandelate Racemase QM/MM study in section 2.3, and its equations are described in the introduction section 1.2.3.1.

The search algorithm works in the following way (the name of the subroutine that performs the task is indicated in italics. The implemented source code is described in the appendix, page ):

1. An initial geometry ${\bf q}_0$ in Cartesian coordinates is chosen. The QM region, the MM region and the corresponding link atoms are defined. Then, three environment atoms are always fixed in order to get rid of the translations and rotations of the whole molecular system. (subr. rdigpr pqqorder)
2. The QM/MM energy and the gradient at ${\bf q}_0$ are calculated using the Roar-cp module. (subr. energy_qmmm)
3. The guess ${\bf B}_0$ Hessian matrix is built up according to figure 3.1 (subr. calchess). To this aim, we implemented numerical second derivatives in the Roar-cp module as required here. Diagonalization of ${\bf B}_0$ allows us to test whether ${\bf q}_0$ lies on the suitable quadratic region of the PES (i.e., zero or one negative eigenvalue for searching a minimum or a transition state, respectively). If this is not the case, the ${\bf B}_0$ matrix is forced to have the convenient number of negative eigenvalues. (subr. corr_eigval)
4. The RFO method plus a procedure to solve the corresponding secular equations through the full diagonalization of the AH matrix are used in order to obtain the displacement vector at each iteration. (subr. ts_rfo: build_ah, diagonalize, choose_eig)

   This vector is scaled by a factor according to the procedure described in section 3.1.1.1 (subr. calc_step). Then, the new geometry ${\bf q}_{k+1}$ is obtained as a result of the current $k$ iteration.

5. If the RMS of the gradient at ${\bf q}_{k+1}$ is smaller than a suitable convergence criterion, it is considered that the corresponding stationary point has been reached and the search ends (subr. eval_step). Otherwise, the algorithm proceeds to step 6.
6. For minima, the modified BFGS formula (m-BFGS) given in equations 3.5, 3.6 and 3.7 is used to update the Hessian matrix. For transition states the Powell formula is employed (equation 3.5). Owing to the update the QM and the MM parts of the resulting ${\bf B}_{k+1}$ Hessian matrix become coupled and the MM part is no longer diagonal (as displayed in the initial Hessian of figure 3.1). (subr. update_hess)

   Diagonalization of ${\bf B}_{k+1}$ allows us to test whether ${\bf q}_{k+1}$ lies on the suitable quadratic region of the PES. If this is not the case, the ${\bf B}_{k+1}$ matrix is forced to have the convenient number of negative eigenvalues (subr. corr_eigval). Then the algorithm proceeds to the step 4 to start a new iteration cycle.

All the calculations, except the evaluation of energies and gradients, were carried out with the algorithm just described. Hereafter this algorithm is called RFO-m-BFGS or RFO-Powell when a minimum or a transition state, respectively, is looked for.

The BFGS and L-BFGS procedure were implemented in a development version of ROAR[235] by Gérald Monard and Mireia Garcia-Viloca.