[gmx-users] Significant differences in Ewald Vs PME for clay sheet

[Gaurav Goel]

A (brief) background:

We have prepared a periodic Na-MMT clay sheet using force-field parameters
as developed by Heinz and co-workers. The box (and clay sheet) are 2.6-2.7
nm in x- and y- directions. Sheet thickness along z is 0.7nm and we add a
vacuum layer to get a box_z=4nm.

On using Ewald summation for electrostatic interactions, we get an exact
correspondence with results of Heinz et al. The equilibrium box-size in
2.6nm in x and y, as reported by Heinz using Ewald summations. Further the
Young's modulus values are in close agreement with other simulation studies
as well as experimental data.

However when we use PME electrostatics, we see the equilibrium box size
increases to 2.7nm. Also, the Young's modulus values are doubled compared
to Ewald summations. This prompted us to dig more into the energy and force
values.

*** ***
As test case, we took a single snapshot of Ewald summation equilibrated
clay sheet configuration (box dimensions: 2.6 x 2.6 x 4 nm) and using
'mdrun -rerun' determined the potential energy, forces and virial
components. Details are below:

Case 1. Ewald (fourierspacing=0.12 ~ 40 wave vectors; ewald_rtol:  1e-5;
rc=1.2nm)
Case 2: PME (fourierspacin=0.08; ewald_rtol: 1e-7; pme_order 8; ewald_3dc;
rc=1.2nm)
We have tried a few other combinations, such as standard Ewald-PME instead
of ewald_3dc, but all the numbers are same as case 2.
                Case 1(Ewald)                case 2 (PME)

Total PE:  -404119                           -404118
Note 1: Individual contributions to the total PE in both cases are
similarly matching.
Note 2: Force versus atom plots as obtained from .trr file is again exactly
matching in the two cases.

However, the various Virial contributions differ significantly:
Vir-xx:   7444.71                             -43303.8
Vir-yy:  -6327.84                             -56858.7

Above numbers clearly indicate a large compressive pressure in case 2
(PME), that is commensurate with observed increase in sheet size along x-
and y- directions in (semi-isotropic) NPT ensemble at zero pressure.

** ** **
Since the virial components are determined from interparticle forces, why
we are seeing this difference in terms of 'force versus atom' data matching
in Ewald versus PME, but the virials showing a significant difference.

We will appreciate any help/suggestions on above issue. Thank you for your
time.

Thanks,
G

-- 
Gaurav Goel, PhD
Assistant Professor
Department of Chemical Engineering
Indian Institute of Technology, Delhi
Hauz Khas, New Delhi 110016
India