Skip to content

Commit

Permalink
add basic rigid body snapshot test
Browse files Browse the repository at this point in the history
  • Loading branch information
chrisjonesBSU committed Sep 23, 2024
1 parent ad603fa commit fdba315
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 94 deletions.
26 changes: 14 additions & 12 deletions gmso/external/convert_hoomd.py
Original file line number Diff line number Diff line change
Expand Up @@ -284,19 +284,19 @@ def _parse_particle_information(
site.name if site.atom_type is None else site.atom_type.name
for site in top.sites
]
masses = list()
charges = list()
for site in top.sites:
masses.append(
masses = np.zeros(len(xyz))
charges = np.zeros(len(xyz))
for idx, site in enumerate(top.sites):
masses[idx] = (
site.mass.to_value(base_units["mass"])
if site.mass
else 1 * base_units["mass"]
)
charges.append(site.charge if site.charge else 0 * u.elementary_charge)
charges[idx] = (site.charge if site.charge else 0 * u.elementary_charge)
# Check for rigid IDs
rigid_ids = [site.rigid_id for site in top.sites]
rigid_ids_set = set(rigid_ids)
if all(rigid_ids): # Need to create rigid bodies
if None not in rigid_ids:
n_rigid = len(rigid_ids_set)
write_rigid = True
else:
Expand All @@ -306,12 +306,14 @@ def _parse_particle_information(
typeids = np.array([unique_types.index(t) for t in types])
if write_rigid:
rigid_masses = np.zeros(n_rigid)
rigid_xyz = np.zeros(n_rigid)
rigid_xyz = np.zeros((n_rigid, 3))
# Rigid particle type defaults to "R"; add to front of list
# TODO: Can we always use "R" here? What if an atom_type is "R"?
unique_types = ["R"] + unique_types
# Rigid particles get type ID 0, move all others up by 1
typeids = np.concatenate(([0] * n_rigid), typeids + 1)
typeids = np.concatenate(
(np.array([0] * n_rigid), typeids + 1)
)
# Update mass list and position list of Frame
for idx, _id in enumerate(rigid_ids_set):
group_indices = np.where(rigid_ids == _id)[0]
Expand All @@ -323,9 +325,9 @@ def _parse_particle_information(
rigid_masses[idx] = sum(group_masses)
rigid_xyz[idx] = com_xyz
# Append rigid center mass and xyz to front
masses = np.concatenate(rigid_masses, masses)
xyz = np.concatenate(rigid_xyz, xyz)
rigid_id_tags = np.concatenate(np.arange(n_rigid), np.array(rigid_ids))
masses = np.concatenate((rigid_masses, masses))
xyz = np.concatenate((rigid_xyz, xyz))
rigid_id_tags = np.concatenate((np.arange(n_rigid), np.array(rigid_ids)))

"""
Permittivity of free space = 2.39725e-4 e^2/((kcal/mol)(angstrom)),
Expand All @@ -349,7 +351,7 @@ def _parse_particle_information(
if write_rigid:
snapshot.particles.body[0:] = rigid_id_tags
elif isinstance(snapshot, gsd.hoomd.Frame):
snapshot.particles.N = top.n_sites
snapshot.particles.N = top.n_sites + n_rigid
snapshot.particles.types = unique_types
snapshot.particles.position = xyz
snapshot.particles.typeid = typeids
Expand Down
105 changes: 23 additions & 82 deletions gmso/tests/test_hoomd.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,9 @@
import numpy as np
import pytest
import unyt as u
from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield

from gmso import ForceField
from gmso.external import from_mbuild
from gmso.external.convert_hoomd import to_hoomd_forcefield, to_hoomd_snapshot
from gmso.external.convert_hoomd import to_hoomd_forcefield, to_hoomd_snapshot, to_gsd_snapshot
from gmso.parameterization import apply
from gmso.tests.base_test import BaseTest
from gmso.tests.utils import get_path
Expand Down Expand Up @@ -53,88 +51,31 @@ def run_hoomd_nvt(snapshot, forces, vhoomd=4):
return sim


@pytest.mark.skipif(not has_hoomd, reason="hoomd is not installed")
@pytest.mark.skipif(not has_mbuild, reason="mbuild not installed")
class TestGsd(BaseTest):
def test_mbuild_comparison(self, oplsaa_forcefield):
compound = mb.load("CCC", smiles=True)
com_box = mb.packing.fill_box(compound, box=[5, 5, 5], n_compounds=2)
base_units = {
"mass": u.g / u.mol,
"length": u.nm,
"energy": u.kJ / u.mol,
}

top = from_mbuild(com_box)
top.identify_connections()
top = apply(top, oplsaa_forcefield, remove_untyped=True)

gmso_snapshot, _ = to_hoomd_snapshot(top, base_units=base_units)
gmso_forces, _ = to_hoomd_forcefield(
top,
r_cut=1.4,
base_units=base_units,
pppm_kwargs={"resolution": (64, 64, 64), "order": 7},
)

integrator_forces = list()
for cat in gmso_forces:
for force in gmso_forces[cat]:
integrator_forces.append(force)

import foyer

oplsaa = foyer.Forcefield(name="oplsaa")
structure = oplsaa.apply(com_box)

d = 10
e = 1 / 4.184
m = 0.9999938574

mb_snapshot, mb_forcefield, _ = create_hoomd_forcefield(
structure,
ref_distance=d,
ref_energy=e,
ref_mass=m,
r_cut=1.4,
init_snap=None,
pppm_kwargs={"Nx": 64, "Ny": 64, "Nz": 64, "order": 7},
)

assert mb_snapshot.particles.N == gmso_snapshot.particles.N
assert np.allclose(
mb_snapshot.particles.position, gmso_snapshot.particles.position
)
assert mb_snapshot.bonds.N == gmso_snapshot.bonds.N
assert mb_snapshot.angles.N == gmso_snapshot.angles.N
assert mb_snapshot.dihedrals.N == gmso_snapshot.dihedrals.N

sorted_gmso_ff = sorted(integrator_forces, key=lambda cls: str(cls.__class__))
sorted_mbuild_ff = sorted(mb_forcefield, key=lambda cls: str(cls.__class__))

for mb_force, gmso_force in zip(sorted_mbuild_ff, sorted_gmso_ff):
if isinstance(
mb_force,
(
hoomd.md.long_range.pppm.Coulomb,
hoomd.md.pair.pair.LJ,
hoomd.md.special_pair.LJ,
hoomd.md.pair.pair.Ewald,
hoomd.md.special_pair.Coulomb,
),
):
continue
keys = mb_force.params.param_dict.keys()
for key in keys:
gmso_key = key.replace("opls_135", "CT")
gmso_key = gmso_key.replace("opls_136", "CT")
gmso_key = gmso_key.replace("opls_140", "HC")
gmso_key = "-".join(sort_connection_strings(gmso_key.split("-")))
mb_params = mb_force.params.param_dict[key]
gmso_params = gmso_force.params.param_dict[gmso_key]
variables = mb_params.keys()
for var in variables:
assert np.isclose(mb_params[var], gmso_params[var])
@pytest.mark.skipif(not has_hoomd, reason="hoomd is not installed")
@pytest.mark.skipif(not has_mbuild, reason="mbuild not installed")
class TestGsd(BaseTest):
def test_rigid_bodies(self):
ethane = mb.lib.molecules.Ethane()
box = mb.fill_box(ethane, n_compounds=2, box=[2, 2, 2])
box_no_rigid = mb.clone(box)
for i, child in enumerate(box.children):
for p in child.particles():
p.rigid_id = i

top = from_mbuild(box)
top_no_rigid = from_mbuild(box_no_rigid)

rigid_ids = [site.rigid_id for site in top.sites]
assert len(rigid_ids) == box.n_particles
assert len(set(rigid_ids)) == 2

snapshot, refs = to_gsd_snapshot(top)
snapshot_no_rigid, refs = to_gsd_snapshot(top_no_rigid)
assert "R" in snapshot.particles.types
assert snapshot.particles.N - 2 == snapshot_no_rigid.particles.N

@pytest.mark.skipif(
int(hoomd_version[0]) < 4, reason="Unsupported features in HOOMD 3"
Expand Down

0 comments on commit fdba315

Please sign in to comment.