mosdef-hub / gmso

Flexible storage of chemical topology for molecular simulation
https://gmso.mosdef.org
MIT License
53 stars 32 forks source link

Add support for `hoomd.md.improper.Periodic` in `to_hoomd_forcefield` #807

Closed chrisjonesBSU closed 2 months ago

chrisjonesBSU commented 4 months ago

This PR starts adding support for periodic improper potentials as they were recently added to Hoomd version 4.5.

I'm running into some issues, I figured I'd get the PR submitted so there are some other eyes on it as well. Here is a gist of the notebook I'm testing it out with..

One note: I made some changes to the environment file because I was having issues with mBuild plugins using mBuild < 17. This is a separate issue, but the pinned versions in environment-dev.yml are changed here to fix that for now.

The contents of the notebook is basically the unit test we have for using GAFF, but with benzene instead of ethanol.

import mbuild as mb
import unyt as u
import gmso
from gmso.external import from_mbuild, to_hoomd_forcefield, to_hoomd_snapshot
from gmso.core.forcefield import ForceField
from gmso.parameterization import apply
from gmso.utils.io import get_fn

benzene = mb.load("c1ccccc1", smiles=True)
benzene.box = mb.Box([5, 5, 5,])
top = from_mbuild(benzene)

base_units = dict(mass=u.g/u.mol, length=u.nm, energy=u.kJ/u.mol)
gaff = ForceField(get_fn("gmso_xmls/test_ffstyles/gaff.xml"))
paramaterized_top = apply(top=top, forcefields=gaff, identify_connections=True)
snap, _ = to_hoomd_snapshot(top, base_units=base_units)
forces, _ = to_hoomd_forcefield(top, r_cut=1.4, base_units=base_units)

It looks like hoomd improper force params are being written with the wild cards?

>>> improper_force = forces["impropers"][0]
>>> for param in improper_force.params:
>>>    print(param)
>>>    print(improper_force.params[param])

ca-*-*-ha
_HOOMDDict{'k': 4.6024, 'n': 2, 'd': 1, 'chi0': 3.141592653589793}

This doesn't match the snapshot improper types

>>> print(snap.impropers.types)
['ca-ca-ca-ha']

I'll keep working on this, but wanted to get it out there for review; maybe I'm missing something obvious at the moment.

chrisjonesBSU commented 3 months ago

So, there is kind of a band-aid fix in here that I think is sort of related to #767, but when the forcefield uses wild cards, these carry over into member_classes, so the force entry to hoomd looks something like [ca, *, *, ha] while the snapshot's particle types are [ca, ca, ca, ha].

Right now, the snapshot writer uses sort_connection_members while the forcefield writer(s) use sort_by_classes

The band-aid fix is to check for * in the forcefield members, and if it's found, use sort_by_types instead.

IMO, I think this is fine for now, but there is still the overall larger issue of having more flexibility in both giving instructions on what to use in the writers (types vs classes) as well as flexibility in combining them (as is discussed in #767)

chrisjonesBSU commented 3 months ago

This should be ready to merge, the only test failing is test_charmm_improper_ff which is unrelated to this PR.

CalCraven commented 2 months ago

Raise and issue about python < 3.12 though once this is merged, since we'll want to fix that when we get to that issue.