Skip to content

Commit

Permalink
Generate dihedral_types in to_hoomd_forces by dihedral, not dihedral_…
Browse files Browse the repository at this point in the history
…types
  • Loading branch information
CalCraven committed Sep 11, 2023
1 parent 37fdf57 commit adebaaa
Show file tree
Hide file tree
Showing 2 changed files with 65 additions and 92 deletions.
93 changes: 63 additions & 30 deletions gmso/external/convert_hoomd.py
Original file line number Diff line number Diff line change
Expand Up @@ -1096,27 +1096,27 @@ def _parse_dihedral_forces(
[site.atom_type.atomclass for site in dihedral.connection_members]
)
unique_dihedrals[unique_members] = dihedral

unique_dtypes = [

Check notice

Code scanning / CodeQL

Unused local variable Note library

Variable unique_dtypes is not used.
dihedral.dihedral_type for dihedral in unique_dihedrals.values()
]
groups = dict()
for dtype in unique_dtypes:
group = potential_types[dtype]
for dihedral in unique_dihedrals.values():
group = potential_types[dihedral.dihedral_type]
if group not in groups:
groups[group] = [dtype]
groups[group] = [dihedral]
else:
groups[group].append(dtype)
groups[group].append(dihedral)

expected_unitsDict = {}
for group in groups:
expected_units_dim = potential_refs[group][
expected_unitsDict[group] = potential_refs[group][
"expected_parameters_dimensions"
]
groups[group] = _convert_params_units(
groups[group],
expected_units_dim,
base_units,
)
# groups[group] = _convert_connection_params_units(
# groups[group],
# expected_units_dim,
# base_units,
# )
dtype_group_map = {
"OPLSTorsionPotential": {
"container": hoomd.md.dihedral.OPLS,
Expand Down Expand Up @@ -1152,19 +1152,24 @@ def _parse_dihedral_forces(
dihedral_forces.append(
dtype_group_map[group]["parser"](
container=dtype_group_map[group]["container"](),
dtypes=groups[group],
dihedrals=groups[group],
expected_units_dim=expected_unitsDict[group],
base_units=base_units,
)
)
return dihedral_forces


def _parse_periodic_dihedral(
container,
dtypes,
container, dihedrals, expected_units_dim, base_units
):
for dtype in dtypes:
member_classes = sort_by_classes(dtype)
container.params["-".join(member_classes)] = {
for dihedral in dihedrals:
dtype = dihedral.dihedral_type
dtype = _convert_single_param_units(

Check warning on line 1168 in gmso/external/convert_hoomd.py

View check run for this annotation

Codecov / codecov/patch

gmso/external/convert_hoomd.py#L1166-L1168

Added lines #L1166 - L1168 were not covered by tests
dtype, expected_units_dim, base_units
)
member_sites = sort_connection_members(dihedral, "atomclass")
container.params["-".join(member_sites)] = {

Check warning on line 1172 in gmso/external/convert_hoomd.py

View check run for this annotation

Codecov / codecov/patch

gmso/external/convert_hoomd.py#L1171-L1172

Added lines #L1171 - L1172 were not covered by tests
"k": dtype.parameters["k"],
"d": 1,
"n": dtype.parameters["n"],
Expand All @@ -1173,15 +1178,17 @@ def _parse_periodic_dihedral(
return container


def _parse_opls_dihedral(
container,
dtypes,
):
for dtype in dtypes:
def _parse_opls_dihedral(container, dihedrals, expected_units_dim, base_units):
for dihedral in dihedrals:
dtype = dihedral.dihedral_type
dtype = _convert_single_param_units(

Check warning on line 1184 in gmso/external/convert_hoomd.py

View check run for this annotation

Codecov / codecov/patch

gmso/external/convert_hoomd.py#L1182-L1184

Added lines #L1182 - L1184 were not covered by tests
dtype, expected_units_dim, base_units
)
member_sites = sort_connection_members(dihedral, "atomclass")

Check warning on line 1187 in gmso/external/convert_hoomd.py

View check run for this annotation

Codecov / codecov/patch

gmso/external/convert_hoomd.py#L1187

Added line #L1187 was not covered by tests
# TODO: The range of ks is mismatched (GMSO go from k0 to k5)
# May need to do a check that k0 == k5 == 0 or raise a warning
member_classes = sort_by_classes(dtype)

Check notice

Code scanning / CodeQL

Unused local variable Note library

Variable member_classes is not used.
container.params["-".join(member_classes)] = {
container.params["-".join(member_sites)] = {

Check warning on line 1191 in gmso/external/convert_hoomd.py

View check run for this annotation

Codecov / codecov/patch

gmso/external/convert_hoomd.py#L1190-L1191

Added lines #L1190 - L1191 were not covered by tests
"k1": dtype.parameters["k1"],
"k2": dtype.parameters["k2"],
"k3": dtype.parameters["k3"],
Expand All @@ -1190,19 +1197,22 @@ def _parse_opls_dihedral(
return container


def _parse_rb_dihedral(
container,
dtypes,
):
def _parse_rb_dihedral(container, dihedrals, expected_units_dim, base_units):
warnings.warn(
"RyckaertBellemansTorsionPotential will be converted to OPLSTorsionPotential."
)
for dtype in dtypes:
for dihedral in dihedrals:
dtype = dihedral.dihedral_type
dtype = _convert_single_param_units(
dtype, expected_units_dim, base_units
)
opls = convert_ryckaert_to_opls(dtype)
member_classes = sort_by_classes(dtype)
member_sites = sort_connection_members(dihedral, "atomclass")
# TODO: The range of ks is mismatched (GMSO go from k0 to k5)
# May need to do a check that k0 == k5 == 0 or raise a warning
container.params["-".join(member_classes)] = {
container.params[
"-".join([site.atom_type.atomclass for site in member_sites])
] = {
"k1": opls.parameters["k1"],
"k2": opls.parameters["k2"],
"k3": opls.parameters["k3"],
Expand Down Expand Up @@ -1420,3 +1430,26 @@ def _convert_params_units(
potential.parameters = converted_params
converted_potentials.append(potential)
return converted_potentials


def _convert_single_param_units(
potential,
expected_units_dim,
base_units,
):
"""Convert parameters' units in the potential to that specified in the base_units."""
converted_params = dict()
for parameter in potential.parameters:
unit_dim = expected_units_dim[parameter]
ind_units = re.sub("[^a-zA-Z]+", " ", unit_dim).split()
for unit in ind_units:
unit_dim = unit_dim.replace(
unit,
f"({str(base_units[unit].value)} * {str(base_units[unit].units)})",
)

converted_params[parameter] = potential.parameters[parameter].to(
unit_dim
)
potential.parameters = converted_params
return potential
64 changes: 2 additions & 62 deletions gmso/tests/test_hoomd.py
Original file line number Diff line number Diff line change
Expand Up @@ -456,64 +456,6 @@ def test_zero_charges(self):
assert not isinstance(force, hoomd.md.special_pair.Coulomb)

def test_forces_connections_match(self):
compound = mb.load("CC", smiles=True)
com_box = mb.packing.fill_box(compound, box=[5, 5, 5], n_compounds=2)
base_units = {
"mass": u.amu,
"length": u.nm,
"energy": u.kJ / u.mol,
}
top = com_box.to_gmso()
top.identify_connections()
ethaneFF = ForceField(get_path("alkanes.xml"))
ethaneFF.atom_types["opls_01"] = ethaneFF.atom_types.pop("opls_140")
ethaneFF.atom_types["opls_01"].name = "opls_01"
ethaneFF.atom_types["opls_01"].atomclass = "opls_01"
ethaneFF.atom_types["opls_1004"] = ethaneFF.atom_types.pop("opls_135")
ethaneFF.atom_types["opls_1004"].name = "opls_1004"
ethaneFF.atom_types["opls_1004"].atomclass = "opls_1004"
xDict = {"bond": {}, "angle": {}, "dihedral": {}}
for dictKey in xDict:
for connection in getattr(ethaneFF, dictKey + "_types"):
newname = connection
for atomclass, atomname in {
"HC": "opls_01",
"CT": "opls_1004",
}.items():
newname = newname.replace(atomclass, atomname)
xDict[dictKey][connection] = newname

for dictKey in xDict:
for oldname, newname in xDict[dictKey].items():
getattr(ethaneFF, dictKey + "_types")[newname] = getattr(
ethaneFF, dictKey + "_types"
).pop(oldname)
getattr(ethaneFF, dictKey + "_types")[newname].name = newname
getattr(ethaneFF, dictKey + "_types")[
newname
].member_types = tuple(newname.split("~"))

# ethaneFF.bond_types["opls_01~opls_1004"] = ethaneFF.bond_types.pop("CT~HC")
# ethaneFF.bond_types["opls_01~opls_01"] = ethaneFF.bond_types.pop("CT~CT")
# ethaneFF.angle_types["opls_01~opls_01~opls_1004"] = ethaneFF.angle_types.pop("CT~CT~HC")
# ethaneFF.angle_types["opls_1004~opls_01~opls_1004"] = ethaneFF.angle_types.pop("HC~CT~HC")
# ethaneFF.dihedral_types["opls_1004~opls_01~opls_01~opls_1004"] = ethaneFF.dihedral_types.pop("HC~CT~CT~HC")
# should sort these opls_01, opls_1004
top = apply(top, ethaneFF, remove_untyped=True)

snapshot, snapshot_base_units = to_hoomd_snapshot(
top, base_units=base_units
)

forces, forces_base_units = to_hoomd_forcefield(
top=top, r_cut=1.4, base_units=base_units
)

assert forces_base_units == snapshot_base_units
for conntype in snapshot.bonds.types:
assert conntype in list(forces["bonds"][0].params.keys())

def test_forces_connections_match2(self):
compound = mb.load("CC", smiles=True)
com_box = mb.packing.fill_box(compound, box=[5, 5, 5], n_compounds=1)
base_units = {
Expand All @@ -527,12 +469,10 @@ def test_forces_connections_match2(self):

top = apply(top, ethaneFF, remove_untyped=True)

snapshot, snapshot_base_units = to_hoomd_snapshot(
top, base_units=base_units
)
snapshot, _ = to_hoomd_snapshot(top, base_units=base_units)
assert "CT-HC" in snapshot.bonds.types

forces, forces_base_units = to_hoomd_forcefield(
forces, _ = to_hoomd_forcefield(
top=top, r_cut=1.4, base_units=base_units
)
assert "CT-HC" in forces["bonds"][0].params.keys()
Expand Down

0 comments on commit adebaaa

Please sign in to comment.