diff --git a/hoomd_organics/base/simulation.py b/hoomd_organics/base/simulation.py index 3e5c7ad2..910b1dcb 100644 --- a/hoomd_organics/base/simulation.py +++ b/hoomd_organics/base/simulation.py @@ -7,7 +7,12 @@ import numpy as np import unyt as u -from hoomd_organics.utils import StdOutLogger, UpdateWalls +from hoomd_organics.utils import ( + StdOutLogger, + UpdateWalls, + calculate_box_length, + validate_ref_value, +) from hoomd_organics.utils.exceptions import ReferenceUnitError @@ -52,6 +57,7 @@ def __init__( self, initial_state, forcefield=None, + reference_values=dict(), r_cut=2.5, dt=0.0001, device=hoomd.device.auto_select(), @@ -83,12 +89,43 @@ def __init__( self.integrator = None self._dt = dt self._reference_values = dict() + self._reference_values = reference_values self._integrate_group = hoomd.filter.All() self._wall_forces = dict() self._create_state(self.initial_state) # Add a gsd and thermo props logger to sim operations self._add_hoomd_writers() + @classmethod + def from_system(cls, system, **kwargs): + """Initialize a simulation from a `hoomd_organics.base.System` + object.""" + + if system.hoomd_forcefield: + return cls( + initial_state=system.hoomd_snapshot, + forcefield=system.hoomd_forcefield, + reference_values=system.reference_values, + **kwargs, + ) + elif kwargs.get("forcefield", None): + return cls( + initial_state=system.hoomd_snapshot, + reference_values=system.reference_values, + **kwargs, + ) + else: + raise ValueError( + "No forcefield provided. Please provide a forcefield " + "or a system with a forcefield." + ) + + @classmethod + def from_snapshot_forces(cls, initial_state, forcefield, **kwargs): + """Initialize a simulation from an initial state object and a + list of HOOMD forces.""" + return cls(initial_state=initial_state, forcefield=forcefield, **kwargs) + @property def forces(self): if self.integrator: @@ -113,50 +150,19 @@ def reference_values(self): return self._reference_values @reference_length.setter - def reference_length(self, length, unit=None): - if isinstance(length, u.array.unyt_quantity): - self._reference_values["length"] = length - elif isinstance(unit, str) and ( - isinstance(length, float) or isinstance(length, int) - ): - self._reference_values["length"] = length * getattr(u, unit) - else: - raise ReferenceUnitError( - f"Invalid reference length input.Please provide reference " - f"length (number) and unit (string) or pass length value as an " - f"{str(u.array.unyt_quantity)}." - ) + def reference_length(self, length): + validated_length = validate_ref_value(length, u.dimensions.length) + self._reference_values["length"] = validated_length @reference_energy.setter - def reference_energy(self, energy, unit=None): - if isinstance(energy, u.array.unyt_quantity): - self._reference_values["energy"] = energy - elif isinstance(unit, str) and ( - isinstance(energy, float) or isinstance(energy, int) - ): - self._reference_values["energy"] = energy * getattr(u, unit) - else: - raise ReferenceUnitError( - f"Invalid reference energy input.Please provide reference " - f"energy (number) and unit (string) or pass energy value as an " - f"{str(u.array.unyt_quantity)}." - ) + def reference_energy(self, energy): + validated_energy = validate_ref_value(energy, u.dimensions.energy) + self._reference_values["energy"] = validated_energy @reference_mass.setter - def reference_mass(self, mass, unit=None): - if isinstance(mass, u.array.unyt_quantity): - self._reference_values["mass"] = mass - elif isinstance(unit, str) and ( - isinstance(mass, float) or isinstance(mass, int) - ): - self._reference_values["mass"] = mass * getattr(u, unit) - else: - raise ReferenceUnitError( - f"Invalid reference mass input.Please provide reference " - f"mass (number) and " - f"unit (string) or pass mass value as an " - f"{str(u.array.unyt_quantity)}." - ) + def reference_mass(self, mass): + validated_mass = validate_ref_value(mass, u.dimensions.mass) + self._reference_values["mass"] = validated_mass @reference_values.setter def reference_values(self, ref_value_dict): @@ -164,12 +170,7 @@ def reference_values(self, ref_value_dict): for k in ref_keys: if k not in ref_value_dict.keys(): raise ValueError(f"Missing reference for {k}.") - if not isinstance(ref_value_dict[k], u.array.unyt_quantity): - raise ReferenceUnitError( - f"{k} reference value must be of type " - f"{str(u.array.unyt_quantity)}" - ) - self._reference_values = ref_value_dict + self.__setattr__(f"reference_{k}", ref_value_dict[k]) @property def box_lengths_reduced(self): @@ -379,7 +380,8 @@ def run_update_volume( period, kT, tau_kt, - final_box_lengths, + final_box_lengths=None, + final_density=None, thermalize_particles=True, ): """Runs an NVT simulation while shrinking or expanding @@ -395,20 +397,54 @@ def run_update_volume( The temperature to use during shrinking. tau_kt : float; required Thermostat coupling period (in simulation time units) - final_box_lengths : np.ndarray, shape=(3,), dtype=float; required + final_box_lengths : np.ndarray, shape=(3,), dtype=float; optional The final box edge lengths in (x, y, z) order + final_density : float; optional + The final density of the simulation """ + if final_box_lengths is None and final_density is None: + raise ValueError( + "Must provide either `final_box_lengths` or `final_density`" + ) + if final_box_lengths is not None and final_density is not None: + raise ValueError( + "Cannot provide both `final_box_lengths` and `final_density`." + ) + if final_box_lengths is not None: + final_box = hoomd.Box( + Lx=final_box_lengths[0], + Ly=final_box_lengths[1], + Lz=final_box_lengths[2], + ) + else: + if not self.reference_values: + raise ReferenceUnitError( + "Missing simulation units. Please " + "provide units for mass, length, and" + " energy." + ) + + if isinstance(final_density, u.unyt_quantity): + density_quantity = final_density.to(u.g / u.cm**3) + else: + density_quantity = u.unyt_quantity( + final_density, u.g / u.cm**3 + ) + mass_g = self.mass.to("g") + L = calculate_box_length(mass_g, density_quantity) + # convert L from cm to reference units + L = ( + L.to(self.reference_length.units) / self.reference_length.value + ).value + final_box = hoomd.Box(Lx=L, Ly=L, Lz=L) + resize_trigger = hoomd.trigger.Periodic(period) box_ramp = hoomd.variant.Ramp( A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps) ) initial_box = self.state.box - final_box = hoomd.Box( - Lx=final_box_lengths[0], - Ly=final_box_lengths[1], - Lz=final_box_lengths[2], - ) + box_resizer = hoomd.update.BoxResize( box1=initial_box, box2=final_box, diff --git a/hoomd_organics/base/system.py b/hoomd_organics/base/system.py index 7a43a6af..8c13173f 100644 --- a/hoomd_organics/base/system.py +++ b/hoomd_organics/base/system.py @@ -12,6 +12,7 @@ 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, @@ -56,7 +57,7 @@ def __init__( self.remove_hydrogens = remove_hydrogens self.remove_charges = remove_charges self.scale_charges = scale_charges - self.target_box = None + self._target_box = None self.all_molecules = [] self._hoomd_snapshot = None self._hoomd_forcefield = [] @@ -68,7 +69,8 @@ def __init__( self.n_mol_types = 0 for mol_item in self._molecules: if isinstance(mol_item, Molecule): - mol_item.assign_mol_name(str(self.n_mol_types)) + if self._force_field: + mol_item.assign_mol_name(str(self.n_mol_types)) self.all_molecules.extend(mol_item.molecules) # if ff is provided in Molecule class if mol_item.force_field: @@ -199,8 +201,18 @@ def hoomd_snapshot(self): @property def hoomd_forcefield(self): - self._hoomd_forcefield = self._create_hoomd_forcefield() - return self._hoomd_forcefield + if self._force_field: + self._hoomd_forcefield = self._create_hoomd_forcefield() + return self._hoomd_forcefield + else: + return self._hoomd_forcefield + + @property + def target_box(self): + if self.reference_length: + return self._target_box / self.reference_length.value + else: + return self._target_box def _remove_hydrogens(self): """Call this method to remove hydrogen atoms from the system. @@ -311,12 +323,7 @@ def _apply_forcefield(self): self._reference_values["energy"] = energy_scale * epsilons[0].unit_array self._reference_values["length"] = length_scale * sigmas[0].unit_array - if self.auto_scale: - self._reference_values["mass"] = mass_scale * masses[ - 0 - ].unit_array.to("amu") - else: - self._reference_values["mass"] = mass_scale * u.g / u.mol + self._reference_values["mass"] = mass_scale * masses[0].unit_array def set_target_box( self, x_constraint=None, y_constraint=None, z_constraint=None @@ -341,14 +348,14 @@ def set_target_box( Lx = Ly = Lz = self._calculate_L() else: constraints = np.array([x_constraint, y_constraint, z_constraint]) - fixed_L = constraints[np.where(constraints is not None)] + fixed_L = constraints[np.not_equal(constraints, None).nonzero()] # Conv from nm to cm for _calculate_L fixed_L *= 1e-7 L = self._calculate_L(fixed_L=fixed_L) - constraints[np.where(constraints is None)] = L + constraints[np.equal(constraints, None).nonzero()] = L Lx, Ly, Lz = constraints - self.target_box = np.array([Lx, Ly, Lz]) + self._target_box = np.array([Lx, Ly, Lz]) def visualize(self): if self.system: @@ -360,7 +367,7 @@ def visualize(self): def _calculate_L(self, fixed_L=None): """Calculates the required box length(s) given the - mass of a sytem and the target density. + mass of a system and the target density. Box edge length constraints can be set by set_target_box(). If constraints are set, this will solve for the required @@ -374,18 +381,14 @@ def _calculate_L(self, fixed_L=None): when solving for L """ - # Convert from amu to grams - M = self.mass * 1.66054e-24 - vol = M / self.density # cm^3 - if fixed_L is None: - L = vol ** (1 / 3) - else: - L = vol / np.prod(fixed_L) - if len(fixed_L) == 1: # L is cm^2 - L = L ** (1 / 2) - # Convert from cm back to nm - L *= 1e7 - return L + mass_quantity = u.unyt_quantity(self.mass, u.g / u.mol).to("g") + density_quantity = u.unyt_quantity(self.density, u.g / u.cm**3) + if fixed_L is not None: + fixed_L = u.unyt_array(fixed_L, u.cm) + L = calculate_box_length( + mass_quantity, density_quantity, fixed_L=fixed_L + ) + return L.to("nm").value class Pack(System): @@ -393,6 +396,21 @@ class Pack(System): The box used for packing is expanded to allow PACKMOL to more easily place all the molecules. + Warnings: + --------- + Note that the default `packing_expand_factor` for pack is 5, which means + that the box density will not be the same as the specified density. This is + because in some cases PACKMOL will not be able to fit all the molecules + into the box if the target box is too small, therefore, we need to expand + the box by a factor (default:5) to allow PACKMOL to fit all the molecules. + + In order to get the + specified density there are two options: + 1) set the `packing_expand_factor` to 1, which will not expand the box, + however, this may result in PACKMOL errors if the box is too small. + 2) Update the box volume after creating the simulation object to the target + box length. This property is called `target_box`. + Parameters ---------- packing_expand_factor : int; optional, default 5 @@ -432,7 +450,7 @@ def _build_system(self): system = mb.packing.fill_box( compound=self.all_molecules, n_compounds=[1 for i in self.all_molecules], - box=list(self.target_box * self.packing_expand_factor), + box=list(self._target_box * self.packing_expand_factor), overlap=0.2, edge=self.edge, ) diff --git a/hoomd_organics/tests/base/sim_data.txt b/hoomd_organics/tests/base/sim_data.txt new file mode 100644 index 00000000..e69de29b diff --git a/hoomd_organics/tests/base/test_simulation.py b/hoomd_organics/tests/base/test_simulation.py index 8b37be39..65fd40c2 100644 --- a/hoomd_organics/tests/base/test_simulation.py +++ b/hoomd_organics/tests/base/test_simulation.py @@ -1,22 +1,44 @@ +import copy import os import pickle import hoomd import numpy as np +import pytest +import unyt as u from hoomd_organics import Simulation from hoomd_organics.tests import BaseTest class TestSimulate(BaseTest): - def test_initialize_from_snap(self, benzene_system): - Simulation( + def test_initialize_from_system(self, benzene_system): + sim = Simulation.from_system(benzene_system) + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) + assert len(sim.forces) == len(benzene_system.hoomd_forcefield) + assert sim.reference_values == benzene_system.reference_values + + def test_initialize_from_system_separate_ff( + self, benzene_cg_system, cg_single_bead_ff + ): + sim = Simulation.from_system( + benzene_cg_system, forcefield=cg_single_bead_ff + ) + sim.run_NVT(kT=0.1, tau_kt=10, n_steps=500) + + def test_initialize_from_system_missing_ff(self, benzene_cg_system): + with pytest.raises(ValueError): + Simulation.from_system(benzene_cg_system) + + def test_initialize_from_state(self, benzene_system): + Simulation.from_snapshot_forces( initial_state=benzene_system.hoomd_snapshot, forcefield=benzene_system.hoomd_forcefield, + reference_values=benzene_system.reference_values, ) def test_no_reference_values(self, benzene_system): - sim = Simulation( + sim = Simulation.from_snapshot_forces( initial_state=benzene_system.hoomd_snapshot, forcefield=benzene_system.hoomd_forcefield, ) @@ -29,24 +51,57 @@ def test_reference_values(self, benzene_system): sim = Simulation( initial_state=benzene_system.hoomd_snapshot, forcefield=benzene_system.hoomd_forcefield, + reference_values=benzene_system.reference_values, ) - sim.reference_values = benzene_system.reference_values assert np.isclose(float(sim.mass.value), benzene_system.mass, atol=1e-4) assert np.allclose(benzene_system.box.lengths, sim.box_lengths.value) - def test_NVT(self, benzene_system): + def test_set_ref_values(self, benzene_system): sim = Simulation( initial_state=benzene_system.hoomd_snapshot, forcefield=benzene_system.hoomd_forcefield, ) - sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) - assert isinstance(sim.method, hoomd.md.methods.NVT) + ref_value_dict = { + "length": 1 * u.angstrom, + "energy": 3.0 * u.kcal / u.mol, + "mass": 1.25 * u.Unit("amu"), + } + sim.reference_values = ref_value_dict + assert sim.reference_length == ref_value_dict["length"] + assert sim.reference_energy == ref_value_dict["energy"] + assert sim.reference_mass == ref_value_dict["mass"] - def test_NPT(self, benzene_system): + def test_set_ref_length(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.reference_length = 1 * u.angstrom + assert sim.reference_length == 1 * u.angstrom + + def test_set_ref_energy(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + ) + sim.reference_energy = 3.0 * u.kcal / u.mol + assert sim.reference_energy == 3.0 * u.kcal / u.mol + + def test_set_ref_mass(self, benzene_system): sim = Simulation( initial_state=benzene_system.hoomd_snapshot, forcefield=benzene_system.hoomd_forcefield, ) + sim.reference_mass = 1.25 * u.amu + assert sim.reference_mass == 1.25 * u.amu + + def test_NVT(self, benzene_system): + sim = Simulation.from_system(benzene_system) + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) + assert isinstance(sim.method, hoomd.md.methods.NVT) + + def test_NPT(self, benzene_system): + sim = Simulation.from_system(benzene_system) sim.run_NPT( kT=1.0, n_steps=500, @@ -57,61 +112,99 @@ def test_NPT(self, benzene_system): assert isinstance(sim.method, hoomd.md.methods.NPT) def test_langevin(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.run_langevin(n_steps=500, kT=1.0, alpha=0.5) assert isinstance(sim.method, hoomd.md.methods.Langevin) def test_NVE(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.run_NVE(n_steps=500) assert isinstance(sim.method, hoomd.md.methods.NVE) def test_displacement_cap(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.run_displacement_cap(n_steps=500, maximum_displacement=1e-4) assert isinstance(sim.method, hoomd.md.methods.DisplacementCapped) - def test_update_volume(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + def test_update_volume_target_box(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_box_lengths=sim.box_lengths * 0.5, + final_box_lengths=sim.box_lengths_reduced * 0.5, ) def test_update_volume_walls(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.add_walls(wall_axis=(1, 0, 0), sigma=1.0, epsilon=1.0, r_cut=1.12) sim.run_update_volume( kT=1.0, tau_kt=0.01, n_steps=500, period=5, - final_box_lengths=sim.box_lengths * 0.5, + final_box_lengths=sim.box_lengths_reduced * 0.5, ) - def test_change_methods(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, + 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 ) + assert np.isclose( + sim.density.to(u.g / u.cm**3).value, + (0.1 * (u.g / u.cm**3)).value, + atol=1e-4, + ) + + def test_update_volume_by_density_factor(self, benzene_system): + sim = Simulation.from_system(benzene_system) + init_density = copy.deepcopy(sim.density) + sim.run_update_volume( + kT=1.0, + tau_kt=0.01, + n_steps=500, + period=1, + final_density=sim.density * 5, + ) + assert np.isclose( + sim.density.value, (init_density * 5).value, atol=1e-4 + ) + + def test_update_volume_missing_values(self, benzene_system): + sim = Simulation.from_system(benzene_system) + with pytest.raises(ValueError): + sim.run_update_volume(kT=1.0, tau_kt=0.01, n_steps=500, period=1) + + def test_update_volume_two_values(self, benzene_system): + sim = Simulation.from_system(benzene_system) + with pytest.raises(ValueError): + sim.run_update_volume( + kT=1.0, + tau_kt=0.01, + n_steps=500, + period=1, + final_box_lengths=sim.box_lengths_reduced * 0.5, + final_density=0.1, + ) + + # def test_update_volume_with_density_no_ref_values(self, benzene_system): + # sim_no_ref = Simulation( + # initial_state=benzene_system.hoomd_snapshot, + # forcefield=benzene_system.hoomd_forcefield, + # ) + # with pytest.raises(ReferenceUnitError): + # sim_no_ref.run_update_volume( + # kT=1.0, + # tau_kt=0.01, + # n_steps=500, + # period=1, + # final_density=0.1, + # ) + + def test_change_methods(self, benzene_system): + sim = Simulation.from_system(benzene_system) sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) assert isinstance(sim.method, hoomd.md.methods.NVT) sim.run_NPT( @@ -120,20 +213,14 @@ def test_change_methods(self, benzene_system): assert isinstance(sim.method, hoomd.md.methods.NPT) def test_change_dt(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) sim.dt = 0.003 sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0) assert sim.dt == 0.003 def test_scale_epsilon(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) epsilons = [] for param in sim._lj_force().params: epsilons.append(sim._lj_force().params[param]["epsilon"]) @@ -145,10 +232,7 @@ def test_scale_epsilon(self, benzene_system): assert np.allclose(i * 0.5, j, atol=1e-3) def test_shift_epsilon(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) epsilons = [] for param in sim._lj_force().params: epsilons.append(sim._lj_force().params[param]["epsilon"]) @@ -160,10 +244,7 @@ def test_shift_epsilon(self, benzene_system): assert np.allclose(i + 1, j, atol=1e-3) def test_scale_sigma(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sigmas = [] for param in sim._lj_force().params: sigmas.append(sim._lj_force().params[param]["sigma"]) @@ -175,10 +256,7 @@ def test_scale_sigma(self, benzene_system): assert np.allclose(i * 0.5, j, atol=1e-3) def test_shift_sigma(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sigmas = [] for param in sim._lj_force().params: sigmas.append(sim._lj_force().params[param]["sigma"]) @@ -190,19 +268,13 @@ def test_shift_sigma(self, benzene_system): assert np.allclose(i + 1, j, atol=1e-3) def test_remove_force(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.remove_force(sim._lj_force()) for i in sim.forces: assert not isinstance(i, hoomd.md.pair.LJ) def test_set_integrate_group(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) assert isinstance(sim.integrate_group, hoomd.filter.All) tag_filter = hoomd.filter.Tags([0, 1, 2, 3]) sim.integrate_group = tag_filter @@ -210,10 +282,7 @@ def test_set_integrate_group(self, benzene_system): sim.run_NVT(n_steps=200, kT=1.0, tau_kt=0.01) def test_pickle_ff(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.pickle_forcefield("forcefield.pickle") assert os.path.isfile("forcefield.pickle") f = open("forcefield.pickle", "rb") @@ -224,15 +293,14 @@ def test_pickle_ff(self, benzene_system): os.remove("forcefield.pickle") def test_save_restart_gsd(self, benzene_system): - sim = Simulation( - initial_state=benzene_system.hoomd_snapshot, - forcefield=benzene_system.hoomd_forcefield, - ) + sim = Simulation.from_system(benzene_system) sim.save_restart_gsd("restart.gsd") assert os.path.isfile("restart.gsd") sim.pickle_forcefield("forcefield.pickle") f = open("forcefield.pickle", "rb") hoomd_ff = pickle.load(f) - Simulation(initial_state="restart.gsd", forcefield=hoomd_ff) + Simulation.from_snapshot_forces( + initial_state="restart.gsd", forcefield=hoomd_ff + ) os.remove("forcefield.pickle") os.remove("restart.gsd") diff --git a/hoomd_organics/tests/base_test.py b/hoomd_organics/tests/base_test.py index a7f1c241..c6666dbe 100644 --- a/hoomd_organics/tests/base_test.py +++ b/hoomd_organics/tests/base_test.py @@ -6,12 +6,16 @@ from gmso.external.convert_mbuild import from_mbuild from hoomd_organics import Molecule, Pack, Polymer, Simulation -from hoomd_organics.library import OPLS_AA +from hoomd_organics.library import OPLS_AA, BeadSpring ASSETS_DIR = os.path.join(os.path.dirname(__file__), "assets") class BaseTest: + @pytest.fixture(autouse=True) + def initdir(self, tmpdir): + tmpdir.chdir() + @pytest.fixture() def benzene_smiles(self): return "c1ccccc1" @@ -209,16 +213,28 @@ def __init__(self, lengths, num_mols, **kwargs): @pytest.fixture() def benzene_system(self, benzene_mb): - benzene = Molecule(num_mols=5, compound=benzene_mb) + benzene = Molecule(num_mols=20, compound=benzene_mb) system = Pack( molecules=[benzene], - density=0.5, + density=0.2, r_cut=2.5, force_field=OPLS_AA(), auto_scale=True, ) return system + @pytest.fixture() + def benzene_cg_system(self, benzene_molecule, benzene_smiles): + benzene_mols = benzene_molecule(n_mols=10) + benzene_mols.coarse_grain(beads={"A": benzene_smiles}) + system = Pack( + molecules=[benzene_mols], + density=0.5, + r_cut=2.5, + auto_scale=False, + ) + return system + @pytest.fixture() def polyethylene_system(self, polyethylene): polyethylene_mol = polyethylene(num_mols=5, lengths=5) @@ -239,3 +255,13 @@ def benzene_simulation(self, benzene_system): forcefield=benzene_system.hoomd_forcefield, ) return sim + + @pytest.fixture() + def cg_single_bead_ff(self): + ff = BeadSpring( + r_cut=2.5, + beads={ + "A": dict(epsilon=1.0, sigma=1.0), + }, + ) + return ff.hoomd_forcefield diff --git a/hoomd_organics/tests/utils/test_utils.py b/hoomd_organics/tests/utils/test_utils.py index b132cda7..16d38462 100644 --- a/hoomd_organics/tests/utils/test_utils.py +++ b/hoomd_organics/tests/utils/test_utils.py @@ -1,7 +1,11 @@ import pytest import unyt as u -from hoomd_organics.utils import check_return_iterable, validate_ref_value +from hoomd_organics.utils import ( + calculate_box_length, + check_return_iterable, + validate_ref_value, +) from hoomd_organics.utils.exceptions import ReferenceUnitError @@ -38,3 +42,23 @@ def test_validate_ref_value(self): with pytest.raises(ValueError): validate_ref_value("test g", u.dimensions.mass) + + def test_calculate_box_length(self): + mass = u.unyt_quantity(4.0, u.g) + density = u.unyt_quantity(0.5, u.g / u.cm**3) + box_length = calculate_box_length(mass, density) + assert box_length == 2.0 * u.cm + + def test_calculate_box_length_fixed_l_1d(self): + mass = u.unyt_quantity(6.0, u.g) + density = u.unyt_quantity(0.5, u.g / u.cm**3) + fixed_L = u.unyt_quantity(3.0, u.cm) + box_length = calculate_box_length(mass, density, fixed_L=fixed_L) + assert box_length == 2.0 * u.cm + + def test_calculate_box_length_fixed_l_2d(self): + mass = u.unyt_quantity(12.0, u.g) + density = u.unyt_quantity(0.5, u.g / u.cm**3) + fixed_L = u.unyt_array([3.0, 2.0], u.cm) + box_length = calculate_box_length(mass, density, fixed_L=fixed_L) + assert box_length == 4.0 * u.cm diff --git a/hoomd_organics/utils/__init__.py b/hoomd_organics/utils/__init__.py index 116f22ac..9558082f 100644 --- a/hoomd_organics/utils/__init__.py +++ b/hoomd_organics/utils/__init__.py @@ -1,4 +1,8 @@ from .actions import * from .base_types import FF_Types from .ff_utils import xml_to_gmso_ff -from .utils import check_return_iterable, validate_ref_value +from .utils import ( + calculate_box_length, + check_return_iterable, + validate_ref_value, +) diff --git a/hoomd_organics/utils/utils.py b/hoomd_organics/utils/utils.py index ead926e0..75826518 100644 --- a/hoomd_organics/utils/utils.py +++ b/hoomd_organics/utils/utils.py @@ -1,3 +1,4 @@ +import numpy as np import unyt as u from hoomd_organics.utils.exceptions import ReferenceUnitError @@ -89,3 +90,39 @@ def _is_float(num): f"a reference value with unit of " f"{dimension} dimension." ) + + +def calculate_box_length(mass, density, fixed_L=None): + """Calculates the required box length(s) given the + mass of a sytem and the target density. + + Box edge length constraints can be set by set_target_box(). + If constraints are set, this will solve for the required + lengths of the remaining non-constrained edges to match + the target density. + + Parameters + ---------- + mass : unyt.unyt_quantity, required + Mass of the system + density : unyt.unyt_quantity, required + Target density of the system + fixed_L : np.array, optional, defualt=None + Array of fixed box lengths to be accounted for + when solving for L + + Returns + ------- + L : float + Box edge length + """ + + vol = mass / density # cm^3 + if fixed_L is None: + L = vol ** (1 / 3) + else: + L = vol / np.prod(fixed_L) + if fixed_L.size == 1: # L is cm^2 + L = L ** (1 / 2) + # L is cm + return L