Calculating Shieldings for Proteins

I'm just going to show some recent work on extending the Polarizable Embedding from being a focused QM/MM model into an approach to generate parameters for entire proteins or DNA. Here is a snapshot of the human protein GB3 (PDB: 2OED) with the central QM-fragment we are interested in - here LYS4 - shown in vdw-spheres. Also in the QM-region are neighbouring residues and stuff (side-chains and water) that binds to the central fragment of interest) all highlighted in cyan in sticks. The rest of the protein along with a lot of water is the MM-region. Only the protein is shown in a classical cartoon rendering.

Enjoy

Oxygen, you little rascal!

Calculations of NMR shielding constants using my recent DFT-PE implementation has been underway for some time and I wanted to share the latest results. DFT-PE is a QM/MM type method where the MM region is treated by a polarizable force field calculated for that particular geometry. Initially, I just went away and calculated the NMR shielding constants for acrolein and a few waters with QM and as large an MM region as I could possibly cut out. Unfortunately, the results looked a bit weird when we compared them to a non-polarizable force field (TIP3P) of lesser quality. Lesser in the sense that the description of the electrostatic potential from TIP3P does not match that of the QM electrostatic potential (magnus paper?).

The problem was that the TIP3P results seemed to converge really fast with respect to the size of the QM region and I, using a high-quality PE potential, was converging slower. I went back to the drawing board and came up with the following analysis

Here we see that for increasing sizes of QM region (red) the NMR shielding constant of ¹⁷O is reduced by 102 ppm going from acrolein in the gas phase to acrolein surrounded by 53 water molecules (follow the diagonal). Going vertically up in each row, we decrease the size of the QM region, gradually replacing QM waters with MM waters (blue). Deviations from the QM result for that particular system size is shown below each illustration.

Acrolein in 53 MM waters deviates by 23.3 ppm compared to the full QM calculation. This deviation is (if the polarizable force field we are using) purely quantum mechanical. If the first solvation shell is included, that is 7 QM waters, the deviation drops to 4.3 ppm and only improves from here.

Needless to say, TIP3P is also doing great here, but the convergence is not nearly as systematic as is observed for PE. So, it looks like the PE potential is doing it wrong for the right reasons, where as TIP3P is right for the wrong reasons ... that is, using DFT with TIP3P, two wrongs does make a right!

Final note: Oxygen is really, really hard to get right. I showed you only the most difficult case I tried – NMR shielding constants for the rest of acrolein are converged much faster using the smaller QM regions shown above. Luckily for later studies on proteins, nobody does that silly Oxygen atom anyways.

This work is licensed under a Creative Commons Attribution 3.0 Unported License.

NMR shielding constants through polarizable embedding

After I finished my Ph.D. in the Jan Jensen group, I've begun working at the University of Southern Denmark with Jacob Kongsted.

Apart from using DALTON instead of GAMESS I've also switched gears and I am now focused on an advanced form of QM/MM which is called polarizable embedding (see this paywalled link for details). Basically your gas-phase Fock operator is extended with interactions from the surroundings. In the language of polarizable embedding (PE) we constructs the effective Fock operator here for Hartree-Fock

$\hat{f}_{eff} = \hat{f}_{HF} + \hat{\nu}_{PE}$

The PE is a sum of interactions from the surroundings onto a molecule of interest. I will deal with the specifics of the PE approach later on.

This is all fine and dandy, but until now the PE model has focused solely on electronic molecular properties. I've extended it so you can calculate nuclear magnetic shielding constants in the PE model using London atomic orbitals (you might know them as Gauge-Including Atomic Orbitals) with contributions from charges, dipoles (induced and static) and quadrupoles. The hope is that we need a small QM region due to a very high-quality potential compared to previous QM/MM studie (see here for an example of needed more than 1000 atoms for converged NMR shieldings).

We'll see how it goes.