cmelab / grits

A toolkit for working with coarse-grain systems
https://grits.readthedocs.io
GNU General Public License v3.0
5 stars 9 forks source link

Mapping of mass from atomistic to cg beads failing because of how pybel smarts matching works #52

Closed chrisjonesBSU closed 1 year ago

chrisjonesBSU commented 1 year ago

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)