From d38aba8038237401324b891bf020a721eaadd389 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 15 Mar 2024 13:02:56 -0600 Subject: [PATCH 01/13] calculate real T from kT --- flowermd/base/simulation.py | 129 +++++++++++++++++++++--------------- 1 file changed, 75 insertions(+), 54 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index c6225937..fe299705 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -54,19 +54,19 @@ class Simulation(hoomd.simulation.Simulation): """ def __init__( - self, - initial_state, - forcefield, - reference_values=dict(), - dt=0.0001, - device=hoomd.device.auto_select(), - seed=42, - gsd_write_freq=1e4, - gsd_file_name="trajectory.gsd", - gsd_max_buffer_size=64 * 1024 * 1024, - log_write_freq=1e3, - log_file_name="sim_data.txt", - thermostat=HOOMDThermostats.MTTK, + self, + initial_state, + forcefield, + reference_values=dict(), + dt=0.0001, + device=hoomd.device.auto_select(), + seed=42, + gsd_write_freq=1e4, + gsd_file_name="trajectory.gsd", + gsd_max_buffer_size=64 * 1024 * 1024, + log_write_freq=1e3, + log_file_name="sim_data.txt", + thermostat=HOOMDThermostats.MTTK, ): if not isinstance(forcefield, Iterable) or isinstance(forcefield, str): raise ValueError( @@ -101,6 +101,7 @@ def __init__( ] self.integrator = None self._dt = dt + self._kT = None self._reference_values = dict() self._reference_values = reference_values self._integrate_group = hoomd.filter.All() @@ -365,10 +366,26 @@ def real_timestep(self): energy = self.reference_energy.to("J") else: energy = 1 * u.J - tau = (mass * (dist**2)) / energy - timestep = self.dt * (tau**0.5) + tau = (mass * (dist ** 2)) / energy + 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 + @property def integrate_group(self): """The group of particles to apply the integrator to. @@ -417,7 +434,7 @@ def thermostat(self, thermostat): The type of thermostat to use. """ if not issubclass( - self._thermostat, hoomd.md.methods.thermostats.Thermostat + self._thermostat, hoomd.md.methods.thermostats.Thermostat ): raise ValueError( f"Invalid thermostat. Please choose from: {HOOMDThermostats}" @@ -594,14 +611,14 @@ def remove_walls(self, wall_axis): self.remove_force(wall_force) def run_update_volume( - self, - final_box_lengths, - n_steps, - period, - kT, - tau_kt, - thermalize_particles=True, - write_at_start=True, + self, + final_box_lengths, + n_steps, + period, + kT, + tau_kt, + thermalize_particles=True, + write_at_start=True, ): """Run an NVT simulation while shrinking or expanding simulation box. @@ -662,6 +679,7 @@ def run_update_volume( ) """ + self._kT = kT 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) @@ -714,14 +732,14 @@ def run_update_volume( self.operations.updaters.remove(box_resizer) def run_langevin( - self, - n_steps, - kT, - tally_reservoir_energy=False, - default_gamma=1.0, - default_gamma_r=(1.0, 1.0, 1.0), - thermalize_particles=True, - write_at_start=True, + self, + n_steps, + kT, + tally_reservoir_energy=False, + default_gamma=1.0, + default_gamma_r=(1.0, 1.0, 1.0), + thermalize_particles=True, + write_at_start=True, ): """Run the simulation using the Langevin dynamics integrator. @@ -746,6 +764,7 @@ def run_langevin( time step. """ + self._kT = kT self.set_integrator_method( integrator_method=hoomd.md.methods.Langevin, method_kwargs={ @@ -768,18 +787,18 @@ def run_langevin( self.operations.updaters.remove(std_out_logger_printer) def run_NPT( - self, - n_steps, - kT, - pressure, - tau_kt, - tau_pressure, - couple="xyz", - box_dof=[True, True, True, False, False, False], - rescale_all=False, - gamma=0.0, - thermalize_particles=True, - write_at_start=True, + self, + n_steps, + kT, + pressure, + tau_kt, + tau_pressure, + couple="xyz", + box_dof=[True, True, True, False, False, False], + rescale_all=False, + gamma=0.0, + thermalize_particles=True, + write_at_start=True, ): """Run the simulation in the NPT ensemble. @@ -812,6 +831,7 @@ def run_NPT( time step. """ + self._kT = kT self.set_integrator_method( integrator_method=hoomd.md.methods.ConstantPressure, method_kwargs={ @@ -839,12 +859,12 @@ def run_NPT( self.operations.updaters.remove(std_out_logger_printer) def run_NVT( - self, - n_steps, - kT, - tau_kt, - thermalize_particles=True, - write_at_start=True, + self, + n_steps, + kT, + tau_kt, + thermalize_particles=True, + write_at_start=True, ): """Run the simulation in the NVT ensemble. @@ -864,6 +884,7 @@ def run_NVT( time step. """ + self._kT = kT self.set_integrator_method( integrator_method=hoomd.md.methods.ConstantVolume, method_kwargs={ @@ -911,10 +932,10 @@ def run_NVE(self, n_steps, write_at_start=True): self.operations.updaters.remove(std_out_logger_printer) def run_displacement_cap( - self, - n_steps, - maximum_displacement=1e-3, - write_at_start=True, + self, + n_steps, + maximum_displacement=1e-3, + write_at_start=True, ): """NVE integrator with a cap on the maximum displacement per time step. From 05d594e4268f66c851b02b1db4a5c5524a18147b Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 15 Mar 2024 13:03:17 -0600 Subject: [PATCH 02/13] unit test for real_temperature --- flowermd/tests/base/test_simulation.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index e0b99c51..43fcb060 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -24,7 +24,7 @@ def test_initialize_from_system(self, benzene_system): assert sim.reference_values == benzene_system.reference_values def test_initialize_from_system_separate_ff( - self, benzene_cg_system, cg_single_bead_ff + self, benzene_cg_system, cg_single_bead_ff ): sim = Simulation.from_system( benzene_cg_system, forcefield=cg_single_bead_ff @@ -183,8 +183,8 @@ def test_update_volume_density(self, benzene_system): final_box_lengths=target_box, ) assert np.isclose( - sim.density.to(u.g / u.cm**3).value, - (0.05 * (u.g / u.cm**3)).value, + sim.density.to(u.g / u.cm ** 3).value, + (0.05 * (u.g / u.cm ** 3)).value, atol=1e-4, ) @@ -382,3 +382,10 @@ 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=500) + assert np.isclose(sim.real_temperature, 35.225, atol=1e-4) From 0c47ee36f53f7dd1f0b75480a3ac73542d24c2b6 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 15 Mar 2024 13:05:58 -0600 Subject: [PATCH 03/13] precommit reformat --- flowermd/base/simulation.py | 108 ++++++++++++------------- flowermd/tests/base/test_simulation.py | 6 +- 2 files changed, 57 insertions(+), 57 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index fe299705..3a6952d5 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -54,19 +54,19 @@ class Simulation(hoomd.simulation.Simulation): """ def __init__( - self, - initial_state, - forcefield, - reference_values=dict(), - dt=0.0001, - device=hoomd.device.auto_select(), - seed=42, - gsd_write_freq=1e4, - gsd_file_name="trajectory.gsd", - gsd_max_buffer_size=64 * 1024 * 1024, - log_write_freq=1e3, - log_file_name="sim_data.txt", - thermostat=HOOMDThermostats.MTTK, + self, + initial_state, + forcefield, + reference_values=dict(), + dt=0.0001, + device=hoomd.device.auto_select(), + seed=42, + gsd_write_freq=1e4, + gsd_file_name="trajectory.gsd", + gsd_max_buffer_size=64 * 1024 * 1024, + log_write_freq=1e3, + log_file_name="sim_data.txt", + thermostat=HOOMDThermostats.MTTK, ): if not isinstance(forcefield, Iterable) or isinstance(forcefield, str): raise ValueError( @@ -366,8 +366,8 @@ def real_timestep(self): energy = self.reference_energy.to("J") else: energy = 1 * u.J - tau = (mass * (dist ** 2)) / energy - timestep = self.dt * (tau ** 0.5) + tau = (mass * (dist**2)) / energy + timestep = self.dt * (tau**0.5) return timestep @property @@ -434,7 +434,7 @@ def thermostat(self, thermostat): The type of thermostat to use. """ if not issubclass( - self._thermostat, hoomd.md.methods.thermostats.Thermostat + self._thermostat, hoomd.md.methods.thermostats.Thermostat ): raise ValueError( f"Invalid thermostat. Please choose from: {HOOMDThermostats}" @@ -611,14 +611,14 @@ def remove_walls(self, wall_axis): self.remove_force(wall_force) def run_update_volume( - self, - final_box_lengths, - n_steps, - period, - kT, - tau_kt, - thermalize_particles=True, - write_at_start=True, + self, + final_box_lengths, + n_steps, + period, + kT, + tau_kt, + thermalize_particles=True, + write_at_start=True, ): """Run an NVT simulation while shrinking or expanding simulation box. @@ -732,14 +732,14 @@ def run_update_volume( self.operations.updaters.remove(box_resizer) def run_langevin( - self, - n_steps, - kT, - tally_reservoir_energy=False, - default_gamma=1.0, - default_gamma_r=(1.0, 1.0, 1.0), - thermalize_particles=True, - write_at_start=True, + self, + n_steps, + kT, + tally_reservoir_energy=False, + default_gamma=1.0, + default_gamma_r=(1.0, 1.0, 1.0), + thermalize_particles=True, + write_at_start=True, ): """Run the simulation using the Langevin dynamics integrator. @@ -787,18 +787,18 @@ def run_langevin( self.operations.updaters.remove(std_out_logger_printer) def run_NPT( - self, - n_steps, - kT, - pressure, - tau_kt, - tau_pressure, - couple="xyz", - box_dof=[True, True, True, False, False, False], - rescale_all=False, - gamma=0.0, - thermalize_particles=True, - write_at_start=True, + self, + n_steps, + kT, + pressure, + tau_kt, + tau_pressure, + couple="xyz", + box_dof=[True, True, True, False, False, False], + rescale_all=False, + gamma=0.0, + thermalize_particles=True, + write_at_start=True, ): """Run the simulation in the NPT ensemble. @@ -859,12 +859,12 @@ def run_NPT( self.operations.updaters.remove(std_out_logger_printer) def run_NVT( - self, - n_steps, - kT, - tau_kt, - thermalize_particles=True, - write_at_start=True, + self, + n_steps, + kT, + tau_kt, + thermalize_particles=True, + write_at_start=True, ): """Run the simulation in the NVT ensemble. @@ -932,10 +932,10 @@ def run_NVE(self, n_steps, write_at_start=True): self.operations.updaters.remove(std_out_logger_printer) def run_displacement_cap( - self, - n_steps, - maximum_displacement=1e-3, - write_at_start=True, + self, + n_steps, + maximum_displacement=1e-3, + write_at_start=True, ): """NVE integrator with a cap on the maximum displacement per time step. diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index 43fcb060..92960f06 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -24,7 +24,7 @@ def test_initialize_from_system(self, benzene_system): assert sim.reference_values == benzene_system.reference_values def test_initialize_from_system_separate_ff( - self, benzene_cg_system, cg_single_bead_ff + self, benzene_cg_system, cg_single_bead_ff ): sim = Simulation.from_system( benzene_cg_system, forcefield=cg_single_bead_ff @@ -183,8 +183,8 @@ def test_update_volume_density(self, benzene_system): final_box_lengths=target_box, ) assert np.isclose( - sim.density.to(u.g / u.cm ** 3).value, - (0.05 * (u.g / u.cm ** 3)).value, + sim.density.to(u.g / u.cm**3).value, + (0.05 * (u.g / u.cm**3)).value, atol=1e-4, ) From d2bad9980ad3e525c648d2c53b51322242dfb063 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 15 Mar 2024 14:40:19 -0600 Subject: [PATCH 04/13] allow temperature to be passed to run methods instead of kT --- flowermd/base/simulation.py | 81 ++++++++++++++++++++++++++++--------- 1 file changed, 62 insertions(+), 19 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 3a6952d5..7da31c08 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -386,6 +386,36 @@ def real_temperature(self): 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 kT.value + + def _setup_temperature(self, kT=None, temperature=None): + if kT and temperature: + raise ValueError( + "Both kT and temperature are provided. Please provide only one." + ) + if not kT and not temperature: + raise ValueError( + "Either kT (unitless) or temperature (with units) must be " + "provided for the simulation." + ) + if kT: + return kT + else: + return self._temperature_to_kT(temperature) + @property def integrate_group(self): """The group of particles to apply the integrator to. @@ -615,8 +645,9 @@ def run_update_volume( final_box_lengths, n_steps, period, - kT, tau_kt, + kT=None, + temperature=None, thermalize_particles=True, write_at_start=True, ): @@ -679,7 +710,7 @@ def run_update_volume( ) """ - self._kT = kT + self._kT = self._setup_temperature(kT, temperature) 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) @@ -707,13 +738,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) @@ -734,7 +765,8 @@ def run_update_volume( def run_langevin( self, n_steps, - kT, + kT=None, + temperature=None, tally_reservoir_energy=False, default_gamma=1.0, default_gamma_r=(1.0, 1.0, 1.0), @@ -764,19 +796,19 @@ def run_langevin( time step. """ - self._kT = kT + self._kT = self._setup_temperature(kT, temperature) 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) + 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), @@ -789,10 +821,11 @@ def run_langevin( def run_NPT( self, n_steps, - kT, pressure, - tau_kt, tau_pressure, + tau_kt, + kT=None, + temperature=None, couple="xyz", box_dof=[True, True, True, False, False, False], rescale_all=False, @@ -831,7 +864,7 @@ def run_NPT( time step. """ - self._kT = kT + self._kT = self._setup_temperature(kT, temperature) self.set_integrator_method( integrator_method=hoomd.md.methods.ConstantPressure, method_kwargs={ @@ -843,12 +876,12 @@ 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) + 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), @@ -861,8 +894,9 @@ def run_NPT( def run_NVT( self, n_steps, - kT, tau_kt, + kT=None, + temperature=None, thermalize_particles=True, write_at_start=True, ): @@ -884,18 +918,18 @@ def run_NVT( time step. """ - self._kT = kT + self._kT = self._setup_temperature(kT, temperature) 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) + 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), @@ -974,7 +1008,14 @@ def run_displacement_cap( 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, + kT_start=None, + temperature_start=None, + kT_final=None, + temperature_final=None, + ): """Create a temperature ramp. Parameters @@ -987,8 +1028,10 @@ def temperature_ramp(self, n_steps, kT_start, kT_final): The final temperature. """ + _kT_start = self._setup_temperature(kT_start, temperature_start) + _kT_final = self._setup_temperature(kT_final, temperature_final) 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"): From ee71094b1a3bd275db68e86d0b2737966c018b2f Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 15 Mar 2024 14:40:44 -0600 Subject: [PATCH 05/13] test temperature for nvt runs --- flowermd/tests/base/test_simulation.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index 92960f06..3e74fa14 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -389,3 +389,18 @@ def test_real_temperature(self, benzene_system): sim.real_temperature sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500) assert np.isclose(sim.real_temperature, 35.225, atol=1e-4) + + 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_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) From 27a2d7483246c73d2ba1b8478d9471cf04fe502a Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 15 Mar 2024 14:58:20 -0600 Subject: [PATCH 06/13] add time_length parameter to run methods --- flowermd/base/simulation.py | 89 +++++++++++++++++++++++++++---------- 1 file changed, 66 insertions(+), 23 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 7da31c08..f4167d9c 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -399,7 +399,7 @@ def _temperature_to_kT(self, temperature): else: temperature = temperature * u.K kT = (temperature * u.boltzmann_constant_mks) / energy - return kT.value + return float(kT) def _setup_temperature(self, kT=None, temperature=None): if kT and temperature: @@ -408,7 +408,7 @@ def _setup_temperature(self, kT=None, temperature=None): ) if not kT and not temperature: raise ValueError( - "Either kT (unitless) or temperature (with units) must be " + "Either kT or temperature must be " "provided for the simulation." ) if kT: @@ -416,6 +416,33 @@ def _setup_temperature(self, kT=None, temperature=None): 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 and time_length: + raise ValueError( + "Both n_steps and time_length are provided. Please provide only" + " one." + ) + if not n_steps and not time_length: + raise ValueError( + "Either n_steps or time_length must be provided for the " + "simulation." + ) + if n_steps: + 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. @@ -643,9 +670,10 @@ def remove_walls(self, wall_axis): def run_update_volume( self, final_box_lengths, - n_steps, period, tau_kt, + n_steps=None, + time_length=None, kT=None, temperature=None, thermalize_particles=True, @@ -711,6 +739,7 @@ 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) @@ -723,7 +752,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 @@ -752,19 +781,20 @@ 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, + n_steps=None, + time_length=None, kT=None, temperature=None, tally_reservoir_energy=False, @@ -797,6 +827,7 @@ def run_langevin( """ 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={ @@ -809,23 +840,24 @@ def run_langevin( ) if thermalize_particles: self._thermalize_system(self._kT) - 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_NPT( self, - n_steps, pressure, 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, @@ -865,6 +897,7 @@ def run_NPT( """ 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={ @@ -882,21 +915,22 @@ def run_NPT( ) if thermalize_particles: self._thermalize_system(self._kT) - 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_NVT( self, - n_steps, tau_kt, kT=None, temperature=None, + n_steps=None, + time_length=None, thermalize_particles=True, write_at_start=True, ): @@ -919,6 +953,7 @@ def run_NVT( """ 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={ @@ -930,16 +965,16 @@ def run_NVT( ) if thermalize_particles: self._thermalize_system(self._kT) - 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_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 @@ -952,22 +987,24 @@ def run_NVE(self, n_steps, write_at_start=True): time step. """ + _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, ): @@ -992,6 +1029,7 @@ def run_displacement_cap( time step. """ + _n_steps = self._setup_n_steps(n_steps, time_length) self.set_integrator_method( integrator_method=hoomd.md.methods.DisplacementCapped, method_kwargs={ @@ -999,18 +1037,19 @@ 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, + n_steps=None, + time_length=None, kT_start=None, temperature_start=None, kT_final=None, @@ -1030,8 +1069,12 @@ def temperature_ramp( """ _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"): From 17b6043d8ec1f7926bc4a3f943668c2d81d2a403 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Mon, 18 Mar 2024 14:34:53 -0600 Subject: [PATCH 07/13] fix None conditions --- flowermd/base/simulation.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index f4167d9c..01b4f2a5 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -402,16 +402,16 @@ def _temperature_to_kT(self, temperature): return float(kT) def _setup_temperature(self, kT=None, temperature=None): - if kT and temperature: + if kT is not None and temperature is not None: raise ValueError( "Both kT and temperature are provided. Please provide only one." ) - if not kT and not temperature: + if kT is None and temperature is None: raise ValueError( "Either kT or temperature must be " "provided for the simulation." ) - if kT: + if kT is not None: return kT else: return self._temperature_to_kT(temperature) @@ -428,17 +428,17 @@ def _time_length_to_n_steps(self, time_length): return int(time_length / real_timestep) def _setup_n_steps(self, n_steps=None, time_length=None): - if n_steps and time_length: + 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 not n_steps and not time_length: + 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: + if n_steps is not None: return n_steps else: return self._time_length_to_n_steps(time_length) From 822861fd43f3bebd252c2457ebbd71cc119fa200 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Wed, 20 Mar 2024 14:46:39 -0600 Subject: [PATCH 08/13] add tests for time length --- flowermd/tests/base/test_simulation.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index 3e74fa14..1cc66000 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -404,3 +404,17 @@ 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_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) From b422a11070bdbccb82746eb09988bc8d7f8fb3a3 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Wed, 20 Mar 2024 15:02:05 -0600 Subject: [PATCH 09/13] update docstrings --- flowermd/base/simulation.py | 125 +++++++++++++++++++++++++++++------- 1 file changed, 103 insertions(+), 22 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 01b4f2a5..a859f4ec 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -692,19 +692,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` @@ -807,10 +820,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. @@ -825,6 +844,13 @@ 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) @@ -869,16 +895,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; @@ -895,6 +927,13 @@ 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) @@ -938,12 +977,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: 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. thermalize_particles: bool, default True When set to True, assigns random velocities to all particles. write_at_start : bool, default True @@ -951,6 +996,13 @@ 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) @@ -979,13 +1031,21 @@ def run_NVE(self, n_steps=None, time_length=None, write_at_start=True): 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( @@ -1018,16 +1078,23 @@ 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( @@ -1059,12 +1126,26 @@ def 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) From e83b3c40eccc8c22c14580ca3fffadd83ba77e4e Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 21 Mar 2024 10:53:26 -0600 Subject: [PATCH 10/13] save reference values to pickle --- flowermd/base/system.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/flowermd/base/system.py b/flowermd/base/system.py index 7bc08545..5dbc1b39 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) From e4ce1c64462e9c10aed967f38bf9b7b72833bf26 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 21 Mar 2024 10:54:00 -0600 Subject: [PATCH 11/13] tests for save ref values --- flowermd/tests/base/test_system.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/flowermd/tests/base/test_system.py b/flowermd/tests/base/test_system.py index aa2accf8..e4ed1669 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 @@ -882,6 +883,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) From 4d91a492e4dc9ee29894730bdb3faf77a9318534 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 21 Mar 2024 10:58:17 -0600 Subject: [PATCH 12/13] added more unit tests for temperature and time length --- flowermd/tests/base/test_simulation.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index 1cc66000..8cc321df 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -395,6 +395,11 @@ def test_NVT_with_temperature(self, 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): @@ -410,6 +415,11 @@ def test_NVT_time_length(self, 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): From 26af373ea693f16e5af4ae4494145ee743925d7b Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 21 Mar 2024 11:33:15 -0600 Subject: [PATCH 13/13] real temperature without energy unit --- flowermd/base/simulation.py | 2 +- flowermd/tests/base/test_simulation.py | 11 ++++++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index a859f4ec..ffe7a927 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -979,7 +979,7 @@ def run_NVT( ---------- tau_kt: float, required Thermostat coupling period (in simulation time units). - kT: int or hoomd.variant.Ramp, optional + 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 diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index 8cc321df..391681f2 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -387,9 +387,18 @@ 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=500) + 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)