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 94.84536% with 5 lines in your changes missing coverage. Please review.

Project coverage is 94.08%. Comparing base (5b3040a) to head (862db26).

❗ Current head 862db26 differs from pull request most recent head bde8b6d

Please upload reports for the commit bde8b6d to get more accurate results.

Files	Patch %	Lines
gmso/abc/abstract_site.py	96.61%	2 Missing ⚠️
gmso/core/element.py	60.00%	2 Missing ⚠️
gmso/core/topology.py	94.11%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #771      +/-   ##
==========================================
- Coverage   94.09%   94.08%   -0.01%     
==========================================
  Files          65       65              
  Lines        6870     6950      +80     
==========================================
+ Hits         6464     6539      +75     
- Misses        406      411       +5     

☔ 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.

[CalCraven]

Looks pretty much good. The only thing I can think to add is maybe a utility to set the molecules to be rigid.

def set_rigid(self, molecule_name=None, molecule_number=None):
"""Set molecule tags to rigid if they match the name or number specified."""
from itertools import compress
all_molecules = all_molecules = [site.molecule for site in ptop.sites] # ptop.unique_site_labels() doesn't work, since each site gets a unique molecule in memory. Might be reasons to make actually only generate one object for each unique Molecule and Residue
if molecule_name:
mask_name = [molecule_name == mol.name for mol in all_molecules]
else:
mask_name = [True] * len(all_molecules)
if molecule_number:
mask_number = [molecule_number == mol.number for mol in all_molecules]
else:
mask_number = [True] * len(all_molecules)
mask_all = np.all((mask_name, mask_number), axis=0)
molecules = list(compress(all_molecules, mask_all))
for mol in molecules:
mol.isrigid = True

[CalCraven]

Suggested change

	else:
	tmp = containers_dict[key](
	name=value[0], number=value[1], isrigid=value[2]
	)
	elif len(value) == 3:
	tmp = containers_dict[key](
	name=value[0], number=value[1], isrigid=value[2]
	)
	else:
	raise ValueError(f"Argument value was passed as {value=}, but should be an indexible iterable of [name, number, isrigid] where name is type string, number is type int, and isrigid is type bool.")

[CalCraven]

Also, do we need to raise an additional error for trying to use isrigid argument in value for Residues, since we're not making residues allowed to be rigid?

[daico007]

Added the code above. Having more check would be good, but it's becoming a bit clumsy. I might need to rethink a better strategy to do these checking (open to suggestion).

[CalCraven]

Also to add, we need a test in test_top.py for the settles section. Doesn't look like it's there.

[CalCraven]

Testing it locally, the settles section is leaving in the unyts. So that needs to get removed I think.
[ settles ] ;Water specific constraint algorithm
; OW_idx funct doh dhh
1 1 0.09903 nm 0.15661 nm

[CalCraven]

One minor typo, then it looks good to me.