From f2fa27b1d591e59a450e1fa518b4e33951211bb7 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Wed, 8 Nov 2023 20:42:09 -0700 Subject: [PATCH 01/37] create EllipsoidChain class --- flowermd/library/__init__.py | 1 + flowermd/library/polymers.py | 57 ++++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+) diff --git a/flowermd/library/__init__.py b/flowermd/library/__init__.py index 3a514de4..4ccfa38b 100644 --- a/flowermd/library/__init__.py +++ b/flowermd/library/__init__.py @@ -16,6 +16,7 @@ PEKK, PPS, LJChain, + EllipsoidChain, PEKK_meta, PEKK_para, PolyEthylene, diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index 746eded0..f35cf9d7 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -3,6 +3,7 @@ import mbuild as mb from mbuild.coordinate_transform import z_axis_transform +from mbuild.lib.recipes import Polymer as mbPolymer from flowermd import CoPolymer, Polymer from flowermd.assets import MON_DIR @@ -261,3 +262,59 @@ def _build(self, length): chain.add_bond([next_bead, last_bead]) last_bead = next_bead return chain + + +class EllipsoidChain(Polymer): + def __init__( + self, + lengths, + num_mols, + bead_length, + bead_name, + bead_mass, + bond_length + ): + self.bead_name = bead_name + self.bead_mass = bead_mass + self.bead_bond_length = bond_length + self.bead_length = bead_length + super(EllipsoidChain, self).__init__(lengths=lengths, num_mols=num_mols) + + def _build(self, length): + # Build up ellipsoid bead + bead = mb.Compound(name="ellipsoid") + head = mb.Compound( + pos=(self.bead_length/2, 0, 0), + name="A", + mass=self.bead_mass/2 + ) + tail = mb.Compound( + pos=(-self.bead_length/2, 0, 0), + name="A", + mass=self.bead_mass/2 + ) + head_mid = mb.Compound( + pos=(self.bead_length/4, 0, 0), + name="B", + mass=0 + ) + tail_mid = mb.Compound( + pos=(-self.bead_length/4, 0, 0), + name="B", + mass=0 + ) + bead.add([head, tail, head_mid, tail_mid]) + # Build the chain + chain = mbPolymer() + chain.add_monomer( + bead, + indices=[0, 1], + orientation=[[1,0,0], [-1,0,0]], + replace=False, + separation=self.bead_bond_length + ) + chain.build(n=length, add_hydrogens=False) + return chain + + + From 0f6238322d4bec1a32dc41f61ca357667a07cd83 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Wed, 8 Nov 2023 20:43:36 -0700 Subject: [PATCH 02/37] remove blank space --- flowermd/library/polymers.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index f35cf9d7..cd171e2e 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -304,7 +304,7 @@ def _build(self, length): mass=0 ) bead.add([head, tail, head_mid, tail_mid]) - # Build the chain + chain = mbPolymer() chain.add_monomer( bead, @@ -315,6 +315,3 @@ def _build(self, length): ) chain.build(n=length, add_hydrogens=False) return chain - - - From f490cfd095c0e100abfd2c5d1d80ec7560452248 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 30 Nov 2023 13:50:20 -0700 Subject: [PATCH 03/37] add class for ellipsoid chain FF --- flowermd/library/forcefields.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/flowermd/library/forcefields.py b/flowermd/library/forcefields.py index 0c119856..b073f4fc 100644 --- a/flowermd/library/forcefields.py +++ b/flowermd/library/forcefields.py @@ -481,3 +481,25 @@ def _check_widths(self): "number of points for table energies and forces." ) return bond_width, angle_width, dih_width + + +class EllipsoidForcefield(BaseHOOMDForcefield): + def __init__(self): + hoomd_forces = self._create_forcefield() + super(TableForcefield, self).__init__(hoomd_forces) + + def _create_forcefield(self): + forces = [] + + # Bonds + + # Angles + + # Gay-Berne Pairs + nlist = hoomd.md.nlist.Cell(buffer=0.40) + + + + + + return forces From 97aff78d939017cfaeac7ae42e8d8db6cb0d6edc Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Fri, 1 Dec 2023 11:15:41 -0700 Subject: [PATCH 04/37] add FF class to library --- flowermd/library/forcefields.py | 56 +++++++++++++++++++++++++++------ 1 file changed, 46 insertions(+), 10 deletions(-) diff --git a/flowermd/library/forcefields.py b/flowermd/library/forcefields.py index b073f4fc..4d45954b 100644 --- a/flowermd/library/forcefields.py +++ b/flowermd/library/forcefields.py @@ -484,22 +484,58 @@ def _check_widths(self): class EllipsoidForcefield(BaseHOOMDForcefield): - def __init__(self): - hoomd_forces = self._create_forcefield() + """A forcefield for modeling anisotropic bead polymers.""" + + def __init__( + self, + epsilon, + lper, + lpar, + bead_length, + r_cut, + bond_k, + bond_r0, + angle_k=None, + angle_theta0=None, + ): + self.epsilon = epsilon + self.lper = lper + self.lpar = lpar + self.bead_length = bead_length + self.r_cut = r_cut + self.bond_k = bond_k + self.bond_r0 = bond_r0 + self.angle_k = angle_k + self.angle_theta0 = angle_theta0 + hoomd_forces = self._create_forcefield() super(TableForcefield, self).__init__(hoomd_forces) def _create_forcefield(self): forces = [] - # Bonds - + bond = hoomd.md.bond.Harmonic() + bond.params["A-A"] = dict(k=self.bond_k, r0=self.bond_r0) + bond.params["B-B"] = dict(k=0, r0=self.bead_length / 2) + forces.append(bond) # Angles - + if all([self.angle_k, self.angle_theta0]): + angle = hoomd.md.angle.Harmonic() + angle.params["B-B-B"] = dict(k=self.angle_k, t0=self.angle_theta0) + forces.append(angle) # Gay-Berne Pairs nlist = hoomd.md.nlist.Cell(buffer=0.40) - - - - - + gb = hoomd.md.pair.aniso.GayBerne(nlist=nlist, default_r_cut=self.r_cut) + gb.params[("R", "R")] = dict( + epsilon=self.epsilon, lperp=self.lperp, lpar=self.lpar + ) + # Add zero pairs + for pair in [ + ("A", "A"), + ("B", "B"), + ("A", "B"), + ("A", "R"), + ("B", "R"), + ]: + gb.params[pair] = dict(epsilon=0.0, lperp=0.0, lpar=0.0) + forces.append(gb) return forces From c6b8326ac2e08036e727dd9d141a224fd2c9c7a6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 1 Dec 2023 20:16:11 +0000 Subject: [PATCH 05/37] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- flowermd/library/__init__.py | 2 +- flowermd/library/polymers.py | 34 ++++++++++------------------------ 2 files changed, 11 insertions(+), 25 deletions(-) diff --git a/flowermd/library/__init__.py b/flowermd/library/__init__.py index 4ccfa38b..0dfcc145 100644 --- a/flowermd/library/__init__.py +++ b/flowermd/library/__init__.py @@ -15,8 +15,8 @@ PEEK, PEKK, PPS, - LJChain, EllipsoidChain, + LJChain, PEKK_meta, PEKK_para, PolyEthylene, diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index cd171e2e..2b24e5d8 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -266,13 +266,7 @@ def _build(self, length): class EllipsoidChain(Polymer): def __init__( - self, - lengths, - num_mols, - bead_length, - bead_name, - bead_mass, - bond_length + self, lengths, num_mols, bead_length, bead_name, bead_mass, bond_length ): self.bead_name = bead_name self.bead_mass = bead_mass @@ -284,34 +278,26 @@ def _build(self, length): # Build up ellipsoid bead bead = mb.Compound(name="ellipsoid") head = mb.Compound( - pos=(self.bead_length/2, 0, 0), - name="A", - mass=self.bead_mass/2 + pos=(self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 2 ) tail = mb.Compound( - pos=(-self.bead_length/2, 0, 0), - name="A", - mass=self.bead_mass/2 + pos=(-self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 2 ) head_mid = mb.Compound( - pos=(self.bead_length/4, 0, 0), - name="B", - mass=0 + pos=(self.bead_length / 4, 0, 0), name="B", mass=0 ) tail_mid = mb.Compound( - pos=(-self.bead_length/4, 0, 0), - name="B", - mass=0 + pos=(-self.bead_length / 4, 0, 0), name="B", mass=0 ) bead.add([head, tail, head_mid, tail_mid]) chain = mbPolymer() chain.add_monomer( - bead, - indices=[0, 1], - orientation=[[1,0,0], [-1,0,0]], - replace=False, - separation=self.bead_bond_length + bead, + indices=[0, 1], + orientation=[[1, 0, 0], [-1, 0, 0]], + replace=False, + separation=self.bead_bond_length, ) chain.build(n=length, add_hydrogens=False) return chain From 4af279527c64047b27bf6236123e9d2332f2a050 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 7 Dec 2023 12:08:24 -0700 Subject: [PATCH 06/37] update particle mass, add bead particle types sequence --- flowermd/library/polymers.py | 10 ++++++---- flowermd/utils/rigid_body.py | 0 2 files changed, 6 insertions(+), 4 deletions(-) create mode 100644 flowermd/utils/rigid_body.py diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index 2b24e5d8..4383f8be 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -273,21 +273,23 @@ def __init__( self.bead_bond_length = bond_length self.bead_length = bead_length super(EllipsoidChain, self).__init__(lengths=lengths, num_mols=num_mols) + # get the indices of the particles in a rigid body + self.bead_constituents_types = ["A", "A", "B", "B"] def _build(self, length): # Build up ellipsoid bead bead = mb.Compound(name="ellipsoid") head = mb.Compound( - pos=(self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 2 + pos=(self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 4 ) tail = mb.Compound( - pos=(-self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 2 + pos=(-self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 4 ) head_mid = mb.Compound( - pos=(self.bead_length / 4, 0, 0), name="B", mass=0 + pos=(self.bead_length / 4, 0, 0), name="B", mass=self.bead_mass / 4 ) tail_mid = mb.Compound( - pos=(-self.bead_length / 4, 0, 0), name="B", mass=0 + pos=(-self.bead_length / 4, 0, 0), name="B", mass=self.bead_mass / 4 ) bead.add([head, tail, head_mid, tail_mid]) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py new file mode 100644 index 00000000..e69de29b From 2460a7fe4967afa6b98ce6d4841aa863f1d45b30 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 7 Dec 2023 12:08:59 -0700 Subject: [PATCH 07/37] update class parent name --- flowermd/library/forcefields.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/flowermd/library/forcefields.py b/flowermd/library/forcefields.py index 4d45954b..aaa33c04 100644 --- a/flowermd/library/forcefields.py +++ b/flowermd/library/forcefields.py @@ -489,7 +489,7 @@ class EllipsoidForcefield(BaseHOOMDForcefield): def __init__( self, epsilon, - lper, + lperp, lpar, bead_length, r_cut, @@ -499,7 +499,7 @@ def __init__( angle_theta0=None, ): self.epsilon = epsilon - self.lper = lper + self.lperp = lperp self.lpar = lpar self.bead_length = bead_length self.r_cut = r_cut @@ -508,7 +508,7 @@ def __init__( self.angle_k = angle_k self.angle_theta0 = angle_theta0 hoomd_forces = self._create_forcefield() - super(TableForcefield, self).__init__(hoomd_forces) + super(EllipsoidForcefield, self).__init__(hoomd_forces) def _create_forcefield(self): forces = [] From fb0708f1ae1e38a5cb4f0f56fe20e7cc7495cc19 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 7 Dec 2023 12:11:02 -0700 Subject: [PATCH 08/37] add func to create rigid frame --- flowermd/utils/rigid_body.py | 85 ++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index e69de29b..924e6240 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -0,0 +1,85 @@ +import gsd +import gsd.hoomd +import numpy as np +from cmeutils.geometry import moit +import hoomd + + +def _get_com_mass_pos_moi(snapshot, rigid_const_idx): + com_mass = [] + com_position = [] + com_moi = [] + for idx in rigid_const_idx: + total_mass = np.sum(snapshot.particles.mass[idx]) + com_mass.append(total_mass) + com = np.sum( + snapshot.particles.position[idx] + * snapshot.particles.mass[idx, np.newaxis], + axis=0, + ) / total_mass + com_position.append(com) + com_moi.append( + moit(points=snapshot.particles.position[idx], + masses=snapshot.particles.mass[idx], + center=com)) + return com_mass, com_position, com_moi + + +def create_rigid_body(snapshot, bead_constituents_types): + """Create rigid bodies from a snapshot. + + Parameters + ---------- + snapshot : gsd.hoomd.Snapshot; required + The snapshot of the system. + bead_constituents_types : list of particle types; required + The list of particle types that make up a rigid body. + + Returns + ------- + rigid_frame : gsd.hoomd.Frame + The snapshot of the rigid bodies. + rigid_constrain : hoomd.md.constrain.Rigid + The rigid body constrain object. + """ + # find typeid sequence of the constituent particles types in a rigid bead + p_types = np.array(snapshot.particles.types) + constituent_type_ids = \ + np.where(p_types[:, None] == bead_constituents_types)[0] + + # find indices that matches the constituent particle types + typeids = snapshot.particles.typeid + bead_len = len(bead_constituents_types) + matches = np.where((typeids.reshape(-1, bead_len) == constituent_type_ids)) + if len(matches[0]) == 0: + raise ValueError("No matches found between particle types in the system" + " and bead constituents particles") + rigid_const_idx = (matches[0] * bead_len + matches[1]).reshape(-1, bead_len) + + n_rigid = rigid_const_idx.shape[0] # number of rigid bodies + + # calculate center of mass and its position for each rigid body + com_mass, com_position, com_moi = _get_com_mass_pos_moi(snapshot, + rigid_const_idx) + + rigid_frame = gsd.hoomd.Frame() + rigid_frame.particles.types = ['R'] + snapshot.particles.types + rigid_frame.particles.N = n_rigid + rigid_frame.particles.typeid = [0] * n_rigid + rigid_frame.particles.mass = com_mass + rigid_frame.particles.position = com_position + rigid_frame.particles.moment_inertia = com_moi + rigid_frame.particles.orientation = [(1.0, 0.0, 0.0, 0.0)] * n_rigid + + # find local coordinates of the particles in the first rigid body + # only need to find the local coordinates for the first rigid body + local_coords = snapshot.particles.position[rigid_const_idx[0]] - \ + com_position[0] + + rigid_constrain = hoomd.md.constrain.Rigid() + rigid_constrain.body['rigid'] = { + "constituent_types": bead_constituents_types, + "positions": local_coords, + "orientations": [(1.0, 0.0, 0.0, 0.0)] * len(local_coords), + } + return rigid_frame, rigid_constrain From dca5a86ca381a520e8ab08485e36fe83a40f82f7 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 7 Dec 2023 12:14:09 -0700 Subject: [PATCH 09/37] fix formatting, remove extra param --- flowermd/library/polymers.py | 22 +++++++++++++--- flowermd/utils/rigid_body.py | 49 ++++++++++++++++++++++-------------- 2 files changed, 48 insertions(+), 23 deletions(-) diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index 4383f8be..42d02ef1 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -265,10 +265,24 @@ def _build(self, length): class EllipsoidChain(Polymer): - def __init__( - self, lengths, num_mols, bead_length, bead_name, bead_mass, bond_length - ): - self.bead_name = bead_name + """Create an ellipsoid polymer chain. + + Parameters + ---------- + lengths : int, required + The number of monomer repeat units in the chain. + num_mols : int, required + The number of chains to create. + bead_length : float, required + The length of the ellipsoid bead. + bead_mass : float, required + The mass of the ellipsoid bead. + bond_length : float, required + The bond length between connected beads. + + """ + + def __init__(self, lengths, num_mols, bead_length, bead_mass, bond_length): self.bead_mass = bead_mass self.bead_bond_length = bond_length self.bead_length = bead_length diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index 924e6240..8e3ada2f 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -1,8 +1,8 @@ import gsd import gsd.hoomd +import hoomd import numpy as np from cmeutils.geometry import moit -import hoomd def _get_com_mass_pos_moi(snapshot, rigid_const_idx): @@ -12,16 +12,22 @@ def _get_com_mass_pos_moi(snapshot, rigid_const_idx): for idx in rigid_const_idx: total_mass = np.sum(snapshot.particles.mass[idx]) com_mass.append(total_mass) - com = np.sum( - snapshot.particles.position[idx] - * snapshot.particles.mass[idx, np.newaxis], - axis=0, - ) / total_mass + com = ( + np.sum( + snapshot.particles.position[idx] + * snapshot.particles.mass[idx, np.newaxis], + axis=0, + ) + / total_mass + ) com_position.append(com) com_moi.append( - moit(points=snapshot.particles.position[idx], - masses=snapshot.particles.mass[idx], - center=com)) + moit( + points=snapshot.particles.position[idx], + masses=snapshot.particles.mass[idx], + center=com, + ) + ) return com_mass, com_position, com_moi @@ -44,26 +50,30 @@ def create_rigid_body(snapshot, bead_constituents_types): """ # find typeid sequence of the constituent particles types in a rigid bead p_types = np.array(snapshot.particles.types) - constituent_type_ids = \ - np.where(p_types[:, None] == bead_constituents_types)[0] + constituent_type_ids = np.where( + p_types[:, None] == bead_constituents_types + )[0] # find indices that matches the constituent particle types typeids = snapshot.particles.typeid bead_len = len(bead_constituents_types) matches = np.where((typeids.reshape(-1, bead_len) == constituent_type_ids)) if len(matches[0]) == 0: - raise ValueError("No matches found between particle types in the system" - " and bead constituents particles") + raise ValueError( + "No matches found between particle types in the system" + " and bead constituents particles" + ) rigid_const_idx = (matches[0] * bead_len + matches[1]).reshape(-1, bead_len) n_rigid = rigid_const_idx.shape[0] # number of rigid bodies # calculate center of mass and its position for each rigid body - com_mass, com_position, com_moi = _get_com_mass_pos_moi(snapshot, - rigid_const_idx) + com_mass, com_position, com_moi = _get_com_mass_pos_moi( + snapshot, rigid_const_idx + ) rigid_frame = gsd.hoomd.Frame() - rigid_frame.particles.types = ['R'] + snapshot.particles.types + rigid_frame.particles.types = ["R"] + snapshot.particles.types rigid_frame.particles.N = n_rigid rigid_frame.particles.typeid = [0] * n_rigid rigid_frame.particles.mass = com_mass @@ -73,11 +83,12 @@ def create_rigid_body(snapshot, bead_constituents_types): # find local coordinates of the particles in the first rigid body # only need to find the local coordinates for the first rigid body - local_coords = snapshot.particles.position[rigid_const_idx[0]] - \ - com_position[0] + local_coords = ( + snapshot.particles.position[rigid_const_idx[0]] - com_position[0] + ) rigid_constrain = hoomd.md.constrain.Rigid() - rigid_constrain.body['rigid'] = { + rigid_constrain.body["rigid"] = { "constituent_types": bead_constituents_types, "positions": local_coords, "orientations": [(1.0, 0.0, 0.0, 0.0)] * len(local_coords), From e8247bcbcafe76a91a30cf1dfbfa0a46b15d8c78 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 7 Dec 2023 12:35:03 -0700 Subject: [PATCH 10/37] fix numpy bug, add box to config --- flowermd/utils/rigid_body.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index 8e3ada2f..de5ffb7e 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -10,12 +10,13 @@ def _get_com_mass_pos_moi(snapshot, rigid_const_idx): com_position = [] com_moi = [] for idx in rigid_const_idx: - total_mass = np.sum(snapshot.particles.mass[idx]) + constituents_mass = np.array(snapshot.particles.mass) + constituents_pos = np.array(snapshot.particles.position) + total_mass = np.sum(constituents_mass[idx]) com_mass.append(total_mass) com = ( np.sum( - snapshot.particles.position[idx] - * snapshot.particles.mass[idx, np.newaxis], + constituents_pos[idx] * constituents_mass[idx, np.newaxis], axis=0, ) / total_mass @@ -23,8 +24,8 @@ def _get_com_mass_pos_moi(snapshot, rigid_const_idx): com_position.append(com) com_moi.append( moit( - points=snapshot.particles.position[idx], - masses=snapshot.particles.mass[idx], + points=constituents_pos[idx], + masses=constituents_mass[idx], center=com, ) ) @@ -80,6 +81,7 @@ def create_rigid_body(snapshot, bead_constituents_types): rigid_frame.particles.position = com_position rigid_frame.particles.moment_inertia = com_moi rigid_frame.particles.orientation = [(1.0, 0.0, 0.0, 0.0)] * n_rigid + rigid_frame.configuration.box = snapshot.configuration.box # find local coordinates of the particles in the first rigid body # only need to find the local coordinates for the first rigid body From 50819ee7b608b6dee29673077c82bb9f8cf96074 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Thu, 7 Dec 2023 13:32:47 -0700 Subject: [PATCH 11/37] fix rigid body sim bugs --- flowermd/base/simulation.py | 30 ++++++++++++++++++++++++++++-- flowermd/utils/rigid_body.py | 2 +- 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 5ac38067..0bc50a00 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -73,6 +73,7 @@ def __init__( log_write_freq=1e3, log_file_name="sim_data.txt", thermostat=HOOMDThermostats.MTTK, + rigid_constraint=None, ): super(Simulation, self).__init__(device, seed) self.initial_state = initial_state @@ -98,9 +99,22 @@ def __init__( self._dt = dt self._reference_values = dict() self._reference_values = reference_values - self._integrate_group = hoomd.filter.All() + if rigid_constraint and not isinstance( + rigid_constraint, hoomd.md.constrain.Rigid + ): + raise ValueError( + "Invalid rigid constraint. Please provide a " + "hoomd.md.constrain.Rigid object." + ) + self._rigid_constraint = rigid_constraint + self._integrate_group = self._create_integrate_group( + rigid=True if rigid_constraint else False + ) self._wall_forces = dict() self._create_state(self.initial_state) + if rigid_constraint: + # place constituent particles in the simulation state + rigid_constraint.create_bodies(self.state) # Add a gsd and thermo props logger to sim operations self._add_hoomd_writers() self._thermostat = thermostat @@ -524,7 +538,14 @@ def set_integrator_method(self, integrator_method, method_kwargs): """ if not self.integrator: # Integrator and method not yet created - self.integrator = hoomd.md.Integrator(dt=self.dt) + self.integrator = hoomd.md.Integrator( + dt=self.dt, + integrate_rotational_dof=True + if self._rigid_constraint + else False, + ) + if self._rigid_constraint: + self.integrator.rigid = self._rigid_constraint self.integrator.forces = self._forcefield self.operations.add(self.integrator) new_method = integrator_method(**method_kwargs) @@ -1088,6 +1109,11 @@ def _lj_force(self): ][0] return lj_force + def _create_integrate_group(self, rigid): + if rigid: + return hoomd.filter.Rigid(("center", "free")) + return hoomd.filter.All() + def _create_state(self, initial_state): """Create the simulation state. diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index de5ffb7e..b634ed99 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -90,7 +90,7 @@ def create_rigid_body(snapshot, bead_constituents_types): ) rigid_constrain = hoomd.md.constrain.Rigid() - rigid_constrain.body["rigid"] = { + rigid_constrain.body["R"] = { "constituent_types": bead_constituents_types, "positions": local_coords, "orientations": [(1.0, 0.0, 0.0, 0.0)] * len(local_coords), From 52cdea3ba3fb6081547f074440f9598d67d444f2 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Fri, 8 Dec 2023 09:57:42 -0700 Subject: [PATCH 12/37] add bonds between B particles after building the chain --- flowermd/library/polymers.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index 42d02ef1..3883c879 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -316,4 +316,12 @@ def _build(self, length): separation=self.bead_bond_length, ) chain.build(n=length, add_hydrogens=False) + # Generate bonds between the mid-particles. + # This is needed to use an angle potential between 2 beads. + chain.freud_generate_bonds( + name_a="B", + name_b="B", + dmin=self.bead_length / 2 - 0.1, + dmax=self.bead_length / 2 + self.bead_bond_length + 0.1, + ) return chain From cb70367f9ba5808a4255b074df39db82cf310016 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 10:03:54 -0700 Subject: [PATCH 13/37] add fix_orientation variable to pack --- flowermd/base/system.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/flowermd/base/system.py b/flowermd/base/system.py index cbe1befd..abd2f9ad 100644 --- a/flowermd/base/system.py +++ b/flowermd/base/system.py @@ -729,9 +729,11 @@ def __init__( base_units=dict(), packing_expand_factor=5, edge=0.2, + fix_orientation=False, ): self.packing_expand_factor = packing_expand_factor self.edge = edge + self.fix_orientation = fix_orientation super(Pack, self).__init__( molecules=molecules, density=density, @@ -746,6 +748,7 @@ def _build_system(self): box=list(self._target_box * self.packing_expand_factor), overlap=0.2, edge=self.edge, + fix_orientation=self.fix_orientation, ) return system From e2aab60d0c5210a7e9fbb16ffa5b2120b3d23fa5 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 10:05:41 -0700 Subject: [PATCH 14/37] add bond, angle and dihedral --- flowermd/utils/rigid_body.py | 52 ++++++++++++++++++++++++++++++------ 1 file changed, 44 insertions(+), 8 deletions(-) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index b634ed99..eea192a5 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -32,7 +32,7 @@ def _get_com_mass_pos_moi(snapshot, rigid_const_idx): return com_mass, com_position, com_moi -def create_rigid_body(snapshot, bead_constituents_types): +def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): """Create rigid bodies from a snapshot. Parameters @@ -74,15 +74,51 @@ def create_rigid_body(snapshot, bead_constituents_types): ) rigid_frame = gsd.hoomd.Frame() - rigid_frame.particles.types = ["R"] + snapshot.particles.types - rigid_frame.particles.N = n_rigid - rigid_frame.particles.typeid = [0] * n_rigid - rigid_frame.particles.mass = com_mass - rigid_frame.particles.position = com_position - rigid_frame.particles.moment_inertia = com_moi - rigid_frame.particles.orientation = [(1.0, 0.0, 0.0, 0.0)] * n_rigid + rigid_frame.particles.types = [bead_name] + snapshot.particles.types + rigid_frame.particles.N = n_rigid + snapshot.particles.N + rigid_frame.particles.typeid = np.concatenate( + (([0] * n_rigid), snapshot.particles.typeid + 1) + ) + rigid_frame.particles.mass = np.concatenate( + (com_mass, snapshot.particles.mass) + ) + rigid_frame.particles.position = np.concatenate( + (com_position, snapshot.particles.position) + ) + rigid_frame.particles.moment_inertia = np.concatenate( + (com_moi, np.zeros((snapshot.particles.N, 3))) + ) + rigid_frame.particles.orientation = [(1.0, 0.0, 0.0, 0.0)] * ( + n_rigid + snapshot.particles.N + ) rigid_frame.configuration.box = snapshot.configuration.box + # set up bonds + if snapshot.bonds.N > 0: + rigid_frame.bonds.N = snapshot.bonds.N + rigid_frame.bonds.types = snapshot.bonds.types + rigid_frame.bonds.typeid = snapshot.bonds.typeid + rigid_frame.bonds.group = [ + list(np.add(g, n_rigid)) for g in snapshot.bonds.group + ] + # set up angles + if snapshot.angles.N > 0: + rigid_frame.angles.N = not snapshot.angles.N + rigid_frame.angles.types = snapshot.angles.types + rigid_frame.angles.typeid = snapshot.angles.typeid + rigid_frame.angles.group = [ + list(np.add(g, n_rigid)) for g in snapshot.angles.group + ] + + # set up dihedrals + if snapshot.dihedrals.N > 0: + rigid_frame.dihedrals.N = snapshot.dihedrals.N + rigid_frame.dihedrals.types = snapshot.dihedrals.types + rigid_frame.dihedrals.typeid = snapshot.dihedrals.typeid + rigid_frame.dihedrals.group = [ + list(np.add(g, n_rigid)) for g in snapshot.dihedrals.group + ] + # find local coordinates of the particles in the first rigid body # only need to find the local coordinates for the first rigid body local_coords = ( From fe25155f7235a789b2d851c42d569e28f51de61e Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 10:16:35 -0700 Subject: [PATCH 15/37] remove create_rigid_body from sim --- flowermd/base/simulation.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 0bc50a00..0ea87608 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -112,9 +112,6 @@ def __init__( ) self._wall_forces = dict() self._create_state(self.initial_state) - if rigid_constraint: - # place constituent particles in the simulation state - rigid_constraint.create_bodies(self.state) # Add a gsd and thermo props logger to sim operations self._add_hoomd_writers() self._thermostat = thermostat From 0c09d8260f5a255b686611b1b5f61cc469ea0917 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 10:46:28 -0700 Subject: [PATCH 16/37] add body tags to rigid frame --- flowermd/utils/rigid_body.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index eea192a5..9a8f7c61 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -91,6 +91,12 @@ def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): rigid_frame.particles.orientation = [(1.0, 0.0, 0.0, 0.0)] * ( n_rigid + snapshot.particles.N ) + rigid_frame.particles.body = np.concatenate( + ( + np.arange(n_rigid), + np.arange(n_rigid).repeat(rigid_const_idx.shape[1]), + ) + ) rigid_frame.configuration.box = snapshot.configuration.box # set up bonds @@ -103,7 +109,7 @@ def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): ] # set up angles if snapshot.angles.N > 0: - rigid_frame.angles.N = not snapshot.angles.N + rigid_frame.angles.N = snapshot.angles.N rigid_frame.angles.types = snapshot.angles.types rigid_frame.angles.typeid = snapshot.angles.typeid rigid_frame.angles.group = [ From 963b8cf625a42f41ec12ebd787b9f3f9a4bfac6c Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 14:43:10 -0700 Subject: [PATCH 17/37] add new functions to __init__ --- flowermd/library/__init__.py | 1 + flowermd/tests/library/test_forcefields.py | 24 +++++++ flowermd/tests/library/test_polymers.py | 14 +++++ flowermd/tests/utils/test_rigid_body.py | 73 ++++++++++++++++++++++ flowermd/utils/__init__.py | 1 + 5 files changed, 113 insertions(+) create mode 100644 flowermd/tests/utils/test_rigid_body.py diff --git a/flowermd/library/__init__.py b/flowermd/library/__init__.py index 0dfcc145..399748eb 100644 --- a/flowermd/library/__init__.py +++ b/flowermd/library/__init__.py @@ -8,6 +8,7 @@ BaseHOOMDForcefield, BaseXMLForcefield, BeadSpring, + EllipsoidForcefield, FF_from_file, TableForcefield, ) diff --git a/flowermd/tests/library/test_forcefields.py b/flowermd/tests/library/test_forcefields.py index b554a906..788b0e2c 100644 --- a/flowermd/tests/library/test_forcefields.py +++ b/flowermd/tests/library/test_forcefields.py @@ -11,6 +11,7 @@ OPLS_AA_DIMETHYLETHER, OPLS_AA_PPS, BeadSpring, + EllipsoidForcefield, FF_from_file, TableForcefield, ) @@ -93,6 +94,29 @@ def test_BeadSpring(self): assert ff.hoomd_forces[3].params[param]["d"] == -1 assert ff.hoomd_forces[3].params[param]["n"] == 1 + def test_ellipsoid_ff(self): + ellipsoid_ff = EllipsoidForcefield( + epsilon=1.0, + lperp=0.5, + lpar=1.0, + bead_length=1, + r_cut=3, + bond_k=500, + bond_r0=0.1, + angle_k=100, + angle_theta0=1.57, + ) + assert len(ellipsoid_ff.hoomd_forces) == 3 + assert isinstance( + ellipsoid_ff.hoomd_forces[-1], hoomd.md.pair.aniso.GayBerne + ) + assert ("R", "R") in list( + dict(ellipsoid_ff.hoomd_forces[-1].params).keys() + ) + assert ("B", "R") in list( + dict(ellipsoid_ff.hoomd_forces[-1].params).keys() + ) + class TestTableForcefield: def test_from_txt_file(self): diff --git a/flowermd/tests/library/test_polymers.py b/flowermd/tests/library/test_polymers.py index d8dc0353..d51d6502 100644 --- a/flowermd/tests/library/test_polymers.py +++ b/flowermd/tests/library/test_polymers.py @@ -5,6 +5,7 @@ PEEK, PEKK, PPS, + EllipsoidChain, LJChain, PEKK_meta, PEKK_para, @@ -105,3 +106,16 @@ def test_lj_chain_sequence_bad_mass(self): def test_copolymer(self): pass + + def test_ellipsoid_chain(self): + ellipsoid_chain = EllipsoidChain( + lengths=4, + num_mols=2, + bead_length=1, + bead_mass=100, + bond_length=0.01, + ) + assert ellipsoid_chain.n_particles == 32 + assert ellipsoid_chain.molecules[0].mass == 400 + assert ellipsoid_chain.molecules[0].n_particles == 16 + assert ellipsoid_chain.molecules[0].n_bonds == 10 diff --git a/flowermd/tests/utils/test_rigid_body.py b/flowermd/tests/utils/test_rigid_body.py new file mode 100644 index 00000000..9f83e24d --- /dev/null +++ b/flowermd/tests/utils/test_rigid_body.py @@ -0,0 +1,73 @@ +import numpy as np + +from flowermd.base import Pack +from flowermd.library.polymers import EllipsoidChain +from flowermd.tests import BaseTest +from flowermd.utils import create_rigid_body + + +class TestRigidBody(BaseTest): + def test_ellipsoid_create_rigid_body(self): + ellipsoid_chain = EllipsoidChain( + lengths=4, + num_mols=2, + bead_length=1, + bead_mass=100, + bond_length=0.01, + ) + system = Pack( + molecules=ellipsoid_chain, density=0.1, fix_orientation=True + ) + + rigid_frame, rigid = create_rigid_body( + system.hoomd_snapshot, + ellipsoid_chain.bead_constituents_types, + bead_name="R", + ) + assert rigid_frame.particles.N == 8 + ellipsoid_chain.n_particles + assert ( + rigid_frame.particles.types + == ["R"] + system.hoomd_snapshot.particles.types + ) + assert rigid_frame.particles.mass[0] == 100 + assert np.all( + np.isclose( + rigid_frame.particles.position[0], + np.mean(system.hoomd_snapshot.particles.position[:4], axis=0), + ) + ) + + points = ( + rigid_frame.particles.position[0] + - system.hoomd_snapshot.particles.position[:4] + ) + + x = points[:, 0] + y = points[:, 1] + z = points[:, 2] + I_xx = np.sum( + (y**2 + z**2) * system.hoomd_snapshot.particles.mass[:4] + ) + I_yy = np.sum( + (x**2 + z**2) * system.hoomd_snapshot.particles.mass[:4] + ) + I_zz = np.sum( + (x**2 + y**2) * system.hoomd_snapshot.particles.mass[:4] + ) + assert np.all( + np.isclose( + rigid_frame.particles.moment_inertia[0], + np.array((I_xx, I_yy, I_zz)), + ) + ) + + assert np.all(rigid_frame.particles.body[:8] == np.arange(8)) + assert np.all(rigid_frame.particles.body[8:12] == [0]) + + assert rigid_frame.bonds.N == 20 + assert rigid_frame.bonds.types == system.hoomd_snapshot.bonds.types + assert rigid_frame.angles.N == 12 + assert rigid_frame.angles.types == system.hoomd_snapshot.angles.types + + assert "R" in list(rigid.body.keys()) + assert rigid.body["R"]["constituent_types"] == ["A", "A", "B", "B"] diff --git a/flowermd/utils/__init__.py b/flowermd/utils/__init__.py index 4200ed64..a869b9e3 100644 --- a/flowermd/utils/__init__.py +++ b/flowermd/utils/__init__.py @@ -1,6 +1,7 @@ from .actions import * from .base_types import HOOMDThermostats from .ff_utils import xml_to_gmso_ff +from .rigid_body import create_rigid_body from .utils import ( calculate_box_length, check_return_iterable, From 42a8f91255d59695b221f5b0a49efd67aed8dd23 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 14:45:11 -0700 Subject: [PATCH 18/37] add unit tests --- flowermd/tests/library/test_forcefields.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/flowermd/tests/library/test_forcefields.py b/flowermd/tests/library/test_forcefields.py index 788b0e2c..eb8622f3 100644 --- a/flowermd/tests/library/test_forcefields.py +++ b/flowermd/tests/library/test_forcefields.py @@ -116,6 +116,9 @@ def test_ellipsoid_ff(self): assert ("B", "R") in list( dict(ellipsoid_ff.hoomd_forces[-1].params).keys() ) + assert ellipsoid_ff.hoomd_forces[-1].params["A", "A"]["epsilon"] == 0.0 + assert ellipsoid_ff.hoomd_forces[-1].params["A", "A"]["lperp"] == 0.0 + assert ellipsoid_ff.hoomd_forces[-1].params["A", "A"]["lpar"] == 0.0 class TestTableForcefield: From ca676610238fc5992546ca14bf6a603997668095 Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 14:49:09 -0700 Subject: [PATCH 19/37] move private func to bottom --- flowermd/utils/rigid_body.py | 54 ++++++++++++++++++------------------ 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index 9a8f7c61..5ed7065c 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -5,33 +5,6 @@ from cmeutils.geometry import moit -def _get_com_mass_pos_moi(snapshot, rigid_const_idx): - com_mass = [] - com_position = [] - com_moi = [] - for idx in rigid_const_idx: - constituents_mass = np.array(snapshot.particles.mass) - constituents_pos = np.array(snapshot.particles.position) - total_mass = np.sum(constituents_mass[idx]) - com_mass.append(total_mass) - com = ( - np.sum( - constituents_pos[idx] * constituents_mass[idx, np.newaxis], - axis=0, - ) - / total_mass - ) - com_position.append(com) - com_moi.append( - moit( - points=constituents_pos[idx], - masses=constituents_mass[idx], - center=com, - ) - ) - return com_mass, com_position, com_moi - - def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): """Create rigid bodies from a snapshot. @@ -138,3 +111,30 @@ def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): "orientations": [(1.0, 0.0, 0.0, 0.0)] * len(local_coords), } return rigid_frame, rigid_constrain + + +def _get_com_mass_pos_moi(snapshot, rigid_const_idx): + com_mass = [] + com_position = [] + com_moi = [] + for idx in rigid_const_idx: + constituents_mass = np.array(snapshot.particles.mass) + constituents_pos = np.array(snapshot.particles.position) + total_mass = np.sum(constituents_mass[idx]) + com_mass.append(total_mass) + com = ( + np.sum( + constituents_pos[idx] * constituents_mass[idx, np.newaxis], + axis=0, + ) + / total_mass + ) + com_position.append(com) + com_moi.append( + moit( + points=constituents_pos[idx], + masses=constituents_mass[idx], + center=com, + ) + ) + return com_mass, com_position, com_moi From 394be0843c580a52970015dc6010cd53423e2b81 Mon Sep 17 00:00:00 2001 From: chrisjonesbsu Date: Fri, 8 Dec 2023 15:25:06 -0700 Subject: [PATCH 20/37] add doc strings --- flowermd/library/forcefields.py | 37 +++++++++++++++++++++++++++++++-- flowermd/library/polymers.py | 17 ++++++++++++++- 2 files changed, 51 insertions(+), 3 deletions(-) diff --git a/flowermd/library/forcefields.py b/flowermd/library/forcefields.py index aaa33c04..e13f3c7e 100644 --- a/flowermd/library/forcefields.py +++ b/flowermd/library/forcefields.py @@ -484,13 +484,46 @@ def _check_widths(self): class EllipsoidForcefield(BaseHOOMDForcefield): - """A forcefield for modeling anisotropic bead polymers.""" + """A forcefield for modeling anisotropic bead polymers. + + Notes + ----- + This is designed to be used with `flowermd.library.polymers.EllipsoidChain` + and uses ghost particles of type "A" and "B" for intra-molecular + interactions of bonds and two-body angles. + Ellipsoid centers (type "R") are used in inter-molecular pair interations. + + The set of interactions are: + 1. `hoomd.md.bond.Harmonic`: Models ellipsoid bonds as tip-to-tip bonds + 2. `hoomd.md.angle.Harmonic`: Models angles of two neighboring ellipsoids. + 3. `hoomd.md.pair.aniso.GayBerne`" Model pair interactions between beads. + + Parameters + ---------- + epsilon : float, required + energy + lpar: float, required + Semi-axis length of the ellipsoid along the major axis. + lperp : float, required + Semi-axis length of the ellipsoid along the minor axis. + r_cut : float, required + Cut off radius for pair interactions + bond_k : float, required + Spring constant in harmonic bond + bond_r0: float, required + Equilibrium tip-to-tip bond length. + angle_k : float, required + Spring constant in harmonic angle. + angle_theta0: float, required + Equilibrium angle between 2 consecutive beads. + + """ def __init__( self, epsilon, - lperp, lpar, + lperp, bead_length, r_cut, bond_k, diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index 3883c879..7f23c3e6 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -267,6 +267,19 @@ def _build(self, length): class EllipsoidChain(Polymer): """Create an ellipsoid polymer chain. + This is a coarse-grained molecule where each monomer is modeled + as an anisotropic bead (i.e. ellipsoid). + + Notes + ----- + In order to form chains of connected ellipsoids, "ghost" + particles of types "A" and "B" are used. + + This is meant to be used with + `flowermd.library.forcefields.EllipsoidForcefield` + and requires using `flowermd.utils.rigid_body` to set up + the rigid bodies correctly in HOOMD-Blue. + Parameters ---------- lengths : int, required @@ -274,11 +287,13 @@ class EllipsoidChain(Polymer): num_mols : int, required The number of chains to create. bead_length : float, required - The length of the ellipsoid bead. + The length of the ellipsoid bead along its major axis. bead_mass : float, required The mass of the ellipsoid bead. bond_length : float, required The bond length between connected beads. + This is used as the bond length between ellipsoid tips + rather than between ellipsoid centers. """ From 2b49384b7dd31bb674f1a1fff808fcb641ba299c Mon Sep 17 00:00:00 2001 From: marjanalbouye Date: Fri, 8 Dec 2023 16:57:26 -0700 Subject: [PATCH 21/37] add example code to docstring --- flowermd/utils/rigid_body.py | 41 ++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index 5ed7065c..954d12c5 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -21,6 +21,47 @@ def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): The snapshot of the rigid bodies. rigid_constrain : hoomd.md.constrain.Rigid The rigid body constrain object. + + Examples + -------- + This example demonstrates how to create a rigid body snapshot from a + snapshot of a system of ellipsoids. The ellipsoids are created using + the EllipsoidChain class from the `flowermd.library.polymers`. + The `rigid_frame` snapshot contains all rigid body centers and the + constituent particles along with information about the positions of + center of mass, orientations and moment of inertia. The + `rigid_constrain` object is used by hoomd in the simulation class to + constrain the rigid bodies. + + :: + + from flowermd.library.polymers import EllipsoidChain + from flowermd.library.forcefields import EllipsoidForcefield + from flowermd.base import Pack + from flowermd.base import Simulation + + ellipsoid_chain = EllipsoidChain(lengths=9, num_mols=20, + bead_length=1, bead_mass=100, + bond_length=0.01) + + system = Pack(molecules=ellipsoid_chain, density=0.1, + fix_orientation=True) + + rigid_frame, rigid = create_rigid_body(system.hoomd_snapshot, + ellipsoid_chain.bead_constituents_types) + + + ellipsoid_ff = EllipsoidForcefield(epsilon= 1.0, lperp=0.5 , + lpar=1.0, bead_length=1, + r_cut=3, bond_k=500, + bond_r0=0.1) + + simulation = Simulation(initial_state=rigid_frame, + forcefield=ellipsoid_ff.hoomd_forces, + rigid_constraint=rigid) + + simulation.run_NVT(n_steps=10000, kT=3.5, tau_kt=1) + """ # find typeid sequence of the constituent particles types in a rigid bead p_types = np.array(snapshot.particles.types) From 418bc293fc157c69e8d4e26eaa2985bd6828b953 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 11 Dec 2023 11:47:30 -0700 Subject: [PATCH 22/37] replace bead_length with lpar in EllipsoidChain class --- flowermd/library/polymers.py | 24 ++++++++++++------------ flowermd/tests/library/test_polymers.py | 2 +- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index 7f23c3e6..4e0defd0 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -286,8 +286,8 @@ class EllipsoidChain(Polymer): The number of monomer repeat units in the chain. num_mols : int, required The number of chains to create. - bead_length : float, required - The length of the ellipsoid bead along its major axis. + lpar : float, required + The semi-axis length of the ellipsoid bead along its major axis. bead_mass : float, required The mass of the ellipsoid bead. bond_length : float, required @@ -297,31 +297,31 @@ class EllipsoidChain(Polymer): """ - def __init__(self, lengths, num_mols, bead_length, bead_mass, bond_length): + def __init__(self, lengths, num_mols, lpar, bead_mass, bond_length): self.bead_mass = bead_mass self.bead_bond_length = bond_length - self.bead_length = bead_length - super(EllipsoidChain, self).__init__(lengths=lengths, num_mols=num_mols) + self.lpar = lpar # get the indices of the particles in a rigid body self.bead_constituents_types = ["A", "A", "B", "B"] + super(EllipsoidChain, self).__init__(lengths=lengths, num_mols=num_mols) def _build(self, length): # Build up ellipsoid bead bead = mb.Compound(name="ellipsoid") head = mb.Compound( - pos=(self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 4 + pos=(self.lpar, 0, 0), name="A", mass=self.bead_mass / 4 ) tail = mb.Compound( - pos=(-self.bead_length / 2, 0, 0), name="A", mass=self.bead_mass / 4 + pos=(-self.lpar, 0, 0), name="A", mass=self.bead_mass / 4 ) head_mid = mb.Compound( - pos=(self.bead_length / 4, 0, 0), name="B", mass=self.bead_mass / 4 + pos=(self.lpar / 2, 0, 0), name="B", mass=self.bead_mass / 4 ) tail_mid = mb.Compound( - pos=(-self.bead_length / 4, 0, 0), name="B", mass=self.bead_mass / 4 + pos=(-self.lpar / 2, 0, 0), name="B", mass=self.bead_mass / 4 ) bead.add([head, tail, head_mid, tail_mid]) - + # Build the bead chain chain = mbPolymer() chain.add_monomer( bead, @@ -336,7 +336,7 @@ def _build(self, length): chain.freud_generate_bonds( name_a="B", name_b="B", - dmin=self.bead_length / 2 - 0.1, - dmax=self.bead_length / 2 + self.bead_bond_length + 0.1, + dmin=self.lpar - 0.1, + dmax=self.lpar + self.bead_bond_length + 0.1, ) return chain diff --git a/flowermd/tests/library/test_polymers.py b/flowermd/tests/library/test_polymers.py index d51d6502..18d85c17 100644 --- a/flowermd/tests/library/test_polymers.py +++ b/flowermd/tests/library/test_polymers.py @@ -111,7 +111,7 @@ def test_ellipsoid_chain(self): ellipsoid_chain = EllipsoidChain( lengths=4, num_mols=2, - bead_length=1, + lpar=0.5, bead_mass=100, bond_length=0.01, ) From e635f137c6ddc4e46d5fb7e0373ebeae88ae2465 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 11 Dec 2023 12:17:16 -0700 Subject: [PATCH 23/37] add rigid body handling of mass --- flowermd/base/simulation.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 0ea87608..97398abd 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -295,7 +295,16 @@ def volume(self): def mass_reduced(self): """The total mass of the system in reduced units.""" with self.state.cpu_local_snapshot as snap: - return sum(snap.particles.mass) + if self._rigid_constraint: + last_body_tag = 0 + for body_tag in snap.particles.body: + if body_tag == last_body_tag: + last_body_tag += 1 + else: + break + return sum(snap.particles.mass[last_body_tag:]) + else: + return sum(snap.particles.mass) @property def mass(self): From 0c719c6afa7346178a04f49ad19b2876c26132f7 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 11 Dec 2023 13:59:19 -0700 Subject: [PATCH 24/37] add overlap to Pack --- flowermd/base/system.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/flowermd/base/system.py b/flowermd/base/system.py index abd2f9ad..721c94d5 100644 --- a/flowermd/base/system.py +++ b/flowermd/base/system.py @@ -701,6 +701,8 @@ class Pack(System): The factor by which to expand the box for packing. edge : float, default 0.2 The space (nm) between the edge of the box and the molecules. + overlap : float, default 0.2 + Minimum separation (nm) between particles of different molecules. .. warning:: @@ -729,10 +731,12 @@ def __init__( base_units=dict(), packing_expand_factor=5, edge=0.2, + overlap=0.2, fix_orientation=False, ): self.packing_expand_factor = packing_expand_factor self.edge = edge + self.overlap = overlap self.fix_orientation = fix_orientation super(Pack, self).__init__( molecules=molecules, @@ -746,7 +750,7 @@ def _build_system(self): compound=self.all_molecules, n_compounds=[1 for i in self.all_molecules], box=list(self._target_box * self.packing_expand_factor), - overlap=0.2, + overlap=self.overlap, edge=self.edge, fix_orientation=self.fix_orientation, ) From 1e9cd6f249a20357e0b3cb20e430c4d305010e56 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 12 Dec 2023 13:42:10 -0700 Subject: [PATCH 25/37] remove bead_length from Ellipsoid FF class, fix unit test --- flowermd/library/forcefields.py | 4 +--- flowermd/tests/library/test_forcefields.py | 1 - 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/flowermd/library/forcefields.py b/flowermd/library/forcefields.py index e13f3c7e..28bda73a 100644 --- a/flowermd/library/forcefields.py +++ b/flowermd/library/forcefields.py @@ -524,7 +524,6 @@ def __init__( epsilon, lpar, lperp, - bead_length, r_cut, bond_k, bond_r0, @@ -534,7 +533,6 @@ def __init__( self.epsilon = epsilon self.lperp = lperp self.lpar = lpar - self.bead_length = bead_length self.r_cut = r_cut self.bond_k = bond_k self.bond_r0 = bond_r0 @@ -548,7 +546,7 @@ def _create_forcefield(self): # Bonds bond = hoomd.md.bond.Harmonic() bond.params["A-A"] = dict(k=self.bond_k, r0=self.bond_r0) - bond.params["B-B"] = dict(k=0, r0=self.bead_length / 2) + bond.params["B-B"] = dict(k=0, r0=self.lpar) forces.append(bond) # Angles if all([self.angle_k, self.angle_theta0]): diff --git a/flowermd/tests/library/test_forcefields.py b/flowermd/tests/library/test_forcefields.py index eb8622f3..7c8eecc4 100644 --- a/flowermd/tests/library/test_forcefields.py +++ b/flowermd/tests/library/test_forcefields.py @@ -99,7 +99,6 @@ def test_ellipsoid_ff(self): epsilon=1.0, lperp=0.5, lpar=1.0, - bead_length=1, r_cut=3, bond_k=500, bond_r0=0.1, From 74eae9d44a3d756296f5ce1a252b716a1aa213d0 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 12 Dec 2023 13:45:50 -0700 Subject: [PATCH 26/37] fix rigid unit test --- flowermd/tests/utils/test_rigid_body.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/flowermd/tests/utils/test_rigid_body.py b/flowermd/tests/utils/test_rigid_body.py index 9f83e24d..081104bc 100644 --- a/flowermd/tests/utils/test_rigid_body.py +++ b/flowermd/tests/utils/test_rigid_body.py @@ -11,12 +11,15 @@ def test_ellipsoid_create_rigid_body(self): ellipsoid_chain = EllipsoidChain( lengths=4, num_mols=2, - bead_length=1, + lpar=0.5, bead_mass=100, bond_length=0.01, ) system = Pack( - molecules=ellipsoid_chain, density=0.1, fix_orientation=True + molecules=ellipsoid_chain, + density=0.1, + base_units=dict(), + fix_orientation=True, ) rigid_frame, rigid = create_rigid_body( From 9d963fe9cca485fb7059d3701e3f01db19b35115 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 12 Dec 2023 14:20:37 -0700 Subject: [PATCH 27/37] add 2 unit tests for rigid simulations --- flowermd/tests/base/test_simulation.py | 45 ++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index 0c8e9810..5c478fee 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -9,7 +9,11 @@ import unyt as u from flowermd import Simulation +from flowermd.base import Pack +from flowermd.library.forcefields import EllipsoidForcefield +from flowermd.library.polymers import EllipsoidChain from flowermd.tests import BaseTest +from flowermd.utils import create_rigid_body class TestSimulate(BaseTest): @@ -293,6 +297,47 @@ def test_pickle_ff(self, benzene_system): assert type(i) is type(j) os.remove("forcefield.pickle") + def test_bad_rigid(self, benzene_system): + with pytest.raises(ValueError): + Simulation.from_system(benzene_system, rigid_constraint="A") + + def test_rigid_sim(self): + ellipsoid_chain = EllipsoidChain( + lengths=4, + num_mols=2, + lpar=0.5, + bead_mass=100, + bond_length=0.01, + ) + system = Pack( + molecules=ellipsoid_chain, + density=0.1, + base_units=dict(), + fix_orientation=True, + ) + ellipsoid_ff = EllipsoidForcefield( + lpar=0.5, + lperp=0.25, + epsilon=1.0, + r_cut=2.0, + bond_k=500, + bond_r0=0.01, + angle_k=250, + angle_theta0=2.2, + ) + rigid_frame, rigid = create_rigid_body( + system.hoomd_snapshot, + ellipsoid_chain.bead_constituents_types, + bead_name="R", + ) + sim = Simulation( + initial_state=rigid_frame, + forcefield=ellipsoid_ff.hoomd_forces, + rigid_constraint=rigid, + ) + sim.run_NVT(n_steps=0, kT=1.0, tau_kt=sim.dt * 100) + assert sim.integrator.integrate_rotational_dof is True + def test_save_restart_gsd(self, benzene_system): sim = Simulation.from_system(benzene_system) sim.save_restart_gsd("restart.gsd") From 375f40c34e1d74ae8ea02904c12aa3432bf5d84e Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 13 Dec 2023 12:20:32 -0700 Subject: [PATCH 28/37] fix check for body tags in mass, add assertion for mass --- flowermd/base/simulation.py | 7 ++++--- flowermd/tests/base/test_simulation.py | 1 + 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 97398abd..ef15b351 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -296,13 +296,14 @@ def mass_reduced(self): """The total mass of the system in reduced units.""" with self.state.cpu_local_snapshot as snap: if self._rigid_constraint: - last_body_tag = 0 + last_body_tag = -1 for body_tag in snap.particles.body: - if body_tag == last_body_tag: + if body_tag > last_body_tag: last_body_tag += 1 else: break - return sum(snap.particles.mass[last_body_tag:]) + cut_index = last_body_tag + 1 + return sum(snap.particles.mass[cut_index:]) else: return sum(snap.particles.mass) diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py index 5c478fee..6c08280f 100644 --- a/flowermd/tests/base/test_simulation.py +++ b/flowermd/tests/base/test_simulation.py @@ -337,6 +337,7 @@ def test_rigid_sim(self): ) sim.run_NVT(n_steps=0, kT=1.0, tau_kt=sim.dt * 100) assert sim.integrator.integrate_rotational_dof is True + assert sim.mass_reduced == 800.0 def test_save_restart_gsd(self, benzene_system): sim = Simulation.from_system(benzene_system) From 72fdfba85221376be9e8927d1698624f5ca9db68 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 13 Dec 2023 12:39:11 -0700 Subject: [PATCH 29/37] ignore E203 in precommit --- .pre-commit-config.yaml | 1 + flowermd/base/simulation.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 779217b1..61b67cb7 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -37,6 +37,7 @@ repos: - id: flake8 args: - --max-line-length=80 + - --extend-ignore=E203 exclude: '__init__.py' - repo: https://github.com/pycqa/pydocstyle diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index ef15b351..9137ad90 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -302,8 +302,8 @@ def mass_reduced(self): last_body_tag += 1 else: break - cut_index = last_body_tag + 1 - return sum(snap.particles.mass[cut_index:]) + # cut_index = last_body_tag + 1 + return sum(snap.particles.mass[last_body_tag + 1 :]) else: return sum(snap.particles.mass) From b37cb0d78032f40573cbbb3e1d0f0fbccbd4207a Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 13 Dec 2023 12:39:28 -0700 Subject: [PATCH 30/37] remove unused var --- flowermd/base/simulation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 9137ad90..288a64bd 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -302,7 +302,6 @@ def mass_reduced(self): last_body_tag += 1 else: break - # cut_index = last_body_tag + 1 return sum(snap.particles.mass[last_body_tag + 1 :]) else: return sum(snap.particles.mass) From 3b382c6c10622d1d3690b738a3587b66e6938447 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Wed, 13 Dec 2023 15:20:16 -0700 Subject: [PATCH 31/37] use >= when checking for ascending next body tag --- flowermd/base/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 288a64bd..562d5604 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -298,7 +298,7 @@ def mass_reduced(self): if self._rigid_constraint: last_body_tag = -1 for body_tag in snap.particles.body: - if body_tag > last_body_tag: + if body_tag >= last_body_tag: last_body_tag += 1 else: break From adfa064bd637003116c05e2b1eda91b38904467f Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Thu, 14 Dec 2023 13:22:02 -0700 Subject: [PATCH 32/37] actaully, we shouldn't use >= in body tag count --- flowermd/base/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py index 562d5604..288a64bd 100644 --- a/flowermd/base/simulation.py +++ b/flowermd/base/simulation.py @@ -298,7 +298,7 @@ def mass_reduced(self): if self._rigid_constraint: last_body_tag = -1 for body_tag in snap.particles.body: - if body_tag >= last_body_tag: + if body_tag > last_body_tag: last_body_tag += 1 else: break From 347b331c255843a8ae971866de7bb545a16553ff Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 19 Feb 2024 18:50:44 +0000 Subject: [PATCH 33/37] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- flowermd/tests/utils/test_rigid_body.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/flowermd/tests/utils/test_rigid_body.py b/flowermd/tests/utils/test_rigid_body.py index 081104bc..08ed20d5 100644 --- a/flowermd/tests/utils/test_rigid_body.py +++ b/flowermd/tests/utils/test_rigid_body.py @@ -48,15 +48,9 @@ def test_ellipsoid_create_rigid_body(self): x = points[:, 0] y = points[:, 1] z = points[:, 2] - I_xx = np.sum( - (y**2 + z**2) * system.hoomd_snapshot.particles.mass[:4] - ) - I_yy = np.sum( - (x**2 + z**2) * system.hoomd_snapshot.particles.mass[:4] - ) - I_zz = np.sum( - (x**2 + y**2) * system.hoomd_snapshot.particles.mass[:4] - ) + I_xx = np.sum((y**2 + z**2) * system.hoomd_snapshot.particles.mass[:4]) + I_yy = np.sum((x**2 + z**2) * system.hoomd_snapshot.particles.mass[:4]) + I_zz = np.sum((x**2 + y**2) * system.hoomd_snapshot.particles.mass[:4]) assert np.all( np.isclose( rigid_frame.particles.moment_inertia[0], From b8279fb24407750b62add6f84ae513188fec6a8d Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Sat, 24 Feb 2024 22:42:18 -0700 Subject: [PATCH 34/37] fix rigid body import --- flowermd/utils/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/flowermd/utils/__init__.py b/flowermd/utils/__init__.py index 65e69050..f3e56cb6 100644 --- a/flowermd/utils/__init__.py +++ b/flowermd/utils/__init__.py @@ -1,5 +1,6 @@ from .actions import * from .base_types import HOOMDThermostats +from .rigid_body import create_rigid_body from .utils import ( _calculate_box_length, get_target_box_mass_density, From b09f86ef4b614a87fab69789cd8278eb5019fc9a Mon Sep 17 00:00:00 2001 From: chrisjones4 Date: Wed, 3 Apr 2024 17:37:58 -0600 Subject: [PATCH 35/37] Switch to z-axis for ellipsoid chain alignment, add param for initial orientaiton to rigid body util --- flowermd/base/system.py | 2 +- flowermd/library/polymers.py | 10 +++++----- flowermd/utils/rigid_body.py | 11 ++++++++--- 3 files changed, 14 insertions(+), 9 deletions(-) diff --git a/flowermd/base/system.py b/flowermd/base/system.py index 0dfad388..5e4d0be8 100644 --- a/flowermd/base/system.py +++ b/flowermd/base/system.py @@ -739,7 +739,7 @@ def __init__( def _build_system(self): for mol in self._molecules: - mol._align_backbones_z_axis(heavy_atoms_only=True) + mol._align_backbones_z_axis(heavy_atoms_only=False) next_idx = 0 system = mb.Compound() for i in range(self.n): diff --git a/flowermd/library/polymers.py b/flowermd/library/polymers.py index 41ac203b..2ac9038f 100644 --- a/flowermd/library/polymers.py +++ b/flowermd/library/polymers.py @@ -320,16 +320,16 @@ def _build(self, length): # Build up ellipsoid bead bead = mb.Compound(name="ellipsoid") head = mb.Compound( - pos=(self.lpar, 0, 0), name="A", mass=self.bead_mass / 4 + pos=(0, 0, self.lpar), name="A", mass=self.bead_mass / 4 ) tail = mb.Compound( - pos=(-self.lpar, 0, 0), name="A", mass=self.bead_mass / 4 + pos=(0, 0, -self.lpar), name="A", mass=self.bead_mass / 4 ) head_mid = mb.Compound( - pos=(self.lpar / 2, 0, 0), name="B", mass=self.bead_mass / 4 + pos=(0, 0, self.lpar / 2), name="B", mass=self.bead_mass / 4 ) tail_mid = mb.Compound( - pos=(-self.lpar / 2, 0, 0), name="B", mass=self.bead_mass / 4 + pos=(0, 0, -self.lpar / 2), name="B", mass=self.bead_mass / 4 ) bead.add([head, tail, head_mid, tail_mid]) # Build the bead chain @@ -337,7 +337,7 @@ def _build(self, length): chain.add_monomer( bead, indices=[0, 1], - orientation=[[1, 0, 0], [-1, 0, 0]], + orientation=[[0, 0, 1], [0, 0, -1]], replace=False, separation=self.bead_bond_length, ) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index 954d12c5..56a0b32e 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -5,7 +5,12 @@ from cmeutils.geometry import moit -def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): +def create_rigid_body( + snapshot, + bead_constituents_types, + bead_name="R", + initial_orientation=[1, 0, 0, 0] +): """Create rigid bodies from a snapshot. Parameters @@ -102,7 +107,7 @@ def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): rigid_frame.particles.moment_inertia = np.concatenate( (com_moi, np.zeros((snapshot.particles.N, 3))) ) - rigid_frame.particles.orientation = [(1.0, 0.0, 0.0, 0.0)] * ( + rigid_frame.particles.orientation = [initial_orientation] * ( n_rigid + snapshot.particles.N ) rigid_frame.particles.body = np.concatenate( @@ -149,7 +154,7 @@ def create_rigid_body(snapshot, bead_constituents_types, bead_name="R"): rigid_constrain.body["R"] = { "constituent_types": bead_constituents_types, "positions": local_coords, - "orientations": [(1.0, 0.0, 0.0, 0.0)] * len(local_coords), + "orientations": [initial_orientation] * len(local_coords), } return rigid_frame, rigid_constrain From ef833c9539eeb88255a356e4f1403c8652622196 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 3 Apr 2024 23:38:36 +0000 Subject: [PATCH 36/37] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- flowermd/utils/rigid_body.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/flowermd/utils/rigid_body.py b/flowermd/utils/rigid_body.py index 56a0b32e..4484b1e6 100644 --- a/flowermd/utils/rigid_body.py +++ b/flowermd/utils/rigid_body.py @@ -6,10 +6,10 @@ def create_rigid_body( - snapshot, - bead_constituents_types, - bead_name="R", - initial_orientation=[1, 0, 0, 0] + snapshot, + bead_constituents_types, + bead_name="R", + initial_orientation=[1, 0, 0, 0], ): """Create rigid bodies from a snapshot. From 6e33305e5d2f14b91537d8c116a97271f1458022 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 8 Apr 2024 10:54:22 -0600 Subject: [PATCH 37/37] fix test --- flowermd/base/system.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flowermd/base/system.py b/flowermd/base/system.py index 5e4d0be8..0dfad388 100644 --- a/flowermd/base/system.py +++ b/flowermd/base/system.py @@ -739,7 +739,7 @@ def __init__( def _build_system(self): for mol in self._molecules: - mol._align_backbones_z_axis(heavy_atoms_only=False) + mol._align_backbones_z_axis(heavy_atoms_only=True) next_idx = 0 system = mb.Compound() for i in range(self.n):