Lammps ewald parameters for water gas density

Hello,

I am trying to understand how to choose appropriate ewald parameters for gas phase spce water system. My system of interest is a box if flexible SPCE water at 600 K, 1 bar, essentially ideal gas conditions. I am using pppm with precision 1e-6 and real space cutoff of 6 A, and running NPT simulations to collect gas phase densities. I have found that the gas density is affected by ewald parameters. Below are specific questions with regard to selecting these parameters:

  • Why does changing ewald parameters (real space cutoff, precision, etc) affect gas density at ideal conditions? While I understand these parameters affect the PPPM calculation, I expected at ideal conditions the density is controlled by the temperature and the density would be precisely ideal, and the potential does not matter. Instead, I am seeing below ideal gas like densities when I use a pppm precision of 1e-6 and and real space cutoff of 6 A. Upon increasing the precision and increasing the real space cutoff, the gas densities go back up to above ideal gas. Also, I have tried removing charges from the system and using only a cut and shifted LJ potential - this resulted in gas densities that at least are above ideal gas conditions also. So it seems as if the long range electrostatics are inducing the system to be less dense than ideal condition density. Do you have suggestions on why this might happen?

  • Similarly, removing the long range vdw tail correction also affects the gas density. Without the long range tail correction, the gas density also dips below the ideal gas density at this state-point. If the system is at ideal conditions, I would think that the potential does not matter, and the density should be ideal. Why does a long range vdw tail affect gas density at these conditions?

Thank you in advance for your time!

Best,

Jo

I cannot give any specific answers on your questions but have two general comments.

  1. it seems to me that you are underestimating how far dispersion interactions reach. a cutoff of 6 \AA is very short and thus you are truncating quite a bit of interactions. you need to keep in mind that those dispersion interactions are always attractive and hence the integral estimate from the tail corrections will be a significant contribution since there is no cancellation.

  2. i would be careful with using the automated cutoffs through setting ewald convergence for non-condensed phase systems. this is based on some heuristics and may not work well for strongly diluted systems. you may be missing some contributions due to underestimating the kspace cutoff or using a non-optimal alpha parameter. i suggest to make manual convergence tests.

Axel.

Hello,

Thanks for the advice. I realized that the mesh points are assigned at the beginning of the simulation, and my simulation changed in volume quite a bit at the beginning, which would explain the inaccuracy. I ended up manually specifying mesh grid points and the issue was resolved.

I now have a slightly different question. I understand in ewald long range calculations, g_ewald is specified analytically from precision and real space cutoff as follows:
g_ewald = (1.35 - 0.15*log(precision))/cutoff

For PPPM however, it seems g_ewald is a function of mesh points in addition to precision and cutoff. I would like to know how the number of mesh grid points relates to g_ewald. My understanding is that g_ewald controls the weight of real/reciprocal space contributions and controls the width of gaussians, which then dictates the number of k-vectors necessary. And number of mesh points control degree of discretization of charge densities. I would like to understand the relationship between these two parameters (g_ewald and mesh grid points). Practically, I assume only one of these parameters should be specified using kspace_modify?

Thanks,

Jo

It is not just the mesh that has an impact, but also the order of the interpolation polynomial.
For the details, you are probably best off looking through the source code and/or read up in the corresponding papers. My knowledge of the theory is not sufficient to give you more advice beyond this. I am copying Stan, who is currently the resident Kspace expert among the LAMMPS developers; perhaps he has some advice for you.

Axel.

Analytic expressions are used to estimate the errors in PPPM. We use a Newton-Raphson solver to find a good g_ewald value based on these expressions. See https://github.com/lammps/lammps/blob/master/src/KSPACE/pppm.cpp#L1285 and

How to mesh up Ewald sums. II. An accurate error estimate for the particle–particle–particle-mesh algorithm

J. Chem. Phys. 109, 7694 (1998); https://doi.org/10.1063/1.477415

Stan