Quickstart example fails

[iGulitch]

Quickstart example fails / is wrong / outdated :

mol = parmed.load("ethane.pdb") - not load, but load_file

Moreover, when using provided ethane.pdb molecule from ./docs/source/files/, got the error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/igor/miniforge3/envs/foyer/lib/python3.12/site-packages/parmed/formats/registry.py", line 194, in load_file
    return cls.parse(filename, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/igor/miniforge3/envs/foyer/lib/python3.12/site-packages/parmed/formats/pdb.py", line 432, in parse
    inst._parse_open_file(fileobj)
  File "/home/igor/miniforge3/envs/foyer/lib/python3.12/site-packages/parmed/formats/pdb.py", line 495, in _parse_open_file
    method_dispatch[rec](line)
  File "/home/igor/miniforge3/envs/foyer/lib/python3.12/site-packages/parmed/formats/pdb.py", line 778, in _end_model
    raise PDBError(f'Coordinate mismatch in model {self._current_model_number}')
parmed.exceptions.PDBError: Coordinate mismatch in model 

I tried benzene ring and styrene molecule instead, and it worked.

ff = foyer.forcefield.load_OPLSAA() - no load_OPLSAA attiribute in module foyer.forcefield, but there is one in foyer.forcefields . Trying out the latter ...

mol_ff = ff.apply(mol)
mol_ff.save("benzene.gro")
mol_ff.save("benzene.top")

worked for benzene. However, for styrene, the line mol_ff = ff.apply(mol) gave :

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/tmp/foyer/foyer/forcefield.py", line 787, in apply
    structure = self.parametrize_system(
                ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/tmp/foyer/foyer/forcefield.py", line 899, in parametrize_system
    _check_dihedrals(
  File "/tmp/foyer/foyer/forcefield.py", line 468, in _check_dihedrals
    _error_or_warn(assert_dihedral_params, msg)
  File "/tmp/foyer/foyer/forcefield.py", line 348, in _error_or_warn
    raise Exception(msg)
Exception: Parameters have not been assigned to all proper dihedrals. Total system dihedrals: 32, Parameterized dihedrals: 30. Note that if your system contains torsions of Ryckaert-Bellemans functional form, all of these torsions are processed as propers.

Thus, no further save of *.gro and *.top files is possible.

Please, help to understand how to make Foyer working.

P.S. I followed the official Installation workflow via git, via conda install and even tried out your Docker container. The results are the same.

[daico007]

Hi, sorry again about the outdated info, we are working on overhauling the docs and supply more up to date examples

To load up the OPLSAA force field:

from foyer import Forcefield

oplsaa = foyer.Forcefield(name="oplsaa")

Regarding the ethane example, that file used to work before, but there might have been changes our dependency's reader that causing the error. You can try to create the compound the smiles string (with mbuild) instead to test that out

import mbuild 
ethane = mbuild.load("CC", smiles=True) 

typed_ethane = oplsaa.apply(ethane)

For the styrene molecule, it appears that the current oplsaa force field XML is missing two dihedral parameters (I am guessing the part C=C sticking out from the aromatic ring). You can either find the parameters from the literature and add those in the XML file (foyer/forcefields/xml/oplsaa.xml under the RBTorsion section). If you're just using them as an example, and try to see that file being written out with the two missing dihedrals, you can turn out the dihedral parameters check:

styrene = mbuild.load("C=Cc1ccccc1,", smiles=True)

typed_styrene = oplsaa.apply(styrene, assert_dihedral_params=False)

[iGulitch]

Hello again @daico007 and thanks for your reply!

I'm looking forward for you to update the Foyer docs, since it's not entirely clear how to work with it except for digging into the source code. For instance, if you didn't provide me with instructions, I would not know how to apply OPLS-AA in Foyer correctly.

I've tried what you suggested and got some more warnings [ working with your latest Foyer container with Python 3.9.19 ] . To be precise, since I don't have a superuser rights on the Linux-driven machine that I'm working with, I create a Singularity container using your Foyer image from Docker Hub as a base one. This requires to manually install foyer module via conda install inside the Singularity container after it's built, otherwise Python command import foyer inside the container leads to mistake ModuleNotFoundError: No module named 'foyer'. Luckily, parmed module is installed automatically at conda install foyer call. Is it the case also for Docker containers?

Given your fixes, I believe that Quickstart example should now be looking like :

import foyer
import parmed
from foyer import Forcefield

mol = parmed.load_file("<name>.pdb")
FF = foyer.Forcefield(name = "oplsaa")

mol_FF = FF.apply(mol)

mol_FF.save("<name>.gro")
mol_FF.save("<name>.top")

Does it look correct?
FF = foyer.Forcefield(name = "oplsaa") leads to the warning : /opt/conda/lib/python3.9/site-packages/foyer/validator.py:165: ValidationWarning: You have empty smart definition(s) warn("You have empty smart definition(s)", ValidationWarning) . Is it OK?
mol_FF = FF.apply(mol) leads to :

/opt/conda/lib/python3.9/site-packages/foyer/forcefield.py:349: UserWarning: Parameters have not been assigned to all impropers. Total system impropers: 6, Parameterized impropers: 0. Note that if your system contains torsions of Ryckaert-Bellemans functional form, all of these torsions are processed as propers
  warnings.warn(msg)

Is it also OK? Nonetheless, both *.top and *.gro files are successfully generated. Thanks! However, warnings worry me, I mean whether the generated files are reliable and useable. Need to check them first.

Once, I also got some error related to from foyer import Forcefield command. I don't know how to reproduce it as I can't remember how I made it to appear [ and scrolling up the terminal is useless ] , but the error message was referring to __init__.py file, in which I found :

from foyer.forcefield import Forcefield
from foyer.forcefields import forcefields

Does it make any difference, which one to call inside the container :
from foyer import Forcefield
from foyer.forcefield import Forcefield
from foyer.forcefield import Forcefield
?

Regarding the ethane example, that file used to work before, but there might have been changes our dependency's reader that causing the error. You can try to create the compound the smiles string (with mbuild) instead to test that out

import mbuild 
ethane = mbuild.load("CC", smiles=True) 
typed_ethane = oplsaa.apply(ethane)

import mbuild requires to manually install mbuild module via conda install inside the Singularity container after it's built, otherwise there is the mistake ModuleNotFoundError: No module named 'mbulid'. Another thing, is that the file ethane.pdb created out of the SMILE sequence doesn't contain the connectivity info, i.e. there are no lines starting with CONNECT. Do you think it might affect the final *.gro and *.top files?

For the styrene molecule, it appears that the current oplsaa force field XML is missing two dihedral parameters (I am guessing the part C=C sticking out from the aromatic ring). You can either find the parameters from the literature and add those in the XML file (foyer/forcefields/xml/oplsaa.xml under the RBTorsion section). If you're just using them as an example, and try to see that file being written out with the two missing dihedrals, you can turn out the dihedral parameters check:

styrene = mbuild.load("C=Cc1ccccc1,", smiles=True)
typed_styrene = oplsaa.apply(styrene, assert_dihedral_params=False)

Given your comment, I build the following work-around example:

import foyer, parmed, mbuild
from foyer import Forcefield

mol = mbuild.load("C=Cc1ccccc1", smiles = True)
FF = foyer.Forcefield(name = "oplsaa")

mol_FF = FF.apply(mol, assert_dihedral_params = False)

which led to the following warnings :

/opt/conda/lib/python3.9/site-packages/foyer/forcefield.py:349: UserWarning: Parameters have not been assigned to all proper dihedrals. Total system dihedrals: 32, Parameterized dihedrals: 30. Note that if your system contains torsions of Ryckaert-Bellemans functional form, all of these torsions are processed as propers.
  warnings.warn(msg)
/opt/conda/lib/python3.9/site-packages/foyer/forcefield.py:349: UserWarning: Parameters have not been assigned to all impropers. Total system impropers: 8, Parameterized impropers: 0. Note that if your system contains torsions of Ryckaert-Bellemans functional form, all of these torsions are processed as propers
  warnings.warn(msg)
/opt/conda/lib/python3.9/site-packages/foyer/forcefield.py:926: UserWarning: Parametrized structure has non-zero charge.Structure's total charge: -0.1150000000000002
  warnings.warn(

Nonetheless, it seemed to work, i.e. mol_FF.save() did the job for both *.gro and *.top files. I've had a brief look at the resulting files, but haven't check them precisely. However, the same warning as above appear after mol_FF = FF.apply(mol, assert_dihedral_params = False) even in that case when I use mol = parmed.load_file("PDBs/styrene.pdb") instead of mol = mbuild.load("C=Cc1ccccc1", smiles = True). I.e. the styrene file might not be corrupted, how do you think? Finally, in foyer/opls_validation/styrene/, I found styrene.gro and styrene.top . Can I find styrene.pdb file those two were generated out of or did you use SMILE sequence, can you remember?

Regarding editing foyer/forcefields/xml/oplsaa.xml file, I gotta ask help from my more experienced colleagues. Moreover, do you suggest me to create a pull request in this GitHub repo, so that line with the two missing FF parms is added to the current OPLS-AA xml file after the check or just to edit my own fork?

[iGulitch]

Hi again, @daico007 !

I got a general Q. Do I understand it correctly that if one applies OPLS-AA forcefield to some complicated molecule and the error is :

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/opt/conda/lib/python3.9/site-packages/foyer/forcefield.py", line 780, in apply
    typemap = self.run_atomtyping(
  File "/opt/conda/lib/python3.9/site-packages/foyer/forcefield.py", line 862, in run_atomtyping
    typemap = find_atomtypes(structure, forcefield=self)
  File "/opt/conda/lib/python3.9/site-packages/foyer/atomtyper.py", line 149, in find_atomtypes
    _resolve_atomtypes(topology_graph, typemap)
  File "/opt/conda/lib/python3.9/site-packages/foyer/atomtyper.py", line 228, in _resolve_atomtypes
    raise FoyerError(
foyer.exceptions.FoyerError: Found no types for atom 4 (6)

then, it means that there is no corresponding info in foyer/forcefields/xml/oplsaa.xml about some atoms?

I got the same error message for both *.pdb file and the molecule created based on the corresponding smiles sequence, having performed the following :

>>> import foyer, parmed, mbuild
>>> from foyer import Forcefield

next command was either :

>>> mol = parmed.load_file(<pdb file>)

or

>>> mol = mbuild.load(<smiles sequence>, smiles = True)

then, FF :

>>> ff = foyer.Forcefield(name = "oplsaa")
/opt/conda/lib/python3.9/site-packages/foyer/validator.py:165: ValidationWarning: You have empty smart definition(s)
  warn("You have empty smart definition(s)", ValidationWarning)

and it's application either via :

>>> mol_ff = ff.apply(mol)

or :

>>> mol_ff = ff.apply(mol, assert_dihedral_params = False)

As a result, I always got the error like written above. The molecule isn't that large and complicated, it's 6PPD basically, i.e. C18--H24--N2.