Write out rigid water for GROMACS top file

[daico007]

Add support to write out the settles section for gromacs top file. Decision still need to be made about what would be considered the appropriate flag/signal/label to indicate that the water is rigid. Current, the file will check if water in the molecule name and if all of it's site.label is `rigid'.

[codecov]

Codecov Report

Attention: Patch coverage is 85.39326% with 13 lines in your changes are missing coverage. Please review.

Project coverage is 93.99%. Comparing base (069170b) to head (d667404).
Report is 4 commits behind head on main.

Files	Patch %	Lines
gmso/formats/top.py	23.07%	10 Missing ⚠️
gmso/abc/abstract_site.py	96.61%	2 Missing ⚠️
gmso/core/topology.py	91.66%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #771      +/-   ##
==========================================
- Coverage   94.02%   93.99%   -0.03%     
==========================================
  Files          68       68              
  Lines        6859     6947      +88     
==========================================
+ Hits         6449     6530      +81     
- Misses        410      417       +7     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[CalCraven]

I worry a bit about setting sites as rigid. Because what if only some sites are rigid and others are not, but they're all part of the same molecule? What happens in this case?

I'm wondering if a better way to leverage rigid molecules are with the molecules tag we've already set up in GMSO. In that way, the rigid nature can be tagged to the molecules as a whole. Something like:
MoleculeType = NamedTuple("Molecule", name=StrictStr, number=StrictInt, isrigid=StrictBool)

gmso/gmso/abc/abstract_site.py

Line 18 in 615fd29

MoleculeType = NamedTuple("Molecule", name=StrictStr, number=StrictInt)

Then, when making unique_molecules in _get_unique_molecules, you could have it grab the tag as well, and then when iterating through the unique_molecules, you can check
if tag.isrigid: write_settles

[CalCraven]

Can we also parse from_mbuild
to check for particles with rigid bodies. I think something like:
molecule_isrigid = molecule.contains_rigid()

site_map[particle]["molecule"] = (
    molecule_tag.name,
    molecule_number,
    molecule_isrigid,
)

[daico007]

Add check for n_sites==3

[daico007]

A quick note when I'm working on this PR. Since we are trying to make rigidity as part of a molecule, which is a namedtuple, this will make the attribute immutable. Is that really the behavior we want? Since the rigidity can be determined for a set of particles, not necessarily the whole molecule (thinking of our surface system with frozen silica base @CalCraven). I am thinking the rigidity probably make more sense to be defined at site level and independent from the molecule info.

[daico007]

I have to walk back on the previous statement a bit, since I realized I got confused between frozen particle with rigid molecule (with go as molecule), so there's still basis to have rigid as an attribute of molecule, but I think we still need some discussion regarding the immutability of that rigid attribute (as it belong to a namedtuple). I am leaning toward having it as a mutable attribute, but then, where it should be stored is another question that we need to answer.

[CalCraven]

I have to walk back on the previous statement a bit, since I realized I got confused between frozen particle with rigid molecule (with go as molecule), so there's still basis to have rigid as an attribute of molecule, but I think we still need some discussion regarding the immutability of that rigid attribute (as it belong to a namedtuple). I am leaning toward having it as a mutable attribute, but then, where it should be stored is another question that we need to answer.

I've bumped into the immutability myself as an issue when trying to modify the molecule.number for using a different numbering system in the LAMMPS writer. A simple workaround is to just assign a new molecule to the site.

from gmso.abc.abstract_site import Molecule
old_molecule = site.molecule # isrigid=False
site.molecule = Molecule(name=old_molecule.name, number=old_molecule.number, isrigid=True)

What do think of that as a viable solution for keeping molecule the way it is mostly?

On a different note, I don't hate the site.isrigid attribute as well. We will just constantly need a method to check a list of sites to all be rigid. I.e.

molList = top.unique_site_labels("molecule", name_only=True)
for mol in molList:
    for site in top.iter_by_molecules(key="name", value="mol"):
        assert site.isrigid

Generally seems feasible, just a bit funky overall.

[CalCraven]

Looks like this is the error.

def test_top(self, spce_water):
  >       spce_water.save("spce_restraint.top", overwrite=True)
  
  /Users/runner/work/gmso/gmso/gmso/tests/test_top.py:24: 
  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
  /Users/runner/work/gmso/gmso/gmso/core/topology.py:1626: in save
      saver(self, filename, **kwargs)
  /Users/runner/work/gmso/gmso/gmso/formats/top.py:181: in write_top
      if "water" in tag.lower() and all(
  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
  
  .0 = <list_iterator object at 0x16028ccd0>
  
      if "water" in tag.lower() and all(
  >       site.molecule.isrigid for site in unique_molecules[tag]["sites"]
      ):
  E   AttributeError: 'Molecule' object has no attribute 'isrigid'
  
  /Users/runner/work/gmso/gmso/gmso/formats/top.py:182: AttributeError

Looks like the test is there, but the namedTuple doesn't have the isrigid attribute, probably due to the back and forth discussed above.

In a standard mosdef workflow, do we want rigidity to be set in GMSO, in the forcefield, or both? IMO, the best practice would be in the forcefield. However, I can see reason to separate it from there for ease of access, and because it would cause further functional API changes to parsing the Forcefield files.

With that in mind, I think we should have a setter that can be called on a set of sites that make up a molecule. Rigidity has to be on a per molecule basis. You can't have a partially rigid molecule, the way you can have a partially constrained molecule. At some point we're going to have to have both. Restraints should be a settable attribute on the connection or connection_type level. Constraints should be set on the molecule level.

In terms of implementing this, it may be time to graduate from the named tuple for molecules and groups. I think we should make a class for them with setters, getters, and the relevant attributes. @daico007 what are your thoughts on this?

[daico007]

Given how much we are using site.molecule, I agree that we should have a more dedicated class for those attribute (site.group, site.molecule, and site.residue). My only hesitant would be how much work will require to implement those.

Regarding the information regarding the rigidity information, my preference would be to have such info stored in the forcefield, but also give the user the ability add in additional constraint on the GMSO side. Also agree on the point that the rigidity should be set to at least the molecule level.

I also want to note that this PR is dealing with a very specific case (Water) which has a very specific handling in GROMACS (writing out its own [ settles ] section.

[bc118]

This issues is directly related to a new issue I opened. #817 .

Therefore please see this for more information, as part of the non-element errors source is found here.