From a10449cc9c7a8f87e7958631e414731c5c2262b2 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 25 Sep 2023 16:02:40 -0600 Subject: [PATCH 01/41] use default auto_scale=False for gsd functions --- hoomd_organics/base/system.py | 28 +++++----------------------- 1 file changed, 5 insertions(+), 23 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index d3428bb5..9d9182f0 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -48,23 +48,9 @@ class System(ABC): target_box attribute. Can be useful when initializing systems at low density and running a shrink simulation to achieve a target density. - r_cut : float, required - The cutoff radius for the Lennard-Jones potential. force_field : hoomd_organics.ForceField or a list of ForceField objects, default=None The force field to be applied to the system for parameterization. - auto_scale : bool, default=False - Set to true to use reduced simulation units. - distance, mass, and energy are scaled by the largest value - present in the system for each. - remove_hydrogens : bool, default False - Set to true to remove hydrogen atoms from the system. - The masses and charges of the hydrogens are absorbed into - the heavy atoms they were bonded to. - remove_charges : bool, default False - Set to true to remove charges from the system. - scale_charges : bool, default False - Set to true to scale charges to net zero. base_units : dict, default {} Dictionary of base units to use for scaling. Dictionary keys are "length", "mass", and "energy". Values should be an @@ -90,7 +76,6 @@ def __init__( molecules, density: float, force_field=None, - auto_scale=False, base_units=dict(), ): self._molecules = check_return_iterable(molecules) @@ -98,7 +83,6 @@ def __init__( if force_field: self._force_field = check_return_iterable(force_field) self.density = density - self.auto_scale = auto_scale self.all_molecules = [] self.gmso_system = None self._reference_values = base_units @@ -110,6 +94,7 @@ def __init__( self._ff_refs = dict() self._snap_refs = dict() self._ff_kwargs = dict() + self.auto_scale = False # Collecting all molecules self.n_mol_types = 0 @@ -383,7 +368,7 @@ def remove_hydrogens(self): site.mass += h.mass site.charge += h.charge self.gmso_system.remove_site(site=h) - # If a snap shot was already made, need to re-create it w/o hydrogens + # If a snapshot was already made, need to re-create it w/o hydrogens if self._hoomd_snapshot: self._create_hoomd_snapshot() @@ -428,7 +413,7 @@ def _create_hoomd_forcefield(self, r_cut, nlist_buffer, pppm_kwargs): r_cut=r_cut, nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs, - auto_scale=self.auto_scale, + auto_scale=False, base_units=self._reference_values if self._reference_values else None, @@ -442,7 +427,7 @@ def _create_hoomd_snapshot(self): """Create a HOOMD snapshot.""" snap, refs = to_gsd_snapshot( top=self.gmso_system, - auto_scale=self.auto_scale, + auto_scale=False, base_units=self._reference_values if self._reference_values else None, @@ -492,6 +477,7 @@ def apply_forcefield( Neighborlist buffer for simulation cell. """ + self.auto_scale = auto_scale if not self._force_field: # TODO: Better erorr message raise ValueError( @@ -648,7 +634,6 @@ def __init__( molecules, density: float, force_field=None, - auto_scale=False, base_units=dict(), packing_expand_factor=5, edge=0.2, @@ -659,7 +644,6 @@ def __init__( molecules=molecules, density=density, force_field=force_field, - auto_scale=auto_scale, base_units=base_units, ) @@ -702,7 +686,6 @@ def __init__( n: int, basis_vector=[0.5, 0.5, 0], force_field=None, - auto_scale=False, base_units=dict(), ): self.x = x @@ -713,7 +696,6 @@ def __init__( molecules=molecules, density=density, force_field=force_field, - auto_scale=auto_scale, base_units=base_units, ) From 7dc97ce143648e43e203385e16657258f3e637af Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 25 Sep 2023 16:27:50 -0600 Subject: [PATCH 02/41] move auto_scale to apply_forcefield --- hoomd_organics/tests/base/test_system.py | 150 ++++++++++------------- 1 file changed, 62 insertions(+), 88 deletions(-) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index 1409a36a..6247fceb 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -96,9 +96,8 @@ def test_system_from_mol2_mol_parameterization(self, benzene_molecule_mol2): molecules=[benzene_mol], density=0.8, force_field=OPLS_AA(), - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.gmso_system.is_typed() assert len(system.hoomd_forcefield) > 0 assert system.n_particles == system.hoomd_snapshot.particles.N @@ -109,9 +108,12 @@ def test_remove_hydrogen(self, benzene_molecule): molecules=[benzene_mol], density=0.8, force_field=OPLS_AA(), + ) + system.apply_forcefield( + r_cut=2.5, + remove_hydrogens=True, auto_scale=True, ) - system.apply_forcefield(r_cut=2.5, remove_hydrogens=True) assert not any( [s.element.atomic_number == 1 for s in system.gmso_system.sites] ) @@ -132,9 +134,10 @@ def test_remove_hydrogen_no_atomic_num(self, benzene_molecule): molecules=[benzene_mol], density=0.8, force_field=OPLS_AA(), - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5, remove_hydrogens=False) + system.apply_forcefield( + r_cut=2.5, remove_hydrogens=False, auto_scale=True + ) for site in system.gmso_system.sites: if site.name == "H": site.element = gmso.core.element.Element( @@ -153,9 +156,10 @@ def test_remove_hydrogen_no_hydrogen(self, benzene_molecule): molecules=[benzene_mol], density=0.8, force_field=OPLS_AA(), - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5, remove_hydrogens=False) + system.apply_forcefield( + r_cut=2.5, remove_hydrogens=False, auto_scale=True + ) hydrogens = [ site for site in system.gmso_system.sites @@ -173,10 +177,12 @@ def test_add_mass_charges(self, benzene_molecule): molecules=[benzene_mols], density=0.8, force_field=OPLS_AA(), - auto_scale=False, ) system.apply_forcefield( - r_cut=2.5, remove_hydrogens=True, scale_charges=False + r_cut=2.5, + remove_hydrogens=True, + scale_charges=False, + auto_scale=False, ) for site in system.gmso_system.sites: assert site.mass.value == (12.011 + 1.008) @@ -194,17 +200,15 @@ def test_target_box(self, benzene_molecule): molecules=[benzene_mol], density=0.1, force_field=OPLS_AA(), - auto_scale=True, ) - low_density_system.apply_forcefield(r_cut=2.5) + low_density_system.apply_forcefield(r_cut=2.5, auto_scale=True) high_density_system = Pack( molecules=[benzene_mol], density=0.9, force_field=OPLS_AA(), - auto_scale=True, ) - high_density_system.apply_forcefield(r_cut=2.5) + high_density_system.apply_forcefield(r_cut=2.5, auto_scale=True) assert all( low_density_system.target_box > high_density_system.target_box ) @@ -222,9 +226,8 @@ def test_ref_length(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert np.allclose( system.reference_length.to("angstrom").value, 3.5, atol=1e-3 @@ -241,9 +244,8 @@ def test_ref_mass(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) total_red_mass = sum(system.hoomd_snapshot.particles.mass) assert np.allclose( system.mass, @@ -257,9 +259,8 @@ def test_ref_energy(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert np.allclose( system.reference_energy.to("kcal/mol").value, 0.066, atol=1e-3 ) @@ -270,9 +271,8 @@ def test_rebuild_snapshot(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) assert system._snap_refs == system.reference_values assert system._ff_refs == system.reference_values init_snap = system.hoomd_snapshot @@ -292,9 +292,8 @@ def test_ref_values_no_autoscale(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_length = 1 * u.angstrom system.reference_energy = 1 * u.kcal / u.mol system.reference_mass = 1 * u.amu @@ -311,12 +310,9 @@ def test_ref_values_no_autoscale(self, polyethylene): def test_set_ref_values(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( - molecules=[polyethylene], - force_field=[OPLS_AA()], - density=1.0, - auto_scale=False, + molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0 ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -333,9 +329,8 @@ def test_set_ref_values_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) ref_value_dict = { "length": "1 angstrom", "energy": "3 kcal/mol", @@ -352,9 +347,8 @@ def test_set_ref_values_missing_key(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -368,9 +362,8 @@ def test_set_ref_values_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -385,9 +378,8 @@ def test_set_ref_length(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_length = 1 * u.angstrom assert system.reference_length == 1 * u.angstrom @@ -409,9 +401,8 @@ def test_ref_length_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_length = "1 angstrom" assert system.reference_length == 1 * u.angstrom @@ -421,9 +412,8 @@ def test_ref_length_invalid_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0" @@ -433,9 +423,8 @@ def test_ref_length_invalid_unit_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0 invalid_unit" @@ -445,9 +434,8 @@ def test_ref_length_invalid_dimension(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_length = 1.0 * u.g @@ -457,9 +445,8 @@ def test_ref_length_invalid_dimension_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0 g" @@ -469,9 +456,8 @@ def test_set_ref_energy(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_energy = 1 * u.kcal / u.mol assert system.reference_energy == 1 * u.kcal / u.mol @@ -481,9 +467,8 @@ def test_set_ref_energy_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 @@ -493,9 +478,8 @@ def test_ref_energy_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_energy = "1 kJ" assert system.reference_energy == 1 * u.kJ @@ -505,9 +489,8 @@ def test_ref_energy_string_comb_units(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_energy = "1 kcal/mol" assert system.reference_energy == 1 * u.kcal / u.mol @@ -517,9 +500,8 @@ def test_ref_energy_invalid_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0" @@ -529,9 +511,8 @@ def test_ref_energy_invalid_unit_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0 invalid_unit" @@ -541,9 +522,8 @@ def test_ref_energy_invalid_dimension(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 * u.g @@ -553,9 +533,8 @@ def test_ref_energy_invalid_dimension_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0 m" @@ -565,9 +544,8 @@ def test_set_ref_mass(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_mass = 1.0 * u.amu assert system.reference_mass == 1.0 * u.amu @@ -578,9 +556,8 @@ def test_set_ref_mass_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_mass = 1.0 @@ -590,9 +567,8 @@ def test_ref_mass_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) system.reference_mass = "1 g" assert system.reference_mass == 1.0 * u.g @@ -602,9 +578,8 @@ def test_ref_mass_invalid_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0" @@ -614,9 +589,8 @@ def test_ref_mass_invalid_unit_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 invalid_unit" @@ -626,9 +600,8 @@ def test_ref_mass_invalid_dimension(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 * u.m @@ -638,9 +611,8 @@ def test_ref_mass_invalid_dimension_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 m" @@ -650,10 +622,9 @@ def test_apply_forcefield_no_forcefield(self, polyethylene): molecules=[polyethylene], force_field=None, density=1.0, - auto_scale=False, ) with pytest.raises(ValueError): - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) def test_forcefield_kwargs_attr(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) @@ -661,10 +632,13 @@ def test_forcefield_kwargs_attr(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) system.apply_forcefield( - r_cut=2.5, nlist_buffer=0.5, pppm_resolution=(4, 4, 4), pppm_order=3 + r_cut=2.5, + auto_scale=False, + nlist_buffer=0.5, + pppm_resolution=(4, 4, 4), + pppm_order=3, ) assert system._ff_kwargs["r_cut"] == 2.5 assert system._ff_kwargs["nlist_buffer"] == 0.5 @@ -680,9 +654,8 @@ def test_lattice_polymer(self, polyethylene): x=1, y=1, n=4, - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(polyethylene.molecules) @@ -700,9 +673,8 @@ def test_lattice_molecule(self, benzene_molecule): x=1, y=1, n=4, - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(benzene_mol.molecules) assert len(system.hoomd_forcefield) > 0 @@ -715,16 +687,18 @@ def test_scale_charges(self, pps): molecules=pps_mol, density=0.5, force_field=OPLS_AA_PPS(), - auto_scale=True, ) - no_scale.apply_forcefield(r_cut=2.5, scale_charges=False) + no_scale.apply_forcefield( + r_cut=2.5, scale_charges=False, auto_scale=True + ) with_scale = Pack( molecules=pps_mol, density=0.5, force_field=OPLS_AA_PPS(), - auto_scale=True, ) - with_scale.apply_forcefield(r_cut=2.5, scale_charges=True) + with_scale.apply_forcefield( + r_cut=2.5, scale_charges=True, auto_scale=True + ) assert abs(no_scale.net_charge.value) > abs(with_scale.net_charge.value) assert np.allclose(0, with_scale.net_charge.value, atol=1e-30) From 4a48597bc98446122e5edfe28273c4b803041a49 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 25 Sep 2023 16:43:19 -0600 Subject: [PATCH 03/41] add unit test for setting ref values with auto_scale=True --- hoomd_organics/tests/base/test_system.py | 28 +++++++++++++++++------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index 6247fceb..b9527aee 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -18,9 +18,8 @@ def test_single_mol_type(self, benzene_molecule): molecules=[benzene_mols], density=0.8, force_field=OPLS_AA(), - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(benzene_mols.molecules) assert system.gmso_system.is_typed() @@ -36,9 +35,8 @@ def test_multiple_mol_types(self, benzene_molecule, ethane_molecule): molecules=[benzene_mol, ethane_mol], density=0.8, force_field=OPLS_AA(), - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 2 assert len(system.all_molecules) == len(benzene_mol.molecules) + len( ethane_mol.molecules @@ -64,9 +62,8 @@ def test_multiple_mol_types_different_ff( molecules=[dimethylether_mol, pps_mol], density=0.8, force_field=[OPLS_AA_DIMETHYLETHER(), OPLS_AA_PPS()], - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 2 assert len(system.all_molecules) == len( dimethylether_mol.molecules @@ -372,6 +369,22 @@ def test_set_ref_values_invalid_type(self, polyethylene): with pytest.raises(ReferenceUnitError): system.reference_values = ref_value_dict + def test_set_ref_values_auto_scale_true(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=5) + system = Pack( + molecules=[polyethylene], + force_field=[OPLS_AA()], + density=1.0, + ) + system.apply_forcefield(r_cut=2.5, auto_scale=True) + ref_value_dict = { + "length": 1 * u.angstrom, + "energy": 3.0 * u.kcal / u.mol, + "mass": 1.25 * u.Unit("amu"), + } + with pytest.warns(): + system.reference_values = ref_value_dict + def test_set_ref_length(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( @@ -389,9 +402,8 @@ def test_set_ref_length_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - auto_scale=False, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=False) with pytest.raises(ReferenceUnitError): system.reference_length = 1.0 From 15fef09ac9c849c00bfbab2437689b192a1b0992 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 25 Sep 2023 16:46:49 -0600 Subject: [PATCH 04/41] add more unit tests --- hoomd_organics/tests/base/test_system.py | 34 ++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index b9527aee..d11dcaa3 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -462,6 +462,17 @@ def test_ref_length_invalid_dimension_string(self, polyethylene): with pytest.raises(ReferenceUnitError): system.reference_length = "1.0 g" + def test_ref_length_auto_scale_true(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=1) + system = Pack( + molecules=[polyethylene], + force_field=[OPLS_AA()], + density=1.0, + ) + system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.reference_length = 1 * u.angstrom + assert system.reference_length == 1 * u.angstrom + def test_set_ref_energy(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( @@ -550,6 +561,17 @@ def test_ref_energy_invalid_dimension_string(self, polyethylene): with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0 m" + def test_set_ref_energy_auto_scale_true(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=1) + system = Pack( + molecules=[polyethylene], + force_field=[OPLS_AA()], + density=1.0, + ) + system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.reference_energy = 1 * u.kcal / u.mol + assert system.reference_energy == 1 * u.kcal / u.mol + def test_set_ref_mass(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( @@ -628,6 +650,18 @@ def test_ref_mass_invalid_dimension_string(self, polyethylene): with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 m" + def test_set_ref_mass_auto_scale_true(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=1) + system = Pack( + molecules=[polyethylene], + force_field=[OPLS_AA()], + density=1.0, + ) + system.apply_forcefield(r_cut=2.5, auto_scale=True) + + system.reference_mass = 1.0 * u.amu + assert system.reference_mass == 1.0 * u.amu + def test_apply_forcefield_no_forcefield(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( From e7bbdfe87ccffb54dd53291ee779cc2a5b7211e9 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 25 Sep 2023 17:01:03 -0600 Subject: [PATCH 05/41] raise error when forcefield type is invalid --- hoomd_organics/base/system.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 9d9182f0..f11df2b0 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -138,7 +138,7 @@ def __init__( else: # if there is only one ff for all molecule types ff_index = 0 - if getattr(self._force_field[ff_index], "gmso_ff"): + if hasattr(self._force_field[ff_index], "gmso_ff"): self._gmso_forcefields_dict[str(i)] = self._force_field[ ff_index ].gmso_ff @@ -242,6 +242,12 @@ def reference_length(self, length): example, unyt.unyt_quantity(1, "nm"). """ + if self.auto_scale: + warnings.warn( + "`auto_scale` was set to True for this system. " + "Setting reference length manually disables auto " + "scaling." + ) validated_length = validate_ref_value(length, u.dimensions.length) self._reference_values["length"] = validated_length @@ -259,6 +265,12 @@ def reference_energy(self, energy): example, unyt.unyt_quantity(1, "kJ/mol"). """ + if self.auto_scale: + warnings.warn( + "`auto_scale` was set to True for this system. " + "Setting reference energy manually disables auto " + "scaling." + ) validated_energy = validate_ref_value(energy, u.dimensions.energy) self._reference_values["energy"] = validated_energy @@ -276,6 +288,12 @@ def reference_mass(self, mass): example, unyt.unyt_quantity(1, "amu"). """ + if self.auto_scale: + warnings.warn( + "`auto_scale` was set to True for this system. " + "Setting reference mass manually disables auto " + "scaling." + ) validated_mass = validate_ref_value(mass, u.dimensions.mass) self._reference_values["mass"] = validated_mass @@ -292,6 +310,13 @@ def reference_values(self, ref_value_dict): length, mass, and energy. """ + if self.auto_scale: + warnings.warn( + "`auto_scale` was set to True for this system. " + "Setting reference values manually disables auto " + "scaling." + ) + ref_keys = ["length", "mass", "energy"] for k in ref_keys: if k not in ref_value_dict.keys(): From 2805233b883bef48013453d5152a964b2449a0ab Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 25 Sep 2023 17:01:25 -0600 Subject: [PATCH 06/41] add unit test for invalid ff --- hoomd_organics/tests/base/test_system.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index d11dcaa3..5a0b5786 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -8,7 +8,7 @@ from hoomd_organics import Lattice, Pack from hoomd_organics.library import OPLS_AA, OPLS_AA_DIMETHYLETHER, OPLS_AA_PPS from hoomd_organics.tests import BaseTest -from hoomd_organics.utils.exceptions import ReferenceUnitError +from hoomd_organics.utils.exceptions import ForceFieldError, ReferenceUnitError class TestSystem(BaseTest): @@ -691,6 +691,24 @@ def test_forcefield_kwargs_attr(self, polyethylene): assert system._ff_kwargs["pppm_kwargs"]["resolution"] == (4, 4, 4) assert system._ff_kwargs["pppm_kwargs"]["order"] == 3 + def test_forcefield_list_hoomd_ff(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=1) + system = Pack( + molecules=[polyethylene], + force_field=[OPLS_AA()], + density=1.0, + ) + system.apply_forcefield(r_cut=2.5, auto_scale=True) + + hoomd_ff = system.hoomd_forcefield + + with pytest.raises(ForceFieldError): + Pack( + molecules=[polyethylene], + force_field=hoomd_ff, + density=1.0, + ) + def test_lattice_polymer(self, polyethylene): polyethylene = polyethylene(lengths=2, num_mols=32) system = Lattice( From 22ec7efc61f0fcdeb0ed73e82b564942a10d2ae0 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 25 Sep 2023 17:26:18 -0600 Subject: [PATCH 07/41] fix unit tests --- hoomd_organics/tests/base_test.py | 7 +++---- hoomd_organics/tests/library/test_tensile.py | 3 +-- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py index 3652dc34..14445056 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -218,7 +218,6 @@ def benzene_system(self, benzene_mb): molecules=[benzene], density=0.2, force_field=OPLS_AA(), - auto_scale=True, ) system.apply_forcefield(r_cut=2.5, auto_scale=True) return system @@ -230,7 +229,6 @@ def benzene_cg_system(self, benzene_molecule, benzene_smiles): system = Pack( molecules=[benzene_mols], density=0.5, - auto_scale=False, ) return system @@ -241,9 +239,10 @@ def polyethylene_system(self, polyethylene): molecules=polyethylene_mol, density=0.5, force_field=OPLS_AA(), - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5, remove_hydrogens=True) + system.apply_forcefield( + r_cut=2.5, remove_hydrogens=True, auto_scale=True + ) return system @pytest.fixture() diff --git a/hoomd_organics/tests/library/test_tensile.py b/hoomd_organics/tests/library/test_tensile.py index 3e48abb9..2c47193d 100644 --- a/hoomd_organics/tests/library/test_tensile.py +++ b/hoomd_organics/tests/library/test_tensile.py @@ -15,9 +15,8 @@ def test_tensile(self): x=1.2, y=1.2, n=4, - auto_scale=True, ) - system.apply_forcefield(r_cut=2.5) + system.apply_forcefield(r_cut=2.5, auto_scale=True) tensile_sim = Tensile( initial_state=system.hoomd_snapshot, From e5fb048533d94ac94f5d5e6b2360b26590b12858 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Tue, 26 Sep 2023 12:53:06 -0600 Subject: [PATCH 08/41] move forcefield logic to apply ff --- hoomd_organics/base/system.py | 57 ++++++++++++++++++----------------- 1 file changed, 29 insertions(+), 28 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index f11df2b0..964fa1a1 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -48,9 +48,6 @@ class System(ABC): target_box attribute. Can be useful when initializing systems at low density and running a shrink simulation to achieve a target density. - force_field : hoomd_organics.ForceField or a list of ForceField objects, - default=None - The force field to be applied to the system for parameterization. base_units : dict, default {} Dictionary of base units to use for scaling. Dictionary keys are "length", "mass", and "energy". Values should be an @@ -75,13 +72,9 @@ def __init__( self, molecules, density: float, - force_field=None, base_units=dict(), ): self._molecules = check_return_iterable(molecules) - self._force_field = None - if force_field: - self._force_field = check_return_iterable(force_field) self.density = density self.all_molecules = [] self.gmso_system = None @@ -128,26 +121,6 @@ def __init__( ) self.n_mol_types += 1 - # Collecting all force-fields only if xml force-field is provided - if self._force_field: - for i in range(self.n_mol_types): - if not self._gmso_forcefields_dict.get(str(i)): - if i < len(self._force_field): - # if there is a ff for each molecule type - ff_index = i - else: - # if there is only one ff for all molecule types - ff_index = 0 - if hasattr(self._force_field[ff_index], "gmso_ff"): - self._gmso_forcefields_dict[str(i)] = self._force_field[ - ff_index - ].gmso_ff - else: - raise ForceFieldError( - msg=f"GMSO Force field in " - f"{self._force_field[ff_index]} is not " - f"provided." - ) # Create mBuild system self.system = self._build_system() # Create GMSO topology @@ -463,6 +436,7 @@ def _create_hoomd_snapshot(self): def apply_forcefield( self, r_cut, + force_field=None, auto_scale=False, scale_charges=False, remove_charges=False, @@ -477,7 +451,11 @@ def apply_forcefield( ---------- r_cut : float The cutoff radius for the Lennard-Jones interactions. - + force_field : hoomd_organics.ForceField or a list of ForceField objects, + default=None + The force field to be applied to the system for parameterization. + If a list of force fields is provided, the length of the list must + be equal to the number of molecule types in the system. auto_scale : bool, default=False Set to true to use reduced simulation units. distance, mass, and energy are scaled by the largest value @@ -502,6 +480,29 @@ def apply_forcefield( Neighborlist buffer for simulation cell. """ + self._force_field = None + if force_field: + self._force_field = check_return_iterable(force_field) + # Collecting all force-fields only if xml force-field is provided + if self._force_field: + for i in range(self.n_mol_types): + if not self._gmso_forcefields_dict.get(str(i)): + if i < len(self._force_field): + # if there is a ff for each molecule type + ff_index = i + else: + # if there is only one ff for all molecule types + ff_index = 0 + if hasattr(self._force_field[ff_index], "gmso_ff"): + self._gmso_forcefields_dict[str(i)] = self._force_field[ + ff_index + ].gmso_ff + else: + raise ForceFieldError( + msg=f"GMSO Force field in " + f"{self._force_field[ff_index]} is not " + f"provided." + ) self.auto_scale = auto_scale if not self._force_field: # TODO: Better erorr message From 006e0a48c8566057be3df703f1eabc912d3deaf3 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Tue, 26 Sep 2023 13:11:41 -0600 Subject: [PATCH 09/41] remove forcefield from pack and lattice --- hoomd_organics/base/system.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 964fa1a1..50acd956 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -433,6 +433,12 @@ def _create_hoomd_snapshot(self): self._snap_refs = self._reference_values.copy() return snap + def _validate_forcefield(self, input_forcefield): + self._force_field = None + if input_forcefield: + self._force_field = check_return_iterable(input_forcefield) + # check if forcefield is xml based. return error if not + def apply_forcefield( self, r_cut, @@ -480,10 +486,13 @@ def apply_forcefield( Neighborlist buffer for simulation cell. """ - self._force_field = None - if force_field: - self._force_field = check_return_iterable(force_field) - # Collecting all force-fields only if xml force-field is provided + # make sure forcefield is xml based + self._validate_forcefield(force_field) + + # check if molecules already had ff in them. If yes, use that or + # return error (can't have two ff per molecule) + + # Collecting all force-fields into a dict if self._force_field: for i in range(self.n_mol_types): if not self._gmso_forcefields_dict.get(str(i)): @@ -659,7 +668,6 @@ def __init__( self, molecules, density: float, - force_field=None, base_units=dict(), packing_expand_factor=5, edge=0.2, @@ -669,7 +677,6 @@ def __init__( super(Pack, self).__init__( molecules=molecules, density=density, - force_field=force_field, base_units=base_units, ) @@ -711,7 +718,6 @@ def __init__( y: float, n: int, basis_vector=[0.5, 0.5, 0], - force_field=None, base_units=dict(), ): self.x = x @@ -721,7 +727,6 @@ def __init__( super(Lattice, self).__init__( molecules=molecules, density=density, - force_field=force_field, base_units=base_units, ) From 17a931ce2f2167d3ebd8c695733cc7f424bf9084 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Tue, 26 Sep 2023 13:12:25 -0600 Subject: [PATCH 10/41] update test units --- hoomd_organics/tests/base/test_system.py | 280 +++++++++++++---------- hoomd_organics/tests/base_test.py | 11 +- 2 files changed, 165 insertions(+), 126 deletions(-) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index 5a0b5786..8d6d58b4 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -14,12 +14,10 @@ class TestSystem(BaseTest): def test_single_mol_type(self, benzene_molecule): benzene_mols = benzene_molecule(n_mols=3) - system = Pack( - molecules=[benzene_mols], - density=0.8, - force_field=OPLS_AA(), + system = Pack(molecules=[benzene_mols], density=0.8) + system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(benzene_mols.molecules) assert system.gmso_system.is_typed() @@ -31,12 +29,10 @@ def test_single_mol_type(self, benzene_molecule): def test_multiple_mol_types(self, benzene_molecule, ethane_molecule): benzene_mol = benzene_molecule(n_mols=3) ethane_mol = ethane_molecule(n_mols=2) - system = Pack( - molecules=[benzene_mol, ethane_mol], - density=0.8, - force_field=OPLS_AA(), + system = Pack(molecules=[benzene_mol, ethane_mol], density=0.8) + system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 2 assert len(system.all_molecules) == len(benzene_mol.molecules) + len( ethane_mol.molecules @@ -58,12 +54,12 @@ def test_multiple_mol_types_different_ff( ): dimethylether_mol = dimethylether_molecule(n_mols=3) pps_mol = pps_molecule(n_mols=2) - system = Pack( - molecules=[dimethylether_mol, pps_mol], - density=0.8, + system = Pack(molecules=[dimethylether_mol, pps_mol], density=0.8) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA_DIMETHYLETHER(), OPLS_AA_PPS()], + auto_scale=True, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.n_mol_types == 2 assert len(system.all_molecules) == len( dimethylether_mol.molecules @@ -89,25 +85,20 @@ def test_multiple_mol_types_different_ff( def test_system_from_mol2_mol_parameterization(self, benzene_molecule_mol2): benzene_mol = benzene_molecule_mol2(n_mols=3) - system = Pack( - molecules=[benzene_mol], - density=0.8, - force_field=OPLS_AA(), + system = Pack(molecules=[benzene_mol], density=0.8) + system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) assert system.gmso_system.is_typed() assert len(system.hoomd_forcefield) > 0 assert system.n_particles == system.hoomd_snapshot.particles.N def test_remove_hydrogen(self, benzene_molecule): benzene_mol = benzene_molecule(n_mols=3) - system = Pack( - molecules=[benzene_mol], - density=0.8, - force_field=OPLS_AA(), - ) + system = Pack(molecules=[benzene_mol], density=0.8) system.apply_forcefield( r_cut=2.5, + force_field=OPLS_AA(), remove_hydrogens=True, auto_scale=True, ) @@ -130,10 +121,12 @@ def test_remove_hydrogen_no_atomic_num(self, benzene_molecule): system = Pack( molecules=[benzene_mol], density=0.8, - force_field=OPLS_AA(), ) system.apply_forcefield( - r_cut=2.5, remove_hydrogens=False, auto_scale=True + r_cut=2.5, + force_field=OPLS_AA(), + remove_hydrogens=False, + auto_scale=True, ) for site in system.gmso_system.sites: if site.name == "H": @@ -152,10 +145,12 @@ def test_remove_hydrogen_no_hydrogen(self, benzene_molecule): system = Pack( molecules=[benzene_mol], density=0.8, - force_field=OPLS_AA(), ) system.apply_forcefield( - r_cut=2.5, remove_hydrogens=False, auto_scale=True + r_cut=2.5, + force_field=OPLS_AA(), + remove_hydrogens=False, + auto_scale=True, ) hydrogens = [ site @@ -173,10 +168,10 @@ def test_add_mass_charges(self, benzene_molecule): system = Pack( molecules=[benzene_mols], density=0.8, - force_field=OPLS_AA(), ) system.apply_forcefield( r_cut=2.5, + force_field=OPLS_AA(), remove_hydrogens=True, scale_charges=False, auto_scale=False, @@ -196,16 +191,15 @@ def test_target_box(self, benzene_molecule): low_density_system = Pack( molecules=[benzene_mol], density=0.1, - force_field=OPLS_AA(), ) - low_density_system.apply_forcefield(r_cut=2.5, auto_scale=True) + low_density_system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True + ) - high_density_system = Pack( - molecules=[benzene_mol], - density=0.9, - force_field=OPLS_AA(), + high_density_system = Pack(molecules=[benzene_mol], density=0.9) + high_density_system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True ) - high_density_system.apply_forcefield(r_cut=2.5, auto_scale=True) assert all( low_density_system.target_box > high_density_system.target_box ) @@ -221,10 +215,11 @@ def test_ref_length(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=5) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) assert np.allclose( system.reference_length.to("angstrom").value, 3.5, atol=1e-3 @@ -239,10 +234,11 @@ def test_ref_mass(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=5) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) total_red_mass = sum(system.hoomd_snapshot.particles.mass) assert np.allclose( system.mass, @@ -254,10 +250,11 @@ def test_ref_energy(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=5) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) assert np.allclose( system.reference_energy.to("kcal/mol").value, 0.066, atol=1e-3 ) @@ -266,10 +263,11 @@ def test_rebuild_snapshot(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) assert system._snap_refs == system.reference_values assert system._ff_refs == system.reference_values init_snap = system.hoomd_snapshot @@ -287,10 +285,11 @@ def test_ref_values_no_autoscale(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_length = 1 * u.angstrom system.reference_energy = 1 * u.kcal / u.mol system.reference_mass = 1 * u.amu @@ -306,10 +305,10 @@ def test_ref_values_no_autoscale(self, polyethylene): def test_set_ref_values(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) - system = Pack( - molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0 + system = Pack(molecules=[polyethylene], density=1.0) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -324,10 +323,11 @@ def test_set_ref_values_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) ref_value_dict = { "length": "1 angstrom", "energy": "3 kcal/mol", @@ -342,10 +342,11 @@ def test_set_ref_values_missing_key(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -357,10 +358,11 @@ def test_set_ref_values_invalid_type(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=5) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -373,10 +375,11 @@ def test_set_ref_values_auto_scale_true(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=5) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -389,10 +392,11 @@ def test_set_ref_length(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_length = 1 * u.angstrom assert system.reference_length == 1 * u.angstrom @@ -400,10 +404,11 @@ def test_set_ref_length_invalid_type(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=5) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_length = 1.0 @@ -411,10 +416,11 @@ def test_ref_length_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_length = "1 angstrom" assert system.reference_length == 1 * u.angstrom @@ -422,10 +428,11 @@ def test_ref_length_invalid_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0" @@ -433,10 +440,11 @@ def test_ref_length_invalid_unit_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0 invalid_unit" @@ -444,10 +452,11 @@ def test_ref_length_invalid_dimension(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_length = 1.0 * u.g @@ -455,10 +464,11 @@ def test_ref_length_invalid_dimension_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0 g" @@ -466,10 +476,11 @@ def test_ref_length_auto_scale_true(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) system.reference_length = 1 * u.angstrom assert system.reference_length == 1 * u.angstrom @@ -477,10 +488,11 @@ def test_set_ref_energy(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_energy = 1 * u.kcal / u.mol assert system.reference_energy == 1 * u.kcal / u.mol @@ -488,10 +500,11 @@ def test_set_ref_energy_invalid_type(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 @@ -499,10 +512,11 @@ def test_ref_energy_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_energy = "1 kJ" assert system.reference_energy == 1 * u.kJ @@ -510,10 +524,11 @@ def test_ref_energy_string_comb_units(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_energy = "1 kcal/mol" assert system.reference_energy == 1 * u.kcal / u.mol @@ -521,10 +536,11 @@ def test_ref_energy_invalid_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0" @@ -532,10 +548,11 @@ def test_ref_energy_invalid_unit_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0 invalid_unit" @@ -543,10 +560,11 @@ def test_ref_energy_invalid_dimension(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 * u.g @@ -554,10 +572,11 @@ def test_ref_energy_invalid_dimension_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0 m" @@ -565,10 +584,11 @@ def test_set_ref_energy_auto_scale_true(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) system.reference_energy = 1 * u.kcal / u.mol assert system.reference_energy == 1 * u.kcal / u.mol @@ -576,10 +596,11 @@ def test_set_ref_mass(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_mass = 1.0 * u.amu assert system.reference_mass == 1.0 * u.amu @@ -588,10 +609,11 @@ def test_set_ref_mass_invalid_type(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_mass = 1.0 @@ -599,10 +621,11 @@ def test_ref_mass_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) system.reference_mass = "1 g" assert system.reference_mass == 1.0 * u.g @@ -610,10 +633,11 @@ def test_ref_mass_invalid_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0" @@ -621,10 +645,11 @@ def test_ref_mass_invalid_unit_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 invalid_unit" @@ -632,10 +657,11 @@ def test_ref_mass_invalid_dimension(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 * u.m @@ -643,10 +669,11 @@ def test_ref_mass_invalid_dimension_string(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False + ) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 m" @@ -654,10 +681,11 @@ def test_set_ref_mass_auto_scale_true(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) system.reference_mass = 1.0 * u.amu assert system.reference_mass == 1.0 * u.amu @@ -666,21 +694,22 @@ def test_apply_forcefield_no_forcefield(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=None, density=1.0, ) with pytest.raises(ValueError): - system.apply_forcefield(r_cut=2.5, auto_scale=False) + system.apply_forcefield( + r_cut=2.5, force_field=None, auto_scale=False + ) def test_forcefield_kwargs_attr(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) system.apply_forcefield( r_cut=2.5, + force_field=[OPLS_AA()], auto_scale=False, nlist_buffer=0.5, pppm_resolution=(4, 4, 4), @@ -695,31 +724,33 @@ def test_forcefield_list_hoomd_ff(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) hoomd_ff = system.hoomd_forcefield with pytest.raises(ForceFieldError): - Pack( + system = Pack( molecules=[polyethylene], - force_field=hoomd_ff, density=1.0, ) + system.apply_forcefield(r_cut=2.5, force_field=hoomd_ff) def test_lattice_polymer(self, polyethylene): polyethylene = polyethylene(lengths=2, num_mols=32) system = Lattice( molecules=[polyethylene], - force_field=[OPLS_AA()], density=1.0, x=1, y=1, n=4, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(polyethylene.molecules) @@ -732,13 +763,14 @@ def test_lattice_molecule(self, benzene_molecule): benzene_mol = benzene_molecule(n_mols=32) system = Lattice( molecules=[benzene_mol], - force_field=OPLS_AA(), density=1.0, x=1, y=1, n=4, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True + ) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(benzene_mol.molecules) assert len(system.hoomd_forcefield) > 0 @@ -750,19 +782,23 @@ def test_scale_charges(self, pps): no_scale = Pack( molecules=pps_mol, density=0.5, - force_field=OPLS_AA_PPS(), ) no_scale.apply_forcefield( - r_cut=2.5, scale_charges=False, auto_scale=True + r_cut=2.5, + force_field=OPLS_AA_PPS(), + scale_charges=False, + auto_scale=True, ) with_scale = Pack( molecules=pps_mol, density=0.5, - force_field=OPLS_AA_PPS(), ) with_scale.apply_forcefield( - r_cut=2.5, scale_charges=True, auto_scale=True + r_cut=2.5, + force_field=OPLS_AA_PPS(), + scale_charges=True, + auto_scale=True, ) assert abs(no_scale.net_charge.value) > abs(with_scale.net_charge.value) assert np.allclose(0, with_scale.net_charge.value, atol=1e-30) diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py index 14445056..d150fc63 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -217,9 +217,10 @@ def benzene_system(self, benzene_mb): system = Pack( molecules=[benzene], density=0.2, - force_field=OPLS_AA(), ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True + ) return system @pytest.fixture() @@ -238,10 +239,12 @@ def polyethylene_system(self, polyethylene): system = Pack( molecules=polyethylene_mol, density=0.5, - force_field=OPLS_AA(), ) system.apply_forcefield( - r_cut=2.5, remove_hydrogens=True, auto_scale=True + r_cut=2.5, + remove_hydrogens=True, + force_field=OPLS_AA(), + auto_scale=True, ) return system From 735eb7558eccc342ab0297adcd97605cb4070c90 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Tue, 26 Sep 2023 13:13:11 -0600 Subject: [PATCH 11/41] update tensile unit test --- hoomd_organics/tests/library/test_tensile.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/hoomd_organics/tests/library/test_tensile.py b/hoomd_organics/tests/library/test_tensile.py index 2c47193d..d7c4bb85 100644 --- a/hoomd_organics/tests/library/test_tensile.py +++ b/hoomd_organics/tests/library/test_tensile.py @@ -10,13 +10,14 @@ def test_tensile(self): pps = PPS(lengths=6, num_mols=32) system = Lattice( molecules=[pps], - force_field=[OPLS_AA_PPS()], density=1.0, x=1.2, y=1.2, n=4, ) - system.apply_forcefield(r_cut=2.5, auto_scale=True) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA_PPS()], auto_scale=True + ) tensile_sim = Tensile( initial_state=system.hoomd_snapshot, From c38968f7f4d9e33755d3a60ce0e52bc29bd02750 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Wed, 27 Sep 2023 15:18:08 -0600 Subject: [PATCH 12/41] assign mol names in apply ff --- hoomd_organics/base/system.py | 40 +++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 50acd956..4379cd3f 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -91,10 +91,14 @@ def __init__( # Collecting all molecules self.n_mol_types = 0 + self.mol_type_idx = [] for mol_item in self._molecules: if isinstance(mol_item, Molecule): - if self._force_field: - mol_item._assign_mol_name(str(self.n_mol_types)) + # extend mol_type_idx by number of particles in mol_item with + # a list with n_mol-types as values + self.mol_type_idx.extend( + [self.n_mol_types] * len(mol_item.n_particles) + ) self.all_molecules.extend(mol_item.molecules) # if ff is provided in Molecule class if mol_item.force_field: @@ -434,11 +438,18 @@ def _create_hoomd_snapshot(self): return snap def _validate_forcefield(self, input_forcefield): - self._force_field = None if input_forcefield: - self._force_field = check_return_iterable(input_forcefield) + return check_return_iterable(input_forcefield) + return None # check if forcefield is xml based. return error if not + def _assign_particle_idx(self): + system_particles = list(self.system.particles()) + # assign particle idx to each particles based on self._mol_type_idx in + # batches of similar idx + for i, mol_idx in enumerate(self._mol_type_idx): + system_particles[i].name = mol_idx + def apply_forcefield( self, r_cut, @@ -487,38 +498,45 @@ def apply_forcefield( """ # make sure forcefield is xml based - self._validate_forcefield(force_field) + _force_field = self._validate_forcefield(force_field) # check if molecules already had ff in them. If yes, use that or # return error (can't have two ff per molecule) # Collecting all force-fields into a dict - if self._force_field: + if _force_field: for i in range(self.n_mol_types): if not self._gmso_forcefields_dict.get(str(i)): - if i < len(self._force_field): + if i < len(_force_field): # if there is a ff for each molecule type ff_index = i else: # if there is only one ff for all molecule types ff_index = 0 - if hasattr(self._force_field[ff_index], "gmso_ff"): - self._gmso_forcefields_dict[str(i)] = self._force_field[ + if hasattr(_force_field[ff_index], "gmso_ff"): + self._gmso_forcefields_dict[str(i)] = _force_field[ ff_index ].gmso_ff else: raise ForceFieldError( msg=f"GMSO Force field in " - f"{self._force_field[ff_index]} is not " + f"{_force_field[ff_index]} is not " f"provided." ) self.auto_scale = auto_scale - if not self._force_field: + if not _force_field: # TODO: Better erorr message raise ValueError( "This method can only be used when the System is " "initialized with an XML type forcefield." ) + if self._gmso_forcefields_dict: + # assign names to all the particles of system based on mol_type to + # match the keys in self._gmso_forcefields_dict and recreate the + # gmso system object + self.gmso_system = self._convert_to_gmso( + self._assign_particle_idx() + ) self.gmso_system = apply( self.gmso_system, self._gmso_forcefields_dict, From f2b916be852ae4c2bac30d5ab062c2ca1ac85b12 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Wed, 27 Sep 2023 15:41:12 -0600 Subject: [PATCH 13/41] change mol name type to str and assign it to parents too --- hoomd_organics/base/system.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 4379cd3f..afb1d7d7 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -91,13 +91,13 @@ def __init__( # Collecting all molecules self.n_mol_types = 0 - self.mol_type_idx = [] + self._mol_type_idx = [] for mol_item in self._molecules: if isinstance(mol_item, Molecule): # extend mol_type_idx by number of particles in mol_item with # a list with n_mol-types as values - self.mol_type_idx.extend( - [self.n_mol_types] * len(mol_item.n_particles) + self._mol_type_idx.extend( + [self.n_mol_types] * mol_item.n_particles ) self.all_molecules.extend(mol_item.molecules) # if ff is provided in Molecule class @@ -448,7 +448,9 @@ def _assign_particle_idx(self): # assign particle idx to each particles based on self._mol_type_idx in # batches of similar idx for i, mol_idx in enumerate(self._mol_type_idx): - system_particles[i].name = mol_idx + system_particles[i].name = str(mol_idx) + if system_particles[i].parent: + system_particles[i].parent.name = str(mol_idx) def apply_forcefield( self, @@ -534,9 +536,8 @@ def apply_forcefield( # assign names to all the particles of system based on mol_type to # match the keys in self._gmso_forcefields_dict and recreate the # gmso system object - self.gmso_system = self._convert_to_gmso( - self._assign_particle_idx() - ) + self._assign_particle_idx() + self.gmso_system = self._convert_to_gmso() self.gmso_system = apply( self.gmso_system, self._gmso_forcefields_dict, From 435f8e3cca76e0b1e00e6589f77efc47d5307e48 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Wed, 27 Sep 2023 16:10:32 -0600 Subject: [PATCH 14/41] assign mol type idx to gmso sites --- hoomd_organics/base/system.py | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index afb1d7d7..0166fb2a 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -444,13 +444,23 @@ def _validate_forcefield(self, input_forcefield): # check if forcefield is xml based. return error if not def _assign_particle_idx(self): - system_particles = list(self.system.particles()) - # assign particle idx to each particles based on self._mol_type_idx in - # batches of similar idx - for i, mol_idx in enumerate(self._mol_type_idx): - system_particles[i].name = str(mol_idx) - if system_particles[i].parent: - system_particles[i].parent.name = str(mol_idx) + # system_particles = list(self.system.particles()) + # # assign particle idx to each particles based on self._mol_type_idx in + # # batches of similar idx + # + # #TODO: This is a hacky way to assign particle idx. Need to find a + # # better way to assign names to + # for i, mol_idx in enumerate(self._mol_type_idx): + # system_particles[i].name = str(mol_idx) + # if system_particles[i].parent: + # system_particles[i].parent.name = str(mol_idx) + # if system_particles[i].root: + # system_particles[i].root.name = str(mol_idx) + # if system_particles[i].root.children: + # for child in system_particles[i].root.children: + # child.name = str(mol_idx) + for i, site in enumerate(self.gmso_system.sites): + site.group = str(self._mol_type_idx[i]) def apply_forcefield( self, @@ -537,10 +547,11 @@ def apply_forcefield( # match the keys in self._gmso_forcefields_dict and recreate the # gmso system object self._assign_particle_idx() - self.gmso_system = self._convert_to_gmso() + # self.gmso_system = self._convert_to_gmso() self.gmso_system = apply( self.gmso_system, self._gmso_forcefields_dict, + match_ff_by="group", identify_connections=True, speedup_by_moltag=True, speedup_by_molgraph=False, From 38a26e833ac392d7433f431d9dbcccdb119642fd Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Wed, 27 Sep 2023 16:16:11 -0600 Subject: [PATCH 15/41] check if _gmso_forcefields_dict exists when getting hoomd forces --- hoomd_organics/base/system.py | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 0166fb2a..d14de545 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -312,7 +312,10 @@ def hoomd_snapshot(self): @property def hoomd_forcefield(self): """List of HOOMD forces.""" - if self._ff_refs != self.reference_values and self._force_field: + if ( + self._ff_refs != self.reference_values + and self._gmso_forcefields_dict + ): self._hoomd_forcefield = self._create_hoomd_forcefield( **self._ff_kwargs ) @@ -443,22 +446,8 @@ def _validate_forcefield(self, input_forcefield): return None # check if forcefield is xml based. return error if not - def _assign_particle_idx(self): - # system_particles = list(self.system.particles()) - # # assign particle idx to each particles based on self._mol_type_idx in - # # batches of similar idx - # - # #TODO: This is a hacky way to assign particle idx. Need to find a - # # better way to assign names to - # for i, mol_idx in enumerate(self._mol_type_idx): - # system_particles[i].name = str(mol_idx) - # if system_particles[i].parent: - # system_particles[i].parent.name = str(mol_idx) - # if system_particles[i].root: - # system_particles[i].root.name = str(mol_idx) - # if system_particles[i].root.children: - # for child in system_particles[i].root.children: - # child.name = str(mol_idx) + def _assign_site_mol_type_idx(self): + """Assign molecule type index to the gmso sites.""" for i, site in enumerate(self.gmso_system.sites): site.group = str(self._mol_type_idx[i]) @@ -546,7 +535,7 @@ def apply_forcefield( # assign names to all the particles of system based on mol_type to # match the keys in self._gmso_forcefields_dict and recreate the # gmso system object - self._assign_particle_idx() + self._assign_site_mol_type_idx() # self.gmso_system = self._convert_to_gmso() self.gmso_system = apply( self.gmso_system, From a47a1a28b50630bbad23a1fa3c9f82dc9b7fed3e Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Wed, 27 Sep 2023 16:22:30 -0600 Subject: [PATCH 16/41] update unit tests based on changes in system --- hoomd_organics/tests/base/test_system.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index 8d6d58b4..e9642aea 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -37,8 +37,8 @@ def test_multiple_mol_types(self, benzene_molecule, ethane_molecule): assert len(system.all_molecules) == len(benzene_mol.molecules) + len( ethane_mol.molecules ) - assert system.all_molecules[0].name == "0" - assert system.all_molecules[-1].name == "1" + assert system.gmso_system.sites[0].group == "0" + assert system.gmso_system.sites[-1].group == "1" assert system.gmso_system.is_typed() assert len(system.hoomd_forcefield) > 0 assert system.n_particles == system.hoomd_snapshot.particles.N @@ -64,8 +64,8 @@ def test_multiple_mol_types_different_ff( assert len(system.all_molecules) == len( dimethylether_mol.molecules ) + len(pps_mol.molecules) - assert system.all_molecules[0].name == "0" - assert system.all_molecules[-1].name == "1" + assert system.gmso_system.sites[0].group == "0" + assert system.gmso_system.sites[-1].group == "1" assert system.gmso_system.is_typed() assert len(system.hoomd_forcefield) > 0 for hoomd_force in system.hoomd_forcefield: @@ -129,6 +129,8 @@ def test_remove_hydrogen_no_atomic_num(self, benzene_molecule): auto_scale=True, ) for site in system.gmso_system.sites: + # this does not work as the name of all sites are changed to reflect + # the molecule type id if site.name == "H": site.element = gmso.core.element.Element( symbol="C", From 90baf1f0a2007b33cfac0c3db0390d19214f5d7d Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 10:58:56 -0600 Subject: [PATCH 17/41] add comments and remove commented lines --- hoomd_organics/base/system.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index d14de545..ab22718a 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -94,8 +94,8 @@ def __init__( self._mol_type_idx = [] for mol_item in self._molecules: if isinstance(mol_item, Molecule): - # extend mol_type_idx by number of particles in mol_item with - # a list with n_mol-types as values + # keep track of molecule types indices to assign to sites + # before applying forcefield self._mol_type_idx.extend( [self.n_mol_types] * mol_item.n_particles ) @@ -498,13 +498,14 @@ def apply_forcefield( Neighborlist buffer for simulation cell. """ + self.auto_scale = auto_scale # make sure forcefield is xml based _force_field = self._validate_forcefield(force_field) # check if molecules already had ff in them. If yes, use that or # return error (can't have two ff per molecule) - # Collecting all force-fields into a dict + # Collecting all force-fields into a dict with mol_type index as key if _force_field: for i in range(self.n_mol_types): if not self._gmso_forcefields_dict.get(str(i)): @@ -524,7 +525,6 @@ def apply_forcefield( f"{_force_field[ff_index]} is not " f"provided." ) - self.auto_scale = auto_scale if not _force_field: # TODO: Better erorr message raise ValueError( @@ -532,11 +532,9 @@ def apply_forcefield( "initialized with an XML type forcefield." ) if self._gmso_forcefields_dict: - # assign names to all the particles of system based on mol_type to - # match the keys in self._gmso_forcefields_dict and recreate the - # gmso system object + # assign names to all the gmso sites based on mol_type to + # match the keys in self._gmso_forcefields_dict before applying ff self._assign_site_mol_type_idx() - # self.gmso_system = self._convert_to_gmso() self.gmso_system = apply( self.gmso_system, self._gmso_forcefields_dict, From 0ee5a6e8b19f7763f7dbc77d35847385d68c3c31 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 11:08:22 -0600 Subject: [PATCH 18/41] use hoomd_organics ff in molecule class --- hoomd_organics/base/molecule.py | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/hoomd_organics/base/molecule.py b/hoomd_organics/base/molecule.py index d544a29f..24bfbb11 100644 --- a/hoomd_organics/base/molecule.py +++ b/hoomd_organics/base/molecule.py @@ -7,17 +7,13 @@ import mbuild as mb from gmso.core.topology import Topology from gmso.external.convert_mbuild import from_mbuild, to_mbuild +from gmso.parameterization import apply from grits import CG_Compound from mbuild.lib.recipes import Polymer as mbPolymer from hoomd_organics.utils import check_return_iterable -from hoomd_organics.utils.base_types import FF_Types from hoomd_organics.utils.exceptions import MoleculeLoadError -from hoomd_organics.utils.ff_utils import ( - _validate_hoomd_ff, - apply_xml_ff, - find_xml_ff, -) +from hoomd_organics.utils.ff_utils import _validate_hoomd_ff class Molecule: @@ -32,11 +28,12 @@ class Molecule: ---------- num_mols : int, required Number of molecules to generate. - force_field : str, default None - The force field to apply to the molecule. - Note that setting force_field will not apply the forcefield to the - molecule. The forcefield in this step is just used for validation - purposes. + force_field : hoomd_organics.ForceField or a list of ForceField objects, + default=None + The force field to be applied to the system for parameterization. + Note that setting `force_field` will not apply the forcefield to the + molecule. The forcefield in this step is just used for validation + purposes. smiles : str, default None The smiles string of the molecule to generate. file : str, default None @@ -391,16 +388,18 @@ def _identify_topology_information(self, gmso_molecule): def _validate_force_field(self): """Validate the force field for the molecule.""" - self.ff_type = None - if isinstance(self.force_field, str): - ff_xml_path, ff_type = find_xml_ff(self.force_field) - self.ff_type = ff_type - self.gmso_molecule = apply_xml_ff(ff_xml_path, self.gmso_molecule) + if hasattr(self.force_field, "gmso_ff"): + self.gmso_molecule = apply( + self.gmso_molecule, + self.force_field.gmso_ff, + identify_connections=True, + speedup_by_moltag=True, + speedup_by_molgraph=False, + ) # Update topology information from typed gmso after applying ff. self._identify_topology_information(self.gmso_molecule) elif isinstance(self.force_field, List): _validate_hoomd_ff(self.force_field, self.topology_information) - self.ff_type = FF_Types.Hoomd def _assign_mol_name(self, name): """Assign a name to the molecule.""" From cff987cba4cbc572541726169291c649bffc502d Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 11:09:00 -0600 Subject: [PATCH 19/41] remove assign_mol_names from molecule --- hoomd_organics/base/molecule.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/hoomd_organics/base/molecule.py b/hoomd_organics/base/molecule.py index 24bfbb11..e8037f87 100644 --- a/hoomd_organics/base/molecule.py +++ b/hoomd_organics/base/molecule.py @@ -401,17 +401,6 @@ def _validate_force_field(self): elif isinstance(self.force_field, List): _validate_hoomd_ff(self.force_field, self.topology_information) - def _assign_mol_name(self, name): - """Assign a name to the molecule.""" - for mol in self.molecules: - mol.name = name - # TODO: This is a hack to make sure that the name of the children is - # also updated, so that when converting to gmso, all the sites have - # the correct name. This needs additional investigation into gmso's - # convert mbuilder to gmso functionality. - for child in mol.children: - child.name = name - class Polymer(Molecule): """Builds a polymer from a monomer. From dfdc8f1d438f37526ed7a4c7166a76335ac953df Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 11:38:53 -0600 Subject: [PATCH 20/41] validate forcefield types in molecule --- hoomd_organics/base/molecule.py | 47 ++++++++++++++++++++++++--------- 1 file changed, 35 insertions(+), 12 deletions(-) diff --git a/hoomd_organics/base/molecule.py b/hoomd_organics/base/molecule.py index e8037f87..0cc53209 100644 --- a/hoomd_organics/base/molecule.py +++ b/hoomd_organics/base/molecule.py @@ -11,8 +11,8 @@ from grits import CG_Compound from mbuild.lib.recipes import Polymer as mbPolymer -from hoomd_organics.utils import check_return_iterable -from hoomd_organics.utils.exceptions import MoleculeLoadError +from hoomd_organics.utils import FF_Types, check_return_iterable +from hoomd_organics.utils.exceptions import ForceFieldError, MoleculeLoadError from hoomd_organics.utils.ff_utils import _validate_hoomd_ff @@ -388,18 +388,41 @@ def _identify_topology_information(self, gmso_molecule): def _validate_force_field(self): """Validate the force field for the molecule.""" - if hasattr(self.force_field, "gmso_ff"): - self.gmso_molecule = apply( - self.gmso_molecule, - self.force_field.gmso_ff, - identify_connections=True, - speedup_by_moltag=True, - speedup_by_molgraph=False, - ) - # Update topology information from typed gmso after applying ff. - self._identify_topology_information(self.gmso_molecule) + if hasattr(self.force_field, "ff_type"): + if self.force_field.ff_type == FF_Types.XML and hasattr( + self.force_field, "gmso_ff" + ): + self.gmso_molecule = apply( + self.gmso_molecule, + self.force_field.gmso_ff, + identify_connections=True, + speedup_by_moltag=True, + speedup_by_molgraph=False, + ) + # Update topology information from typed gmso after applying ff. + self._identify_topology_information(self.gmso_molecule) + elif self.force_field.ff_type == FF_Types.HOOMD and hasattr( + self.force_field, "hoomd_forcefield" + ): + _validate_hoomd_ff( + self.force_field.hoomd_forcefield, self.topology_information + ) + else: + raise ForceFieldError( + "Unsupported forcefield type. Please " + "check forcefield classes from " + "`hoomd_orgnaics.library.forcefields` " + "for supported forcefields." + ) elif isinstance(self.force_field, List): _validate_hoomd_ff(self.force_field, self.topology_information) + else: + raise ForceFieldError( + "Unsupported forcefield type. Forcefields " + "should be of type " + "`hoomd_organics.ForceField` or a list of" + " `hoomd.md.force.Force` objects." + ) class Polymer(Molecule): From a7635e70008ea9f39e9e6e8f27143e30a969bd47 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 11:39:32 -0600 Subject: [PATCH 21/41] assign types to forcefields in library --- hoomd_organics/library/forcefields.py | 8 ++++++++ hoomd_organics/utils/base_types.py | 8 ++------ 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/hoomd_organics/library/forcefields.py b/hoomd_organics/library/forcefields.py index 28645358..f07a41bf 100644 --- a/hoomd_organics/library/forcefields.py +++ b/hoomd_organics/library/forcefields.py @@ -6,6 +6,7 @@ import hoomd from hoomd_organics.assets import FF_DIR +from hoomd_organics.utils import FF_Types class GAFF(foyer.Forcefield): @@ -18,6 +19,7 @@ def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml"): "The XML file was obtained from the antefoyer package: " "https://github.com/rsdefever/antefoyer/tree/master/antefoyer" ) + self.ff_type = FF_Types.XML self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() @@ -27,6 +29,7 @@ class OPLS_AA(foyer.Forcefield): def __init__(self, name="oplsaa"): super(OPLS_AA, self).__init__(name=name) self.description = "opls-aa forcefield found in the Foyer package." + self.ff_type = FF_Types.XML self.gmso_ff = ffutils.FoyerFFs().load(name).to_gmso_ff() @@ -44,6 +47,7 @@ def __init__(self, forcefield_files=f"{FF_DIR}/pps_opls.xml"): "experimental PPS papers. The spring constant taken " "from the equivalent angle in GAFF." ) + self.ff_type = FF_Types.XML self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() @@ -56,6 +60,7 @@ def __init__(self, forcefield_files=f"{FF_DIR}/benzene_opls.xml"): "Based on hoomd_organics.forcefields.OPLS_AA. " "Trimmed down to include only benzene parameters." ) + self.ff_type = FF_Types.XML self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() @@ -70,6 +75,7 @@ def __init__(self, forcefield_files=f"{FF_DIR}/dimethylether_opls.xml"): "Based on hoomd_organics.forcefields.OPLS_AA. " "Trimmed down to include only dimethyl ether parameters." ) + self.ff_type = FF_Types.XML self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() @@ -78,6 +84,7 @@ class FF_from_file(foyer.Forcefield): def __init__(self, xml_file): super(FF_from_file, self).__init__(forcefield_files=xml_file) + self.ff_type = FF_Types.XML self.gmso_ff = ffutils.FoyerFFs().load(xml_file).to_gmso_ff() @@ -147,6 +154,7 @@ def __init__( self.r_cut = r_cut self.exclusions = exclusions self.hoomd_forcefield = self._create_forcefield() + self.ff_type = FF_Types.HOOMD def _create_forcefield(self): """Create the hoomd force objects.""" diff --git a/hoomd_organics/utils/base_types.py b/hoomd_organics/utils/base_types.py index 8bef2646..291f91a1 100644 --- a/hoomd_organics/utils/base_types.py +++ b/hoomd_organics/utils/base_types.py @@ -2,12 +2,8 @@ class FF_Types: - opls = "opls" - pps_opls = "pps_opls" - oplsaa = "oplsaa" - gaff = "gaff" - custom = "custom" - Hoomd = "Hoomd" + XML = "XML" + HOOMD = "HOOMD" class HOOMDThermostats: From 078ddf984740828f48f9efd4539cb0b64e125241 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 11:42:54 -0600 Subject: [PATCH 22/41] update molecule unit tests --- hoomd_organics/tests/base/test_molecule.py | 46 ++++++++-------------- 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/hoomd_organics/tests/base/test_molecule.py b/hoomd_organics/tests/base/test_molecule.py index c25f6c27..4c181ea3 100644 --- a/hoomd_organics/tests/base/test_molecule.py +++ b/hoomd_organics/tests/base/test_molecule.py @@ -1,6 +1,7 @@ import pytest from hoomd_organics import CoPolymer, Molecule, Polymer +from hoomd_organics.library import OPLS_AA, FF_from_file from hoomd_organics.tests import BaseTest from hoomd_organics.utils import FF_Types, exceptions @@ -41,20 +42,10 @@ def test_molecule_topology_benzene(self, benzene_mb): def test_validate_force_field_oplsaa(self, benzene_mb): molecule = Molecule( - num_mols=2, force_field="oplsaa", compound=benzene_mb + num_mols=2, force_field=OPLS_AA(), compound=benzene_mb ) - assert molecule.ff_type == FF_Types.oplsaa - assert set(molecule.topology_information["particle_types"]) == { - "opls_145", - "opls_146", - } - assert any(molecule.topology_information["particle_charge"]) - - def test_validate_force_field_xml_file(self, benzene_mb): - molecule = Molecule( - num_mols=2, force_field="oplsaa.xml", compound=benzene_mb - ) - assert molecule.ff_type == FF_Types.oplsaa + assert molecule.force_field.ff_type == FF_Types.XML + assert molecule.gmso_molecule.is_typed() assert set(molecule.topology_information["particle_types"]) == { "opls_145", "opls_146", @@ -63,27 +54,18 @@ def test_validate_force_field_xml_file(self, benzene_mb): def test_validate_force_field_xml_file_path(self, benzene_mb, benzene_xml): molecule = Molecule( - num_mols=2, force_field=benzene_xml, compound=benzene_mb + num_mols=2, + force_field=FF_from_file(benzene_xml), + compound=benzene_mb, ) - assert molecule.ff_type == FF_Types.custom + assert molecule.force_field.ff_type == FF_Types.XML + assert molecule.gmso_molecule.is_typed() assert set(molecule.topology_information["particle_types"]) == { "opls_145", "opls_146", } assert any(molecule.topology_information["particle_charge"]) - def test_validate_force_field_not_xml_file(self, benzene_mb): - with pytest.raises(ValueError): - Molecule(num_mols=2, force_field="oplsaa.txt", compound=benzene_mb) - - def test_validate_force_field_not_supported(self, benzene_mb): - with pytest.raises(ValueError): - Molecule(num_mols=2, force_field="oplsaa2", compound=benzene_mb) - - def test_validate_force_field_invalid_xml_file(self, benzene_mb): - with pytest.raises(ValueError): - Molecule(num_mols=2, force_field="oplsaa2.xml", compound=benzene_mb) - def test_validate_force_field_hoomd_ff_aa( self, benzene_mb, benzene_hoomd_ff ): @@ -91,7 +73,8 @@ def test_validate_force_field_hoomd_ff_aa( molecule = Molecule( num_mols=2, force_field=hoomd_ff, compound=benzene_mb ) - assert molecule.ff_type == FF_Types.Hoomd + assert molecule.force_field == hoomd_ff + assert not molecule.gmso_molecule.is_typed() def test_validate_fore_field_hoomd_ff_ua( self, benzene_mb, benzene_hoomd_ff @@ -100,7 +83,8 @@ def test_validate_fore_field_hoomd_ff_ua( molecule = Molecule( num_mols=2, force_field=hoomd_ff, compound=benzene_mb ) - assert molecule.ff_type == FF_Types.Hoomd + assert molecule.force_field == hoomd_ff + assert not molecule.gmso_molecule.is_typed() def test_validate_force_field_hoomd_ff_missing_pair( self, benzene_mb, benzene_hoomd_ff @@ -130,7 +114,9 @@ def test_validate_force_field_hoomd_ff_missing_Coulomb( ): hoomd_ff = benzene_hoomd_ff(include_hydrogen=True) typed_molecule = Molecule( - num_mols=2, force_field=benzene_xml, compound=benzene_mb + num_mols=2, + force_field=FF_from_file(benzene_xml), + compound=benzene_mb, ) with pytest.raises(exceptions.MissingCoulombPotentialError): Molecule( From efa2c309b06a01b605a88d6e5cd87a5bae89d545 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 11:46:54 -0600 Subject: [PATCH 23/41] add unit test for invalid ff type --- hoomd_organics/tests/base/test_molecule.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/hoomd_organics/tests/base/test_molecule.py b/hoomd_organics/tests/base/test_molecule.py index 4c181ea3..62f4fea1 100644 --- a/hoomd_organics/tests/base/test_molecule.py +++ b/hoomd_organics/tests/base/test_molecule.py @@ -125,6 +125,14 @@ def test_validate_force_field_hoomd_ff_missing_Coulomb( compound=typed_molecule.gmso_molecule, ) + def test_validate_force_field_invalid_ff_type(self, benzene_mb): + with pytest.raises(exceptions.ForceFieldError): + Molecule( + num_mols=2, + force_field="invalid_ff.xml", + compound=benzene_mb, + ) + def test_coarse_grain_with_single_beads(self, benzene_smiles): molecule = Molecule(num_mols=2, smiles=benzene_smiles) molecule.coarse_grain(beads={"A": benzene_smiles}) From 1cc6f549d891cf7e10d9118399752175eb140665 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 11:55:14 -0600 Subject: [PATCH 24/41] define a base class for xml forcefields --- hoomd_organics/library/forcefields.py | 38 ++++++++++++++------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/hoomd_organics/library/forcefields.py b/hoomd_organics/library/forcefields.py index f07a41bf..60562089 100644 --- a/hoomd_organics/library/forcefields.py +++ b/hoomd_organics/library/forcefields.py @@ -9,7 +9,18 @@ from hoomd_organics.utils import FF_Types -class GAFF(foyer.Forcefield): +class BaseXMLForcefield(foyer.Forcefield): + """Base XML forcefield class.""" + + def __init__(self, forcefield_files=None, name=None): + super(BaseXMLForcefield, self).__init__( + forcefield_files=forcefield_files, name=name + ) + self.ff_type = FF_Types.XML + self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() + + +class GAFF(BaseXMLForcefield): """GAFF forcefield class.""" def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml"): @@ -19,21 +30,17 @@ def __init__(self, forcefield_files=f"{FF_DIR}/gaff.xml"): "The XML file was obtained from the antefoyer package: " "https://github.com/rsdefever/antefoyer/tree/master/antefoyer" ) - self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() -class OPLS_AA(foyer.Forcefield): +class OPLS_AA(BaseXMLForcefield): """OPLS All Atom forcefield class.""" def __init__(self, name="oplsaa"): super(OPLS_AA, self).__init__(name=name) self.description = "opls-aa forcefield found in the Foyer package." - self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(name).to_gmso_ff() -class OPLS_AA_PPS(foyer.Forcefield): +class OPLS_AA_PPS(BaseXMLForcefield): """OPLS All Atom for PPS molecule forcefield class.""" def __init__(self, forcefield_files=f"{FF_DIR}/pps_opls.xml"): @@ -47,11 +54,9 @@ def __init__(self, forcefield_files=f"{FF_DIR}/pps_opls.xml"): "experimental PPS papers. The spring constant taken " "from the equivalent angle in GAFF." ) - self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() -class OPLS_AA_BENZENE(foyer.Forcefield): +class OPLS_AA_BENZENE(BaseXMLForcefield): """OPLS All Atom for benzene molecule forcefield class.""" def __init__(self, forcefield_files=f"{FF_DIR}/benzene_opls.xml"): @@ -60,11 +65,9 @@ def __init__(self, forcefield_files=f"{FF_DIR}/benzene_opls.xml"): "Based on hoomd_organics.forcefields.OPLS_AA. " "Trimmed down to include only benzene parameters." ) - self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() -class OPLS_AA_DIMETHYLETHER(foyer.Forcefield): +class OPLS_AA_DIMETHYLETHER(BaseXMLForcefield): """OPLS All Atom for dimethyl ether molecule forcefield class.""" def __init__(self, forcefield_files=f"{FF_DIR}/dimethylether_opls.xml"): @@ -79,13 +82,12 @@ def __init__(self, forcefield_files=f"{FF_DIR}/dimethylether_opls.xml"): self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() -class FF_from_file(foyer.Forcefield): +class FF_from_file(BaseXMLForcefield): """Forcefield class for loading a forcefield from an XML file.""" - def __init__(self, xml_file): - super(FF_from_file, self).__init__(forcefield_files=xml_file) - self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(xml_file).to_gmso_ff() + def __init__(self, forcefield_files): + super(FF_from_file, self).__init__(forcefield_files=forcefield_files) + self.description = "Forcefield loaded from an XML file. " class BeadSpring: From c27df87f0c963768a158990ca24ea1b63a03f00c Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 12:07:40 -0600 Subject: [PATCH 25/41] fix minor bugs in forcefields --- hoomd_organics/library/forcefields.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/hoomd_organics/library/forcefields.py b/hoomd_organics/library/forcefields.py index 60562089..e7e45e99 100644 --- a/hoomd_organics/library/forcefields.py +++ b/hoomd_organics/library/forcefields.py @@ -20,6 +20,15 @@ def __init__(self, forcefield_files=None, name=None): self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() +class BaseHOOMDForcefield: + """Base HOOMD forcefield class.""" + + def __init__(self, hoomd_forces=None): + super(BaseHOOMDForcefield, self).__init__() + self.hoomd_forces = hoomd_forces + self.ff_type = FF_Types.HOOMD + + class GAFF(BaseXMLForcefield): """GAFF forcefield class.""" @@ -90,7 +99,7 @@ def __init__(self, forcefield_files): self.description = "Forcefield loaded from an XML file. " -class BeadSpring: +class BeadSpring(BaseHOOMDForcefield): """Bead-spring forcefield class. Given a dictionary of bead types, this class creates a list @@ -155,8 +164,7 @@ def __init__( self.dihedrals = dihedrals self.r_cut = r_cut self.exclusions = exclusions - self.hoomd_forcefield = self._create_forcefield() - self.ff_type = FF_Types.HOOMD + self.hoomd_forces = self._create_forcefield() def _create_forcefield(self): """Create the hoomd force objects.""" From a7620b585667a2c43aa458f19a53f53a4dfafc37 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 12:45:37 -0600 Subject: [PATCH 26/41] define base forcefield classes --- hoomd_organics/base/__init__.py | 1 + hoomd_organics/base/forcefield.py | 30 ++++++++++++++++++++++++++++++ hoomd_organics/library/__init__.py | 2 ++ 3 files changed, 33 insertions(+) create mode 100644 hoomd_organics/base/forcefield.py diff --git a/hoomd_organics/base/__init__.py b/hoomd_organics/base/__init__.py index 2c588be7..91dd661e 100644 --- a/hoomd_organics/base/__init__.py +++ b/hoomd_organics/base/__init__.py @@ -1,4 +1,5 @@ """Base classes for hoomd_organics.""" +from .forcefield import BaseHOOMDForcefield, BaseXMLForcefield from .molecule import CoPolymer, Molecule, Polymer from .simulation import Simulation from .system import Lattice, Pack, System diff --git a/hoomd_organics/base/forcefield.py b/hoomd_organics/base/forcefield.py new file mode 100644 index 00000000..69dd0863 --- /dev/null +++ b/hoomd_organics/base/forcefield.py @@ -0,0 +1,30 @@ +"""Base forcefield classes.""" +import forcefield_utilities as ffutils +import foyer + +from hoomd_organics.utils import FF_Types + + +class BaseXMLForcefield(foyer.Forcefield): + """Base XML forcefield class.""" + + def __init__(self, forcefield_files=None, name=None): + super(BaseXMLForcefield, self).__init__( + forcefield_files=forcefield_files, name=name + ) + self.ff_type = FF_Types.XML + self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() + + +class BaseHOOMDForcefield: + """Base HOOMD forcefield class.""" + + def __init__(self, hoomd_forces): + self.ff_type = FF_Types.HOOMD + self.hoomd_forces = hoomd_forces + if hoomd_forces is None: + raise NotImplementedError( + "`hoomd_forces` must be defined in the subclass." + ) + if not isinstance(hoomd_forces, list): + raise TypeError("`hoomd_forces` must be a list.") diff --git a/hoomd_organics/library/__init__.py b/hoomd_organics/library/__init__.py index dac0acd3..4f730899 100644 --- a/hoomd_organics/library/__init__.py +++ b/hoomd_organics/library/__init__.py @@ -5,6 +5,8 @@ OPLS_AA_BENZENE, OPLS_AA_DIMETHYLETHER, OPLS_AA_PPS, + BaseHOOMDForcefield, + BaseXMLForcefield, BeadSpring, FF_from_file, ) From 12445fa8116d69d432a0fb04e74ace6c2f2a4bc4 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 12:56:35 -0600 Subject: [PATCH 27/41] use name or forcefield_files when creating gmso_ff --- hoomd_organics/base/forcefield.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/hoomd_organics/base/forcefield.py b/hoomd_organics/base/forcefield.py index 69dd0863..752fe3a0 100644 --- a/hoomd_organics/base/forcefield.py +++ b/hoomd_organics/base/forcefield.py @@ -13,7 +13,9 @@ def __init__(self, forcefield_files=None, name=None): forcefield_files=forcefield_files, name=name ) self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() + self.gmso_ff = ( + ffutils.FoyerFFs().load(forcefield_files or name).to_gmso_ff() + ) class BaseHOOMDForcefield: From 5910ddb57bb8c3861ba0ef634fd67335bd387fc8 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 12:57:15 -0600 Subject: [PATCH 28/41] update validate_ff --- hoomd_organics/base/molecule.py | 50 +++++++++++++-------------------- 1 file changed, 20 insertions(+), 30 deletions(-) diff --git a/hoomd_organics/base/molecule.py b/hoomd_organics/base/molecule.py index 0cc53209..63e19886 100644 --- a/hoomd_organics/base/molecule.py +++ b/hoomd_organics/base/molecule.py @@ -11,7 +11,8 @@ from grits import CG_Compound from mbuild.lib.recipes import Polymer as mbPolymer -from hoomd_organics.utils import FF_Types, check_return_iterable +from hoomd_organics.base import BaseHOOMDForcefield, BaseXMLForcefield +from hoomd_organics.utils import check_return_iterable from hoomd_organics.utils.exceptions import ForceFieldError, MoleculeLoadError from hoomd_organics.utils.ff_utils import _validate_hoomd_ff @@ -388,40 +389,29 @@ def _identify_topology_information(self, gmso_molecule): def _validate_force_field(self): """Validate the force field for the molecule.""" - if hasattr(self.force_field, "ff_type"): - if self.force_field.ff_type == FF_Types.XML and hasattr( - self.force_field, "gmso_ff" - ): - self.gmso_molecule = apply( - self.gmso_molecule, - self.force_field.gmso_ff, - identify_connections=True, - speedup_by_moltag=True, - speedup_by_molgraph=False, - ) - # Update topology information from typed gmso after applying ff. - self._identify_topology_information(self.gmso_molecule) - elif self.force_field.ff_type == FF_Types.HOOMD and hasattr( - self.force_field, "hoomd_forcefield" - ): - _validate_hoomd_ff( - self.force_field.hoomd_forcefield, self.topology_information - ) - else: - raise ForceFieldError( - "Unsupported forcefield type. Please " - "check forcefield classes from " - "`hoomd_orgnaics.library.forcefields` " - "for supported forcefields." - ) + if isinstance(self.force_field, BaseXMLForcefield): + self.gmso_molecule = apply( + self.gmso_molecule, + self.force_field.gmso_ff, + identify_connections=True, + speedup_by_moltag=True, + speedup_by_molgraph=False, + ) + # Update topology information from typed gmso after applying ff. + self._identify_topology_information(self.gmso_molecule) + elif isinstance(self.force_field, BaseHOOMDForcefield): + _validate_hoomd_ff(self.force_field, self.topology_information) elif isinstance(self.force_field, List): _validate_hoomd_ff(self.force_field, self.topology_information) else: raise ForceFieldError( "Unsupported forcefield type. Forcefields " - "should be of type " - "`hoomd_organics.ForceField` or a list of" - " `hoomd.md.force.Force` objects." + "should be a subclass of " + "`hoomd_organics.base.forcefield.BaseXMLForcefield` or " + "`hoomd_organics.base.forcefield.BaseHOOMDForcefield` or a " + "list of `hoomd.md.force.Force` objects. \n" + "Please check `hoomd_organics.library.forcefields` for " + "examples of supported forcefields." ) From a59856cb5fb459874073bbcb32dde15b9e22027e Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 12:57:58 -0600 Subject: [PATCH 29/41] change to hoomd_forces variable --- .../tests/library/test_forcefields.py | 38 +++++++++---------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/hoomd_organics/tests/library/test_forcefields.py b/hoomd_organics/tests/library/test_forcefields.py index 7537a032..5f5c64c2 100644 --- a/hoomd_organics/tests/library/test_forcefields.py +++ b/hoomd_organics/tests/library/test_forcefields.py @@ -55,37 +55,37 @@ def test_BeadSpring(self): dihedrals={"A-A-A-A": dict(phi0=0.0, k=100, d=-1, n=1)}, ) - assert isinstance(ff.hoomd_forcefield[0], hoomd.md.pair.pair.LJ) - assert isinstance(ff.hoomd_forcefield[1], hoomd.md.bond.Harmonic) - assert isinstance(ff.hoomd_forcefield[2], hoomd.md.angle.Harmonic) - assert isinstance(ff.hoomd_forcefield[3], hoomd.md.dihedral.Periodic) + assert isinstance(ff.hoomd_forces[0], hoomd.md.pair.pair.LJ) + assert isinstance(ff.hoomd_forces[1], hoomd.md.bond.Harmonic) + assert isinstance(ff.hoomd_forces[2], hoomd.md.angle.Harmonic) + assert isinstance(ff.hoomd_forces[3], hoomd.md.dihedral.Periodic) pair_types = [("A", "A"), ("A", "B"), ("B", "B")] - for param in ff.hoomd_forcefield[0].params: + for param in ff.hoomd_forces[0].params: assert param in pair_types if param == ("A", "A"): - assert ff.hoomd_forcefield[0].params[param]["sigma"] == 1.0 + assert ff.hoomd_forces[0].params[param]["sigma"] == 1.0 if param == ("B", "B"): - assert ff.hoomd_forcefield[0].params[param]["epsilon"] == 2.0 + assert ff.hoomd_forces[0].params[param]["epsilon"] == 2.0 if param == ("A", "B"): - assert ff.hoomd_forcefield[0].params[param]["epsilon"] == 1.5 + assert ff.hoomd_forces[0].params[param]["epsilon"] == 1.5 bond_types = [("A-A"), ("A-B")] - for param in ff.hoomd_forcefield[1].params: + for param in ff.hoomd_forces[1].params: assert param in bond_types - assert ff.hoomd_forcefield[1].params[param]["r0"] == 1.1 - assert ff.hoomd_forcefield[1].params[param]["k"] == 300 + assert ff.hoomd_forces[1].params[param]["r0"] == 1.1 + assert ff.hoomd_forces[1].params[param]["k"] == 300 angle_types = [("A-A-A"), ("A-B-A")] - for param in ff.hoomd_forcefield[2].params: + for param in ff.hoomd_forces[2].params: assert param in angle_types - assert ff.hoomd_forcefield[2].params[param]["t0"] == 2.0 - assert ff.hoomd_forcefield[2].params[param]["k"] == 200 + assert ff.hoomd_forces[2].params[param]["t0"] == 2.0 + assert ff.hoomd_forces[2].params[param]["k"] == 200 dihedral_types = [("A-A-A-A")] - for param in ff.hoomd_forcefield[3].params: + for param in ff.hoomd_forces[3].params: assert param in dihedral_types - assert ff.hoomd_forcefield[3].params[param]["phi0"] == 0.0 - assert ff.hoomd_forcefield[3].params[param]["k"] == 100 - assert ff.hoomd_forcefield[3].params[param]["d"] == -1 - assert ff.hoomd_forcefield[3].params[param]["n"] == 1 + assert ff.hoomd_forces[3].params[param]["phi0"] == 0.0 + assert ff.hoomd_forces[3].params[param]["k"] == 100 + assert ff.hoomd_forces[3].params[param]["d"] == -1 + assert ff.hoomd_forces[3].params[param]["n"] == 1 From 39dce43b202c9b0eb034f952a05bf4d759b99d28 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 13:02:17 -0600 Subject: [PATCH 30/41] check molecule forcefield in system init --- hoomd_organics/base/system.py | 18 ++++++++++++------ hoomd_organics/library/forcefields.py | 25 +++---------------------- 2 files changed, 15 insertions(+), 28 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index ab22718a..c10397e3 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -10,13 +10,15 @@ from gmso.external import from_mbuild, to_gsd_snapshot, to_hoomd_forcefield from gmso.parameterization import apply +from hoomd_organics.base.forcefield import ( + BaseHOOMDForcefield, + BaseXMLForcefield, +) from hoomd_organics.base.molecule import Molecule from hoomd_organics.utils import ( - FF_Types, calculate_box_length, check_return_iterable, validate_ref_value, - xml_to_gmso_ff, ) from hoomd_organics.utils.exceptions import ForceFieldError, MoleculeLoadError @@ -102,12 +104,16 @@ def __init__( self.all_molecules.extend(mol_item.molecules) # if ff is provided in Molecule class if mol_item.force_field: - if mol_item.ff_type == FF_Types.Hoomd: - self._hoomd_forcefield.extend(mol_item.force_field) - else: + if isinstance(mol_item.force_field, BaseHOOMDForcefield): + self._hoomd_forcefield.extend( + mol_item.force_field.hoomd_forces + ) + elif isinstance(mol_item.force_field, BaseXMLForcefield): self._gmso_forcefields_dict[ str(self.n_mol_types) - ] = xml_to_gmso_ff(mol_item.force_field) + ] = mol_item.force_field.gmso_ff + elif isinstance(mol_item.force_field, list): + self._hoomd_forcefield.extend(mol_item.force_field) self.n_mol_types += 1 elif isinstance(mol_item, mb.Compound): mol_item.name = str(self.n_mol_types) diff --git a/hoomd_organics/library/forcefields.py b/hoomd_organics/library/forcefields.py index e7e45e99..0f6432fe 100644 --- a/hoomd_organics/library/forcefields.py +++ b/hoomd_organics/library/forcefields.py @@ -2,33 +2,13 @@ import itertools import forcefield_utilities as ffutils -import foyer import hoomd from hoomd_organics.assets import FF_DIR +from hoomd_organics.base import BaseHOOMDForcefield, BaseXMLForcefield from hoomd_organics.utils import FF_Types -class BaseXMLForcefield(foyer.Forcefield): - """Base XML forcefield class.""" - - def __init__(self, forcefield_files=None, name=None): - super(BaseXMLForcefield, self).__init__( - forcefield_files=forcefield_files, name=name - ) - self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() - - -class BaseHOOMDForcefield: - """Base HOOMD forcefield class.""" - - def __init__(self, hoomd_forces=None): - super(BaseHOOMDForcefield, self).__init__() - self.hoomd_forces = hoomd_forces - self.ff_type = FF_Types.HOOMD - - class GAFF(BaseXMLForcefield): """GAFF forcefield class.""" @@ -164,7 +144,8 @@ def __init__( self.dihedrals = dihedrals self.r_cut = r_cut self.exclusions = exclusions - self.hoomd_forces = self._create_forcefield() + hoomd_forces = self._create_forcefield() + super(BeadSpring, self).__init__(hoomd_forces) def _create_forcefield(self): """Create the hoomd force objects.""" From 1205c2b8a28d5b9dcbc66de1f43ede5a1476e8c8 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 13:17:48 -0600 Subject: [PATCH 31/41] add more checks when validating ff in system --- hoomd_organics/base/system.py | 79 +++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 35 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index c10397e3..4c70e0f5 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -102,7 +102,7 @@ def __init__( [self.n_mol_types] * mol_item.n_particles ) self.all_molecules.extend(mol_item.molecules) - # if ff is provided in Molecule class + # if ff is provided in the Molecule class use that as the ff if mol_item.force_field: if isinstance(mol_item.force_field, BaseHOOMDForcefield): self._hoomd_forcefield.extend( @@ -447,10 +447,49 @@ def _create_hoomd_snapshot(self): return snap def _validate_forcefield(self, input_forcefield): + if input_forcefield is None and not self._gmso_forcefields_dict: + raise ForceFieldError( + "Forcefield is not provided. Valid forcefield " + "must be provided either during Molecule " + "initialization or when calling the " + "`apply_forcefield` method of the System " + "class." + ) + + if input_forcefield and self._gmso_forcefields_dict: + raise ForceFieldError( + "Forcefield is provided both during Molecule " + "initialization and when calling the " + "`apply_forcefield` method of the System " + "class. Please provide the forcefield only " + "once." + ) + + if input_forcefield and not isinstance( + input_forcefield, (BaseHOOMDForcefield, BaseXMLForcefield) + ): + raise ForceFieldError( + "Forcefield must be an instance of either " + " `BaseHOOMDForcefield` or " + "`BaseXMLForcefield`. \n" + "Please check " + "`hoomd_organics.library.forcefields` for " + "examples of supported forcefields." + ) if input_forcefield: - return check_return_iterable(input_forcefield) - return None - # check if forcefield is xml based. return error if not + _force_field = check_return_iterable(input_forcefield) + # Collecting all force-fields into a dict with mol_type index as key + for i in range(self.n_mol_types): + if not self._gmso_forcefields_dict.get(str(i)): + if i < len(_force_field): + # if there is a ff for each molecule type + ff_index = i + else: + # if there is only one ff for all molecule types + ff_index = 0 + self._gmso_forcefields_dict[str(i)] = _force_field[ + ff_index + ].gmso_ff def _assign_site_mol_type_idx(self): """Assign molecule type index to the gmso sites.""" @@ -505,38 +544,8 @@ def apply_forcefield( """ self.auto_scale = auto_scale - # make sure forcefield is xml based - _force_field = self._validate_forcefield(force_field) + self._validate_forcefield(force_field) - # check if molecules already had ff in them. If yes, use that or - # return error (can't have two ff per molecule) - - # Collecting all force-fields into a dict with mol_type index as key - if _force_field: - for i in range(self.n_mol_types): - if not self._gmso_forcefields_dict.get(str(i)): - if i < len(_force_field): - # if there is a ff for each molecule type - ff_index = i - else: - # if there is only one ff for all molecule types - ff_index = 0 - if hasattr(_force_field[ff_index], "gmso_ff"): - self._gmso_forcefields_dict[str(i)] = _force_field[ - ff_index - ].gmso_ff - else: - raise ForceFieldError( - msg=f"GMSO Force field in " - f"{_force_field[ff_index]} is not " - f"provided." - ) - if not _force_field: - # TODO: Better erorr message - raise ValueError( - "This method can only be used when the System is " - "initialized with an XML type forcefield." - ) if self._gmso_forcefields_dict: # assign names to all the gmso sites based on mol_type to # match the keys in self._gmso_forcefields_dict before applying ff From dd2bd9b124327313aa10f5826ae1746236045352 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 13:40:14 -0600 Subject: [PATCH 32/41] check type of all items in ff list --- hoomd_organics/base/system.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 4c70e0f5..4e77d29d 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -465,19 +465,20 @@ def _validate_forcefield(self, input_forcefield): "once." ) - if input_forcefield and not isinstance( - input_forcefield, (BaseHOOMDForcefield, BaseXMLForcefield) - ): - raise ForceFieldError( - "Forcefield must be an instance of either " - " `BaseHOOMDForcefield` or " - "`BaseXMLForcefield`. \n" - "Please check " - "`hoomd_organics.library.forcefields` for " - "examples of supported forcefields." - ) if input_forcefield: _force_field = check_return_iterable(input_forcefield) + if not all( + isinstance(ff, (BaseHOOMDForcefield, BaseXMLForcefield)) + for ff in _force_field + ): + raise ForceFieldError( + "Forcefield must be an instance of either " + " `BaseHOOMDForcefield` or " + "`BaseXMLForcefield`. \n" + "Please check " + "`hoomd_organics.library.forcefields` for " + "examples of supported forcefields." + ) # Collecting all force-fields into a dict with mol_type index as key for i in range(self.n_mol_types): if not self._gmso_forcefields_dict.get(str(i)): From 42cfc3f1916618ddf5ddf5f80a89c067ab1f5048 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 14:14:57 -0600 Subject: [PATCH 33/41] fix bugs in unit tests --- hoomd_organics/tests/base/test_system.py | 2 +- hoomd_organics/tests/base_test.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index e9642aea..7f9a3803 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -698,7 +698,7 @@ def test_apply_forcefield_no_forcefield(self, polyethylene): molecules=[polyethylene], density=1.0, ) - with pytest.raises(ValueError): + with pytest.raises(ForceFieldError): system.apply_forcefield( r_cut=2.5, force_field=None, auto_scale=False ) diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py index d150fc63..aa76484f 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -264,4 +264,4 @@ def cg_single_bead_ff(self): "A": dict(epsilon=1.0, sigma=1.0), }, ) - return ff.hoomd_forcefield + return ff.hoomd_forces From e5534f474bbc3ee5b37b4811e09488c9c9db9eec Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 14:15:18 -0600 Subject: [PATCH 34/41] use correct ff_type --- hoomd_organics/utils/ff_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hoomd_organics/utils/ff_utils.py b/hoomd_organics/utils/ff_utils.py index 6f80c73a..1845939e 100644 --- a/hoomd_organics/utils/ff_utils.py +++ b/hoomd_organics/utils/ff_utils.py @@ -33,7 +33,7 @@ def find_xml_ff(ff_source): if not ff_source.endswith(".xml"): raise ValueError("ForceField file type must be XML.") ff_xml_path = ff_source - ff_type = FF_Types.custom + ff_type = FF_Types.XML elif not xml_directory.get(ff_source.split(".xml")[0]): raise ValueError( "{} forcefield is not supported. Supported XML forcefields " @@ -42,7 +42,7 @@ def find_xml_ff(ff_source): else: ff_key = ff_source.split(".xml")[0] ff_xml_path = xml_directory.get(ff_key) - ff_type = getattr(FF_Types, ff_key) + ff_type = FF_Types.XML return ff_xml_path, ff_type From 17ea9aac229ed72651b8755baa7cb878b900cd72 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 14:28:06 -0600 Subject: [PATCH 35/41] add tests for base forcefield classes --- hoomd_organics/tests/base/test_forcefield.py | 56 ++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 hoomd_organics/tests/base/test_forcefield.py diff --git a/hoomd_organics/tests/base/test_forcefield.py b/hoomd_organics/tests/base/test_forcefield.py new file mode 100644 index 00000000..aa8613d6 --- /dev/null +++ b/hoomd_organics/tests/base/test_forcefield.py @@ -0,0 +1,56 @@ +import pytest + +from hoomd_organics.base import BaseHOOMDForcefield, BaseXMLForcefield +from hoomd_organics.tests import BaseTest +from hoomd_organics.utils import FF_Types + + +class TestBaseForcefield(BaseTest): + def test_base_xml_forcefield(self, benzene_xml): + base_xml_ff = BaseXMLForcefield(forcefield_files=benzene_xml) + assert base_xml_ff.gmso_ff is not None + assert base_xml_ff.ff_type == FF_Types.XML + + def test_base_xml_forcefield_no_files(self): + with pytest.raises(TypeError): + BaseXMLForcefield(forcefield_files=None) + + def test_base_xml_forcefield_name(self): + base_xml_ff = BaseXMLForcefield(name="oplsaa") + assert base_xml_ff.gmso_ff is not None + assert base_xml_ff.ff_type == FF_Types.XML + + def test_base_xml_forcefield_invalid_name(self): + with pytest.raises(Exception): + BaseXMLForcefield(name="invalid_name") + + def test_base_xml_forcefield_invalid_files(self): + with pytest.raises(Exception): + BaseXMLForcefield(forcefield_files="invalid_files") + + def test_base_hoomd_forcefield(self): + class TestHOOMDFF(BaseHOOMDForcefield): + def __init__(self): + hoomd_forces = [] + super().__init__(hoomd_forces=hoomd_forces) + + test_hoomd_ff = TestHOOMDFF() + assert test_hoomd_ff.hoomd_forces == [] + assert test_hoomd_ff.ff_type == FF_Types.HOOMD + + def test_base_hoomd_forcefield_no_forces(self): + class TestHOOMDFF(BaseHOOMDForcefield): + def __init__(self): + super().__init__(hoomd_forces=None) + + with pytest.raises(NotImplementedError): + TestHOOMDFF() + + def test_base_hoomd_forcefield_invalid_forces_type(self): + class TestHOOMDFF(BaseHOOMDForcefield): + def __init__(self): + hoomd_forces = "invalid_type" + super().__init__(hoomd_forces=hoomd_forces) + + with pytest.raises(TypeError): + TestHOOMDFF() From 93127ca993f8abe4e39f666eb87768ad7a2cddda Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 14:42:41 -0600 Subject: [PATCH 36/41] add new unit tests for validate_forcefield --- hoomd_organics/tests/base/test_system.py | 46 ++++++++++++++++++++++++ hoomd_organics/tests/base_test.py | 6 ++-- 2 files changed, 50 insertions(+), 2 deletions(-) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index 7f9a3803..ef772363 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -703,6 +703,52 @@ def test_apply_forcefield_no_forcefield(self, polyethylene): r_cut=2.5, force_field=None, auto_scale=False ) + def test_apply_forcefield_no_forcefield_w_mol_ff(self, benzene_molecule): + benzene_mol = benzene_molecule(n_mols=3, force_field=OPLS_AA()) + system = Pack( + molecules=[benzene_mol], + density=1.0, + ) + system.apply_forcefield(r_cut=2.5, auto_scale=True) + assert system.gmso_system.is_typed() + assert len(system.hoomd_forcefield) > 0 + assert system.hoomd_snapshot.particles.N == benzene_mol.n_particles + + def test_apply_forcefield_w_mol_ff(self, benzene_molecule): + benzene_mol = benzene_molecule(n_mols=3, force_field=OPLS_AA()) + system = Pack( + molecules=[benzene_mol], + density=1.0, + ) + with pytest.raises(ForceFieldError): + system.apply_forcefield( + r_cut=2.5, force_field=OPLS_AA(), auto_scale=True + ) + + def test_validate_forcefield_invalid_ff_type(self, benzene_molecule): + benzene_mol = benzene_molecule(n_mols=1) + system = Pack( + molecules=[benzene_mol], + density=1.0, + ) + with pytest.raises(ForceFieldError): + system.apply_forcefield( + r_cut=2.5, force_field="invalid_ff.xml", auto_scale=True + ) + + def test_validate_forcefield_mult_ff_invalid_type( + self, dimethylether_molecule, pps_molecule + ): + dimethylether_mol = dimethylether_molecule(n_mols=3) + pps_mol = pps_molecule(n_mols=2) + system = Pack(molecules=[dimethylether_mol, pps_mol], density=0.8) + with pytest.raises(ForceFieldError): + system.apply_forcefield( + r_cut=2.5, + force_field=["invalid_ff", OPLS_AA_PPS()], + auto_scale=True, + ) + def test_forcefield_kwargs_attr(self, polyethylene): polyethylene = polyethylene(lengths=5, num_mols=1) system = Pack( diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py index aa76484f..87c3659d 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -111,8 +111,10 @@ def _hoomd_ff(include_hydrogen, invalid_pair=False): @pytest.fixture() def benzene_molecule(self, benzene_smiles): - def _benzene_molecule(n_mols): - benzene = Molecule(num_mols=n_mols, smiles=benzene_smiles) + def _benzene_molecule(n_mols, force_field=None): + benzene = Molecule( + num_mols=n_mols, smiles=benzene_smiles, force_field=force_field + ) return benzene return _benzene_molecule From 7f0d87f26d05ef0e994f6aece96ef749ca28e405 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 15:02:35 -0600 Subject: [PATCH 37/41] add unit tests for ff_utils --- hoomd_organics/tests/utils/test_ff_utils.py | 33 ++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/hoomd_organics/tests/utils/test_ff_utils.py b/hoomd_organics/tests/utils/test_ff_utils.py index 24b296ce..8f0dbafa 100644 --- a/hoomd_organics/tests/utils/test_ff_utils.py +++ b/hoomd_organics/tests/utils/test_ff_utils.py @@ -1,10 +1,41 @@ +import os + +import pytest from gmso.core.forcefield import ForceField from hoomd_organics.tests import BaseTest -from hoomd_organics.utils import xml_to_gmso_ff +from hoomd_organics.utils.base_types import FF_Types +from hoomd_organics.utils.ff_utils import find_xml_ff, xml_to_gmso_ff class TestFFUtils(BaseTest): + def test_find_xml_ff(self): + ff_xml_path, ff_type = find_xml_ff("oplsaa.xml") + assert ff_type == FF_Types.XML + assert os.path.exists(ff_xml_path) + + def test_find_xml_only_file_name(self): + ff_xml_path, ff_type = find_xml_ff("oplsaa") + assert ff_type == FF_Types.XML + assert os.path.exists(ff_xml_path) + + def test_find_xml_ff_path(self, benzene_xml): + ff_xml_path, ff_type = find_xml_ff(benzene_xml) + assert ff_type == FF_Types.XML + assert os.path.exists(ff_xml_path) + + def test_find_xml_invalid_extension(self): + with pytest.raises(ValueError): + find_xml_ff("oplsaa.txt") + + def test_find_xml_not_supported_name(self): + with pytest.raises(ValueError): + find_xml_ff("oplsaa2") + + def test_find_xml_not_supported_path(self): + with pytest.raises(ValueError): + find_xml_ff("oplsaa2.xml") + def test_xml_to_gmso_ff(self, benzene_xml): gmso_ff = xml_to_gmso_ff(benzene_xml) assert isinstance(gmso_ff, ForceField) From e46fc921298f03f854c0544bb7b3997eb8a4cb57 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 15:36:43 -0600 Subject: [PATCH 38/41] add test for initiating molecule with beadspring ff --- hoomd_organics/base/molecule.py | 4 +++- hoomd_organics/tests/base/test_molecule.py | 19 ++++++++++++++++++- 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/hoomd_organics/base/molecule.py b/hoomd_organics/base/molecule.py index 63e19886..bad4f20e 100644 --- a/hoomd_organics/base/molecule.py +++ b/hoomd_organics/base/molecule.py @@ -400,7 +400,9 @@ def _validate_force_field(self): # Update topology information from typed gmso after applying ff. self._identify_topology_information(self.gmso_molecule) elif isinstance(self.force_field, BaseHOOMDForcefield): - _validate_hoomd_ff(self.force_field, self.topology_information) + _validate_hoomd_ff( + self.force_field.hoomd_forces, self.topology_information + ) elif isinstance(self.force_field, List): _validate_hoomd_ff(self.force_field, self.topology_information) else: diff --git a/hoomd_organics/tests/base/test_molecule.py b/hoomd_organics/tests/base/test_molecule.py index 62f4fea1..975d3cc6 100644 --- a/hoomd_organics/tests/base/test_molecule.py +++ b/hoomd_organics/tests/base/test_molecule.py @@ -1,7 +1,7 @@ import pytest from hoomd_organics import CoPolymer, Molecule, Polymer -from hoomd_organics.library import OPLS_AA, FF_from_file +from hoomd_organics.library import OPLS_AA, BeadSpring, FF_from_file from hoomd_organics.tests import BaseTest from hoomd_organics.utils import FF_Types, exceptions @@ -133,6 +133,23 @@ def test_validate_force_field_invalid_ff_type(self, benzene_mb): compound=benzene_mb, ) + def test_validate_forcefield_hoomd_ff(self, benzene_smiles): + molecule = Molecule(num_mols=1, smiles=benzene_smiles) + molecule.coarse_grain(beads={"A": benzene_smiles}) + beadspring_ff = BeadSpring( + r_cut=2.5, + beads={ + "A": dict(epsilon=1.0, sigma=1.0), + }, + ) + cg_molecule = Molecule( + num_mols=1, + force_field=beadspring_ff, + compound=molecule.gmso_molecule, + ) + assert cg_molecule.force_field == beadspring_ff + assert not molecule.gmso_molecule.is_typed() + def test_coarse_grain_with_single_beads(self, benzene_smiles): molecule = Molecule(num_mols=2, smiles=benzene_smiles) molecule.coarse_grain(beads={"A": benzene_smiles}) From 831e7f498705115edde0984ba820b8929dc1c828 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 15:37:01 -0600 Subject: [PATCH 39/41] add unit test for ff_utils --- hoomd_organics/tests/utils/test_ff_utils.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/hoomd_organics/tests/utils/test_ff_utils.py b/hoomd_organics/tests/utils/test_ff_utils.py index 8f0dbafa..035d0276 100644 --- a/hoomd_organics/tests/utils/test_ff_utils.py +++ b/hoomd_organics/tests/utils/test_ff_utils.py @@ -28,6 +28,10 @@ def test_find_xml_invalid_extension(self): with pytest.raises(ValueError): find_xml_ff("oplsaa.txt") + def test_find_xml_invalid_file(self, benzene_mol2): + with pytest.raises(ValueError): + find_xml_ff(benzene_mol2) + def test_find_xml_not_supported_name(self): with pytest.raises(ValueError): find_xml_ff("oplsaa2") From 94a65a7d93c5808f42f1c8ad08f6e362b81f1170 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 28 Sep 2023 15:45:26 -0600 Subject: [PATCH 40/41] remove FF_type class --- hoomd_organics/base/forcefield.py | 4 ---- hoomd_organics/library/forcefields.py | 4 ---- hoomd_organics/tests/base/test_forcefield.py | 4 ---- hoomd_organics/tests/base/test_molecule.py | 4 +--- hoomd_organics/tests/utils/test_ff_utils.py | 10 +++------- hoomd_organics/utils/__init__.py | 2 +- hoomd_organics/utils/base_types.py | 5 ----- hoomd_organics/utils/ff_utils.py | 7 ++----- 8 files changed, 7 insertions(+), 33 deletions(-) diff --git a/hoomd_organics/base/forcefield.py b/hoomd_organics/base/forcefield.py index 752fe3a0..c416ea4c 100644 --- a/hoomd_organics/base/forcefield.py +++ b/hoomd_organics/base/forcefield.py @@ -2,8 +2,6 @@ import forcefield_utilities as ffutils import foyer -from hoomd_organics.utils import FF_Types - class BaseXMLForcefield(foyer.Forcefield): """Base XML forcefield class.""" @@ -12,7 +10,6 @@ def __init__(self, forcefield_files=None, name=None): super(BaseXMLForcefield, self).__init__( forcefield_files=forcefield_files, name=name ) - self.ff_type = FF_Types.XML self.gmso_ff = ( ffutils.FoyerFFs().load(forcefield_files or name).to_gmso_ff() ) @@ -22,7 +19,6 @@ class BaseHOOMDForcefield: """Base HOOMD forcefield class.""" def __init__(self, hoomd_forces): - self.ff_type = FF_Types.HOOMD self.hoomd_forces = hoomd_forces if hoomd_forces is None: raise NotImplementedError( diff --git a/hoomd_organics/library/forcefields.py b/hoomd_organics/library/forcefields.py index 0f6432fe..196628ec 100644 --- a/hoomd_organics/library/forcefields.py +++ b/hoomd_organics/library/forcefields.py @@ -1,12 +1,10 @@ """All pre-defined forcefield classes for use in hoomd_organics.""" import itertools -import forcefield_utilities as ffutils import hoomd from hoomd_organics.assets import FF_DIR from hoomd_organics.base import BaseHOOMDForcefield, BaseXMLForcefield -from hoomd_organics.utils import FF_Types class GAFF(BaseXMLForcefield): @@ -67,8 +65,6 @@ def __init__(self, forcefield_files=f"{FF_DIR}/dimethylether_opls.xml"): "Based on hoomd_organics.forcefields.OPLS_AA. " "Trimmed down to include only dimethyl ether parameters." ) - self.ff_type = FF_Types.XML - self.gmso_ff = ffutils.FoyerFFs().load(forcefield_files).to_gmso_ff() class FF_from_file(BaseXMLForcefield): diff --git a/hoomd_organics/tests/base/test_forcefield.py b/hoomd_organics/tests/base/test_forcefield.py index aa8613d6..0c0d4d2b 100644 --- a/hoomd_organics/tests/base/test_forcefield.py +++ b/hoomd_organics/tests/base/test_forcefield.py @@ -2,14 +2,12 @@ from hoomd_organics.base import BaseHOOMDForcefield, BaseXMLForcefield from hoomd_organics.tests import BaseTest -from hoomd_organics.utils import FF_Types class TestBaseForcefield(BaseTest): def test_base_xml_forcefield(self, benzene_xml): base_xml_ff = BaseXMLForcefield(forcefield_files=benzene_xml) assert base_xml_ff.gmso_ff is not None - assert base_xml_ff.ff_type == FF_Types.XML def test_base_xml_forcefield_no_files(self): with pytest.raises(TypeError): @@ -18,7 +16,6 @@ def test_base_xml_forcefield_no_files(self): def test_base_xml_forcefield_name(self): base_xml_ff = BaseXMLForcefield(name="oplsaa") assert base_xml_ff.gmso_ff is not None - assert base_xml_ff.ff_type == FF_Types.XML def test_base_xml_forcefield_invalid_name(self): with pytest.raises(Exception): @@ -36,7 +33,6 @@ def __init__(self): test_hoomd_ff = TestHOOMDFF() assert test_hoomd_ff.hoomd_forces == [] - assert test_hoomd_ff.ff_type == FF_Types.HOOMD def test_base_hoomd_forcefield_no_forces(self): class TestHOOMDFF(BaseHOOMDForcefield): diff --git a/hoomd_organics/tests/base/test_molecule.py b/hoomd_organics/tests/base/test_molecule.py index 975d3cc6..04de1d1b 100644 --- a/hoomd_organics/tests/base/test_molecule.py +++ b/hoomd_organics/tests/base/test_molecule.py @@ -3,7 +3,7 @@ from hoomd_organics import CoPolymer, Molecule, Polymer from hoomd_organics.library import OPLS_AA, BeadSpring, FF_from_file from hoomd_organics.tests import BaseTest -from hoomd_organics.utils import FF_Types, exceptions +from hoomd_organics.utils import exceptions class TestMolecule(BaseTest): @@ -44,7 +44,6 @@ def test_validate_force_field_oplsaa(self, benzene_mb): molecule = Molecule( num_mols=2, force_field=OPLS_AA(), compound=benzene_mb ) - assert molecule.force_field.ff_type == FF_Types.XML assert molecule.gmso_molecule.is_typed() assert set(molecule.topology_information["particle_types"]) == { "opls_145", @@ -58,7 +57,6 @@ def test_validate_force_field_xml_file_path(self, benzene_mb, benzene_xml): force_field=FF_from_file(benzene_xml), compound=benzene_mb, ) - assert molecule.force_field.ff_type == FF_Types.XML assert molecule.gmso_molecule.is_typed() assert set(molecule.topology_information["particle_types"]) == { "opls_145", diff --git a/hoomd_organics/tests/utils/test_ff_utils.py b/hoomd_organics/tests/utils/test_ff_utils.py index 035d0276..49064b18 100644 --- a/hoomd_organics/tests/utils/test_ff_utils.py +++ b/hoomd_organics/tests/utils/test_ff_utils.py @@ -4,24 +4,20 @@ from gmso.core.forcefield import ForceField from hoomd_organics.tests import BaseTest -from hoomd_organics.utils.base_types import FF_Types from hoomd_organics.utils.ff_utils import find_xml_ff, xml_to_gmso_ff class TestFFUtils(BaseTest): def test_find_xml_ff(self): - ff_xml_path, ff_type = find_xml_ff("oplsaa.xml") - assert ff_type == FF_Types.XML + ff_xml_path = find_xml_ff("oplsaa.xml") assert os.path.exists(ff_xml_path) def test_find_xml_only_file_name(self): - ff_xml_path, ff_type = find_xml_ff("oplsaa") - assert ff_type == FF_Types.XML + ff_xml_path = find_xml_ff("oplsaa") assert os.path.exists(ff_xml_path) def test_find_xml_ff_path(self, benzene_xml): - ff_xml_path, ff_type = find_xml_ff(benzene_xml) - assert ff_type == FF_Types.XML + ff_xml_path = find_xml_ff(benzene_xml) assert os.path.exists(ff_xml_path) def test_find_xml_invalid_extension(self): diff --git a/hoomd_organics/utils/__init__.py b/hoomd_organics/utils/__init__.py index c2f8fe48..4200ed64 100644 --- a/hoomd_organics/utils/__init__.py +++ b/hoomd_organics/utils/__init__.py @@ -1,5 +1,5 @@ from .actions import * -from .base_types import FF_Types, HOOMDThermostats +from .base_types import HOOMDThermostats from .ff_utils import xml_to_gmso_ff from .utils import ( calculate_box_length, diff --git a/hoomd_organics/utils/base_types.py b/hoomd_organics/utils/base_types.py index 291f91a1..0b29a68b 100644 --- a/hoomd_organics/utils/base_types.py +++ b/hoomd_organics/utils/base_types.py @@ -1,11 +1,6 @@ import hoomd -class FF_Types: - XML = "XML" - HOOMD = "HOOMD" - - class HOOMDThermostats: BERENDSEN = hoomd.md.methods.thermostats.Berendsen BUSSI = hoomd.md.methods.thermostats.Bussi diff --git a/hoomd_organics/utils/ff_utils.py b/hoomd_organics/utils/ff_utils.py index 1845939e..bb4c4cfe 100644 --- a/hoomd_organics/utils/ff_utils.py +++ b/hoomd_organics/utils/ff_utils.py @@ -6,7 +6,6 @@ from hoomd_organics.assets import FF_DIR -from .base_types import FF_Types from .exceptions import ( MissingAnglePotentialError, MissingBondPotentialError, @@ -33,7 +32,6 @@ def find_xml_ff(ff_source): if not ff_source.endswith(".xml"): raise ValueError("ForceField file type must be XML.") ff_xml_path = ff_source - ff_type = FF_Types.XML elif not xml_directory.get(ff_source.split(".xml")[0]): raise ValueError( "{} forcefield is not supported. Supported XML forcefields " @@ -42,12 +40,11 @@ def find_xml_ff(ff_source): else: ff_key = ff_source.split(".xml")[0] ff_xml_path = xml_directory.get(ff_key) - ff_type = FF_Types.XML - return ff_xml_path, ff_type + return ff_xml_path def xml_to_gmso_ff(ff_xml): - ff_xml_path, ff_type = find_xml_ff(ff_xml) + ff_xml_path = find_xml_ff(ff_xml) gmso_ff = ffutils.FoyerFFs().load(ff_xml_path).to_gmso_ff() return gmso_ff From a0da90545d7dc61219e17e6f4626c42018470c20 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 29 Sep 2023 11:46:39 -0600 Subject: [PATCH 41/41] fix docstring --- hoomd_organics/base/molecule.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/hoomd_organics/base/molecule.py b/hoomd_organics/base/molecule.py index bad4f20e..2a7322dc 100644 --- a/hoomd_organics/base/molecule.py +++ b/hoomd_organics/base/molecule.py @@ -29,12 +29,12 @@ class Molecule: ---------- num_mols : int, required Number of molecules to generate. - force_field : hoomd_organics.ForceField or a list of ForceField objects, - default=None - The force field to be applied to the system for parameterization. - Note that setting `force_field` will not apply the forcefield to the - molecule. The forcefield in this step is just used for validation - purposes. + force_field : hoomd_organics.ForceField or a list of + `hoomd.md.force.Force` objects, default=None + The force field to be applied to the molecule for parameterization. + Note that setting `force_field` does not actually apply the + forcefield to the molecule. The forcefield in this step is mainly + used for validation purposes. smiles : str, default None The smiles string of the molecule to generate. file : str, default None