cmelab / grits Goto Github PK
View Code? Open in Web Editor NEWA toolkit for working with coarse-grain systems
Home Page: https://grits.readthedocs.io
License: GNU General Public License v3.0
A toolkit for working with coarse-grain systems
Home Page: https://grits.readthedocs.io
License: GNU General Public License v3.0
Looks like the quay container is stalled, let's move this to Github Container Repo.
If we want to use the CG_System gsd in MSIBI, the beads should have a mass.
grits is painfully slow on real simulation morphologies...
shout out to @chrisjonesBSU for helpful discussion and alerting me to this issue.
In the instance that you want to map a PEEK monomer to a single bead, you may do the following as in Chris's code (needs mosdef-hub/mbuild#851):
peek_para = mb.load("Oc1ccc(Oc2ccc(C(=O)c3ccccc3)cc2)cc1",smiles=True)
chain = Polymer()
chain.add_monomer(compound=peek_para,
indices = [22, 29],
separation = 0.1376,
replace=True
)
chain.build(n=2, add_hydrogens = True)
cg_beads = {"_A": "Oc1ccc(Oc2ccc(C(=O)c3ccccc3)cc2)cc1"}
with warnings.catch_warnings():
warnings.simplefilter("ignore")
cg_chain = CG_Compound(chain, cg_beads)
cg_chain.visualize(show_atomistic=True).show()
cg_chain.visualize().show()
This is due to the way ring structures are handled.
It could be fixed by adding an allow_overlap kwarg that could bypass this logic:
cg_chain = CG_Compound(chain, cg_beads, allow_overlap=False)
In any workflow that might actually use grits, we would need to be able to save and reload the mapping:
From a conversation with @chrisjonesBSU
chris:
Both backmap and CG_compound use some kind of mapping scheme, but it's formatted differently for each module. backmap takes a dictionary of {bead: smiles} but CG_compound takes a list of [bead, SMARTS]. I think it would be great if the same format was used in both. In this case, I think a dictionary would be best.
me:
great point! the mapping operators are not what I want them to be yet. π€ Basically I want a CG_Compound to be able to be created using SMARTS or a mapping operator. The mapping operator would be solely based on particle indices and would be saved as an attribute to the CG_Compound when the compound is created from SMARTS. Allowing a CG_Compound to be created from a mapping operator would be helpful for using it on simulation outputs that are not in optimal conformations. e.g., We can get the mapping from the initial frame (ideal positions--recognizable to smarts matching) using smarts and use the mapping on the last frame (maybe weird conformation). Another update I want to make is that we keep track of the anchors and bonds in the CG_Compound as theyβre created so that no mapping is necessary to the backmap function. For instance, the only arg to backmap would be the CG compound!
This might be more difficult than expected as I'm unsure if RDKit can infer the bond order: rdkit/rdkit#4082
I could do something like:
import mbuild as mb
from rdkit import Chem
p3ht = mb.load("../grits/tests/assets/P3HT_16.mol2")
mol = p3ht.to_pybel()
mol.OBMol.PerceiveBondOrders()
mol.write(format="sdf", filename="tmp.sdf", overwrite=True)
m = Chem.MolFromMolFile("tmp.sdf")
Chem.SanitizeMol(m)
patt = Chem.MolFromSmarts('c1sccc1')
m.GetSubstructMatches(patt)
But at that point we might as well just use openbabel π
Line 744 in 62b7d93
At this line, if mass_scale is set incorrectly such that we get an empty array, we should a throw a user warning regarding mass_scale being correctly set. Something like "no heavy atoms found using default value"
Ever since we added mapping from atomistic to CG bead masses, the unit test that checks for equivalence in mass has failed. I think I figured out why. openbabel's smarts matching only returns atom indices for the heavy atoms, so when mapping back to the mbuild compounds, the xyz mapping is mostly correct, but the masses are off (hydrogens aren't included).
For example:
SMARTS matching on benzene only returns a list of length 6 even though the pybel molecule has 12 atoms
>>> mol = pybel.readstring("smi","c1ccccc1")
>>> mol.addh()
>>> print(len(mol.atoms))
>>> smarts = pybel.Smarts("c1ccccc1")
>>> smarts.findall(mol)
12
[(1, 6, 5, 4, 3, 2)]
So, where the mapping actually happens, hydrogens are never accounted for
173 def _cg_particles(self):
174 """Set the beads in the coarse-structure."""
175 for key, inds in self.mapping.items():
176 name, smarts = key.split("...")
177 for group in inds:
178 mass = sum([self.atomistic[i].mass for i in group])
179 bead_xyz = self.atomistic.xyz[group, :]
180 avg_xyz = np.mean(bead_xyz, axis=0)
181 bead = Bead(name=name, pos=avg_xyz, smarts=smarts, mass=mass)
182 self.add(bead)
If we knew mapping between detailed and coarse representations were the same 3 atoms (determined by index order. For example, the ath, bth, and cth atoms in the middle ring of perylene) we could use that mapping to get around hard-to-come-up-with SMILES strings
Using SMARTS matching on an United Atom system without the add_hydrogens flag gives a confusing result.
This works
pps = mb.load("c1ccc(S)cc1", smiles=True)
pps.save("pps_box_H.gsd", overwrite=True)
gsdfile = 'pps_box_H.gsd'
system = grits.CG_System(
gsdfile,
beads={"_A": "c1ccc(S)cc1"},
)
system._compounds[0].visualize(show_atomistic = True)
This gives a warning
pps = mb.load("c1ccc(S)cc1", smiles=True)
pps.remove(pps.particles_by_element('H'))
pps.save("pps_box_noH.gsd", overwrite=True)
gsdfile = 'pps_box_noH.gsd'
system = grits.CG_System(
gsdfile,
beads={"_A": "c1ccc(S)cc1"},
#add_hydrogens=True
)
system._compounds[0].visualize(show_atomistic = True)
/Users/noah/miniconda3/envs/grits/lib/python3.8/site-packages/grits/coarsegrain.py:139: UserWarning: c1ccc(S)cc1 not found in compound!
warn(f"{smart_str} not found in compound!")
/Users/noah/miniconda3/envs/grits/lib/python3.8/site-packages/grits/coarsegrain.py:169: UserWarning: Some atoms have been left out of coarse-graining!
warn("Some atoms have been left out of coarse-graining!")
This warning is confusing when reading from a United Atom gsd and using the corresponding SMARTS string. It's not obvious that the hydrogens need to be added back in this case.
We should either add to the warning specifying that we might need the add_hydrogens=True flag or make the behavior automatically add hydrogens back to United Atom systems.
It looks like snap.bonds.groups
is being populated correctly, but bonds.types
and bonds.typeid
are not when saving the coarse-grained GSD file. This is an issue if you're wanting to use GRiTS to generate coarse-grained trajectories for use in MSIBI, or other coarse-graining workflows that need bond and angle information between specific types of beads.
I think enough information is there to finish adding the needed information to snap.bonds
but populating angles (and maybe dihedrals) would require setting up the infrastructure similarly to how it's being done for bonds.
From discussion with @chrisjonesBSU:
So, for [msibi] to work though, the trajectory passed in when creating a state has to already be coarse grained, right? It was ran as an atomistic trajectory, but each frame of the .gsd file would have to be mapped to whatever coarse grain scheme is being used.
Let's work up a workflow that goes from gsd->CG_compound->gsd (usable by msibi)--a good candidate is the propane simulations.
Grits has a conda-forge package! We should mention this in the installation section of the documentation.
It would be nice to be able to use a network-x to categorize a non-atomistic system (e.g., bonded A and B beads) and design a smarts-esque grammar that could group ABB beads, for example.
Fine-graining is not implemented for system sized morphologies. It would be a good proof of concept to show that we can take a UA trajectory, add hydrogens, and run the relaxation workflow Matty described in Fig 8 of http://dx.doi.org/10.1080/08927022.2017.1296958
We could use a similar logic to the new add_hydrogens
kwarg to CG_System
.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
π Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. πππ
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google β€οΈ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.