Oxygen, you little rascal! Part 2.

Background
In the ongoing saga (read part I) with acrolein solvated in water, we were having problems with our high quality embedding potential just not converging nicely when increasing system sizes, QM-region sizes. We were using a functional to get good shielding constants (the KT3 functional) with a basis set designed to get good shielding constants (the apcS-1 basis set) but no matter what tried, the widely adopted TIP3P potential for water just did it better than we did. And we were just confused!

Solution
The problem turned out to be the basis set that with the diffuse functions leaked electron density out into the MM-region when using the PE potential and this was causing some disturbance with the induced dipoles. On the contrary, our belief is that TIP3P is a much too soft potential and there is a large degree of error cancellation because of that. Both are effects I would think to be very small and have no real impact but I had also forgotten that we were talking about nuclear shielding constants which are sensitive as f*ck the wave function.

I switched to the smaller pcS-1 basis set and all that leaking-trouble went away. Everything that we did not expect was gone. TIP3P was converging rather slowly and the PE-model was quickly converging as is shown below for an average of 120 snapshots.

The PE numbers (blue) are converged at QM size 3 but even size 1 and 2 show some good indications. For size 1, acrolein is solvated in purely classical water. As we would expect, the improved potential means that we are converging (fast) towards a value for the absolute shielding of acrolein in water valued -223.2 ppm.

Outlook
We will return to the use of diffuse basis functions and a proper quantum mechanical fix in a later paper as this work here sparked my co-worker +Jógvan Magnus Haugaard Olsen into actually picking up work he had done on repulsion during his Ph.D. (pdf of his thesis)

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

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.