From 3a27d99976fd137238b053fc5a5111eb69253ab2 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 17 Sep 2024 13:40:32 -0600 Subject: [PATCH 01/19] add rigid_body field to atom, populate in from_mbuild --- gmso/core/atom.py | 11 +++++++++++ gmso/external/convert_mbuild.py | 1 + 2 files changed, 12 insertions(+) diff --git a/gmso/core/atom.py b/gmso/core/atom.py index dfdcc2de..46d4e559 100644 --- a/gmso/core/atom.py +++ b/gmso/core/atom.py @@ -73,6 +73,12 @@ class Atom(Site): alias="restraint", ) + rigid_id_: Optional[int] = Field( + None, + description="Rigid body ID used to assign rigid body contraints in HOOMD-Blue topologies.", + alias="rigid_id", + ) + model_config = ConfigDict( alias_to_fields=dict( **Site.model_config["alias_to_fields"], @@ -125,6 +131,11 @@ def restraint(self): """Return the restraint of this atom.""" return self.__dict__.get("restraint_") + @property + def rigid_id(self): + """Return the rigid body ID for this atom.""" + return self.__dict__.get("rigid_id_") + @field_serializer("charge_") def serialize_charge(self, charge_: Union[u.unyt_quantity, None]): if charge_ is None: diff --git a/gmso/external/convert_mbuild.py b/gmso/external/convert_mbuild.py index 5ac724ce..fc6323a9 100644 --- a/gmso/external/convert_mbuild.py +++ b/gmso/external/convert_mbuild.py @@ -272,6 +272,7 @@ def _parse_site(site_map, particle, search_method, infer_element=False): element=ele, charge=charge, mass=mass, + rigid_id=particle.rigid_id, molecule=site_map[particle]["molecule"], residue=site_map[particle]["residue"], group=site_map[particle]["group"], From bbb65c3eacdc2177fe834225aca7d84aee739c02 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 17 Sep 2024 14:06:15 -0600 Subject: [PATCH 02/19] move rigid_id property to abstract_site from atom, write out rigid_id in to_mbuild --- gmso/abc/abstract_site.py | 11 +++++++++++ gmso/core/atom.py | 11 ----------- gmso/external/convert_mbuild.py | 1 + 3 files changed, 12 insertions(+), 11 deletions(-) diff --git a/gmso/abc/abstract_site.py b/gmso/abc/abstract_site.py index 399a2b6d..abbc6862 100644 --- a/gmso/abc/abstract_site.py +++ b/gmso/abc/abstract_site.py @@ -209,6 +209,12 @@ class Site(GMSOBase): alias="position", ) + rigid_id_: Optional[int] = Field( + None, + description="Rigid body ID used to assign rigid body contraints in HOOMD-Blue topologies.", + alias="rigid_id", + ) + model_config = ConfigDict( alias_to_fields={ "name": "name_", @@ -250,6 +256,11 @@ def residue(self): """Return the residue assigned to the site.""" return self.__dict__.get("residue_") + @property + def rigid_id(self): + """Return the rigid body ID for this atom.""" + return self.__dict__.get("rigid_id_") + @field_serializer("position_") def serialize_position(self, position_: PositionType): return unyt_to_dict(position_) diff --git a/gmso/core/atom.py b/gmso/core/atom.py index 46d4e559..dfdcc2de 100644 --- a/gmso/core/atom.py +++ b/gmso/core/atom.py @@ -73,12 +73,6 @@ class Atom(Site): alias="restraint", ) - rigid_id_: Optional[int] = Field( - None, - description="Rigid body ID used to assign rigid body contraints in HOOMD-Blue topologies.", - alias="rigid_id", - ) - model_config = ConfigDict( alias_to_fields=dict( **Site.model_config["alias_to_fields"], @@ -131,11 +125,6 @@ def restraint(self): """Return the restraint of this atom.""" return self.__dict__.get("restraint_") - @property - def rigid_id(self): - """Return the rigid body ID for this atom.""" - return self.__dict__.get("rigid_id_") - @field_serializer("charge_") def serialize_charge(self, charge_: Union[u.unyt_quantity, None]): if charge_ is None: diff --git a/gmso/external/convert_mbuild.py b/gmso/external/convert_mbuild.py index fc6323a9..c8ca1e50 100644 --- a/gmso/external/convert_mbuild.py +++ b/gmso/external/convert_mbuild.py @@ -251,6 +251,7 @@ def _parse_particle(particle_map, site): charge=charge, mass=mass, ) + particle.rigid_id = site.rigid_id particle_map[site] = particle return particle From 92534dbe5d0b33098c4a7252d257afd8e00d50c3 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 17 Sep 2024 15:32:52 -0600 Subject: [PATCH 03/19] add rigid types and typeids --- gmso/external/convert_hoomd.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 1eb7b8ce..8a817177 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -284,8 +284,6 @@ def _parse_particle_information( site.name if site.atom_type is None else site.atom_type.name for site in top.sites ] - unique_types = sorted(list(set(types))) - typeids = np.array([unique_types.index(t) for t in types]) masses = list() charges = list() for site in top.sites: @@ -295,6 +293,23 @@ def _parse_particle_information( else 1 * base_units["mass"] ) charges.append(site.charge if site.charge else 0 * u.elementary_charge) + # Check for rigid IDs + rigid_ids = [site.rigid_id for site in top.sites] + if all(rigid_ids): # Need to create rigid bodies + n_rigid = len(set(rigid_ids)) + write_rigid = True + else: + write_rigid = False + n_rigid = 0 + unique_types = sorted(list(set(types))) + typeids = np.array([unique_types.index(t) for t in types]) + if write_rigid: + # Rigid particle type defaults to "R"; add to front of list + 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) + # Update mass list + # Update position list """ Permittivity of free space = 2.39725e-4 e^2/((kcal/mol)(angstrom)), @@ -309,7 +324,7 @@ def _parse_particle_information( ) ** 0.5 if isinstance(snapshot, hoomd.Snapshot): - snapshot.particles.N = top.n_sites + snapshot.particles.N = top.n_sites + n_rigid snapshot.particles.types = unique_types snapshot.particles.position[0:] = xyz snapshot.particles.typeid[0:] = typeids From ab004935cc48d819b0d32ce4dc4af28811fe9e12 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 18 Sep 2024 15:53:14 -0600 Subject: [PATCH 04/19] add test for rigid_id assignment --- gmso/tests/test_convert_mbuild.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/gmso/tests/test_convert_mbuild.py b/gmso/tests/test_convert_mbuild.py index 7b7a9c44..97620756 100644 --- a/gmso/tests/test_convert_mbuild.py +++ b/gmso/tests/test_convert_mbuild.py @@ -259,3 +259,14 @@ def test_nontop_level_compound(self, mb_ethane): cpd.add(mb_ethane) with pytest.raises(AssertionError): from_mbuild(mb_ethane) + + def test_rigid_ids(self, mb_ethane): + box = mb.fill_box(mb_ethane, n_compounds=3, box=[2, 2, 2]) + for i, child in enumerate(box.children): + for p in child.particles(): + p.rigid_id = i + + top = from_mbuild(box) + assert set([site.rigid_id for site in top.sites]) == {0, 1, 2} + for particle, site in zip(box.particles(), top.sites): + assert particle.rigid_id == site.rigid_id From f6480bcabafefbafa2040827e7faede2cb18583f Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 23 Sep 2024 11:27:00 -0600 Subject: [PATCH 05/19] write out rigid xyz, com, tags --- gmso/external/convert_hoomd.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 8a817177..11d9e287 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -295,21 +295,38 @@ def _parse_particle_information( charges.append(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 - n_rigid = len(set(rigid_ids)) + n_rigid = len(rigid_ids_set) write_rigid = True else: write_rigid = False n_rigid = 0 unique_types = sorted(list(set(types))) typeids = np.array([unique_types.index(t) for t in types]) + rigid_id_tags = None if write_rigid: + rigid_masses = np.zeros(n_rigid) + rigid_xyz = np.zeros(n_rigid) # 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) - # Update mass list - # Update position list + # Update mass list and position list of Frame + for idx, _id in enumerate(rigid_ids_set): + group_indices = np.where(rigid_ids == _id)[0] + group_positions = xyz[group_indices] + group_masses = masses[group_indices] + com_xyz = np.sum(group_positions.T * group_masses, axis=1) / sum( + group_masses + ) + 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)) """ Permittivity of free space = 2.39725e-4 e^2/((kcal/mol)(angstrom)), @@ -330,6 +347,7 @@ def _parse_particle_information( snapshot.particles.typeid[0:] = typeids snapshot.particles.mass[0:] = masses snapshot.particles.charge[0:] = charges / charge_factor + snapshot.particles.body = rigid_id_tags elif isinstance(snapshot, gsd.hoomd.Frame): snapshot.particles.N = top.n_sites snapshot.particles.types = unique_types @@ -337,11 +355,7 @@ def _parse_particle_information( snapshot.particles.typeid = typeids snapshot.particles.mass = masses snapshot.particles.charge = charges / charge_factor - if rigid_bodies: - warnings.warn( - "Rigid bodies detected, but not yet implemented for GSD", - NotYetImplementedWarning, - ) + snapshot.particles.body = rigid_id_tags def _parse_pairs_information( From ad603fa7d2238bbdd5174084e1815b152d1f5cf7 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 23 Sep 2024 11:37:21 -0600 Subject: [PATCH 06/19] add if statement before writing rigid tags --- gmso/external/convert_hoomd.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 11d9e287..f5ece04c 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -304,7 +304,6 @@ def _parse_particle_information( n_rigid = 0 unique_types = sorted(list(set(types))) typeids = np.array([unique_types.index(t) for t in types]) - rigid_id_tags = None if write_rigid: rigid_masses = np.zeros(n_rigid) rigid_xyz = np.zeros(n_rigid) @@ -347,7 +346,8 @@ def _parse_particle_information( snapshot.particles.typeid[0:] = typeids snapshot.particles.mass[0:] = masses snapshot.particles.charge[0:] = charges / charge_factor - snapshot.particles.body = rigid_id_tags + if write_rigid: + snapshot.particles.body[0:] = rigid_id_tags elif isinstance(snapshot, gsd.hoomd.Frame): snapshot.particles.N = top.n_sites snapshot.particles.types = unique_types @@ -355,7 +355,8 @@ def _parse_particle_information( snapshot.particles.typeid = typeids snapshot.particles.mass = masses snapshot.particles.charge = charges / charge_factor - snapshot.particles.body = rigid_id_tags + if write_rigid: + snapshot.particles.body = rigid_id_tags def _parse_pairs_information( From fdba3156fc7d401592ceab81a120e7b0dbab1bde Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 23 Sep 2024 12:40:43 -0600 Subject: [PATCH 07/19] add basic rigid body snapshot test --- gmso/external/convert_hoomd.py | 26 ++++---- gmso/tests/test_hoomd.py | 105 ++++++++------------------------- 2 files changed, 37 insertions(+), 94 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index f5ece04c..19f287fd 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -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: @@ -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] @@ -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)), @@ -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 diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index f3032871..5c2dbc3e 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -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 @@ -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" From 3a06f245275ef92884790a56aed06739934f175c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 23 Sep 2024 18:40:56 +0000 Subject: [PATCH 08/19] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gmso/external/convert_hoomd.py | 6 ++---- gmso/tests/test_hoomd.py | 12 ++++++------ 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 19f287fd..62cfaffa 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -292,7 +292,7 @@ def _parse_particle_information( if site.mass else 1 * base_units["mass"] ) - charges[idx] = (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) @@ -311,9 +311,7 @@ def _parse_particle_information( # 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( - (np.array([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] diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index 5c2dbc3e..8990202d 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -1,16 +1,19 @@ import forcefield_utilities as ffutils import hoomd -import numpy as np import pytest import unyt as u + from gmso import ForceField from gmso.external import from_mbuild -from gmso.external.convert_hoomd import to_hoomd_forcefield, to_hoomd_snapshot, to_gsd_snapshot +from gmso.external.convert_hoomd import ( + to_gsd_snapshot, + to_hoomd_forcefield, + to_hoomd_snapshot, +) from gmso.parameterization import apply from gmso.tests.base_test import BaseTest from gmso.tests.utils import get_path from gmso.utils.io import has_hoomd, has_mbuild, import_ -from gmso.utils.sorting import sort_connection_strings if has_hoomd: hoomd = import_("hoomd") @@ -51,9 +54,6 @@ 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): From 88429736e969d8d7d06e0ed86ba6f3880ceea7d5 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 25 Sep 2024 16:00:02 -0600 Subject: [PATCH 09/19] remove rigid id from convert mbuild, remove rigid id prop from abstract site --- gmso/abc/abstract_site.py | 11 ----------- gmso/external/convert_mbuild.py | 2 -- gmso/tests/test_convert_mbuild.py | 11 ----------- 3 files changed, 24 deletions(-) diff --git a/gmso/abc/abstract_site.py b/gmso/abc/abstract_site.py index d9678a8e..fec774af 100644 --- a/gmso/abc/abstract_site.py +++ b/gmso/abc/abstract_site.py @@ -220,12 +220,6 @@ def __repr__(self): alias="position", ) - rigid_id_: Optional[int] = Field( - None, - description="Rigid body ID used to assign rigid body contraints in HOOMD-Blue topologies.", - alias="rigid_id", - ) - model_config = ConfigDict( alias_to_fields={ "name": "name_", @@ -267,11 +261,6 @@ def residue(self): """Return the residue assigned to the site.""" return self.__dict__.get("residue_") - @property - def rigid_id(self): - """Return the rigid body ID for this atom.""" - return self.__dict__.get("rigid_id_") - @field_serializer("position_") def serialize_position(self, position_: PositionType): return unyt_to_dict(position_) diff --git a/gmso/external/convert_mbuild.py b/gmso/external/convert_mbuild.py index c8ca1e50..5ac724ce 100644 --- a/gmso/external/convert_mbuild.py +++ b/gmso/external/convert_mbuild.py @@ -251,7 +251,6 @@ def _parse_particle(particle_map, site): charge=charge, mass=mass, ) - particle.rigid_id = site.rigid_id particle_map[site] = particle return particle @@ -273,7 +272,6 @@ def _parse_site(site_map, particle, search_method, infer_element=False): element=ele, charge=charge, mass=mass, - rigid_id=particle.rigid_id, molecule=site_map[particle]["molecule"], residue=site_map[particle]["residue"], group=site_map[particle]["group"], diff --git a/gmso/tests/test_convert_mbuild.py b/gmso/tests/test_convert_mbuild.py index 97620756..7b7a9c44 100644 --- a/gmso/tests/test_convert_mbuild.py +++ b/gmso/tests/test_convert_mbuild.py @@ -259,14 +259,3 @@ def test_nontop_level_compound(self, mb_ethane): cpd.add(mb_ethane) with pytest.raises(AssertionError): from_mbuild(mb_ethane) - - def test_rigid_ids(self, mb_ethane): - box = mb.fill_box(mb_ethane, n_compounds=3, box=[2, 2, 2]) - for i, child in enumerate(box.children): - for p in child.particles(): - p.rigid_id = i - - top = from_mbuild(box) - assert set([site.rigid_id for site in top.sites]) == {0, 1, 2} - for particle, site in zip(box.particles(), top.sites): - assert particle.rigid_id == site.rigid_id From 6a383c61774cbe5e8ce7e89c3e5f1c902fe69933 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 25 Sep 2024 23:11:40 -0600 Subject: [PATCH 10/19] parse rigid sites from site.molecule, update test --- gmso/external/convert_hoomd.py | 21 +++++++++++---------- gmso/tests/test_hoomd.py | 16 +++++++--------- 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 62cfaffa..88e1639e 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -293,18 +293,16 @@ def _parse_particle_information( else 1 * base_units["mass"] ) 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 None not in rigid_ids: - n_rigid = len(rigid_ids_set) - write_rigid = True - else: - write_rigid = False - n_rigid = 0 + unique_types = sorted(list(set(types))) typeids = np.array([unique_types.index(t) for t in types]) - if write_rigid: + # Check for rigid molecules + rigid_mols = any([site.molecule.isrigid for site in top.sites]) + if rigid_mols: + rigid_ids = [site.molecule.number for site in top.sites] + rigid_ids_set = set(rigid_ids) + n_rigid = len(rigid_ids_set) + write_rigid = True rigid_masses = np.zeros(n_rigid) rigid_xyz = np.zeros((n_rigid, 3)) # Rigid particle type defaults to "R"; add to front of list @@ -326,6 +324,9 @@ def _parse_particle_information( 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))) + else: + write_rigid = False + n_rigid = 0 """ Permittivity of free space = 2.39725e-4 e^2/((kcal/mol)(angstrom)), diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index 8990202d..75452331 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -60,21 +60,19 @@ 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) + for site in top.sites: + site.molecule.isrigid = True + + top_no_rigid = from_mbuild(box) - rigid_ids = [site.rigid_id for site in top.sites] - assert len(rigid_ids) == box.n_particles - assert len(set(rigid_ids)) == 2 + rigid_ids = [site.molecule.number for site in top.sites] + assert set(rigid_ids) == {0, 1} snapshot, refs = to_gsd_snapshot(top) snapshot_no_rigid, refs = to_gsd_snapshot(top_no_rigid) assert "R" in snapshot.particles.types + assert "R" not in snapshot_no_rigid.particles.types assert snapshot.particles.N - 2 == snapshot_no_rigid.particles.N @pytest.mark.skipif( From 7c656243c273746c12f59bf42e8c964175df0b97 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Sun, 29 Sep 2024 21:39:58 -0600 Subject: [PATCH 11/19] remove un-needed var --- gmso/external/convert_hoomd.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 88e1639e..8859459c 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -302,7 +302,6 @@ def _parse_particle_information( rigid_ids = [site.molecule.number for site in top.sites] rigid_ids_set = set(rigid_ids) n_rigid = len(rigid_ids_set) - write_rigid = True rigid_masses = np.zeros(n_rigid) rigid_xyz = np.zeros((n_rigid, 3)) # Rigid particle type defaults to "R"; add to front of list @@ -325,7 +324,6 @@ def _parse_particle_information( xyz = np.concatenate((rigid_xyz, xyz)) rigid_id_tags = np.concatenate((np.arange(n_rigid), np.array(rigid_ids))) else: - write_rigid = False n_rigid = 0 """ @@ -347,7 +345,7 @@ def _parse_particle_information( snapshot.particles.typeid[0:] = typeids snapshot.particles.mass[0:] = masses snapshot.particles.charge[0:] = charges / charge_factor - if write_rigid: + if n_rigid: snapshot.particles.body[0:] = rigid_id_tags elif isinstance(snapshot, gsd.hoomd.Frame): snapshot.particles.N = top.n_sites + n_rigid @@ -356,7 +354,7 @@ def _parse_particle_information( snapshot.particles.typeid = typeids snapshot.particles.mass = masses snapshot.particles.charge = charges / charge_factor - if write_rigid: + if n_rigid: snapshot.particles.body = rigid_id_tags From 8da6fe571a0058dff07bfe4db9021d1521b80fb9 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 16 Oct 2024 13:55:49 -0600 Subject: [PATCH 12/19] update tests --- gmso/external/convert_hoomd.py | 9 ++++++--- gmso/tests/test_hoomd.py | 11 ++++++++++- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 8859459c..0d0d7780 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -284,15 +284,17 @@ def _parse_particle_information( site.name if site.atom_type is None else site.atom_type.name for site in top.sites ] - masses = np.zeros(len(xyz)) - charges = np.zeros(len(xyz)) + masses = np.zeros(top.n_sites) + # charges = np.zeros(top.n_sites) + charges = list() 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[idx] = site.charge if site.charge else 0 * u.elementary_charge + # charges[idx] = site.charge if site.charge else 0 * u.elementary_charge + charges.append(site.charge if site.charge else 0 * u.elementary_charge) unique_types = sorted(list(set(types))) typeids = np.array([unique_types.index(t) for t in types]) @@ -322,6 +324,7 @@ def _parse_particle_information( # Append rigid center mass and xyz to front masses = np.concatenate((rigid_masses, masses)) xyz = np.concatenate((rigid_xyz, xyz)) + charges = np.concatenate((np.zeros(n_rigid), np.array(charges))) rigid_id_tags = np.concatenate((np.arange(n_rigid), np.array(rigid_ids))) else: n_rigid = 0 diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index 75452331..058cdfc6 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -1,5 +1,6 @@ import forcefield_utilities as ffutils import hoomd +import numpy as np import pytest import unyt as u @@ -28,7 +29,7 @@ def run_hoomd_nvt(snapshot, forces, vhoomd=4): sim = hoomd.Simulation(device=cpu) sim.create_state_from_snapshot(snapshot) - integrator = hoomd.md.Integrator(dt=0.001) + integrator = hoomd.md.Integrator(dt=0.0001) integrator.forces = list(set().union(*forces.values())) temp = 300 * u.K @@ -74,6 +75,14 @@ def test_rigid_bodies(self): assert "R" in snapshot.particles.types assert "R" not in snapshot_no_rigid.particles.types assert snapshot.particles.N - 2 == snapshot_no_rigid.particles.N + assert np.array_equal(snapshot.particles.typeid[:2], np.array([0, 0])) + assert np.array_equal( + snapshot.particles.body[2:10], np.array([0] * ethane.n_particles) + ) + assert np.array_equal( + snapshot.particles.body[10:], np.array([1] * ethane.n_particles) + ) + assert snapshot.particles.mass[0] == ethane.mass @pytest.mark.skipif( int(hoomd_version[0]) < 4, reason="Unsupported features in HOOMD 3" From 75289bb75bc7c1f35848a071ca91911aaa1aa9be Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 17 Oct 2024 09:56:11 -0600 Subject: [PATCH 13/19] fix charge array in np.concatenate call --- gmso/external/convert_hoomd.py | 10 ++++------ gmso/tests/test_hoomd.py | 2 +- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 0d0d7780..11bef706 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -285,16 +285,14 @@ def _parse_particle_information( for site in top.sites ] masses = np.zeros(top.n_sites) - # charges = np.zeros(top.n_sites) - charges = list() + charges = np.zeros(top.n_sites) 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[idx] = site.charge if site.charge else 0 * u.elementary_charge - charges.append(site.charge if site.charge else 0 * u.elementary_charge) + charges[idx] = site.charge if site.charge else 0 * u.elementary_charge unique_types = sorted(list(set(types))) typeids = np.array([unique_types.index(t) for t in types]) @@ -313,7 +311,7 @@ def _parse_particle_information( 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] + group_indices = np.where(np.array(rigid_ids) == _id)[0] group_positions = xyz[group_indices] group_masses = masses[group_indices] com_xyz = np.sum(group_positions.T * group_masses, axis=1) / sum( @@ -324,7 +322,7 @@ def _parse_particle_information( # Append rigid center mass and xyz to front masses = np.concatenate((rigid_masses, masses)) xyz = np.concatenate((rigid_xyz, xyz)) - charges = np.concatenate((np.zeros(n_rigid), np.array(charges))) + charges = np.concatenate((np.zeros(n_rigid), charges)) rigid_id_tags = np.concatenate((np.arange(n_rigid), np.array(rigid_ids))) else: n_rigid = 0 diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index ad984924..8035d429 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -82,7 +82,7 @@ def test_rigid_bodies(self): assert np.array_equal( snapshot.particles.body[10:], np.array([1] * ethane.n_particles) ) - assert snapshot.particles.mass[0] == ethane.mass + assert np.allclose(snapshot.particles.mass[0], ethane.mass, atol=1e-2) @pytest.mark.skipif( int(hoomd_version[0]) < 4, reason="Unsupported features in HOOMD 3" From 0d6e7d256059b7307d9de954f7bbc950349aa10c Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 17 Oct 2024 16:48:12 -0600 Subject: [PATCH 14/19] adjust snapshot groups by num rigid bodies --- gmso/external/convert_hoomd.py | 81 +++++++++++++++++++++++----------- gmso/tests/test_hoomd.py | 11 +++++ 2 files changed, 66 insertions(+), 26 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 11bef706..737afa82 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -126,7 +126,7 @@ def to_gsd_snapshot( NotYetImplementedWarning, ) - _parse_particle_information( + n_rigid = _parse_particle_information( gsd_snapshot, top, base_units, @@ -135,15 +135,15 @@ def to_gsd_snapshot( u.unyt_array([lx, ly, lz]), ) if parse_special_pairs: - _parse_pairs_information(gsd_snapshot, top) + _parse_pairs_information(gsd_snapshot, top, n_rigid) if top.n_bonds > 0: - _parse_bond_information(gsd_snapshot, top) + _parse_bond_information(gsd_snapshot, top, n_rigid) if top.n_angles > 0: - _parse_angle_information(gsd_snapshot, top) + _parse_angle_information(gsd_snapshot, top, n_rigid) if top.n_dihedrals > 0: - _parse_dihedral_information(gsd_snapshot, top) + _parse_dihedral_information(gsd_snapshot, top, n_rigid) if top.n_impropers > 0: - _parse_improper_information(gsd_snapshot, top) + _parse_improper_information(gsd_snapshot, top, n_rigid) return gsd_snapshot, base_units @@ -224,7 +224,7 @@ def to_hoomd_snapshot( NotYetImplementedWarning, ) - _parse_particle_information( + n_rigid = _parse_particle_information( hoomd_snapshot, top, base_units, @@ -233,15 +233,15 @@ def to_hoomd_snapshot( u.unyt_array([lx, ly, lz]), ) if parse_special_pairs: - _parse_pairs_information(hoomd_snapshot, top) + _parse_pairs_information(hoomd_snapshot, top, n_rigid) if top.n_bonds > 0: - _parse_bond_information(hoomd_snapshot, top) + _parse_bond_information(hoomd_snapshot, top, n_rigid) if top.n_angles > 0: - _parse_angle_information(hoomd_snapshot, top) + _parse_angle_information(hoomd_snapshot, top, n_rigid) if top.n_dihedrals > 0: - _parse_dihedral_information(hoomd_snapshot, top) + _parse_dihedral_information(hoomd_snapshot, top, n_rigid) if top.n_impropers > 0: - _parse_improper_information(hoomd_snapshot, top) + _parse_improper_information(hoomd_snapshot, top, n_rigid) hoomd_snapshot.wrap() return hoomd_snapshot, base_units @@ -357,13 +357,23 @@ def _parse_particle_information( snapshot.particles.charge = charges / charge_factor if n_rigid: snapshot.particles.body = rigid_id_tags + return n_rigid -def _parse_pairs_information( - snapshot, - top, -): - """Parse scaled pair types.""" +def _parse_pairs_information(snapshot, top, n_rigid=0): + """Parse sacled pair types. + + Parameters + ---------- + snapshot : gsd.hoomd.Frame or hoomd.Snapshot + The target Snapshot object. + top : gmso.Topology + Topology object holding system information + n_rigid : int + The number of rigid bodies found in `_parse_particle_information()` + Used to adjust pair group indices. + + """ pair_types = list() pair_typeids = list() pairs = list() @@ -383,21 +393,27 @@ def _parse_pairs_information( if pair_type not in pair_types: pair_types.append(pair_type) pair_typeids.append(pair_types.index(pair_type)) - pairs.append((top.get_index(pair[0]), top.get_index(pair[1]))) + pairs.append( + (top.get_index(pair[0] + n_rigid), top.get_index(pair[1] + n_rigid)) + ) if isinstance(snapshot, hoomd.Snapshot): snapshot.pairs.N = len(pairs) snapshot.pairs.group[:] = np.reshape(pairs, (-1, 2)) snapshot.pairs.types = pair_types snapshot.pairs.typeid[:] = pair_typeids + if n_rigid: + snapshot.pairs.group += n_rigid elif isinstance(snapshot, gsd.hoomd.Frame): snapshot.pairs.N = len(pairs) snapshot.pairs.group = np.reshape(pairs, (-1, 2)) snapshot.pairs.types = pair_types snapshot.pairs.typeid = pair_typeids + if n_rigid: + snapshot.pairs.group += n_rigid -def _parse_bond_information(snapshot, top): +def _parse_bond_information(snapshot, top, n_rigid=0): """Parse bonds information from topology. Parameters @@ -406,6 +422,9 @@ def _parse_bond_information(snapshot, top): The target Snapshot object. top : gmso.Topology Topology object holding system information + n_rigid : int + The number of rigid bodies found in `_parse_particle_information()` + Used to adjust bond group indices. """ snapshot.bonds.N = top.n_bonds @@ -426,7 +445,7 @@ def _parse_bond_information(snapshot, top): bond_types.append(bond_type) bond_groups.append( - sorted(tuple(top.get_index(site) for site in connection_members)) + sorted(tuple(top.get_index(site) + n_rigid for site in connection_members)) ) unique_bond_types = list(set(bond_types)) @@ -444,7 +463,7 @@ def _parse_bond_information(snapshot, top): warnings.warn(f"{len(unique_bond_types)} unique bond types detected") -def _parse_angle_information(snapshot, top): +def _parse_angle_information(snapshot, top, n_rigid=0): """Parse angles information from topology. Parameters @@ -453,6 +472,9 @@ def _parse_angle_information(snapshot, top): The target Snapshot object. top : gmso.Topology Topology object holding system information + n_rigid : int + The number of rigid bodies found in `_parse_particle_information()` + Used to adjust angle group indices. """ snapshot.angles.N = top.n_angles @@ -472,7 +494,9 @@ def _parse_angle_information(snapshot, top): angle_type = "-".join([site.name for site in connection_members]) angle_types.append(angle_type) - angle_groups.append(tuple(top.get_index(site) for site in connection_members)) + angle_groups.append( + tuple(top.get_index(site) + n_rigid for site in connection_members) + ) unique_angle_types = list(set(angle_types)) angle_typeids = [unique_angle_types.index(i) for i in angle_types] @@ -490,7 +514,7 @@ def _parse_angle_information(snapshot, top): warnings.warn(f"{len(unique_angle_types)} unique angle types detected") -def _parse_dihedral_information(snapshot, top): +def _parse_dihedral_information(snapshot, top, n_rigid=0): """Parse dihedral information from topology. Parameters @@ -499,6 +523,9 @@ def _parse_dihedral_information(snapshot, top): The target Snapshot object. top : gmso.Topology Topology object holding system information + n_rigid : int + The number of rigid bodies found in `_parse_particle_information()` + Used to adjust dihedral group indices. """ snapshot.dihedrals.N = top.n_dihedrals @@ -517,7 +544,7 @@ def _parse_dihedral_information(snapshot, top): dihedral_types.append(dihedral_type) dihedral_groups.append( - tuple(top.get_index(site) for site in connection_members) + tuple(top.get_index(site) + n_rigid for site in connection_members) ) unique_dihedral_types = list(set(dihedral_types)) @@ -536,7 +563,7 @@ def _parse_dihedral_information(snapshot, top): warnings.warn(f"{len(unique_dihedral_types)} unique dihedral types detected") -def _parse_improper_information(snapshot, top): +def _parse_improper_information(snapshot, top, n_rigid=0): """Parse impropers information from topology. Parameters @@ -545,6 +572,8 @@ def _parse_improper_information(snapshot, top): The target Snapshot object. top : gmso.Topology Topology object holding system information + n_rigid : int + The number of rigid bodies found in `_parse_particle_information()` """ snapshot.impropers.N = top.n_impropers @@ -563,7 +592,7 @@ def _parse_improper_information(snapshot, top): improper_types.append(improper_type) improper_groups.append( - tuple(top.get_index(site) for site in connection_members) + tuple(top.get_index(site) + n_rigid for site in connection_members) ) unique_improper_types = list(set(improper_types)) diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index 8035d429..c941dd70 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -61,6 +61,10 @@ 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]) + cg_comp = mb.Compound() + for child in box.children: + cg_comp.add(mb.Compound(name="_A", pos=child.center, mass=child.mass)) + cg_top = from_mbuild(cg_comp) top = from_mbuild(box) for site in top.sites: site.molecule.isrigid = True @@ -71,6 +75,7 @@ def test_rigid_bodies(self): assert set(rigid_ids) == {0, 1} snapshot, refs = to_gsd_snapshot(top) + cg_snapshot, refs = to_gsd_snapshot(cg_top) snapshot_no_rigid, refs = to_gsd_snapshot(top_no_rigid) assert "R" in snapshot.particles.types assert "R" not in snapshot_no_rigid.particles.types @@ -83,6 +88,12 @@ def test_rigid_bodies(self): snapshot.particles.body[10:], np.array([1] * ethane.n_particles) ) assert np.allclose(snapshot.particles.mass[0], ethane.mass, atol=1e-2) + for i in range(2): + assert np.allclose( + snapshot.particles.position[i], + cg_snapshot.particles.position[i], + atol=1e-3, + ) @pytest.mark.skipif( int(hoomd_version[0]) < 4, reason="Unsupported features in HOOMD 3" From a72d898ede99db1cfca1def0bbafeff353b40f6d Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 17 Oct 2024 17:32:00 -0600 Subject: [PATCH 15/19] update tests --- gmso/external/convert_hoomd.py | 2 +- gmso/tests/test_hoomd.py | 33 ++++++++++++++++++++++----------- 2 files changed, 23 insertions(+), 12 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 737afa82..ec6ddb80 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -394,7 +394,7 @@ def _parse_pairs_information(snapshot, top, n_rigid=0): pair_types.append(pair_type) pair_typeids.append(pair_types.index(pair_type)) pairs.append( - (top.get_index(pair[0] + n_rigid), top.get_index(pair[1] + n_rigid)) + (top.get_index(pair[0]) + n_rigid, top.get_index(pair[1]) + n_rigid) ) if isinstance(snapshot, hoomd.Snapshot): diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index c941dd70..c411a593 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -61,22 +61,20 @@ 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]) - cg_comp = mb.Compound() - for child in box.children: - cg_comp.add(mb.Compound(name="_A", pos=child.center, mass=child.mass)) - cg_top = from_mbuild(cg_comp) top = from_mbuild(box) for site in top.sites: site.molecule.isrigid = True + top.identify_connections() top_no_rigid = from_mbuild(box) + top_no_rigid.identify_connections() rigid_ids = [site.molecule.number for site in top.sites] assert set(rigid_ids) == {0, 1} snapshot, refs = to_gsd_snapshot(top) - cg_snapshot, refs = to_gsd_snapshot(cg_top) snapshot_no_rigid, refs = to_gsd_snapshot(top_no_rigid) + # Check that snapshot has rigid particles added assert "R" in snapshot.particles.types assert "R" not in snapshot_no_rigid.particles.types assert snapshot.particles.N - 2 == snapshot_no_rigid.particles.N @@ -88,12 +86,25 @@ def test_rigid_bodies(self): snapshot.particles.body[10:], np.array([1] * ethane.n_particles) ) assert np.allclose(snapshot.particles.mass[0], ethane.mass, atol=1e-2) - for i in range(2): - assert np.allclose( - snapshot.particles.position[i], - cg_snapshot.particles.position[i], - atol=1e-3, - ) + # Check that topology isn't changed by using rigid bodies + assert snapshot.bonds.N == snapshot_no_rigid.bonds.N + assert snapshot.angles.N == snapshot_no_rigid.angles.N + assert snapshot.dihedrals.N == snapshot_no_rigid.dihedrals.N + # Check that particle indices in snapshot groups are adjust by num rigid bodies + for group1, group2 in zip(snapshot.bonds.group, snapshot_no_rigid.bonds.group): + assert np.array_equal(np.array(group1), np.array(group2) + 2) + for group1, group2 in zip( + snapshot.angles.group, snapshot_no_rigid.angles.group + ): + assert np.array_equal(np.array(group1), np.array(group2) + 2) + for group1, group2 in zip( + snapshot.dihedrals.group, snapshot_no_rigid.dihedrals.group + ): + assert np.array_equal(np.array(group1), np.array(group2) + 2) + for group1, group2 in zip( + snapshot.impropers.group, snapshot_no_rigid.impropers.group + ): + assert np.array_equal(np.array(group1), np.array(group2) + 2) @pytest.mark.skipif( int(hoomd_version[0]) < 4, reason="Unsupported features in HOOMD 3" From 5e8635f3a9a0f7cbd992e0c9a238a636866f8db1 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 17 Oct 2024 17:36:14 -0600 Subject: [PATCH 16/19] fix _parse_pairs --- gmso/external/convert_hoomd.py | 4 ---- gmso/tests/test_hoomd.py | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index ec6ddb80..cb777f13 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -402,15 +402,11 @@ def _parse_pairs_information(snapshot, top, n_rigid=0): snapshot.pairs.group[:] = np.reshape(pairs, (-1, 2)) snapshot.pairs.types = pair_types snapshot.pairs.typeid[:] = pair_typeids - if n_rigid: - snapshot.pairs.group += n_rigid elif isinstance(snapshot, gsd.hoomd.Frame): snapshot.pairs.N = len(pairs) snapshot.pairs.group = np.reshape(pairs, (-1, 2)) snapshot.pairs.types = pair_types snapshot.pairs.typeid = pair_typeids - if n_rigid: - snapshot.pairs.group += n_rigid def _parse_bond_information(snapshot, top, n_rigid=0): diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index c411a593..49e1c12d 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -90,7 +90,7 @@ def test_rigid_bodies(self): assert snapshot.bonds.N == snapshot_no_rigid.bonds.N assert snapshot.angles.N == snapshot_no_rigid.angles.N assert snapshot.dihedrals.N == snapshot_no_rigid.dihedrals.N - # Check that particle indices in snapshot groups are adjust by num rigid bodies + # Check if particle indices in snapshot groups are adjusted by N rigid bodies for group1, group2 in zip(snapshot.bonds.group, snapshot_no_rigid.bonds.group): assert np.array_equal(np.array(group1), np.array(group2) + 2) for group1, group2 in zip( From b96a83c1655593f7b1f25177f493d2bc4e3ed366 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Fri, 18 Oct 2024 15:58:32 -0600 Subject: [PATCH 17/19] use unyt_arrays for charges and mass instead of lists --- gmso/external/convert_hoomd.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index cb777f13..41b51186 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -284,15 +284,17 @@ def _parse_particle_information( site.name if site.atom_type is None else site.atom_type.name for site in top.sites ] - masses = np.zeros(top.n_sites) - charges = np.zeros(top.n_sites) - for idx, site in enumerate(top.sites): - masses[idx] = ( + masses = u.unyt_array( + [ site.mass.to_value(base_units["mass"]) if site.mass else 1 * base_units["mass"] - ) - charges[idx] = site.charge if site.charge else 0 * u.elementary_charge + for site in top.sites + ] + ) + charges = u.unyt_array( + [site.charge if site.charge else 0 * u.elementary_charge for site in top.sites] + ) unique_types = sorted(list(set(types))) typeids = np.array([unique_types.index(t) for t in types]) @@ -302,8 +304,9 @@ def _parse_particle_information( rigid_ids = [site.molecule.number for site in top.sites] rigid_ids_set = set(rigid_ids) n_rigid = len(rigid_ids_set) - rigid_masses = np.zeros(n_rigid) - rigid_xyz = np.zeros((n_rigid, 3)) + rigid_charges = np.zeros(n_rigid) * charges.units + rigid_masses = np.zeros(n_rigid) * masses.units + rigid_xyz = np.zeros((n_rigid, 3)) * xyz.units # 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 @@ -322,8 +325,8 @@ def _parse_particle_information( # Append rigid center mass and xyz to front masses = np.concatenate((rigid_masses, masses)) xyz = np.concatenate((rigid_xyz, xyz)) - charges = np.concatenate((np.zeros(n_rigid), charges)) - rigid_id_tags = np.concatenate((np.arange(n_rigid), np.array(rigid_ids))) + charges = np.concatenate((rigid_charges, charges)) + rigid_id_tags = np.concatenate((np.arange(n_rigid), rigid_ids)) else: n_rigid = 0 From ffe639c4e370f7f69749cf024d7b55b34a0aed52 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 29 Oct 2024 22:11:46 -0600 Subject: [PATCH 18/19] handle mix of rigid and non rigid --- gmso/external/convert_hoomd.py | 8 ++++++-- gmso/tests/test_hoomd.py | 23 +++++++++++++++++++++++ 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 41b51186..76a071c1 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -301,8 +301,10 @@ def _parse_particle_information( # Check for rigid molecules rigid_mols = any([site.molecule.isrigid for site in top.sites]) if rigid_mols: - rigid_ids = [site.molecule.number for site in top.sites] - rigid_ids_set = set(rigid_ids) + rigid_ids = [ + site.molecule.number if site.molecule.isrigid else -1 for site in top.sites + ] + rigid_ids_set = set([i for i in rigid_ids if i != -1]) n_rigid = len(rigid_ids_set) rigid_charges = np.zeros(n_rigid) * charges.units rigid_masses = np.zeros(n_rigid) * masses.units @@ -314,6 +316,8 @@ def _parse_particle_information( 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): + if _id == 1: + continue group_indices = np.where(np.array(rigid_ids) == _id)[0] group_positions = xyz[group_indices] group_masses = masses[group_indices] diff --git a/gmso/tests/test_hoomd.py b/gmso/tests/test_hoomd.py index 49e1c12d..161448fd 100644 --- a/gmso/tests/test_hoomd.py +++ b/gmso/tests/test_hoomd.py @@ -106,6 +106,29 @@ def test_rigid_bodies(self): ): assert np.array_equal(np.array(group1), np.array(group2) + 2) + def test_rigid_bodies_mix(self): + ethane = mb.lib.molecules.Ethane() + ethane.name = "ethane" + benzene = mb.load("c1ccccc1", smiles=True) + benzene.name = "benzene" + box = mb.fill_box([benzene, ethane], n_compounds=[1, 1], box=[2, 2, 2]) + top = from_mbuild(box) + for site in top.sites: + if site.molecule.name == "benzene": + site.molecule.isrigid = True + + snapshot, refs = to_gsd_snapshot(top) + assert snapshot.particles.typeid[0] == 0 + assert snapshot.particles.N == top.n_sites + 1 + assert np.array_equal( + snapshot.particles.body[-ethane.n_particles :], + np.array([-1] * ethane.n_particles), + ) + assert np.array_equal( + snapshot.particles.body[1 : benzene.n_particles + 1], + np.array([0] * benzene.n_particles), + ) + @pytest.mark.skipif( int(hoomd_version[0]) < 4, reason="Unsupported features in HOOMD 3" ) From 4eaf697508499cc9496101f57aeb10dafbf29b89 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 30 Oct 2024 09:44:46 -0600 Subject: [PATCH 19/19] remove if statement to filter -1 ids --- gmso/external/convert_hoomd.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/gmso/external/convert_hoomd.py b/gmso/external/convert_hoomd.py index 76a071c1..81632286 100644 --- a/gmso/external/convert_hoomd.py +++ b/gmso/external/convert_hoomd.py @@ -316,8 +316,6 @@ def _parse_particle_information( 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): - if _id == 1: - continue group_indices = np.where(np.array(rigid_ids) == _id)[0] group_positions = xyz[group_indices] group_masses = masses[group_indices]