diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 8efa4ec8..49c21610 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -39,12 +39,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(), ): self._molecules = check_return_iterable(molecules) @@ -52,21 +48,18 @@ 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.gmso_system = None + self._reference_values = base_units self._hoomd_snapshot = None self._hoomd_forcefield = [] - self._reference_values = base_units 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() - self.gmso_system = None + self._ff_kwargs = dict() # Collecting all molecules self.n_mol_types = 0 @@ -120,16 +113,10 @@ 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() - 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 +147,7 @@ def net_charge(self): @property def box(self): + # TODO: Use gmso system here? return self.system.box @property @@ -205,12 +193,16 @@ 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 def hoomd_forcefield(self): if self._ff_refs != self.reference_values and self._force_field: - self._hoomd_forcefield = self._create_hoomd_forcefield() + self._hoomd_forcefield = self._create_hoomd_forcefield( + **self._ff_kwargs + ) return self._hoomd_forcefield @property @@ -220,7 +212,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. @@ -252,6 +244,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): """""" @@ -278,11 +273,13 @@ def _convert_to_gmso(self): topology.identify_connections() return topology - def _create_hoomd_forcefield(self): + 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, base_units=self._reference_values if self._reference_values @@ -304,7 +301,23 @@ 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, + ): + 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, @@ -312,11 +325,14 @@ 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 +349,20 @@ def _apply_forcefield(self): self._reference_values["length"] = length_scale * sigmas[0].unit_array self._reference_values["mass"] = mass_scale * masses[0].unit_array + if remove_hydrogens: + 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 + ) + self._hoomd_snapshot = self._create_hoomd_snapshot() + def set_target_box( self, x_constraint=None, y_constraint=None, z_constraint=None ): @@ -429,12 +459,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, @@ -445,11 +471,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, ) @@ -485,16 +507,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 @@ -505,11 +523,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_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_system.py b/hoomd_organics/tests/base/test_system.py index def50057..1409a36a 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( @@ -146,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): @@ -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 @@ -168,18 +165,18 @@ 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) 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,24 +193,25 @@ 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 ) 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 ) @@ -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, @@ -259,9 +257,9 @@ 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) assert np.allclose( system.reference_energy.to("kcal/mol").value, 0.066, atol=1e-3 ) @@ -272,9 +270,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 +292,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 +314,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 +333,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 +352,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 +368,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 +385,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 +397,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 +409,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 +421,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 +433,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 +445,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 +457,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 +469,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 +481,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 +493,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 +505,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 +517,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 +529,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 +541,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 +553,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 +565,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 +578,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 +590,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 +602,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 +614,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 +626,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,24 +638,51 @@ 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" + 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( 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 +697,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 +714,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..3652dc34 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -217,10 +217,10 @@ 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, auto_scale=True) return system @pytest.fixture() @@ -230,7 +230,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 +240,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() diff --git a/hoomd_organics/tests/library/test_tensile.py b/hoomd_organics/tests/library/test_tensile.py index 45d4aed6..90aadd6a 100644 --- a/hoomd_organics/tests/library/test_tensile.py +++ b/hoomd_organics/tests/library/test_tensile.py @@ -12,12 +12,13 @@ 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, auto_scale=True, ) + system.apply_forcefield(r_cut=2.5) + tensile_sim = Tensile( initial_state=system.hoomd_snapshot, forcefield=system.hoomd_forcefield,