Gromacs's virtual_sitesn

[keisuke-yanagisawa]

Hi everyone,

Now I'm searching a way to generate gromacs topology file with virtual_sitesn direction (extra point at center-of-mass of selected atoms).

I have written below code:

import parmed as pmd
from parmed.gromacs import gromacstop, gromacsgro

amber = pmd.load_file("amber.parm7", xyz="amber.rst7") 
# input data does not have any EPs

# add EPs for specific residue ("ZZZ" in this case)
atomtype_VIR = pmd.AtomType("VIR", None, 0, 0)
atomtype_VIR.set_lj_params(0,0)
atomtype_VIR.add_nbfix("VIR", rmin, epsilon) # assumed that rmin and epsilon have been defined
for res in amber.residues:
  if res.name == "ZZZ":
    atom_VIR = pmd.ExtraPoint(name="VIR", type="VIR")
    atom_VIR.atom_type = atomtype_VIR
    amber.add_atom_to_residue(atom_VIR, res)

topdata = gromacstop.GromacsTopologyFile.from_structure(amber)
topdata.write("gromacs.top") # case 1
topdata.write("gromacs.top", combine="all") # case 2

In the case 1, the code can output topology file but res ZZZ does not have VIR atom. On the other hand, case 2 raised an AttributeError:

AttributeError                            Traceback (most recent call last)
<ipython-input-87-70cbaad2419e> in <module>
      2 topdata = gromacstop.GromacsTopologyFile.from_structure(amber)
      3 print(topdata)
----> 4 topdata.write("gromacs.top", combine="all")

~/.local/pyenv/versions/anaconda3-2018.12/envs/md_analysis/lib/python3.7/site-packages/parmed/gromacs/gromacstop.py in write(self, dest, combine, parameters, molfile, itp)
   1659                 GromacsTopologyFile._write_molecule(self, _molfile, 'system',
   1660                                                     params,
-> 1661                                                     parameters == 'inline')
   1662                 if not itp :
   1663                     dest.write('[ system ]\n; Name\n')

~/.local/pyenv/versions/anaconda3-2018.12/envs/md_analysis/lib/python3.7/site-packages/parmed/gromacs/gromacstop.py in _write_molecule(struct, dest, title, params, writeparams)
   2018         # Virtual sites
   2019         if EPs:
-> 2020             ftypes = set(type(a.frame_type) for a in EPs)
   2021             for ftype in ftypes:
   2022                 if ftype is TwoParticleExtraPointFrame:

~/.local/pyenv/versions/anaconda3-2018.12/envs/md_analysis/lib/python3.7/site-packages/parmed/gromacs/gromacstop.py in <genexpr>(.0)
   2018         # Virtual sites
   2019         if EPs:
-> 2020             ftypes = set(type(a.frame_type) for a in EPs)
   2021             for ftype in ftypes:
   2022                 if ftype is TwoParticleExtraPointFrame:

~/.local/pyenv/versions/anaconda3-2018.12/envs/md_analysis/lib/python3.7/site-packages/parmed/topologyobjects.py in frame_type(self)
   1107         if self._frame_type is not None:
   1108             return self._frame_type
-> 1109         if len(self.parent.bonds) == 2:
   1110             mybond = None
   1111             otherbond = None

AttributeError: 'NoneType' object has no attribute 'bonds'

I realized the error says I have to determine the type of frame, but I could not find corresponding frame type to virtual_sitesn in Gromacs.

virtual_sitesn is used by some known tools (such as MARTINI), thus I believe it is possibly helpful enhancement while there may be difficult point to implement. I would be glad if admins discuss it, and furthermore, it implements.

[swails]

I'd be happy to accept any contribution that adds this, but it is beyond my experience to do it quickly.

The only virtual site geometries that I've supported are the ones supported by Amber and used in common all-atom force fields -- specifically those used by TIP[45]P water models.

[keisuke-yanagisawa]

Dear swails,

Thank you for your comments. According to the Gromacs documentation, three types of virtual_sitesn can be defined:

• COG (center-of-geometry): atoms are equally weighted regardless the mass of atoms
• COM (center-of-mass): atoms are weighted by mass of atoms
• COW (center-of-weights): atoms are weighted by manually defined values

MARTINI uses COM and I haven't seen COW yet, thus I believe the priority order is COM > COG >>>> COW if you kindly will implement it.

Probably below documentation pages are helpful: Types of virtual_sitesn: http://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#tab-topfile2 Notation of virtual_sitesn: http://manual.gromacs.org/current/reference-manual/topologies/particle-type.html?highlight=virtual_site