Does increasing Lennard-Jones cutoff means higher accuracy?

Question

I am currently trying to simulate a large box of QM/MM water molecules with SPC/Fw model and GGA DFT functional (a cube of 40 angstrom sides). To validate the QM and MM simulations separately, I am testing a pure MM box of same dimensions with SPC/Fw model, and I was trying to set the LJ and coulomb cutoff length. I am using CP2K software and pairwise-style cutoff. By the conventions I should set the cutoff to <1/2 of smallest box dimension (~19 angstroms). But I encountered these two papers about this issue for liquid water:

• Liquid water simulation: A critical examination of cutoff length

• The Origin of Layer Structure Artifacts in Simulations of Liquid Water

both published in 2006, which suggests that larger cutoff lengths can lead to "layering" and other artifacts in pure water simulations. Elsewhere, I can find papers which suggest to keep the cutoff length to 9-12 angstroms (as discussed here: Rule-of-thumb for Morse potential cutoff in molecular dynamics?). So my question is: Does it makes sense then to limit cutoffs till 12 angstroms? And are there any significant issues particularly for QM/MM simulations?

For more details on the QM/MM setup: I am using an additive QM/MM system with electrostatic embedding.

Answer 1

Very interesting question!

The LJ potential is quite benign in that it has a $r^{-6}$ decay, which is pretty rapid. However, MD codes typically go a bit farther than a straightaway truncation at $r_{\rm max}$, but instead also apply an analytical correction for the region $[r_{\rm max},\infty)$. The analytical correction for the cut-off can be justified by long-range statistics: things look more and more uniform the farther you go from a point.

However, the cutoff needs to be large enough so that things really become uniform at that stage, that is, the pair correlation function $g(r)=1$ for $r \ge r_{\rm max}$; otherwise you end up affecting the behavior of the system.

While in principle increasing the cutoff should make simulations more accurate, the findings of the two papers you linked would seem to insinuate that this is not true.

However, one should not think that it is the cutoff that is to blame, but rather the water model. It looks like the water models have either been trained for systems that are too small, and/or the training model calculations have not been fully converged with respect to all computational parameters. What one is seeing is emergent behavior from an incorrect water model: when the calculations are performed numerically accurately, the system stops behaving in a physical manner. A smaller cutoff implies forcing an uniform structure at long range; if you instead use a larger cutoff, the systems starts to exhibit long-range structures.

So, while my general advice is to converge the cutoffs so that the system's properties do not change anymore, if you have a case where it has been reported that too large values for the cutoff result in inaccurate behavior, then you need to limit the cutoffs to smaller values.

Answer 2

I am afraid I can't be of much help regarding the QM part of your simulation, but I can give some thoughts on the MM part. In short: you should use the cutoff that was used to validate your force field. While a proper ab initio method is expected to improve with a longer cutoff, this is not necessarily true for force fields. The reason for that is that force fields are inevitably overfitted to either QM data or existing experimental data, meaning that the further you stray from the validated conditions, the less certain you can be in your results. In the original SPC/Fw paper, the cutoff they use for validation is 9 Å, so this is the cutoff I would use for the MM part.

Note that this may not be necessarily true with a complicating QM method on top, and I am not sure whether a standard practice for a "best cutoff" exists in this case. To me, this sounds like the type of question that warrants its own research project.