Mismatch in number of surface atoms

[ParthGoyal1508]

Hi Marcello,

We were trying to find the surface atoms of a water slab of size 555 nm^3 inside a vacuum box of dimension 5560 nm^3. Using different configuration files for the same simulation time step (i.e using gro, xtc, trr) we are getting different number of interfacial (i.e water-vapour interface) atoms.

Using .gro and .xtc configuration file we are getting this output : [array([[<AtomGroup with 1032 atoms>], [<AtomGroup with 1068 atoms>]], dtype=object)]

while for .trr configuration file we are getting this output: [array([[<AtomGroup with 1036 atoms>], [<AtomGroup with 1068 atoms>]], dtype=object)]

---

Code Used:

import MDAnalysis as mda import pytim u = mda.Universe("t100000.tpr","11.trr") inter = pytim.ITIM(u,max_layers=1, cluster_cut=3.5) print (inter.layers)

---

similar for rest of the configuration files.

Files: https://github.com/ParthGoyal1508/Files

[Marcello-Sega]

Hi,

thanks for reporting this, I will have a look at it soon

[Marcello-Sega]

Hi,

this doesn't seem to be neither a pytim nor a mdanalysis problem, but just a difference in your gromacs files.

I get the same output as you do:

import MDAnalysis as mda
import pytim
import numpy as np
u1 = mda.Universe("t100000.tpr","11.trr")
inter1 = pytim.ITIM(u1,max_layers=1, cluster_cut=3.5)
print (inter1.layers)

u2 = mda.Universe("t100000.tpr","./11.gro")
inter2 = pytim.ITIM(u2,max_layers=1, cluster_cut=3.5)
print (inter2.layers)

[[<AtomGroup with 1036 atoms>]
 [<AtomGroup with 1068 atoms>]]
[[<AtomGroup with 1032 atoms>]
 [<AtomGroup with 1068 atoms>]]

And the culprit is one water molecule (first atom id 2801, or index 2800)

cond = np.isin(inter1.layers[0,0].atoms,inter2.layers[0,0].atoms)
diff = inter1.layers[0,0].atoms[~cond]
print(diff)

<AtomGroup [<Atom 2801: OW of type opls_113 of resname SOL, resid 700 and segid seg_0_SOL and altLoc  >, <Atom 2802: HW1 of type opls_114 of resname SOL, resid 700 and segid seg_0_SOL and altLoc  >, <Atom 2803: HW2 of type opls_114 of resname SOL, resid 700 and segid seg_0_SOL and altLoc  >, <Atom 2804: MW of type opls_115 of resname SOL, resid 700 and segid seg_0_SOL and altLoc  >]>

In fact, the positions read from the trr are different than from those read from the gro

print(diff.positions)

[[  7.435667   48.917595  319.70837  ]
 [  6.6767006  49.341118  319.30737  ]
 [  7.586691   48.138733  319.17285  ]
 [  7.357843   48.87211   319.5885   ]]

u2.atoms[diff.atoms.ids].positions
array([[  7.44    ,  48.920002, 319.71002 ],
       [  6.68    ,  49.34    , 319.31    ],
       [  7.59    ,  48.14    , 319.16998 ],
       [  7.36    ,  48.870003, 319.59    ]], dtype=float32)

The first things that comes to mind is the precision difference between .gro (and .xtc) and the full float precision of the .trr

Inspecting directly the .gro file one sees that indeed the coordinates are read correctly

  701SOL     OW 2801   0.744   4.892  31.971  0.7084 -0.8023  0.6354
  701SOL    HW1 2802   0.668   4.934  31.931  3.0810  2.0216 -0.8727
  701SOL    HW2 2803   0.759   4.814  31.917 -0.3298 -0.6751  0.1576
  701SOL     MW 2804   0.736   4.887  31.959  0.0000 -0.0000 -0.0000

and a quick run of gmx dump on both the .trr and the .xtc:

      x[ 2800]={ 7.43567e-01,  4.89176e+00,  3.19708e+01}
      x[ 2801]={ 6.67670e-01,  4.93411e+00,  3.19307e+01}
      x[ 2802]={ 7.58669e-01,  4.81387e+00,  3.19173e+01}
      x[ 2803]={ 7.35784e-01,  4.88721e+00,  3.19589e+01}

      x[ 2800]={ 7.44000e-01,  4.89200e+00,  3.19710e+01}
      x[ 2801]={ 6.68000e-01,  4.93400e+00,  3.19310e+01}
      x[ 2802]={ 7.59000e-01,  4.81400e+00,  3.19170e+01}
      x[ 2803]={ 7.36000e-01,  4.88700e+00,  3.19590e+01}

shows that the precision with which the positions are stored (1e-3) is the issue.

[Marcello-Sega]

My previous (now deleted) comment about the constraints was just wrong, it's a plain precision issue.

Of course, needless to say, different atomic positions will lead to different surface atoms

[ParthGoyal1508]

Yeah! The issue is with the precision only. Thanks for the help.