Skip to content

Conversation

@zhang-ivy
Copy link
Contributor

This PR includes the following changes that fix external bond perception for glycans:

  • Add exception cases (OME, OLP) to external bond perception function
  • Only designate an atom as containing an external bond if the number of bonds is less than the expected valence
  • Split original external bond perception code from glycan code in _write_omm_residues()

@zhang-ivy
Copy link
Contributor Author

@jchodera @peastman : This PR is necessary for openmm/openmmforcefields#180 (comment) (which also needs to be reviewed and merged)
Could you review when possible?

@zhang-ivy
Copy link
Contributor Author

zhang-ivy commented Nov 8, 2021

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

@jchodera
Copy link
Contributor

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

@zhang-ivy
Copy link
Contributor Author

zhang-ivy commented Apr 11, 2022

@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

Copy link
Contributor

@swails swails left a 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?

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
Copy link
Contributor

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?

Copy link
Contributor Author

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.

@swails
Copy link
Contributor

swails commented Apr 25, 2022

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

@zhang-ivy
Copy link
Contributor Author

@swails : Yes, understandable -- I need to refresh my memory on the glycam conversion pipeline, but I'm working on this.

@zhang-ivy
Copy link
Contributor Author

@swails : I've written the following test function (based on openmmforcefields code) that uses the attached input files:

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 (glycam/GLYCAM_06j-1.xml). What's the easiest way to do this? Is there an easy way in parmed to load in 2 ffxmls and check that they contain the same parameters? Alternatively, I could look into ways for pythonic ways for checking whether 2 xml files are the same.
glycam.zip

@swails
Copy link
Contributor

swails commented May 9, 2022

I'll add this. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants