-
Notifications
You must be signed in to change notification settings - Fork 164
Add fixes to external bond perception for glycans #1206
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: master
Are you sure you want to change the base?
Conversation
|
@jchodera @peastman : This PR is necessary for openmm/openmmforcefields#180 (comment) (which also needs to be reviewed and merged) |
|
I'd like to note that in order to preserve atom ordering in openmm/openmmforcefields#180 (comment), I had to add the following line (after line 728) to this branch (which I don't think should be included/merged): external_bonds = [external_bonds[0], external_bonds[-1]] + external_bonds[1:-1] if len(external_bonds) > 2 and len(external_bonds[-1].name) <= 2 else external_bonds # Hack to get rid of atom ordering diffs |
|
@zhang-ivy : Can you remind us if this is still relevant and what we need to do to review this? We can get a bugfix out quickly if needed. |
|
@jchodera : Yes this PR is necessary to properly convert amber glycam forcefields to openmm forcefield xmls. Just needs to be reviewed and merged, but no rush on getting a bug fix out quickly, since we already converted the forcefields and merged PRs
|
swails
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we get a test case added here?
parmed/openmm/parameters.py
Outdated
| else: | ||
| element_name = atom.element_name | ||
|
|
||
| # If the number of bonds within the residue does not equal the number of bonds the atom should have, the atom has an external bond |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This comment looks out of date now? Why this change?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, it's out of date. I changed this because it was a mistake in my original code -- we are assuming that any extra bonds are external bonds. If there aren't extra bonds, there are no external bonds.
|
I'm a little concerned that GLYCAM support may risk being broken by refactors if we don't have a test case for it. Even a regression test is enough (take a GLYCAM file and check that the generated parameter set matches something computed with this version of the code). |
|
@swails : Yes, understandable -- I need to refresh my memory on the glycam conversion pipeline, but I'm working on this. |
|
@swails : I've written the following test function (based on def check_glycam_conversion():
"""
Takes leaprc and lib files and uses parmed to write out a ffxml.
Simplified from openmmtools/amber/convert_amber.py's `convert_leaprc()`
Note that the output of this function is not the final GLYCAM_06j-1.xml in openmmforcefields/amber/ffxml nor
openmm/app/data/amber14/. Additional post processing steps are necessary to obtain the final ffxml version.
See `modify_glycan_ffxml()` in openmmforcefields/amber/convert_amber.py.
"""
import os
import parmed
from io import StringIO
import re
_loadoffre = re.compile(r'loadOff (\S*)', re.I)
_sourcere = re.compile(r'source (\S*)', re.I)
ignore = {'solvents.lib', 'atomic_ions.lib', 'ions94.lib', 'ions91.lib',
'phosphoaa10.lib'}
input_file = "glycam/leaprc.GLYCAM_06j-1"
output_file = "glycam/GLYCAM_06j-1.xml"
# Pre-process leaprc before converting
with open(input_file, 'rt') as f:
lines = [ line if '#' not in line else line[:line.index('#')] for line in f.readlines() ]
new_lines = []
for line in lines:
if (ignore is not None and _loadoffre.findall(line) and
_loadoffre.findall(line)[0] in ignore):
continue
new_lines += line
leaprc = StringIO(''.join(new_lines))
# Convert leaprc into ffxml
override_level = 1
params = parmed.amber.AmberParameterSet.from_leaprc(leaprc)
params = parmed.openmm.OpenMMParameterSet.from_parameterset(params, remediate_residues=(False))
if override_level:
for name, residue in params.residues.items():
residue.override_level = override_level
params.write(output_file, provenance=None, write_unused=True, improper_dihedrals_ordering='amber', skip_duplicates=False, is_glycam=True)
We'll want to check that the ffxml written by this function is the same as the attached file ( |
|
I'll add this. Thanks! |
This PR includes the following changes that fix external bond perception for glycans: