cmelab / grits Goto Github PK

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)

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!

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.

grits is painfully slow on real simulation morphologies...

• use freud clustering? save doing the same work on a system of identical molecules
• move away from mbuild compound? hierarchical nature not needed, we could use a simpler faster data structure

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 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.

In any workflow that might actually use grits, we would need to be able to save and reload the mapping:

1. Create CG_Compound
2. Save to simulation input format (e.g., gsd)
3. Run simulation
4. Load new structure from simulation output
5. Load mapping
6. Backmap!

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"

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.

• make package dir
• make pip installable
• reorganize/refactor?
• run pre-commit

Grits has a conda-forge package! We should mention this in the installation section of the documentation.

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()

The result is not ideal:

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)

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

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 😂

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.

• create environment.yml
• build container
• write tests

If we want to use the CG_System gsd in MSIBI, the beads should have a mass.