Dihedrals not found - Githubissues

kierannp commented 1 year ago

Bug summary I am trying to type a molecule using proper dihedrals but am receiving an error saying that the dihedrals can't be found:

Exception Traceback (most recent call last) /Users/knp/projects/silk/messing.ipynb Cell 2 in <cell line: 26>() 22 gly = remove_both_terminals(gly) 24 charmm = Forcefield(forcefield_files='ffs/charmm_foyer/glycine.xml',debug=True, validation=True) ---> 26 typed_gly = charmm.apply(gly,verbose=True) 28 typed_gly.save('gly_test.top') 29 typed_gly.save('gly_test.gro')

File ~/miniforge3/envs/mosdef/lib/python3.8/site-packages/foyer/forcefield.py:786, in Forcefield.apply(self, structure, references_file, use_residue_map, assert_bond_params, assert_angle_params, assert_dihedral_params, assert_improper_params, verbose, args, kwargs) 780 typemap = self.run_atomtyping( 781 structure, use_residue_map=use_residue_map, kwargs 782 ) 784 self._apply_typemap(structure, typemap) --> 786 structure = self.parametrize_system( 787 structure=structure, 788 references_file=references_file, 789 assert_bond_params=assert_bond_params, 790 assert_angle_params=assert_angle_params, 791 assert_dihedral_params=assert_dihedral_params, 792 assert_improper_params=assert_improper_params, 793 verbose=verbose, 794 args, 795 **kwargs, 796 ) ... --> 347 raise Exception(msg) 348 else: 349 warnings.warn(msg)

Exception: Parameters have not been assigned to all proper dihedrals. Total system dihedrals: 6, Parameterized dihedrals: 0. Note that if your system contains torsions of Ryckaert-Bellemans functional form, all of these torsions are processed as propers. Code to reproduce the behavior

Please include a code snippet that can be used to reproduce this bug.

import mbuild as mb
from foyer import Forcefield
gly_smiles = 'C(C(=O)O)N'

class Glycine(mb.Compound):
    def __init__(self):
        super(Glycine, self).__init__()
        glycine = mb.load(gly_smiles,smiles=True)
        self.add(glycine)
        self.name = "Glycine"
        self.amine = glycine[4]
        self.carboxyl = glycine[1]
        amine_h = [glycine[8], glycine[9]]
        carboxyl_o = [glycine[3], glycine[7]]
        self.indices = [amine_h, carboxyl_o]

def remove_both_terminals(mol):
    amine_h, carboxyl_o = mol.indices
    mol.remove(amine_h[0])
    mol.remove(carboxyl_o)
    return mol

gly = Glycine()
gly = remove_both_terminals(gly)

charmm = Forcefield(forcefield_files='ffs/charmm_foyer/glycine.xml',debug=True, validation=True)

typed_gly = charmm.apply(gly,verbose=True)

typed_gly.save('gly_test.top')
typed_gly.save('gly_test.gro')

Software versions

Which version of Foyer are you using? (python -c "import foyer; print(foyer.__version__)") 0.11.1
Which version of Python (python --version)? Python 3.8.13
Which operating system? macOS Monterey 12.4

chrisiacovella commented 1 year ago

This isn't an issue with the Foyer code itself, rather that there are 6 dihedrals found in the molecule that haven't been found in the forcefield xml file. I.e., they need to be defined by haven't yet.

I'll note that in some cases you may have quartets that don't require dihedral parameters (especially for CHARMM as the dihedral parameters are corrections to the 1-4 interactions not a replacement). If that is the case, then you'd just want to modify your apply call to tell it to ignore this (i.e., raise a warning, but not throw an error):

typed_gly = charmm.apply(gly,verbose=True, assert_dihedral_params=False)

Bonds, angles, dihedrals, and impropers each have an associated assert flag that can be set (i.e., assert_bond_params, assert_angle_params, assert_dihedral_params, assert_improper_params). This is especially useful as, e.g., a CG forcefield might not need to define dihedral or angle parameters or one might use fictitious bonds in mbuild to orient building blocks in the system, but not want them to be actually chemically bound. See the documentation about this

I'm going to close this issue as it appears the code is functioning as expected, but send me a message if you have more questions.

CalCraven commented 1 year ago

@chrisiacovella I think there is more going on here than meets the eye. All of the missing dihedrals seem to be found in the forcefield. Additionally, the atomtyping works using the GMSO methodology:

from forcefield_utilities import FoyerFFs
loader = FoyerFFs()
gmso_ff = loader.load('ffs/charmm_foyer/glycine.xml').to_gmso_ff()
from gmso.parameterization import apply
from gmso.external.convert_mbuild import from_mbuild
gmso_top = from_mbuild(gly) #object created using above recipe
apply(gmso_top, gmso_ff, identify_connections=True)
print(gmso_top.is_typed)

Since this seems to work properly, there seems to be something more than just missing parameters causing the openmm issues seen here. I can spend a little bit of time better showing the dihedrals in the forcefield, and the identities of the missing dihedrals.

daico007 commented 1 year ago

I checked this XML manually and the dihedral is actually there (for example: <Proper class1="O" class2="C" class3="CT2" class4="NH1" periodicity1="1" phase1="0.0" k1="0.0" />), but was just not being found. This might be on the foyer side (or our use of openmm/parmed).

kierannp commented 1 year ago

glycine.txt

kierannp commented 1 year ago

Heres the xml for this molecule