Tuesday, February 15, 2011

QM/MM course: notes from week 2

QM/MM Yang Paper
View more presentations from molmodbasics.

I presented the method used in the paper by Parks, Hu, Rudolph, and Yang.
Links to the papers can be found here.

The QM region is defined on slide 2 and is treated at the B3LYP/6-31G(d) level if theory while the rest of the protein and solvent is treated by the Charmm22 and TIP3P force field, respectively.  Gaussian is used for the QM calculations and Sigma is used for the rest.

The covalent boundary is treated by Yang's pseudobond method, where the link atom is described by a pseudopotential that gives the atom a valence of 1, but an equilibrium bond length close to that of a C-C single bond.

An explicit solvation model is used, with a 12-Å shell of frozen water to handle the water-vacuum boundary. No other constraints are imposed.

An approximate TS  is located with the Quadratic String Method (QSM) - a method used to interpolate between reactants and products.  It is not clear from the paper whether the approximate TS is refined further with QST3.

Protein/solvent dynamic is handled with Yang's Minimum Free Energy Path method, where the "QM atoms are optimized in the environment of the fluctuating MM subsystem.":
1. A QM/MM energy and gradient is computed for a reference QM and MM geometry along the QSM path.
2. The charges of the QM atoms are computed, and are used to run all and all MM MD (keeping the QM geometry frozen).
3. The energy and gradients are corrected using a Boltzmann weighting scheme and approximating the QM/MM electrostatic energy as point charge interactions.
4. The gradient is then used to minimize the QM region (updating the Boltzmann correction at each step), using the trajectory geometries from 1.  This is apparently done with Gaussian using the "external" keyword (slide 10)
5. Steps 1-4 are repeated so self consistency for each point along the QSM path.
6. Upon self-consistency, the vibrational frequencies of the QM region are computed.

It is not clear from the paper what charge scheme (slide 13) the study uses.

From our discussion
This is perhaps the most sophisticated treatment of dynamics in QM/MM to date, but also the potentially most complicated and error-prone.

From a "dynamics perspective" it is best to keep the QM region as small as possible, since its conformational space is not explored extensively, but the accuracy of the QM/MM interaction energy may suffer.

The parameters in the pseudobond method is basis set dependent (and optimized for 6-31G*).  It is not clear how transferable the parameters are to larger basis sets in the QM/MM region.  One could also mix the basis and use 6-31G* near the pseudobond, but this might cause problems for short side chains such as cysteine (slide 2)

No comments: