ERROR for assigning angles.

[laowang-os]

Hello, I'm a beginner of Foyer. Forgive my boldness in initiating this corrspondence without prior acquaintance.
Here is the situations I met.

Firstly I ran the following codes for generating a '.pdb' stucture file by mbuild successfully.
import mbuild as mb
comp_1 = mb.load('[H]c1[nH]c([H])c2c1CCC2', smiles=True)
comp_1.save('comp_1.pdb')

Secondly I ran the following codes for assigning parameters to my '.pdb' file.
import foyer
import parmed
mol = parmed.load_file("comp_1.pdb")
ff = foyer.Forcefield(forcefield_files=['oplsaa.xml'])
mol_ff = ff.apply(mol)
Then I got the error as following:
FoyerError: Found multiple types for atom 3 (6): ['opls_544', 'opls_141'].
So I used the 'overrides' command to change the line from
<Type name="opls_544" class="CS" element="C" mass="12.01100" def="[C;X3;r5]([C;%opls_543])" desc="C3 in pyrrole" overrides="opls_142" doi="10.1021/jp981200o"/>
to
<Type name="opls_544" class="CS" element="C" mass="12.01100" def="[C;X3;r5]([C;%opls_543])" desc="C3 in pyrrole" overrides="opls_142,opls_141" doi="10.1021/jp981200o"/>
at the same time I created a new file called 'my_oplsaa.xml'.

Thirdly I ran the codes
import foyer
import parmed
mol = parmed.load_file("comp_1.pdb")
ff = foyer.Forcefield(forcefield_files=['my_oplsaa.xml'])
mol_ff = ff.apply(mol)
Then I got the following issue.
Exception: Parameters have not been assigned to all angles. Total system angles: 33, Parameterized angles: 31

For now I have two questions: 1. How could I know which two angles didn't parameterized? 2. Once I know the exact angles how could I add the information or modify the '.xml' file and where is the references?

Additionally, There is also other questions about the Foyer website. In the subsection 'Validation of force field files' of reference section, the link 'XML schema definition' and 'syntax' are invalid nowadays.

I appreciate your time for these issues and regret any disruptiopn my massage may cause.

[chrisjonesBSU]

No need to apologize! The documentation for the mosdef packages needs to be updated, and raising and issue with a question is totally fine.

We created the General Molecular Simulation Object (GMSO) package that is essentially replacing foyer as the user-facing package (it still uses foyer in the background). I'd suggest using GMSO in your workflow. It has a more useful API for seeing and editing your topology and forcefield parameters.

In my foyer installation, I edited the xml file to fix the override issue, just as you did, then running this script:

import mbuild as mb
import gmso
from gmso.core.forcefield import ForceField
from gmso.external import from_mbuild
from gmso.parameterization import apply


gmso_oplsaa = ForceField("oplsaa") # load the oplsaa from the installed version of foyer into a GMSO forcefield
comp_1 = mb.load('[H]c1[nH]c([H])c2c1CCC2', smiles=True) # Create the combound in mBuild
gmso_topology = from_mbuild(comp_1) # Convert mBuild Compound into a GMSO Toplogy
apply(top=gmso_topology, forcefields=gmso_oplsaa, identify_connections=True) # Use the apply method from GMSO

This raises this error message:

ParameterizationError: No parameters found for connection <Angle Angle, id: 140190552781472> , group: angle_type, identifiers: [['opls_544', 'opls_544', 'opls_136'], ['CS', 'CS', 'CT']] in the Forcefield.

You can make it so missing types are ignored, so you can proceed without an error, then using the API to begin inspecting what might be missing:

gmso_oplsaa = ForceField("oplsaa") # load the oplsaa from the installed version of foyer into a GMSO forcefield
comp_1 = mb.load('[H]c1[nH]c([H])c2c1CCC2', smiles=True) # Create the combound in mBuild
gmso_topology = from_mbuild(comp_1) # Convert mBuild Compound into a GMSO Toplogy
apply(top=gmso_topology, forcefields=gmso_oplsaa, identify_connections=True, ignore_params=["angles", "improper"], remove_untyped=False) # Use the apply method from GMSO

Iterate through the angles types, you can see some are None

for angle in gmso_topology.angles:
    print(angle.angle_type)

Look at the atom types for any angle types that are None

for angle in gmso_topology.angles:
    if angle.angle_type is None:
        atom_types = [a.atom_type.name for a in angle.connection_members]
        atom_classes = [a.atom_type.atomclass for a in angle.connection_members]
        print(atom_types)
        print(atom_classes)
        print()


['opls_544', 'opls_544', 'opls_136']
['CS', 'CS', 'CT']

['opls_544', 'opls_544', 'opls_136']
['CS', 'CS', 'CT']

[chrisjonesBSU]

One more thing to add. Bonds, angles and dihedrals use the atom classes (instead of atom types), so if you were to add the missing angle (and its parameters) to your fork of foyer, I would suggest using CS, CS, CT for that angle entry.

[laowang-os]

One more thing to add. Bonds, angles and dihedrals use the atom classes (instead of atom types), so if you were to add the missing angle (and its parameters) to your fork of foyer, I would suggest using CS, CS, CT for that angle entry.

THANKS!!!😊

[laowang-os]

No need to apologize! The documentation for the mosdef packages needs to be updated, and raising and issue with a question is totally fine.

We created the General Molecular Simulation Object (GMSO) package that is essentially replacing foyer as the user-facing package (it still uses foyer in the background). I'd suggest using GMSO in your workflow. It has a more useful API for seeing and editing your topology and forcefield parameters.

In my foyer installation, I edited the xml file to fix the override issue, just as you did, then running this script:

import mbuild as mb
import gmso
from gmso.core.forcefield import ForceField
from gmso.external import from_mbuild
from gmso.parameterization import apply


gmso_oplsaa = ForceField("oplsaa") # load the oplsaa from the installed version of foyer into a GMSO forcefield
comp_1 = mb.load('[H]c1[nH]c([H])c2c1CCC2', smiles=True) # Create the combound in mBuild
gmso_topology = from_mbuild(comp_1) # Convert mBuild Compound into a GMSO Toplogy
apply(top=gmso_topology, forcefields=gmso_oplsaa, identify_connections=True) # Use the apply method from GMSO

This raises this error message:

ParameterizationError: No parameters found for connection <Angle Angle, id: 140190552781472> , group: angle_type, identifiers: [['opls_544', 'opls_544', 'opls_136'], ['CS', 'CS', 'CT']] in the Forcefield.

You can make it so missing types are ignored, so you can proceed without an error, then using the API to begin inspecting what might be missing:

gmso_oplsaa = ForceField("oplsaa") # load the oplsaa from the installed version of foyer into a GMSO forcefield
comp_1 = mb.load('[H]c1[nH]c([H])c2c1CCC2', smiles=True) # Create the combound in mBuild
gmso_topology = from_mbuild(comp_1) # Convert mBuild Compound into a GMSO Toplogy
apply(top=gmso_topology, forcefields=gmso_oplsaa, identify_connections=True, ignore_params=["angles", "improper"], remove_untyped=False) # Use the apply method from GMSO

Iterate through the angles types, you can see some are None

for angle in gmso_topology.angles:
    print(angle.angle_type)

Look at the atom types for any angle types that are None

for angle in gmso_topology.angles:
    if angle.angle_type is None:
        atom_types = [a.atom_type.name for a in angle.connection_members]
        atom_classes = [a.atom_type.atomclass for a in angle.connection_members]
        print(atom_types)
        print(atom_classes)
        print()


['opls_544', 'opls_544', 'opls_136']
['CS', 'CS', 'CT']

['opls_544', 'opls_544', 'opls_136']
['CS', 'CS', 'CT']

Hello, Sir! @chrisjonesBSU
Another question is how to modify the xml file according to the missing information？ What should be my reference? Is there any convenient way？

Intuitively, the oplsaa.xml file in foyer should be able to handle most of the molecular structure information, and my molecular structure is not weird. So it should be recognizable with great probability. But there is no recognition.
Thus, I suspect the problem that the angle cannot be defined may not be the oplsaa.xml file, but it may be the problem when mbuild builds the structure or some other problems. Because when I built this repeating unit structure into a molecular chain, there was another problem that the xml file could not identify the first atom.
The codes I ran are showing below.
import gmso
import numpy as np
import mbuild as mb
from gmso.core.forcefield import ForceField
from gmso.external import from_mbuild
from gmso.parameterization import apply
from mbuild.lib.recipes.polymer import Polymer

comp_1 = mb.load('[H]c1[nH]c([H])c2c1CCC2', smiles=True)
chain.add_monomer(compound=comp_1,
indices=[10, 8],
separation=.15,
replace=True)
chain.add_end_groups(mb.load('[H]c1[nH]c([H])c2c1CCC2',smiles=True),
index=10,
separation=0.15, label="head", duplicate=False)
chain.add_end_groups(mb.load('[H]c1[nH]c([H])c2c1CCC2',smiles=True),
index=8,
separation=0.18, label="tail", duplicate=False)
chain.build(n=20, sequence='A')
chain.visualize(show_ports=True)

gmso_oplsaa = ForceField("my_oplsaa.xml")
gmso_topology = from_mbuild(chain)
apply(top=gmso_topology, forcefields=gmso_oplsaa, identify_connections=True)
Here is the ERROR
FoyerError: Found no types for atom numbered 0 which is atomic number 6. Forcefield file is missing this atomtype, so try to add SMARTS definitions to account for this atom.

How can I figure the problem out and correct it?
THANKS for your time and patience!!!

[iGulitch]

@laowang-os , the error FoyerError: Found no types for atom numbered 0 which is atomic number 6. Forcefield file is missing this atomtype, so try to add SMARTS definitions to account for this atom. is self-explanatory : in your "my_oplsaa.xml" file, an atom type is missing for a particular atom. Basically, it means there is no line like :
<Type name="<name>" class="<class>" element="<element>" mass="<mass>" def="<SMARTS definition>">
in the <AtomTypes> section, where SMARTS sequence "def" satisfies the chemical definition of your atom numbered 0 which is atomic number 6. Atomic number 6 means carbon, number 0 means an atom labeled with an index 0.
To fix the issue, you need to write this missing line manually into your "my_oplsaa.xml" file, but prior to that your have to figure out the SMARTS definition of the problematic atom. Creating of a SMARTS definition is a whole different story.

[laowang-os]

@laowang-os , the error FoyerError: Found no types for atom numbered 0 which is atomic number 6. Forcefield file is missing this atomtype, so try to add SMARTS definitions to account for this atom. is self-explanatory : in your "my_oplsaa.xml" file, an atom type is missing for a particular atom. Basically, it means there is no line like : <Type name="<name>" class="<class>" element="<element>" mass="<mass>" def="<SMARTS definition>"> in the <AtomTypes> section, where SMARTS sequence "def" satisfies the chemical definition of your atom numbered 0 which is atomic number 6. Atomic number 6 means carbon, number 0 means an atom labeled with an index 0. To fix the issue, you need to write this missing line manually into your "my_oplsaa.xml" file, but prior to that your have to figure out the SMARTS definition of the problematic atom. Creating of a SMARTS definition is a whole different story.

Thanks!!! How to write this missing line manually into my "my_oplsaa.xml" file conveniently and quickly? Would you like to share your experience or your way to handle the missing info with me? I would appreciate it!!!

[chrisjonesBSU]

@laowang-os you'll have to learn a bit about SMARTS strings.

I ran your code, and changed the particle with index 0 (the one indicated by the error message) to show up as pink. So, you can see the carbon atom that does not have a definition in your .xml file.

Like @iGulitch said, the SMARTS definitions are found in the atom types section.

Here is an example of opls_135 where the def= field defines A carbon with 4 bonds (X4), and that one of these bonds is another carbon, and the other 3 are hydrogens.

  <Type name="opls_135" class="CT" element="C" mass="12.01100" def="[C;X4](C)(H)(H)H" desc="alkane CH3" doi="10.1021/ja9621760"/>

SMARTS definitions get a lot more complicated than this example, but, you can look through the XML to see how some of these things are handled. For your atom type, you should think about a SMARTS def for a carbon that is in a 5 membered ring, has 3 bonds (2 other carbons and 1 nitrogen). You might have to indicate that the carbon atoms it is bonded to are also part of a 5 membered ring.

The current oplsaa.xml file doesn't have the definition for all opls atom types, or SMARTS definitions for all possible chemical environments.

[laowang-os]

@laowang-os you'll have to learn a bit about SMARTS strings.

I ran your code, and changed the particle with index 0 (the one indicated by the error message) to show up as pink. So, you can see the carbon atom that does not have a definition in your .xml file.

Like [@iGulitch](https://github.com/iGulitch) said, the SMARTS definitions are found in the atom types section.

Here is an example of opls_135 where the def= field defines A carbon with 4 bonds (X4), and that one of these bonds is another carbon, and the other 3 are hydrogens.

  <Type name="opls_135" class="CT" element="C" mass="12.01100" def="[C;X4](C)(H)(H)H" desc="alkane CH3" doi="10.1021/ja9621760"/>

SMARTS definitions get a lot more complicated than this example, but, you can look through the XML to see how some of these things are handled. For your atom type, you should think about a SMARTS def for a carbon that is in a 5 membered ring, has 3 bonds (2 other carbons and 1 nitrogen). You might have to indicate that the carbon atoms it is bonded to are also part of a 5 membered ring.

The current oplsaa.xml file doesn't have the definition for all opls atom types, or SMARTS definitions for all possible chemical environments.

@chrisjonesBSU Sir, Thank you very much !!! And I'd like to accept your suggestions!
According to your suggestions now I can know the atoms for a missing angle, dihedral,improper. But for angles-like type how can I know the information related to it except atom classes like force constant and equilibrium angle in the case of I didn't found the ralated information from the ´gromacs/top/oplsaa.ff´?
​My sincere apologies for taking up your time. As a beginner, I've been quite troubled by these questions. Thank you for your time and patience！

[CalCraven]

Hi @laowang-os,

Thanks for this discussion and for getting into the nitty gritty of how to handle missing parameters. The reality is that almost everyone who has spent time using MoSDeF (or any automated forcefielding tool for that matter) has run into some of these not so subtle issues.

Here are some points I wanted to throw into the mix.

1. We have plans to add a forcefielding workshop module to GMSO, such that users that are looking to modify or aggregate their own custom forcefields have some guidance on identifying missing attributes, and writing proper SMARTS strings that do not interfere with those already in their existing forcefields, whether that be OPLSAA or other forcefields they have on hand. For now, the best guiding example of this procedures can be found in a workshop we hosted last year at the FOMMS Conference.
2. For writing SMARTS, I highly recommend using this website to check the logic for any given string https://smarts.plus/. However, not all SMARTS logic is identifiable by Foyer, so keep that in mind. Here's a link of some of the syntax we're missing.
3. While your molecule looks chemically reasonable, we all have to be aware that chemical space is vast. So while we might expect to be able to simulate a given molecule, using pre-packaged forcefields rarely works out of the box. Forcefielding is more of an art than a science. The grim reality is that there's always a trade-off for selecting any given parameter set, and each forcefield is application specific. OPLS, for instance, is always evolving and adding new parameters, but is by no means exhaustive. It's up to the user to determine if extrapolating a given parameter to their molecule of interest is reasonable. In this vein, IT'S IMPERATIVE that a user goes to the original publication for their given parameters, and verifies that the chemical context is reasonable for their use case. That's a core tenet of TRUE simulations, where the original sources for any used parameters must be accessible. There is a list of common forcefields available at this site for anyone to look at, and see exactly what molecules those parameters were generated from.

Please keep reaching out about forcefields you're interested in using, as the developer team is always willing to provide feedback and advice covering some of these subtleties and how to. best support them in MoSDeF. Especially forcefields that have unique potential forms or have increased complexity.

[CalCraven]

But for angles-like type how can I know the information related to it except atom classes like force constant and equilibrium angle in the case of I didn't found the ralated information from the ´gromacs/top/oplsaa.ff´?
​My sincere apologies for taking up your time. As a beginner, I've been quite troubled by these questions. Thank you for your time and patience！

I would follow a three step process generally:

1. Identify the missing atomtype
2. Write proper SMARTS for the atomtype, and get correct parameters from either modifying a chemically similar atom, or from another similar (preferably OPLS) forcefield.
3. If you got the atom_type from an OPLS forcefield, it should come with an associated atom_class. If not, again, you have to use your judgement by looking through the literature at Professor Jorgensen's work and find something similar.
4. Try to atomtype it in GMSO to see if the all of the dihedral, angle, and bond_types exists in your forcefield. This paper has a more recent set of parameters that may be more extensive than the Foyer OPLS, that is converted from GROMACS.
5. Run a set of control simulations with known quantities using your new parameters to validate that they are reasonable for your use case.

[laowang-os]

Hi @laowang-os,

Thanks for this discussion and for getting into the nitty gritty of how to handle missing parameters. The reality is that almost everyone who has spent time using MoSDeF (or any automated forcefielding tool for that matter) has run into some of these not so subtle issues.

Here are some points I wanted to throw into the mix.

1. We have plans to add a forcefielding workshop module to GMSO, such that users that are looking to modify or aggregate their own custom forcefields have some guidance on identifying missing attributes, and writing proper SMARTS strings that do not interfere with those already in their existing forcefields, whether that be OPLSAA or other forcefields they have on hand. For now, the best guiding example of this procedures can be found in a workshop we hosted last year at the FOMMS Conference.
2. For writing SMARTS, I highly recommend using this website to check the logic for any given string https://smarts.plus/. However, not all SMARTS logic is identifiable by Foyer, so keep that in mind. Here's a link of some of the syntax we're missing.
3. While your molecule looks chemically reasonable, we all have to be aware that chemical space is vast. So while we might expect to be able to simulate a given molecule, using pre-packaged forcefields rarely works out of the box. Forcefielding is more of an art than a science. The grim reality is that there's always a trade-off for selecting any given parameter set, and each forcefield is application specific. OPLS, for instance, is always evolving and adding new parameters, but is by no means exhaustive. It's up to the user to determine if extrapolating a given parameter to their molecule of interest is reasonable. In this vein, IT'S IMPERATIVE that a user goes to the original publication for their given parameters, and verifies that the chemical context is reasonable for their use case. That's a core tenet of TRUE simulations, where the original sources for any used parameters must be accessible. There is a list of common forcefields available at this site for anyone to look at, and see exactly what molecules those parameters were generated from.

Please keep reaching out about forcefields you're interested in using, as the developer team is always willing to provide feedback and advice covering some of these subtleties and how to. best support them in MoSDeF. Especially forcefields that have unique potential forms or have increased complexity.

Thank you for your wonderful explanation!!!😊That really helps me a lot!!!

[laowang-os]

But for angles-like type how can I know the information related to it except atom classes like force constant and equilibrium angle in the case of I didn't found the ralated information from the ´gromacs/top/oplsaa.ff´?
​My sincere apologies for taking up your time. As a beginner, I've been quite troubled by these questions. Thank you for your time and patience！

I would follow a three step process generally:

1. Identify the missing atomtype
2. Write proper SMARTS for the atomtype, and get correct parameters from either modifying a chemically similar atom, or from another similar (preferably OPLS) forcefield.
3. If you got the atom_type from an OPLS forcefield, it should come with an associated atom_class. If not, again, you have to use your judgement by looking through the literature at Professor Jorgensen's work and find something similar.
4. Try to atomtype it in GMSO to see if the all of the dihedral, angle, and bond_types exists in your forcefield. This paper has a more recent set of parameters that may be more extensive than the Foyer OPLS, that is converted from GROMACS.
5. Run a set of control simulations with known quantities using your new parameters to validate that they are reasonable for your use case.

Hi @CalCraven ! Thanks!!! But SORRY, this paper you mentioned in your explanation cannot be opened, may I ask for a help for updating this link? I would appreciate it!!! 😊

[CalCraven]

Try this link for the doi, sorry about that.
10.1021/acs.jpcb.3c06602