Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Refactor apply_forcefield #61

Merged
merged 10 commits into from
Sep 25, 2023
Merged

Refactor apply_forcefield #61

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
33fa89c
move parameters to apply_forcefield
chrisjonesBSU Sep 22, 2023
f407bc9
change remove_hydrogens to public
chrisjonesBSU Sep 22, 2023
4caec72
add error message when calling apply_forcefield
chrisjonesBSU Sep 22, 2023
b4aee94
fix r_cut param
chrisjonesBSU Sep 22, 2023
666ec7d
fix tests
chrisjonesBSU Sep 22, 2023
16d1676
add check for ref values back into self.hoomd_forcefield
chrisjonesBSU Sep 23, 2023
63a9245
fix unit test
chrisjonesBSU Sep 23, 2023
b93a70f
fix kwarg
chrisjonesBSU Sep 25, 2023
b512cca
fix unit test
chrisjonesBSU Sep 25, 2023
637efa1
add a couple more unit tests
chrisjonesBSU Sep 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 56 additions & 42 deletions hoomd_organics/base/system.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,34 +39,27 @@
self,
molecules,
density: float,
r_cut: float,
force_field=None,
auto_scale=False,
remove_hydrogens=False,
remove_charges=False,
scale_charges=False,
base_units=dict(),
):
self._molecules = check_return_iterable(molecules)
self._force_field = None
if force_field:
self._force_field = check_return_iterable(force_field)
self.density = density
self.r_cut = r_cut
self.auto_scale = auto_scale
self.remove_hydrogens = remove_hydrogens
self.remove_charges = remove_charges
self.scale_charges = scale_charges
self._target_box = None
self.all_molecules = []
self.gmso_system = None
self._reference_values = base_units
self._hoomd_snapshot = None
self._hoomd_forcefield = []
self._reference_values = base_units
self._gmso_forcefields_dict = dict()
self._target_box = None
# Reference values used when last writing snapshot and forcefields
self._ff_refs = dict()
self._snap_refs = dict()
self.gmso_system = None
self._ff_kwargs = dict()

# Collecting all molecules
self.n_mol_types = 0
Expand Down Expand Up @@ -120,16 +113,10 @@
f"{self._force_field[ff_index]} is not "
f"provided."
)
# Create mBuild system
self.system = self._build_system()
# Create GMSO topology
self.gmso_system = self._convert_to_gmso()
if self._force_field:
self._apply_forcefield()
if self.remove_hydrogens:
self._remove_hydrogens()
self._hoomd_forcefield = (
self._create_hoomd_forcefield() if self._force_field else []
)
self._hoomd_snapshot = self._create_hoomd_snapshot()

@abstractmethod
def _build_system(self):
Expand Down Expand Up @@ -160,6 +147,7 @@

@property
def box(self):
# TODO: Use gmso system here?
return self.system.box

@property
Expand Down Expand Up @@ -205,12 +193,16 @@
def hoomd_snapshot(self):
if self._snap_refs != self.reference_values:
self._hoomd_snapshot = self._create_hoomd_snapshot()
if self._hoomd_snapshot is None: # Hasn't been created yet
self._hoomd_snapshot = self._create_hoomd_snapshot()

Check warning on line 197 in hoomd_organics/base/system.py

View check run for this annotation

Codecov / codecov/patch

hoomd_organics/base/system.py#L197 

Added line #L197 was not covered by tests
return self._hoomd_snapshot

@property
def hoomd_forcefield(self):
if self._ff_refs != self.reference_values and self._force_field:
self._hoomd_forcefield = self._create_hoomd_forcefield()
self._hoomd_forcefield = self._create_hoomd_forcefield(
**self._ff_kwargs
)
return self._hoomd_forcefield

@property
Expand All @@ -220,7 +212,7 @@
else:
return self._target_box

def _remove_hydrogens(self):
def remove_hydrogens(self):
"""Call this method to remove hydrogen atoms from the system.
The masses and charges of the hydrogens are absorbed into
the heavy atoms they were bonded to.
Expand Down Expand Up @@ -252,6 +244,9 @@
site.mass += h.mass
site.charge += h.charge
self.gmso_system.remove_site(site=h)
# If a snap shot was already made, need to re-create it w/o hydrogens
if self._hoomd_snapshot:
self._create_hoomd_snapshot()

def _scale_charges(self):
""""""
Expand All @@ -278,11 +273,13 @@
topology.identify_connections()
return topology

def _create_hoomd_forcefield(self):
def _create_hoomd_forcefield(self, r_cut, nlist_buffer, pppm_kwargs):
force_list = []
ff, refs = to_hoomd_forcefield(
top=self.gmso_system,
r_cut=self.r_cut,
r_cut=r_cut,
nlist_buffer=nlist_buffer,
pppm_kwargs=pppm_kwargs,
auto_scale=self.auto_scale,
base_units=self._reference_values
if self._reference_values
Expand All @@ -304,19 +301,38 @@
self._snap_refs = self._reference_values.copy()
return snap

def _apply_forcefield(self):
def apply_forcefield(
self,
r_cut,
auto_scale=False,
scale_charges=False,
remove_charges=False,
remove_hydrogens=False,
pppm_resolution=(8, 8, 8),
pppm_order=4,
nlist_buffer=0.4,
):
if not self._force_field:
# TODO: Better erorr message
raise ValueError(
"This method can only be used when the System is "
"initialized with an XML type forcefield."
)
self.gmso_system = apply(
self.gmso_system,
self._gmso_forcefields_dict,
identify_connections=True,
speedup_by_moltag=True,
speedup_by_molgraph=False,
)
if self.remove_charges:

if remove_charges:
for site in self.gmso_system.sites:
site.charge = 0
if self.scale_charges and not self.remove_charges:

if scale_charges and not remove_charges:
self._scale_charges()

epsilons = [
s.atom_type.parameters["epsilon"] for s in self.gmso_system.sites
]
Expand All @@ -333,6 +349,20 @@
self._reference_values["length"] = length_scale * sigmas[0].unit_array
self._reference_values["mass"] = mass_scale * masses[0].unit_array

if remove_hydrogens:
self.remove_hydrogens()

pppm_kwargs = {"resolution": pppm_resolution, "order": pppm_order}
self._ff_kwargs = {
"r_cut": r_cut,
"nlist_buffer": nlist_buffer,
"pppm_kwargs": pppm_kwargs,
}
self._hoomd_forcefield = self._create_hoomd_forcefield(
r_cut=r_cut, nlist_buffer=nlist_buffer, pppm_kwargs=pppm_kwargs
)
self._hoomd_snapshot = self._create_hoomd_snapshot()

def set_target_box(
self, x_constraint=None, y_constraint=None, z_constraint=None
):
Expand Down Expand Up @@ -429,12 +459,8 @@
self,
molecules,
density: float,
r_cut: float,
force_field=None,
auto_scale=False,
remove_hydrogens=False,
remove_charges=False,
scale_charges=False,
base_units=dict(),
packing_expand_factor=5,
edge=0.2,
Expand All @@ -445,11 +471,7 @@
molecules=molecules,
density=density,
force_field=force_field,
r_cut=r_cut,
auto_scale=auto_scale,
remove_hydrogens=remove_hydrogens,
remove_charges=remove_charges,
scale_charges=scale_charges,
base_units=base_units,
)

Expand Down Expand Up @@ -485,16 +507,12 @@
self,
molecules,
density: float,
r_cut: float,
x: float,
y: float,
n: int,
basis_vector=[0.5, 0.5, 0],
force_field=None,
auto_scale=False,
remove_hydrogens=False,
remove_charges=False,
scale_charges=False,
base_units=dict(),
):
self.x = x
Expand All @@ -505,11 +523,7 @@
molecules=molecules,
density=density,
force_field=force_field,
r_cut=r_cut,
auto_scale=auto_scale,
remove_hydrogens=remove_hydrogens,
remove_charges=remove_charges,
scale_charges=scale_charges,
base_units=base_units,
)

Expand Down
4 changes: 2 additions & 2 deletions hoomd_organics/tests/base/test_simulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,11 +151,11 @@ def test_update_volume_walls(self, benzene_system):
def test_update_volume_density(self, benzene_system):
sim = Simulation.from_system(benzene_system)
sim.run_update_volume(
kT=1.0, tau_kt=0.01, n_steps=500, period=1, final_density=0.1
kT=1.0, tau_kt=0.01, n_steps=500, period=1, final_density=0.05
)
assert np.isclose(
sim.density.to(u.g / u.cm**3).value,
(0.1 * (u.g / u.cm**3)).value,
(0.05 * (u.g / u.cm**3)).value,
atol=1e-4,
)

Expand Down
Loading
Loading