FoyerError: Found no types for atom 0 (carbon).

[mikemhenry]

I'm working on a 'hello world' example for using foyer + mbuild + hoomd to run a hoomd simulation. I ran into an issue when trying to create a forcefield file. I created a minimal example below:

import mbuild as mb

class CH2(mb.Compound):
    def __init__(self):
        super(CH2, self).__init__()
        self.add(mb.Particle(name='C', pos=[0,0,0]), label='C[$]')
        self.add(mb.Particle(name='H', pos=[-0.109, 0, 0.0]), label='HC[$]')    
        self.add(mb.Particle(name='H', pos=[0.109, 0, 0.0]), label='HC[$]')
ch2 = CH2()
ch2.save("poly_gen.hoomdxml", overwrite=True, forcefield_name="oplsaa")

Traceback (most recent call last):
  File "test.py", line 13, in <module>
    ch2.save("poly_gen.hoomdxml", overwrite=True, forcefield_name="oplsaa")
  File "/home/mike/Projects/mbuild/mbuild/compound.py", line 1300, in save
    structure = ff.apply(structure)
  File "/home/mike/Projects/foyer/foyer/forcefield.py", line 239, in apply
    system = self.createSystem(topology, *args, **kwargs)
  File "/home/mike/Projects/foyer/foyer/forcefield.py", line 301, in createSystem
    find_atomtypes(atoms=list(topology.atoms()), forcefield=self)
  File "/home/mike/Projects/foyer/foyer/atomtyper.py", line 42, in find_atomtypes
    _resolve_atomtypes(atoms)
  File "/home/mike/Projects/foyer/foyer/atomtyper.py", line 95, in _resolve_atomtypes
    n, atom.element.name))
foyer.exceptions.FoyerError: Found no types for atom 0 (carbon).

[mikemhenry]

I tried checking out a few older releases and got the same error so it is not an obvious regression.

[summeraz]

The CH2 defined in your class is under-coordinated which is why Foyer isn't finding any atom types. The alkane CH2 in the OPLS force field in Foyer is defined as being bonded to both two hydrogens and two carbons.

[mikemhenry]

So the issue is that the atoms are not bonded?

[mikemhenry]

So this is the actual problem I'm working with:

class P3HT_Mon(mb.Compound):
    def __init__(self):
        super(P3HT_Mon, self).__init__()

        mb.load('p3ht.pdb', compound=self)
        # I used VMD to figure out what index to remove
        self.remove([self[8],
                 self[26],
                 self[27],
                 self[28],
                 self[12],
                 self[30],
                 self[31],
                 self[32]])
        # Now we need to add the ports
        # Add ports anchored to the correct carbons
        self.add(mb.Port(anchor=self[9]), label='up')
        self.add(mb.Port(anchor=self[7]), label='down')
        
        # Move the ports approximately half a C-C bond length away from the carbon
        mb.translate(self['up'], [0, 0.154/2, 0]) 
        mb.translate(self['down'], [0, -0.154/2, 0]) 
        

p3ht_mon = P3HT_Mon()
#p3ht_mon.visualize(show_ports=True)

p3ht_mon.save("p3ht_mon.hoomdxml", overwrite=True, forcefield_name="oplsaa")

(I had to add .txt to the pdb file)
If I save p3ht_mon as a hoomdxml without applying the force field the file looks good. Bonds are defined and hoomd can load it. When I try and apply the force field I get this error:

FoyerError: Found no types for atom 7 (carbon).
p3ht.pdb.txt

[summeraz]

Oh right, I hadn't even noticed the molecule didn't have bonds. So then right, Foyer will look for a carbon without bonds and since there aren't any carbons of that type defined in the OPLS force field you will get that error.

In the case of your P3HT_Mon class, SMARTS for atoms in thiophenes are not yet included in the OPLS force field XML in Foyer (we've been adding SMARTS for various atom types as we go, but the force field remains incomplete...). A quick glance at the atom type descriptions for OPLS shipped with Gromacs suggests those parameters might not exist in OPLS, but I'm not positive about that. If you have/find those OPLS parameters, adding SMARTS to the OPLS force field in Foyer would likely be the best option.

[mikemhenry]

Thanks! I'm sure I will be back with some more questions. I can track down what we use for thiophenes and see if I can add the to the force field correctly.

[mikemhenry]

I just want to make sure I'm on the right track,
I found this website that will help type atoms using OPLSAA
http://erg.biophys.msu.ru/tpp/

It created this file
output000.itp.txt
which has some opls_t* types for the thiophene.

So do I just need to add missing opls types and the bond, angle, dihedral parameters?

I'll have to add the SMARTS stuff as well so foyer can type things correctly

[summeraz]

Ah okay, so those atom types are not included in the OPLS force field in Foyer (which features the same atom types as the force field shipped with Gromacs). This is a similar issue to what I have run into with the OPLS parameters for silica, in that the parameters do exist but are not part of the standard OPLS force field. I'm not sure if we want to add in these additional parameters (i.e. those for silica and thiophenes) to the OPLS force field in Foyer or just have these be cases where the user creates their own force field XML file. @ctk3b any ideas on this?

[ctk3b]

If we do add them, we should definitely add the appropriate DOIs for where the parameters came from. Depending on where they come from we can discuss adding them to the main force field or not.

[mikemhenry]

Is there an automated way of adding them? Even if it doesn't go in the standard oplsaa force field that shipped with foyer, do you guys have some tool to help generate the openMM.xml files?

…

On Thu, Mar 23, 2017 at 11:26 AM, Christoph Klein ***@***.***> wrote: If we do add them, we should definitely add the appropriate DOIs for where the parameters came from. Depending on where they come from we can discuss adding them to the main force field or not. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#104 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ALOI3iu9vmQFm7eswQAL_DEKoY7fForyks5roqs9gaJpZM4MmBka> .

[ctk3b]

The .xml files were generated with ParmEd. Something along the lines of the following should work:

import parmed as pmd
from parmed.openmm.parameters import OpenMMParameterSet

oplsaa = pmd.load_file('my_gromacs_forcefield.itp')
xml = OpenMMParameterSet.from_parameterset(oplsaa.parameterset)
xml.write('foo.xml')

The caveat here is that RB-torsions aren't being written right now - we need to commit our modification to ParmEd that fixes this.

[summeraz]

There were a few edits we had to make to the parameters.py file in the parmed/openmm directory in the parmed source code to get this to work. I've attached below the edited parameters.py file we used as well as a file showing the changes that were made.

parameters.py.txt
git-diff.txt

[mikemhenry]

Depending how long it takes for those edits to get incorporated it might be worth forking parmed and just using git submodules

[mikemhenry]

I feel like I'm close but still missing something.
If I use the code posted by @ctk3b with this itp file, I get an error:
withimporpoutput000.itp.txt

ValueError: invalid literal for int() with base 10: 'improper_Z_CA_X_Y'

Since I don't care about impropers, I deleted them creating this file:
output000.itp.txt

Which when I used the code above worked, except this is the contents of foo.xml

<ForceField>
 <Info>
  <DateGenerated>2017-03-23</DateGenerated>
 </Info>
</ForceField>

I did use the file provided by @summeraz to patch ParmEd. While I'm sure that I could edit the oplsaa.xml file by hand, I'd rather exhaust automatic methods first.

[ctk3b]

Do you have a corresponding .gro or other structure file that would be used with this?

[mikemhenry]

Yes, I've got the PDB and TRP file: Input PDB file: input000.pdb <http://erg.biophys.msu.ru/erg/tpp/tppmktop.py?polyp3htpd_AEYQeh#> Output ITP file: output000.itp <http://erg.biophys.msu.ru/erg/tpp/tppmktop.py?polyp3htpd_AEYQeh#> Output RTP file: output000.rtp <http://erg.biophys.msu.ru/erg/tpp/tppmktop.py?polyp3htpd_AEYQeh#>

…

On Thu, Mar 23, 2017 at 2:39 PM, Christoph Klein ***@***.***> wrote: Do you have a corresponding .gro or other structure file that would be used with this? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#104 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ALOI3ip6LMdD6bZCXOxhTcONkT6kLhFeks5rotiYgaJpZM4MmBka> .

[ctk3b]

Looks like the links broke. Can you attach just the .pdb? What I think the problem is is that the .itp file produced by this service defines a [ moleculetype ]. In order for parmed to know what to do with it, it needs a concrete instance of that moleculetype. So you need a .gro file and a [ system ] and [ molecules ] section in the .itp file.

So far when converting force fields we had always been converting files that defined [ atomtypes ], [ bondtypes ] etc. which parmed reads and creates types for.

[mikemhenry]

here is the pdb file:
p3ht.pdb.txt

[ctk3b]

Do you have a full, runnable set of gromacs input files? I've created an overarching .top file which should help but it's missing all of the opls_t* types which I haven't seen before.

output.top:

#include "output.itp"
#include "oplsaa.ff/forcefield.itp"

[ defaults ]
1 3 yes 0.5 0.5

[ system ]
p3ht

[ molecules ]
RES 1

[mikemhenry]

I don't have have a gromacs file. What I'm working with is a pdb file for P3HT that I've loaded into mbuild. Since foyer didn't have those oplsaa types for the thiophenes I went looking for something that could help type them and generate force field parameters.

[ctk3b]

Right I was just wondering if that service you were trying out generated a full set of gromacs input files (I just tried it out and the answer appears to be no). It was unclear to me at this point where some of the parameter values are defined and apparently the answer is here but I can't immediately see if those parameters were derived by this group or are from some paper. Just from reading the comments, most of them are just furan parameters but I don't know about opls_t566.

[mikemhenry]

I'm not sure where they come from which is an issue for publication, but my primary goal is a 'hello world' example that we can use to show foyer and mbuild in action. As a lab we simulate a lot of organic compounds with sulfur so we will need to contribute some oplsaa parameters/citations that add those parmenters. Right now though I'm more interested in a path of least resistance that will let me have foyer generate a forcefeild for P3HT. It doesn't have to be precise, just not cause a simulation to blow up.

…

On Mar 23, 2017 15:48, "Christoph Klein" ***@***.***> wrote: Right I was just wondering if that service you were trying out generated a full set of gromacs input files (I just tried it out and the answer appears to be no). It was unclear to me at this point where some of the parameter values are defined and apparently the answer is here <https://bitbucket.org/comcon1/oplsaa-erg_ff/src/22632ba8b7c7f0c9932ea337efc9577bdaf0cf94/ffoplsaa-add-nb.itp?at=default&fileviewer=file-view-default> but I can't immediately see if those parameters were derived by this group or are from some paper. Just from reading the comments, most of them are just furan parameters but I don't know about opls_t566. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#104 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ALOI3s3CGBul3fJIWGKeBJSWJ3T4av6sks5roujDgaJpZM4MmBka> .

[ctk3b]

Right now though I'm more interested in a path of least resistance that
will let me have foyer generate a forcefeild for P3HT

If opls_t566 is the only custom atomtype and all of the bonded parameters can just come from the furan types (which I would be surprised by), then I would just manually add that atomtype to a copy of our oplsaa.xml. Otherwise to avoid doing too much manual copy pasting:

• download the ERG group's gromacs formatted OPLS files
• convert them to .xml using parmed with our patch
• add the relevant SMARTS strings for the atom types you need
  • you can probably copy almost all of these from our oplsaa.xml
• use your build/save script as is with your customized .xml file

I'm happy to help smooth out any pain points that you encounter here as we're in the process of working on this type of documentation #102

Longer term, I would like to collect the proper publication history of the various parameters before making any decisions about including them with foyer or as a separate "opls extension".

[mikemhenry]

I'm on board with the long term plan! I think curating the foyer oplsaa parameters is the right way to go. In the readme or some wiki, people could link their own forcefield parameter files that might not be general or documented well enough to make it into the foyer mainline but could be useful for others.

Going back to the issue at hand, there 8 opls_t* atom types that the service suggests using so I'll go the route of taking the gromacs file. I might not be able to touch this until next week but I'll tweak the example I'm working on to approach it as a user who has a gromacs + pdb file and they want to use mbuild to initialize their system and foyer to make the forcefield. That way it's useful for our lab and can help with #102 for documentation.

[mikemhenry]

So I was was able to create a usable forcefield.xml file by feeding this file into parmed:

#define _FF_OPLS
#define _FF_OPLSAA

; This force field uses a format that requires Gromacs 3.1.4 or later.
;
; References for the OPLS-AA force field: 
;
; W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives,
; J. Am. Chem. Soc. 118, 11225-11236 (1996).
; W. L. Jorgensen and N. A. McDonald, Theochem 424, 145-155 (1998).
; W. L. Jorgensen and N. A. McDonald, J. Phys. Chem. B 102, 8049-8059 (1998).
; R. C. Rizzo and W. L. Jorgensen, J. Am. Chem. Soc. 121, 4827-4836 (1999).
; M. L. Price, D. Ostrovsky, and W. L. Jorgensen, J. Comp. Chem. (2001).
; E. K. Watkins and W. L. Jorgensen, J. Phys. Chem. A 105, 4118-4125 (2001).
; G. A. Kaminski, R.A. Friesner, J.Tirado-Rives and W.L. Jorgensen, J. Phys. Chem. B 105, 6474 (2001).
;
[ defaults ]
; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
1		3		yes		0.5	0.5

#include "ffnonbonded.itp"
#include "lps-nb.itp"
#include "lighthavestingcofactors-nb.itp"
#include "ffoplsaa-add-nb.itp"
#include "ffbonded.itp"
#include "ffoplsaa-add-bon.itp"
#include "gbsa.itp"
#include "lighthavestingcofactors-bon.itp"

I'll then pull out the parameters I need and merge them into the oplsaa parameters from foyer (to take advantage of the atom typing work that's already done).

[mikemhenry]

Thanks everyone for helping me out with this issue. I was able to get things working by just adding 2 new atom types.