Skip to content

Commit

Permalink
Merge pull request #42 from marjanAlbouye/update-sim-interface
Browse files Browse the repository at this point in the history 
Update simulation interface and add density to update_volume
  • Loading branch information
@marjanalbooyeh
marjanalbooyeh authored Sep 15, 2023
2 parents 4a626f9 + 332b46a commit 15fd092
Show file tree
Hide file tree
Showing 8 changed files with 371 additions and 158 deletions.
144 changes: 90 additions & 54 deletions hoomd_organics/base/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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(),
Expand Down Expand Up @@ -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:
Expand All @@ -113,63 +150,27 @@ 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):
ref_keys = ["length", "mass", "energy"]
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):
Expand Down Expand Up @@ -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
Expand All @@ -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,
Expand Down
72 changes: 45 additions & 27 deletions hoomd_organics/base/system.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 = []
Expand All @@ -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:
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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:
Expand All @@ -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
Expand All @@ -374,25 +381,36 @@ 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):
"""Uses PACKMOL via mbuild.packing.fill_box.
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
Expand Down Expand Up @@ -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,
)
Expand Down
Empty file added hoomd_organics/tests/base/sim_data.txt
View file
Empty file.
Loading

0 comments on commit 15fd092

Please sign in to comment.