From 3c1a20eb1511ca276fb71b103821c1d6c41dcfd4 Mon Sep 17 00:00:00 2001 From: Jenny Fothergill Date: Fri, 10 Jun 2022 11:02:59 -0600 Subject: [PATCH 1/4] Start adding function for finegraining --- grits/finegrain.py | 42 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/grits/finegrain.py b/grits/finegrain.py index 7c4e3ec..04e2427 100644 --- a/grits/finegrain.py +++ b/grits/finegrain.py @@ -6,8 +6,46 @@ from mbuild import Compound, Particle, load -from grits.utils import align, get_hydrogen, get_index - +from grits.utils import align, get_hydrogen, get_index, snap_molecules + + + def _get_compounds( + gsdfile, scale, conversion_dict, frame + ): + """Get compounds for each molecule in the gsd trajectory.""" + with gsd.hoomd.open(gsdfile) as t: + snap = t[frame] + + # Use the conversion dictionary to map particle type to element symbol + if conversion_dict is not None: + snap.particles.types = [ + conversion_dict[i].symbol for i in snap.particles.types + ] + # Break apart the snapshot into separate molecules + molecules = snap_molecules(snap) + mol_inds = [] + for i in range(max(molecules) + 1): + mol_inds.append(np.where(molecules == i)[0]) + + # Convert each unique molecule to a compound + system = Compound() + for inds in mol_inds: + l = len(inds) + compound = comp_from_snapshot(snap, inds, scale=scale), + mol = compound.to_pybel() + mol.OBMol.PerceiveBondOrders() + + # Add hydrogens + n_atoms = mol.OBMol.NumAtoms() + # This is a goofy work around necessary for the aromaticity + # to be set correctly. + with tempfile.NamedTemporaryFile() as f: + mol.write(format="mol2", filename=f.name, overwrite=True) + mol = list(pybel.readfile("mol2", f.name))[0] + + mol.addh() + n_atoms2 = mol.OBMol.NumAtoms() + system.add(Compound().from_pybel(mol)) def backmap(cg_compound): """Backmap a fine-grained representation onto a coarse one. From 2e1f1409afaddb7d92bc2a230a6a96f5872069ca Mon Sep 17 00:00:00 2001 From: Jenny Fothergill Date: Tue, 5 Jul 2022 17:05:09 -0600 Subject: [PATCH 2/4] Fix indentation and imports --- grits/finegrain.py | 87 +++++++++++++++++++++++++--------------------- 1 file changed, 47 insertions(+), 40 deletions(-) diff --git a/grits/finegrain.py b/grits/finegrain.py index 04e2427..0cced2f 100644 --- a/grits/finegrain.py +++ b/grits/finegrain.py @@ -2,50 +2,57 @@ __all__ = ["backmap"] import itertools as it +import tempfile from collections import defaultdict +import numpy as np from mbuild import Compound, Particle, load +from openbabel import pybel + +from grits.utils import ( + align, + comp_from_snapshot, + get_hydrogen, + get_index, + snap_molecules, +) + + +def _get_compounds(snap, conversion_dict=None, scale=1.0): + """Get compounds for each molecule in the gsd snapshot.""" + # Use the conversion dictionary to map particle type to element symbol + if conversion_dict is not None: + snap.particles.types = [ + conversion_dict[i].symbol for i in snap.particles.types + ] + # Break apart the snapshot into separate molecules + molecules = snap_molecules(snap) + mol_inds = [] + for i in range(max(molecules) + 1): + mol_inds.append(np.where(molecules == i)[0]) + + # Convert each unique molecule to a compound + system = Compound() + for inds in mol_inds: + l = len(inds) + compound = comp_from_snapshot(snap, inds, scale=scale) + mol = compound.to_pybel() + mol.OBMol.PerceiveBondOrders() + + # Add hydrogens + n_atoms = mol.OBMol.NumAtoms() + # This is a goofy work around necessary for the aromaticity + # to be set correctly. + with tempfile.NamedTemporaryFile() as f: + mol.write(format="mol2", filename=f.name, overwrite=True) + mol = list(pybel.readfile("mol2", f.name))[0] + + mol.addh() + n_atoms2 = mol.OBMol.NumAtoms() + return mol + system.add(Compound().from_pybel(mol)) + return system -from grits.utils import align, get_hydrogen, get_index, snap_molecules - - - def _get_compounds( - gsdfile, scale, conversion_dict, frame - ): - """Get compounds for each molecule in the gsd trajectory.""" - with gsd.hoomd.open(gsdfile) as t: - snap = t[frame] - - # Use the conversion dictionary to map particle type to element symbol - if conversion_dict is not None: - snap.particles.types = [ - conversion_dict[i].symbol for i in snap.particles.types - ] - # Break apart the snapshot into separate molecules - molecules = snap_molecules(snap) - mol_inds = [] - for i in range(max(molecules) + 1): - mol_inds.append(np.where(molecules == i)[0]) - - # Convert each unique molecule to a compound - system = Compound() - for inds in mol_inds: - l = len(inds) - compound = comp_from_snapshot(snap, inds, scale=scale), - mol = compound.to_pybel() - mol.OBMol.PerceiveBondOrders() - - # Add hydrogens - n_atoms = mol.OBMol.NumAtoms() - # This is a goofy work around necessary for the aromaticity - # to be set correctly. - with tempfile.NamedTemporaryFile() as f: - mol.write(format="mol2", filename=f.name, overwrite=True) - mol = list(pybel.readfile("mol2", f.name))[0] - - mol.addh() - n_atoms2 = mol.OBMol.NumAtoms() - system.add(Compound().from_pybel(mol)) def backmap(cg_compound): """Backmap a fine-grained representation onto a coarse one. From c23c70d44b053c550375e4589456c5a5a64c7de1 Mon Sep 17 00:00:00 2001 From: Jenny Fothergill Date: Tue, 5 Jul 2022 18:56:59 -0600 Subject: [PATCH 3/4] Add option whether to shift coordinates --- grits/utils.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/grits/utils.py b/grits/utils.py index 656afe6..08da0e5 100644 --- a/grits/utils.py +++ b/grits/utils.py @@ -25,7 +25,7 @@ def default(self, obj): return super(NumpyEncoder, self).default(obj) -def comp_from_snapshot(snapshot, indices, scale=1.0): +def comp_from_snapshot(snapshot, indices, scale=1.0, shift_coords=True): """Convert particles by indices from a Snapshot to a Compound. Parameters @@ -36,6 +36,8 @@ def comp_from_snapshot(snapshot, indices, scale=1.0): Indices of the particles to be added to the compound. scale : float, default 1.0 Value by which to scale the length values + shift_coords : bool, default True + Whether to shift the coords by half the box lengths Returns ------- @@ -56,8 +58,11 @@ def comp_from_snapshot(snapshot, indices, scale=1.0): lengths=box[:3] * scale, tilt_factors=box[3:] ) - # GSD and HOOMD snapshots center their boxes on the origin (0,0,0) - shift = np.array(comp.box.lengths) / 2 + if shift_coords: + # GSD and HOOMD snapshots center their boxes on the origin (0,0,0) + shift = np.array(comp.box.lengths) / 2 + else: + shift = np.zeros(3) particle_dict = {} # Add particles for i in range(n_atoms): From a6ec58d24802009c0f502eac940a1ae946775cda Mon Sep 17 00:00:00 2001 From: Jenny Fothergill Date: Tue, 5 Jul 2022 18:58:08 -0600 Subject: [PATCH 4/4] Change function name, write directly from pybel.mol to gsd mbuild compounds are very slow... --- grits/finegrain.py | 55 ++++++++++++++++++++++++++++++++++++---------- 1 file changed, 44 insertions(+), 11 deletions(-) diff --git a/grits/finegrain.py b/grits/finegrain.py index 0cced2f..f2de616 100644 --- a/grits/finegrain.py +++ b/grits/finegrain.py @@ -5,9 +5,11 @@ import tempfile from collections import defaultdict +import gsd.hoomd import numpy as np +from ele import element_from_atomic_number from mbuild import Compound, Particle, load -from openbabel import pybel +from openbabel import openbabel, pybel from grits.utils import ( align, @@ -18,8 +20,22 @@ ) -def _get_compounds(snap, conversion_dict=None, scale=1.0): - """Get compounds for each molecule in the gsd snapshot.""" +def _add_hydrogens(snap, conversion_dict=None, scale=1.0): + """Add hydrogens to a united atom gsd snapshot. + + Parameters + ---------- + snap : gsd.hoomd.Snapshot + The united-atom gsd snapshot + conversion_dict : dict, default None + Dictionary to map particle types to their element. + scale : float, default 1.0 + Factor by which to scale length values. + + Returns + ------- + gsd.hoomd.Snapshot + """ # Use the conversion dictionary to map particle type to element symbol if conversion_dict is not None: snap.particles.types = [ @@ -31,16 +47,20 @@ def _get_compounds(snap, conversion_dict=None, scale=1.0): for i in range(max(molecules) + 1): mol_inds.append(np.where(molecules == i)[0]) - # Convert each unique molecule to a compound - system = Compound() + types = [] + coords = [] + bonds = [] + # Read each molecule into a compound + # Then convert to pybel to add hydrogens for inds in mol_inds: l = len(inds) - compound = comp_from_snapshot(snap, inds, scale=scale) + compound = comp_from_snapshot( + snap, inds, scale=scale, shift_coords=False + ) mol = compound.to_pybel() mol.OBMol.PerceiveBondOrders() # Add hydrogens - n_atoms = mol.OBMol.NumAtoms() # This is a goofy work around necessary for the aromaticity # to be set correctly. with tempfile.NamedTemporaryFile() as f: @@ -48,10 +68,23 @@ def _get_compounds(snap, conversion_dict=None, scale=1.0): mol = list(pybel.readfile("mol2", f.name))[0] mol.addh() - n_atoms2 = mol.OBMol.NumAtoms() - return mol - system.add(Compound().from_pybel(mol)) - return system + for atom in mol.atoms: + coords.append(np.array(atom.coords) / 10) + types.append(element_from_atomic_number(atom.atomicnum).symbol) + for bond in openbabel.OBMolBondIter(mol.OBMol): + bonds.append([bond.GetBeginAtomIdx() - 1, bond.GetEndAtomIdx() - 1]) + alltypes = list(set(types)) + typeids = [alltypes.index(i) for i in types] + new_snap = gsd.hoomd.Snapshot() + new_snap.configuration.box = snap.configuration.box + new_snap.particles.N = len(types) + new_snap.particles.types = alltypes + new_snap.particles.typeid = typeids + new_snap.particles.position = np.array(coords) + new_snap.bonds.N = len(bonds) + new_snap.bonds.group = np.array(bonds) + new_snap.validate() + return new_snap def backmap(cg_compound):