diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 04ba2d44..a3daa130 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -102,6 +102,7 @@ def __init__( ] self.integrator = None self._dt = dt + self._kT = None self._reference_values = dict() self._reference_values = reference_values if rigid_constraint and not isinstance( @@ -389,6 +390,79 @@ def real_timestep(self): timestep = self.dt * (tau**0.5) return timestep + @property + def real_temperature(self): + """The temperature of the simulation in Kelvin.""" + if not self._kT: + raise ValueError( + "Temperature is not set. Please specify the temperature when " + "running the simulation, using one of the following run" + " methods: `run_nvt`, `run_npt`, `run_update_volume`." + ) + if self._reference_values.get("energy"): + energy = self.reference_energy.to("J") + else: + energy = 1 * u.J + temperature = (self._kT * energy) / u.boltzmann_constant_mks + return temperature + + def _temperature_to_kT(self, temperature): + """Convert temperature to kT.""" + if self._reference_values.get("energy"): + energy = self.reference_energy.to("J") + else: + energy = 1 * u.J + if isinstance(temperature, u.unyt_array) or isinstance( + temperature, u.unyt_quantity + ): + temperature = temperature.to("K") + else: + temperature = temperature * u.K + kT = (temperature * u.boltzmann_constant_mks) / energy + return float(kT) + + def _setup_temperature(self, kT=None, temperature=None): + if kT is not None and temperature is not None: + raise ValueError( + "Both kT and temperature are provided. Please provide only one." + ) + if kT is None and temperature is None: + raise ValueError( + "Either kT or temperature must be " + "provided for the simulation." + ) + if kT is not None: + return kT + else: + return self._temperature_to_kT(temperature) + + def _time_length_to_n_steps(self, time_length): + """Convert time length to number of steps.""" + if isinstance(time_length, u.unyt_array) or isinstance( + time_length, u.unyt_quantity + ): + time_length = time_length.to("s") + else: + time_length = time_length * u.s + real_timestep = self.real_timestep.to("s") + return int(time_length / real_timestep) + + def _setup_n_steps(self, n_steps=None, time_length=None): + if n_steps is not None and time_length is not None: + raise ValueError( + "Both n_steps and time_length are provided. Please provide only" + " one." + ) + if n_steps is None and time_length is None: + raise ValueError( + "Either n_steps or time_length must be provided for the " + "simulation." + ) + if n_steps is not None: + return n_steps + else: + return self._time_length_to_n_steps(time_length) + @property def integrate_group(self): """The group of particles to apply the integrator to. @@ -623,10 +697,12 @@ def remove_walls(self, wall_axis): def run_update_volume( self, final_box_lengths, - n_steps, period, - kT, tau_kt, + n_steps=None, + time_length=None, + kT=None, + temperature=None, thermalize_particles=True, write_at_start=True, ): @@ -643,19 +719,32 @@ def run_update_volume( ---------- final_box_lengths : np.ndarray or unyt.array.unyt_array, shape=(3,), required # noqa: E501 The final box edge lengths in (x, y, z) order. - n_steps : int, required - Number of steps to run during volume update. period : int, required The number of steps ran between each box update iteration. - kT : float or hoomd.variant.Ramp, required - The temperature to use during volume update. tau_kt : float, required Thermostat coupling period (in simulation time units). + n_steps : int, optional + Number of steps to run during volume update. + time_length : unyt.unyt_quantity or float, optional + The length of time to run the simulation. If no unit is provided, + the time is assumed to be in seconds. + kT : float or hoomd.variant.Ramp, optional + The temperature to use during volume update. + temperature : unyt.unyt_quantity or float, optional + The temperature to use during volume update. If no unit is provided, + Kelvin is assumed. write_at_start : bool, default True When set to True, triggers writers that evaluate to True for the initial step to execute before the next simulation time step. + Notes + ----- + For the temperature, either `kT` or `temperature` must be provided. + If both are provided, an error will be raised. And for the number of + steps, either `n_steps` or `time_length` must be provided. If both are + provided, an error will be raised. + Examples -------- In this example, a low density system is initialized with `Pack` @@ -689,6 +778,8 @@ def run_update_volume( ) """ + self._kT = self._setup_temperature(kT, temperature) + _n_steps = self._setup_n_steps(n_steps, time_length) if self.reference_length and hasattr(final_box_lengths, "to"): ref_unit = self.reference_length.units final_box_lengths = final_box_lengths.to(ref_unit) @@ -701,7 +792,7 @@ def run_update_volume( ) resize_trigger = hoomd.trigger.Periodic(period) box_ramp = hoomd.variant.Ramp( - A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps) + A=0, B=1, t_start=self.timestep, t_ramp=int(_n_steps) ) initial_box = self.state.box @@ -716,13 +807,13 @@ def run_update_volume( integrator_method=hoomd.md.methods.ConstantVolume, method_kwargs={ "thermostat": self._initialize_thermostat( - {"kT": kT, "tau": tau_kt} + {"kT": self._kT, "tau": tau_kt} ), "filter": self.integrate_group, }, ) if thermalize_particles: - self._thermalize_system(kT) + self._thermalize_system(self._kT) if self._wall_forces: wall_update = UpdateWalls(sim=self) @@ -730,20 +821,22 @@ def run_update_volume( trigger=resize_trigger, action=wall_update ) self.operations.updaters.append(wall_updater) - std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) + std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( trigger=hoomd.trigger.Periodic(self._std_out_freq), action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) - self.run(steps=n_steps + 1, write_at_start=write_at_start) + self.run(steps=_n_steps + 1, write_at_start=write_at_start) self.operations.updaters.remove(std_out_logger_printer) self.operations.updaters.remove(box_resizer) def run_langevin( self, - n_steps, - kT, + n_steps=None, + time_length=None, + kT=None, + temperature=None, tally_reservoir_energy=False, default_gamma=1.0, default_gamma_r=(1.0, 1.0, 1.0), @@ -754,10 +847,16 @@ def run_langevin( Parameters ---------- - n_steps : int, required + n_steps : int, optional Number of steps to run the simulation. - kT : int or hoomd.variant.Ramp, required + time_length : unyt.unyt_quantity or float, optional + The length of time to run the simulation. If no unit is provided, + the time is assumed to be in seconds. + kT : int or hoomd.variant.Ramp, optional The temperature to use during the simulation. + temperature : unyt.unyt_quantity or float, optional + The temperature to use during the simulation. If no unit is + provided, Kelvin is assumed. tally_reservoir_energy : bool, default False When set to True, energy exchange between the thermal reservoir and the particles is tracked. @@ -772,35 +871,46 @@ def run_langevin( for the initial step to execute before the next simulation time step. + Notes + ----- + For the temperature, either `kT` or `temperature` must be provided. + If both are provided, an error will be raised. And for the number of + steps, either `n_steps` or `time_length` must be provided. If both are + provided, an error will be raised. + """ + self._kT = self._setup_temperature(kT, temperature) + _n_steps = self._setup_n_steps(n_steps, time_length) self.set_integrator_method( integrator_method=hoomd.md.methods.Langevin, method_kwargs={ "filter": self.integrate_group, - "kT": kT, + "kT": self._kT, "tally_reservoir_energy": tally_reservoir_energy, "default_gamma": default_gamma, "default_gamma_r": default_gamma_r, }, ) if thermalize_particles: - self._thermalize_system(kT) - std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) + self._thermalize_system(self._kT) + std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( trigger=hoomd.trigger.Periodic(self._std_out_freq), action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) - self.run(steps=n_steps, write_at_start=write_at_start) + self.run(steps=_n_steps, write_at_start=write_at_start) self.operations.updaters.remove(std_out_logger_printer) def run_NPT( self, - n_steps, - kT, pressure, - tau_kt, tau_pressure, + tau_kt, + kT=None, + temperature=None, + n_steps=None, + time_length=None, couple="xyz", box_dof=[True, True, True, False, False, False], rescale_all=False, @@ -812,16 +922,22 @@ def run_NPT( Parameters ---------- - n_steps: int, required - Number of steps to run the simulation. - kT: int or hoomd.variant.Ramp, required - The temperature to use during the simulation. pressure: int or hoomd.variant.Ramp, required The pressure to use during the simulation. - tau_kt: float, required - Thermostat coupling period (in simulation time units). tau_pressure: float, required Barostat coupling period. + tau_kt: float, required + Thermostat coupling period (in simulation time units). + kT: int or hoomd.variant.Ramp, optional + The temperature to use during the simulation. + temperature: unyt.unyt_quantity or float, optional + The temperature to use during the simulation. If no unit is + provided, Kelvin is assumed. + n_steps: int, optional + Number of steps to run the simulation. + time_length: unyt.unyt_quantity or float, optional + The length of time to run the simulation. If no unit is provided, + the time is assumed to be in seconds. couple: str, default "xyz" Couplings of diagonal elements of the stress tensor/ box_dof: list of bool; @@ -838,7 +954,16 @@ def run_NPT( for the initial step to execute before the next simulation time step. + Notes + ----- + For the temperature, either `kT` or `temperature` must be provided. + If both are provided, an error will be raised. And for the number of + steps, either `n_steps` or `time_length` must be provided. If both are + provided, an error will be raised. + """ + self._kT = self._setup_temperature(kT, temperature) + _n_steps = self._setup_n_steps(n_steps, time_length) self.set_integrator_method( integrator_method=hoomd.md.methods.ConstantPressure, method_kwargs={ @@ -850,26 +975,28 @@ def run_NPT( "gamma": gamma, "filter": self.integrate_group, "thermostat": self._initialize_thermostat( - {"kT": kT, "tau": tau_kt} + {"kT": self._kT, "tau": tau_kt} ), }, ) if thermalize_particles: - self._thermalize_system(kT) - std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) + self._thermalize_system(self._kT) + std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( trigger=hoomd.trigger.Periodic(self._std_out_freq), action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) - self.run(steps=n_steps, write_at_start=write_at_start) + self.run(steps=_n_steps, write_at_start=write_at_start) self.operations.updaters.remove(std_out_logger_printer) def run_NVT( self, - n_steps, - kT, tau_kt, + kT=None, + temperature=None, + n_steps=None, + time_length=None, thermalize_particles=True, write_at_start=True, ): @@ -877,12 +1004,18 @@ def run_NVT( Parameters ---------- - n_steps: int, required - Number of steps to run the simulation. - kT: int or hoomd.variant.Ramp, required - The temperature to use during the simulation. tau_kt: float, required Thermostat coupling period (in simulation time units). + kT: float or hoomd.variant.Ramp, optional + The temperature to use during the simulation. + temperature: unyt.unyt_quantity or float, optional + The temperature to use during the simulation. If no unit is + provided, Kelvin is assumed. + n_steps: int, optional + Number of steps to run the simulation. + time_length: unyt.unyt_quantity or float, optional + The length of time to run the simulation. If no unit is provided, + the time is assumed to be in seconds. thermalize_particles: bool, default True When set to True, assigns random velocities to all particles. write_at_start : bool, default True @@ -890,56 +1023,75 @@ def run_NVT( for the initial step to execute before the next simulation time step. + Notes + ----- + For the temperature, either `kT` or `temperature` must be provided. + If both are provided, an error will be raised. And for the number of + steps, either `n_steps` or `time_length` must be provided. If both are + provided, an error will be raised. + """ + self._kT = self._setup_temperature(kT, temperature) + _n_steps = self._setup_n_steps(n_steps, time_length) self.set_integrator_method( integrator_method=hoomd.md.methods.ConstantVolume, method_kwargs={ "thermostat": self._initialize_thermostat( - {"kT": kT, "tau": tau_kt} + {"kT": self._kT, "tau": tau_kt} ), "filter": self.integrate_group, }, ) if thermalize_particles: - self._thermalize_system(kT) - std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) + self._thermalize_system(self._kT) + std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( trigger=hoomd.trigger.Periodic(self._std_out_freq), action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) - self.run(steps=n_steps, write_at_start=write_at_start) + self.run(steps=_n_steps, write_at_start=write_at_start) self.operations.updaters.remove(std_out_logger_printer) - def run_NVE(self, n_steps, write_at_start=True): + def run_NVE(self, n_steps=None, time_length=None, write_at_start=True): """Run the simulation in the NVE ensemble. Parameters ---------- - n_steps: int, required + n_steps: int, optional Number of steps to run the simulation. + time_length: unyt.unyt_quantity or float, optional + The length of time to run the simulation. If no unit is provided, + the time is assumed to be in seconds. write_at_start : bool, default True When set to True, triggers writers that evaluate to True for the initial step to execute before the next simulation time step. + Notes + ----- + For the number of steps, either `n_steps` or `time_length` must be + provided. If both are provided, an error will be raised. + """ + _n_steps = self._setup_n_steps(n_steps, time_length) self.set_integrator_method( integrator_method=hoomd.md.methods.ConstantVolume, method_kwargs={"filter": self.integrate_group}, ) - std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) + std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( trigger=hoomd.trigger.Periodic(self._std_out_freq), action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) - self.run(steps=n_steps, write_at_start=write_at_start) + self.run(steps=_n_steps, write_at_start=write_at_start) self.operations.updaters.remove(std_out_logger_printer) def run_displacement_cap( self, - n_steps, + n_steps=None, + time_length=None, maximum_displacement=1e-3, write_at_start=True, ): @@ -953,17 +1105,25 @@ def run_displacement_cap( Parameters ---------- - n_steps : int, required + n_steps : int, optional Number of steps to run the simulation. + time_length : unyt.unyt_quantity or float, optional + The length of time to run the simulation. If no unit is provided, + the time is assumed to be in seconds. maximum_displacement : float, default 1e-3 Maximum displacement per step (length) - write_at_start : bool, default True When set to True, triggers writers that evaluate to True for the initial step to execute before the next simulation time step. + Notes + ----- + For the number of steps, either `n_steps` or `time_length` must be + provided. If both are provided, an error will be raised. + """ + _n_steps = self._setup_n_steps(n_steps, time_length) self.set_integrator_method( integrator_method=hoomd.md.methods.DisplacementCapped, method_kwargs={ @@ -971,30 +1131,58 @@ def run_displacement_cap( "maximum_displacement": maximum_displacement, }, ) - std_out_logger = StdOutLogger(n_steps=n_steps, sim=self) + std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self) std_out_logger_printer = hoomd.update.CustomUpdater( trigger=hoomd.trigger.Periodic(self._std_out_freq), action=std_out_logger, ) self.operations.updaters.append(std_out_logger_printer) - self.run(steps=n_steps, write_at_start=write_at_start) + self.run(steps=_n_steps, write_at_start=write_at_start) self.operations.updaters.remove(std_out_logger_printer) - def temperature_ramp(self, n_steps, kT_start, kT_final): + def temperature_ramp( + self, + n_steps=None, + time_length=None, + kT_start=None, + temperature_start=None, + kT_final=None, + temperature_final=None, + ): """Create a temperature ramp. Parameters ---------- - n_steps : int, required + n_steps : int, optional The number of steps to ramp the temperature over. - kT_start : float, required + time_length : unyt.unyt_quantity or float, optional + The length of time to ramp the temperature over. If no unit is + provided, the time is assumed to be in seconds. + kT_start : float, optional The starting temperature. - kT_final : float, required + temperature_start : unyt.unyt_quantity or float, optional + The starting temperature. If no unit is provided, Kelvin is assumed. + kT_final : float, optional The final temperature. + temperature_final : unyt.unyt_quantity or float, optional + The final temperature. If no unit is provided, Kelvin is assumed. + + Notes + ----- + For the temperature, either `kT` or `temperature` must be provided. + If both are provided, an error will be raised. And for the number of + steps, either `n_steps` or `time_length` must be provided. If both are + provided, an error will be raised. """ + _kT_start = self._setup_temperature(kT_start, temperature_start) + _kT_final = self._setup_temperature(kT_final, temperature_final) + _n_steps = self._setup_n_steps(n_steps, time_length) return hoomd.variant.Ramp( - A=kT_start, B=kT_final, t_start=self.timestep, t_ramp=int(n_steps) + A=_kT_start, + B=_kT_final, + t_start=self.timestep, + t_ramp=int(_n_steps), ) def pickle_forcefield(self, file_path="forcefield.pickle"): diff --git a/flowermd/base/system.py b/flowermd/base/system.py index 6416be53..22e3b626 100644 --- a/flowermd/base/system.py +++ b/flowermd/base/system.py @@ -404,6 +404,23 @@ def to_gsd(self, file_name): with gsd.hoomd.open(file_name, "w") as traj: traj.append(self.hoomd_snapshot) + def save_reference_values(self, file_path="reference_values.pickle"): + """Save the reference values of the system to a pickle file. + + Parameters + ---------- + file_path : str, default "reference_values.pickle" + The path to save the pickle file to. + + """ + if not self.reference_values: + raise ValueError( + "Reference values have not been set. " + "See System.reference_values" + ) + f = open(file_path, "wb") + pickle.dump(self.reference_values, f) + def _convert_to_gmso(self): """Convert the mbuild system to a gmso system.""" topology = from_mbuild(self.system) diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index a756f808..1bec5a01 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -425,3 +425,58 @@ def test_flush(self, benzene_system): with gsd.hoomd.open("trajectory.gsd") as traj: assert len(traj) > 0 os.remove("trajectory.gsd") + + def test_real_temperature(self, benzene_system): + sim = Simulation.from_system(benzene_system) + with pytest.raises(ValueError): + sim.real_temperature + sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=100) + assert np.isclose(sim.real_temperature, 35.225, atol=1e-4) + + def test_real_temperature_no_energy_units(self, benzene_system): + sim = Simulation( + initial_state=benzene_system.hoomd_snapshot, + forcefield=benzene_system.hoomd_forcefield, + reference_values=dict(), + ) + sim.run_NVT(kT=1e-10, tau_kt=0.01, n_steps=100) + assert np.isclose(sim.real_temperature, 7.2429e12) + + def test_NVT_with_temperature(self, benzene_system): + sim = Simulation.from_system(benzene_system) + sim.run_NVT(tau_kt=0.01, n_steps=500, temperature=35.225) + assert np.isclose(sim._kT, 1.0, atol=1e-4) + + def test_NVT_with_temperature_units(self, benzene_system): + sim = Simulation.from_system(benzene_system) + sim.run_NVT(tau_kt=0.01, n_steps=500, temperature=35.225 * u.K) + assert np.isclose(sim._kT, 1.0, atol=1e-4) + + def test_NVT_with_temperature_and_kT(self, benzene_system): + sim = Simulation.from_system(benzene_system) + with pytest.raises(ValueError): + sim.run_NVT(tau_kt=0.01, n_steps=500, temperature=35.225, kT=1.0) + + def test_NVT_no_temperature(self, benzene_system): + sim = Simulation.from_system(benzene_system) + with pytest.raises(ValueError): + sim.run_NVT(tau_kt=0.01, n_steps=500) + + def test_NVT_time_length(self, benzene_system): + sim = Simulation.from_system(benzene_system) + sim.run_NVT(kT=1.0, tau_kt=0.01, time_length=2 * u.Unit("ps")) + assert sim.timestep == int(2 * u.Unit("ps") / sim.real_timestep) + + def test_NVT_time_length_no_units(self, benzene_system): + sim = Simulation.from_system(benzene_system) + sim.run_NVT(kT=1.0, tau_kt=0.01, time_length=2e-12) + assert sim.timestep == int(2e-12 / sim.real_timestep) + + def test_NVT_time_length_n_steps(self, benzene_system): + sim = Simulation.from_system(benzene_system) + with pytest.raises(ValueError): + sim.run_NVT( + kT=1.0, tau_kt=0.01, time_length=2 * u.Unit("ps"), n_steps=1000 + ) + with pytest.raises(ValueError): + sim.run_NVT(kT=1.0, tau_kt=0.01) diff --git a/flowermd/tests/base/test_system.py b/flowermd/tests/base/test_system.py index cade5f96..a40e45bb 100644 --- a/flowermd/tests/base/test_system.py +++ b/flowermd/tests/base/test_system.py @@ -1,4 +1,5 @@ import os +import pickle import gmso import hoomd @@ -891,6 +892,35 @@ def test_pickle_ff_not_ff(self, polyethylene): with pytest.raises(ValueError): system.pickle_forcefield("forcefield.pickle") + def test_save_reference_values(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=1) + system = Pack( + molecules=[polyethylene], + density=1.0, + ) + system.apply_forcefield( + r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True + ) + system.save_reference_values() + assert os.path.isfile( + os.path.join(os.getcwd(), "reference_values.pickle") + ) + ref_values = pickle.load( + open(os.path.join(os.getcwd(), "reference_values.pickle"), "rb") + ) + assert ref_values == system.reference_values + os.remove(os.path.join(os.getcwd(), "reference_values.pickle")) + + def test_save_empty_reference_values(self, polyethylene): + polyethylene = polyethylene(lengths=5, num_mols=1) + system = Pack( + molecules=[polyethylene], + density=1.0, + base_units=dict(), + ) + with pytest.raises(ValueError): + system.save_reference_values() + def test_mass_density(self, benzene_molecule): benzene_mol = benzene_molecule(n_mols=100) system = Pack(molecules=[benzene_mol], density=0.8 * u.g / u.cm**3)