Is it possible to dump the coordinates of M sites along with coordinates of water molecules while using tip4p water model?

LAMMPS LAMMPS General Discussion
,
1

Along with the coordinates of O, H and H atoms of water molecules, is there a way of dumping the coordinates of M site as well in the dump file? Seems there isn’t, just wanted to confirm.

Thanks,
Chaitanya

2

I suppose that if you need the M site to be saved in the dump file so very badly to avoid recreating its position at each microstate based on the position of the other atoms in the molecule, you can always create an initial microstate where the dummy atom exists together with the rest of the atoms of the water molecules and naturally assign it a given atom type. Then you set up the force-field by yourself manually and just include the atom type of the dummy atom in the group-ID of the dump command.

3

That would be a very nice way to do it actually. Thanks a lot!

Chaitanya

4

Indeed, if you use the tip4p pair styles in LAMMPS, the M point does not really exist during the simulation. Its position is only calculated when needed based on the positions of H-O-H.

Therefore if what you want is the position of the M point, you can either recalculate it yourself based on H-O-H during post-mortem analysis (for instance using Python-MDAnalysis), or create an explicit 4 points water model as suggested by Cecilia, but you will probably loose in computational efficiency. In that case I think that you will need to assign a small (e.g. 1e-10) mass to the M point.

Best,
Simon

5

I mean, it’s possible, and here is an elaborately over-engineered solution for doing so:

  1. Read your initial data file in with extra/atom/types 1 which we will assign to the charge sites, and set its pair_coeffs to zero and mass to whatever
  2. Use create_atoms random ${typeM} ${nH2O} NULL to register those sites, set type ${typeM} mol v_molnum to register them with the TIP4P molecules, and neigh_modify to remove them from neighbour lists, where molnum is crafted to give them the right molecule numbers
  3. Chunk atoms by molecule using compute chunk/atom. Get the oxygen position for each molecule by masking (such as variable xO atom (type==${typeO})*x and the orientation of each molecule by using compute dipole/chunk. Then propagate all six numbers throughout each molecule using compute spread/chunk/atom, and now each atom can calculate the position of its molecule’s M site (which is the oxygen position plus molecular dipole times some scaling factor).
  4. Create variables masked to the M position for type M particles and to actual position for non-M particles (such as variable xM atom (type==${typeM})*v_xM+(type!=${typeM})*x). Then you can use dump custom to dump type-M and non-type-M positions together, and then use something like sed afterwards to rewrite the column names.

Or, I guess, you can just be boring and use fix rigid with explicit sites. :wink: