From 33fa89cc35df7cdd3a3a8d1e435d1946f3e17362 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Fri, 22 Sep 2023 13:40:42 -0600 Subject: [PATCH 01/10] move parameters to apply_forcefield --- hoomd_organics/base/system.py | 57 +++++++++++++++++++++-------------- 1 file changed, 35 insertions(+), 22 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 8efa4ec8..4b5b57ee 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -39,12 +39,11 @@ def __init__( self, molecules, density: float, - r_cut: float, force_field=None, auto_scale=False, - remove_hydrogens=False, - remove_charges=False, - scale_charges=False, + # remove_hydrogens=False, + # remove_charges=False, + # scale_charges=False, base_units=dict(), ): self._molecules = check_return_iterable(molecules) @@ -52,17 +51,14 @@ def __init__( if force_field: self._force_field = check_return_iterable(force_field) self.density = density - self.r_cut = r_cut self.auto_scale = auto_scale - self.remove_hydrogens = remove_hydrogens - self.remove_charges = remove_charges - self.scale_charges = scale_charges - self._target_box = None self.all_molecules = [] + self._reference_values = base_units self._hoomd_snapshot = None self._hoomd_forcefield = [] - self._reference_values = base_units + self.r_cut = None self._gmso_forcefields_dict = dict() + self._target_box = None # Reference values used when last writing snapshot and forcefields self._ff_refs = dict() self._snap_refs = dict() @@ -122,14 +118,6 @@ def __init__( ) self.system = self._build_system() self.gmso_system = self._convert_to_gmso() - if self._force_field: - self._apply_forcefield() - if self.remove_hydrogens: - self._remove_hydrogens() - self._hoomd_forcefield = ( - self._create_hoomd_forcefield() if self._force_field else [] - ) - self._hoomd_snapshot = self._create_hoomd_snapshot() @abstractmethod def _build_system(self): @@ -160,6 +148,7 @@ def net_charge(self): @property def box(self): + # TODO: Use gmso system here? return self.system.box @property @@ -205,6 +194,8 @@ def reference_values(self, ref_value_dict): def hoomd_snapshot(self): if self._snap_refs != self.reference_values: self._hoomd_snapshot = self._create_hoomd_snapshot() + if self._hoomd_snapshot is None: # Hasn't been created yet + self._hoomd_snapshot = self._create_hoomd_snapshot() return self._hoomd_snapshot @property @@ -278,11 +269,13 @@ def _convert_to_gmso(self): topology.identify_connections() return topology - def _create_hoomd_forcefield(self): + def _create_hoomd_forcefield(self, nlist_buffer, pppm_kwargs): force_list = [] ff, refs = to_hoomd_forcefield( top=self.gmso_system, r_cut=self.r_cut, + nlist_buffer=nlist_buffer, + pppm_kwargs=pppm_kwargs, auto_scale=self.auto_scale, base_units=self._reference_values if self._reference_values @@ -304,7 +297,17 @@ def _create_hoomd_snapshot(self): self._snap_refs = self._reference_values.copy() return snap - def _apply_forcefield(self): + def apply_forcefield( + self, + r_cut, + auto_scale=False, + scale_charges=False, + remove_charges=False, + remove_hydrogens=False, + pppm_resolution=(8, 8, 8), + pppm_order=4, + nlist_buffer=0.4, + ): self.gmso_system = apply( self.gmso_system, self._gmso_forcefields_dict, @@ -312,10 +315,10 @@ def _apply_forcefield(self): speedup_by_moltag=True, speedup_by_molgraph=False, ) - if self.remove_charges: + if remove_charges: for site in self.gmso_system.sites: site.charge = 0 - if self.scale_charges and not self.remove_charges: + if scale_charges and not remove_charges: self._scale_charges() epsilons = [ s.atom_type.parameters["epsilon"] for s in self.gmso_system.sites @@ -333,6 +336,16 @@ def _apply_forcefield(self): self._reference_values["length"] = length_scale * sigmas[0].unit_array self._reference_values["mass"] = mass_scale * masses[0].unit_array + pppm_kwargs = {"resolution": pppm_resolution, "order": pppm_order} + self._hoomd_forcefield = ( + self._create_hoomd_forcefield( + nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs + ) + if self._force_field + else [] + ) + self._hoomd_snapshot = self._create_hoomd_snapshot() + def set_target_box( self, x_constraint=None, y_constraint=None, z_constraint=None ): From f407bc9aeb322e450325eea7b1c568c143bb60fe Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Fri, 22 Sep 2023 14:30:35 -0600 Subject: [PATCH 02/10] change remove_hydrogens to public --- hoomd_organics/base/system.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 4b5b57ee..2e305e31 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -41,9 +41,6 @@ def __init__( density: float, force_field=None, auto_scale=False, - # remove_hydrogens=False, - # remove_charges=False, - # scale_charges=False, base_units=dict(), ): self._molecules = check_return_iterable(molecules) @@ -116,7 +113,9 @@ def __init__( f"{self._force_field[ff_index]} is not " f"provided." ) + # Create mBuild system self.system = self._build_system() + # Create GMSO topology self.gmso_system = self._convert_to_gmso() @abstractmethod @@ -211,7 +210,7 @@ def target_box(self): else: return self._target_box - def _remove_hydrogens(self): + def remove_hydrogens(self): """Call this method to remove hydrogen atoms from the system. The masses and charges of the hydrogens are absorbed into the heavy atoms they were bonded to. @@ -243,6 +242,9 @@ 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 self._hoomd_snapshot: + self._create_hoomd_snapshot() def _scale_charges(self): """""" @@ -344,6 +346,8 @@ def apply_forcefield( if self._force_field else [] ) + if remove_hydrogens: + self.remove_hydrogens() self._hoomd_snapshot = self._create_hoomd_snapshot() def set_target_box( From 4caec72a340cfb23dde2956af83146c80b4c3c76 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Fri, 22 Sep 2023 14:41:31 -0600 Subject: [PATCH 03/10] add error message when calling apply_forcefield --- hoomd_organics/base/system.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 2e305e31..6b36f6d6 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -310,6 +310,12 @@ def apply_forcefield( pppm_order=4, nlist_buffer=0.4, ): + if not self._force_field: + # TODO: Better erorr message + raise ValueError( + "This method can only be used when the System is " + "initialized with an XML type forcefield." + ) self.gmso_system = apply( self.gmso_system, self._gmso_forcefields_dict, @@ -317,11 +323,14 @@ def apply_forcefield( speedup_by_moltag=True, speedup_by_molgraph=False, ) + if remove_charges: for site in self.gmso_system.sites: site.charge = 0 + if scale_charges and not remove_charges: self._scale_charges() + epsilons = [ s.atom_type.parameters["epsilon"] for s in self.gmso_system.sites ] @@ -338,16 +347,13 @@ def apply_forcefield( self._reference_values["length"] = length_scale * sigmas[0].unit_array self._reference_values["mass"] = mass_scale * masses[0].unit_array - pppm_kwargs = {"resolution": pppm_resolution, "order": pppm_order} - self._hoomd_forcefield = ( - self._create_hoomd_forcefield( - nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs - ) - if self._force_field - else [] - ) if remove_hydrogens: self.remove_hydrogens() + + pppm_kwargs = {"resolution": pppm_resolution, "order": pppm_order} + self._hoomd_forcefield = self._create_hoomd_forcefield( + nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs + ) self._hoomd_snapshot = self._create_hoomd_snapshot() def set_target_box( From b4aee940de96c4caebf011ab9393fb284f31abdd Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Fri, 22 Sep 2023 15:05:01 -0600 Subject: [PATCH 04/10] fix r_cut param --- hoomd_organics/base/system.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 6b36f6d6..786d0562 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -53,7 +53,6 @@ def __init__( self._reference_values = base_units self._hoomd_snapshot = None self._hoomd_forcefield = [] - self.r_cut = None self._gmso_forcefields_dict = dict() self._target_box = None # Reference values used when last writing snapshot and forcefields @@ -271,11 +270,11 @@ def _convert_to_gmso(self): topology.identify_connections() return topology - def _create_hoomd_forcefield(self, nlist_buffer, pppm_kwargs): + def _create_hoomd_forcefield(self, r_cut, nlist_buffer, pppm_kwargs): force_list = [] ff, refs = to_hoomd_forcefield( top=self.gmso_system, - r_cut=self.r_cut, + r_cut=r_cut, nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs, auto_scale=self.auto_scale, @@ -352,7 +351,7 @@ def apply_forcefield( pppm_kwargs = {"resolution": pppm_resolution, "order": pppm_order} self._hoomd_forcefield = self._create_hoomd_forcefield( - nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs + r_cut=r_cut, nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs ) self._hoomd_snapshot = self._create_hoomd_snapshot() @@ -452,12 +451,8 @@ def __init__( self, molecules, density: float, - r_cut: float, force_field=None, auto_scale=False, - remove_hydrogens=False, - remove_charges=False, - scale_charges=False, base_units=dict(), packing_expand_factor=5, edge=0.2, @@ -468,11 +463,7 @@ def __init__( molecules=molecules, density=density, force_field=force_field, - r_cut=r_cut, auto_scale=auto_scale, - remove_hydrogens=remove_hydrogens, - remove_charges=remove_charges, - scale_charges=scale_charges, base_units=base_units, ) From 666ec7dc375e54c1781fcdad7ecb21a0c87ab862 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Fri, 22 Sep 2023 15:47:11 -0600 Subject: [PATCH 05/10] fix tests --- hoomd_organics/base/system.py | 8 -- hoomd_organics/tests/base/test_system.py | 99 ++++++++++++------------ hoomd_organics/tests/base_test.py | 7 +- 3 files changed, 50 insertions(+), 64 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 786d0562..89199526 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -499,16 +499,12 @@ def __init__( self, molecules, density: float, - r_cut: float, x: float, y: float, n: int, basis_vector=[0.5, 0.5, 0], force_field=None, auto_scale=False, - remove_hydrogens=False, - remove_charges=False, - scale_charges=False, base_units=dict(), ): self.x = x @@ -519,11 +515,7 @@ def __init__( molecules=molecules, density=density, force_field=force_field, - r_cut=r_cut, auto_scale=auto_scale, - remove_hydrogens=remove_hydrogens, - remove_charges=remove_charges, - scale_charges=scale_charges, base_units=base_units, ) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index def50057..2e056cd8 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -17,10 +17,10 @@ def test_single_mol_type(self, benzene_molecule): system = Pack( molecules=[benzene_mols], density=0.8, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(benzene_mols.molecules) assert system.gmso_system.is_typed() @@ -35,10 +35,10 @@ def test_multiple_mol_types(self, benzene_molecule, ethane_molecule): system = Pack( molecules=[benzene_mol, ethane_mol], density=0.8, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert system.n_mol_types == 2 assert len(system.all_molecules) == len(benzene_mol.molecules) + len( ethane_mol.molecules @@ -63,10 +63,10 @@ def test_multiple_mol_types_different_ff( system = Pack( molecules=[dimethylether_mol, pps_mol], density=0.8, - r_cut=2.5, force_field=[OPLS_AA_DIMETHYLETHER(), OPLS_AA_PPS()], auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert system.n_mol_types == 2 assert len(system.all_molecules) == len( dimethylether_mol.molecules @@ -95,10 +95,10 @@ def test_system_from_mol2_mol_parameterization(self, benzene_molecule_mol2): system = Pack( molecules=[benzene_mol], density=0.8, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert system.gmso_system.is_typed() assert len(system.hoomd_forcefield) > 0 assert system.n_particles == system.hoomd_snapshot.particles.N @@ -108,11 +108,10 @@ def test_remove_hydrogen(self, benzene_molecule): system = Pack( molecules=[benzene_mol], density=0.8, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, - remove_hydrogens=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,11 +131,10 @@ def test_remove_hydrogen_no_atomic_num(self, benzene_molecule): system = Pack( molecules=[benzene_mol], density=0.8, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, - remove_hydrogens=False, ) + system.apply_forcefield(r_cut=2.5, remove_hydrogens=False) for site in system.gmso_system.sites: if site.name == "H": site.element = gmso.core.element.Element( @@ -154,11 +152,10 @@ def test_remove_hydrogen_no_hydrogen(self, benzene_molecule): system = Pack( molecules=[benzene_mol], density=0.8, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, - remove_hydrogens=False, ) + system.apply_forcefield(r_cut=2.5, remove_hydrogens=False) hydrogens = [ site for site in system.gmso_system.sites @@ -175,11 +172,11 @@ def test_add_mass_charges(self, benzene_molecule): system = Pack( molecules=[benzene_mols], density=0.8, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=False, - remove_hydrogens=True, - scale_charges=False, + ) + system.apply_forcefield( + r_cut=2.5, remove_hydrogens=True, scale_charges=False ) for site in system.gmso_system.sites: assert site.mass.value == (12.011 + 1.008) @@ -196,17 +193,18 @@ def test_target_box(self, benzene_molecule): low_density_system = Pack( molecules=[benzene_mol], density=0.1, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, ) + low_density_system.apply_forcefield(r_cut=2.5) + high_density_system = Pack( molecules=[benzene_mol], density=0.9, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, ) + high_density_system.apply_forcefield(r_cut=2.5) assert all( low_density_system.target_box > high_density_system.target_box ) @@ -224,9 +222,9 @@ def test_ref_length(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert np.allclose( system.reference_length.to("angstrom").value, 3.5, atol=1e-3 @@ -243,9 +241,9 @@ def test_ref_mass(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) total_red_mass = sum(system.hoomd_snapshot.particles.mass) assert np.allclose( system.mass, @@ -262,6 +260,7 @@ def test_ref_energy(self, polyethylene): r_cut=2.5, auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert np.allclose( system.reference_energy.to("kcal/mol").value, 0.066, atol=1e-3 ) @@ -272,9 +271,9 @@ def test_rebuild_snapshot(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) assert system._snap_refs == system.reference_values assert system._ff_refs == system.reference_values init_snap = system.hoomd_snapshot @@ -294,9 +293,9 @@ def test_ref_values_no_autoscale(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_length = 1 * u.angstrom system.reference_energy = 1 * u.kcal / u.mol system.reference_mass = 1 * u.amu @@ -316,9 +315,9 @@ def test_set_ref_values(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -335,9 +334,9 @@ def test_set_ref_values_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) ref_value_dict = { "length": "1 angstrom", "energy": "3 kcal/mol", @@ -354,9 +353,9 @@ def test_set_ref_values_missing_key(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -370,9 +369,9 @@ def test_set_ref_values_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) ref_value_dict = { "length": 1 * u.angstrom, "energy": 3.0 * u.kcal / u.mol, @@ -387,9 +386,9 @@ def test_set_ref_length(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_length = 1 * u.angstrom assert system.reference_length == 1 * u.angstrom @@ -399,9 +398,9 @@ def test_set_ref_length_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_length = 1.0 @@ -411,9 +410,9 @@ def test_ref_length_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_length = "1 angstrom" assert system.reference_length == 1 * u.angstrom @@ -423,9 +422,9 @@ def test_ref_length_invalid_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0" @@ -435,9 +434,9 @@ def test_ref_length_invalid_unit_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0 invalid_unit" @@ -447,9 +446,9 @@ def test_ref_length_invalid_dimension(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_length = 1.0 * u.g @@ -459,9 +458,9 @@ def test_ref_length_invalid_dimension_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_length = "1.0 g" @@ -471,9 +470,9 @@ def test_set_ref_energy(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_energy = 1 * u.kcal / u.mol assert system.reference_energy == 1 * u.kcal / u.mol @@ -483,9 +482,9 @@ def test_set_ref_energy_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 @@ -495,9 +494,9 @@ def test_ref_energy_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_energy = "1 kJ" assert system.reference_energy == 1 * u.kJ @@ -507,9 +506,9 @@ def test_ref_energy_string_comb_units(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_energy = "1 kcal/mol" assert system.reference_energy == 1 * u.kcal / u.mol @@ -519,9 +518,9 @@ def test_ref_energy_invalid_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0" @@ -531,9 +530,9 @@ def test_ref_energy_invalid_unit_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0 invalid_unit" @@ -543,9 +542,9 @@ def test_ref_energy_invalid_dimension(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 * u.g @@ -555,9 +554,9 @@ def test_ref_energy_invalid_dimension_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_energy = "1.0 m" @@ -567,9 +566,9 @@ def test_set_ref_mass(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_mass = 1.0 * u.amu assert system.reference_mass == 1.0 * u.amu @@ -580,9 +579,9 @@ def test_set_ref_mass_invalid_type(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_mass = 1.0 @@ -592,9 +591,9 @@ def test_ref_mass_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) system.reference_mass = "1 g" assert system.reference_mass == 1.0 * u.g @@ -604,9 +603,9 @@ def test_ref_mass_invalid_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0" @@ -616,9 +615,9 @@ def test_ref_mass_invalid_unit_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 invalid_unit" @@ -628,9 +627,9 @@ def test_ref_mass_invalid_dimension(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_energy = 1.0 * u.m @@ -640,9 +639,9 @@ def test_ref_mass_invalid_dimension_string(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=False, ) + system.apply_forcefield(r_cut=2.5) with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 m" @@ -652,12 +651,12 @@ def test_lattice_polymer(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, x=1, y=1, n=4, auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(polyethylene.molecules) @@ -672,12 +671,12 @@ def test_lattice_molecule(self, benzene_molecule): molecules=[benzene_mol], force_field=OPLS_AA(), density=1.0, - r_cut=2.5, x=1, y=1, n=4, auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) assert system.n_mol_types == 1 assert len(system.all_molecules) == len(benzene_mol.molecules) assert len(system.hoomd_forcefield) > 0 @@ -689,19 +688,17 @@ def test_scale_charges(self, pps): no_scale = Pack( molecules=pps_mol, density=0.5, - r_cut=2.4, force_field=OPLS_AA_PPS(), auto_scale=True, - scale_charges=False, ) + no_scale.apply_forcefield(r_cut=2.5, scale_charges=False) with_scale = Pack( molecules=pps_mol, density=0.5, - r_cut=2.4, force_field=OPLS_AA_PPS(), auto_scale=True, - scale_charges=True, ) + with_scale.apply_forcefield(r_cut=2.5, scale_charges=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 c6666dbe..1f4fe89c 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -217,10 +217,9 @@ def benzene_system(self, benzene_mb): system = Pack( molecules=[benzene], density=0.2, - r_cut=2.5, force_field=OPLS_AA(), - auto_scale=True, ) + system.apply_forcefield(r_cut=2.5, remove_hydrogens=True) return system @pytest.fixture() @@ -230,7 +229,6 @@ def benzene_cg_system(self, benzene_molecule, benzene_smiles): system = Pack( molecules=[benzene_mols], density=0.5, - r_cut=2.5, auto_scale=False, ) return system @@ -241,11 +239,10 @@ def polyethylene_system(self, polyethylene): system = Pack( molecules=polyethylene_mol, density=0.5, - r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, - remove_hydrogens=True, ) + system.apply_forcefield(r_cut=2.5, remove_hydrogens=True) return system @pytest.fixture() From 16d1676649ce577d9c98750b963802d554c617bf Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Fri, 22 Sep 2023 22:51:02 -0600 Subject: [PATCH 06/10] add check for ref values back into self.hoomd_forcefield --- hoomd_organics/base/system.py | 14 +++++++++++--- hoomd_organics/tests/base/test_system.py | 7 +++---- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 89199526..f07586e5 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -50,6 +50,7 @@ def __init__( self.density = density self.auto_scale = auto_scale self.all_molecules = [] + self.gmso_system = None self._reference_values = base_units self._hoomd_snapshot = None self._hoomd_forcefield = [] @@ -58,7 +59,7 @@ def __init__( # Reference values used when last writing snapshot and forcefields self._ff_refs = dict() self._snap_refs = dict() - self.gmso_system = None + self._ff_kwargs = dict() # Collecting all molecules self.n_mol_types = 0 @@ -198,8 +199,10 @@ def hoomd_snapshot(self): @property def hoomd_forcefield(self): - if self._ff_refs != self.reference_values and self._force_field: - self._hoomd_forcefield = self._create_hoomd_forcefield() + if self._ff_refs != self.reference_values: + self._hoomd_forcefield = self._create_hoomd_forcefield( + **self._ff_kwargs + ) return self._hoomd_forcefield @property @@ -350,6 +353,11 @@ def apply_forcefield( self.remove_hydrogens() pppm_kwargs = {"resolution": pppm_resolution, "order": pppm_order} + self._ff_kwargs = { + "r_cut": r_cut, + "nlist_buffer": nlist_buffer, + "pppm_kwargs": pppm_kwargs, + } self._hoomd_forcefield = self._create_hoomd_forcefield( r_cut=r_cut, nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs ) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index 2e056cd8..51d9b984 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -144,7 +144,7 @@ def test_remove_hydrogen_no_atomic_num(self, benzene_molecule): mass=1.008 * Unit("amu"), ) - system._remove_hydrogens() + system.remove_hydrogens() assert system.gmso_system.n_sites == 6 def test_remove_hydrogen_no_hydrogen(self, benzene_molecule): @@ -165,7 +165,7 @@ def test_remove_hydrogen_no_hydrogen(self, benzene_molecule): system.gmso_system.remove_site(h_site) with pytest.warns(): - system._remove_hydrogens() + system.remove_hydrogens() def test_add_mass_charges(self, benzene_molecule): benzene_mols = benzene_molecule(n_mols=1) @@ -211,7 +211,7 @@ def test_target_box(self, benzene_molecule): def test_mass(self, pps_molecule): pps_mol = pps_molecule(n_mols=20) - system = Pack(molecules=[pps_mol], density=1.0, r_cut=2.5) + system = Pack(molecules=[pps_mol], density=1.0) assert np.allclose( system.mass, ((12.011 * 6) + (1.008 * 6) + 32.06) * 20, atol=1e-4 ) @@ -257,7 +257,6 @@ def test_ref_energy(self, polyethylene): molecules=[polyethylene], force_field=[OPLS_AA()], density=1.0, - r_cut=2.5, auto_scale=True, ) system.apply_forcefield(r_cut=2.5) From 63a9245b562a3f2191eb9d6bb5286789c3331468 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Fri, 22 Sep 2023 23:17:14 -0600 Subject: [PATCH 07/10] fix unit test --- hoomd_organics/base/system.py | 2 +- hoomd_organics/tests/base/test_simulation.py | 4 ++-- hoomd_organics/tests/base_test.py | 3 ++- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index f07586e5..49c21610 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -199,7 +199,7 @@ def hoomd_snapshot(self): @property def hoomd_forcefield(self): - if self._ff_refs != self.reference_values: + if self._ff_refs != self.reference_values and self._force_field: self._hoomd_forcefield = self._create_hoomd_forcefield( **self._ff_kwargs ) diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py index 76fa7112..948121dd 100644 --- a/hoomd_organics/tests/base/test_simulation.py +++ b/hoomd_organics/tests/base/test_simulation.py @@ -151,11 +151,11 @@ def test_update_volume_walls(self, benzene_system): def test_update_volume_density(self, benzene_system): sim = Simulation.from_system(benzene_system) sim.run_update_volume( - kT=1.0, tau_kt=0.01, n_steps=500, period=1, final_density=0.1 + kT=1.0, tau_kt=0.01, n_steps=500, period=1, final_density=0.05 ) assert np.isclose( sim.density.to(u.g / u.cm**3).value, - (0.1 * (u.g / u.cm**3)).value, + (0.05 * (u.g / u.cm**3)).value, atol=1e-4, ) diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py index 1f4fe89c..3652dc34 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -218,8 +218,9 @@ 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, remove_hydrogens=True) + system.apply_forcefield(r_cut=2.5, auto_scale=True) return system @pytest.fixture() From b93a70fc8c4878464e58d5542a6a467a6380e01b Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 25 Sep 2023 10:44:49 -0600 Subject: [PATCH 08/10] fix kwarg --- hoomd_organics/tests/library/test_tensile.py | 1 - 1 file changed, 1 deletion(-) diff --git a/hoomd_organics/tests/library/test_tensile.py b/hoomd_organics/tests/library/test_tensile.py index 45d4aed6..0c758777 100644 --- a/hoomd_organics/tests/library/test_tensile.py +++ b/hoomd_organics/tests/library/test_tensile.py @@ -12,7 +12,6 @@ def test_tensile(self): molecules=[pps], force_field=[OPLS_AA_PPS()], density=1.0, - r_cut=2.5, x=1.2, y=1.2, n=4, From b512cca3647ef522ffb1be61eaa66754d64ba6c8 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 25 Sep 2023 11:22:50 -0600 Subject: [PATCH 09/10] fix unit test --- hoomd_organics/tests/library/test_tensile.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/hoomd_organics/tests/library/test_tensile.py b/hoomd_organics/tests/library/test_tensile.py index 0c758777..90aadd6a 100644 --- a/hoomd_organics/tests/library/test_tensile.py +++ b/hoomd_organics/tests/library/test_tensile.py @@ -17,6 +17,8 @@ def test_tensile(self): n=4, auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) + tensile_sim = Tensile( initial_state=system.hoomd_snapshot, forcefield=system.hoomd_forcefield, From 637efa1b1b70becde4eee7dcd4de8da7b7176682 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 25 Sep 2023 12:49:54 -0600 Subject: [PATCH 10/10] add a couple more unit tests --- hoomd_organics/tests/base/test_system.py | 27 ++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/hoomd_organics/tests/base/test_system.py b/hoomd_organics/tests/base/test_system.py index 51d9b984..1409a36a 100644 --- a/hoomd_organics/tests/base/test_system.py +++ b/hoomd_organics/tests/base/test_system.py @@ -644,6 +644,33 @@ def test_ref_mass_invalid_dimension_string(self, polyethylene): with pytest.raises(ReferenceUnitError): system.reference_mass = "1.0 m" + 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, + auto_scale=False, + ) + with pytest.raises(ValueError): + system.apply_forcefield(r_cut=2.5) + + 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, + auto_scale=False, + ) + system.apply_forcefield( + r_cut=2.5, 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 + assert system._ff_kwargs["pppm_kwargs"]["resolution"] == (4, 4, 4) + assert system._ff_kwargs["pppm_kwargs"]["order"] == 3 + def test_lattice_polymer(self, polyethylene): polyethylene = polyethylene(lengths=2, num_mols=32) system = Lattice(