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.

FragIt and Metals

FragIt is a tool that helps you setup a fragment calculation, more specifically for use in the Fragment Molecular Orbital (FMO) method and its Testosterone induced cousin the Effective Fragment Molecular Orbital (EFMO) method.

Since my switch to the University of Southern Denmark, FragIt's capabilities have grown a bit. Now it supports Molecular Fragmentation with Conjugate Caps so it can dump coordinate files and you can estimate interaction energies quickly and I am trying to switch to a more flexible configuration system, but that is a topic for another post on a related blog.

However, my proudest moment regarding FragIt must have been today when I finally fixed an issue with FragIt that has been haunting me for some time. Metal ions. The problem is that FragIt uses the excellent Open Babel toolkit as its underlying work horse, but the charge models that are implemented are misbehaving when metals are present (through the Python bindings). The partial charges are not present for the metals which I guess makes sense

My solution was rather elegant:

while metals are present

• store a copy of them in a separate list
• delete them from the system
• calculate the partial charges of the system without the metal ions.

Finally, reinsert the metal ions and use their formal charge (which you are certainly most interested in) as the partial charge of the ion. From here, everything is as usual and you can do all your usual tricks with FragIt as you are used to.

The only catch is that whatever metal ions you want to include has to be last in whatever structure file you are using and I should mention that the web interface is not yet updated with this feature.

That's what conferences are for! Stay tuned for a 1.0.7 point release sometime in the near future.

btw - metal ions can be mispelled to say meta lions which I at this hour find rather amusing.