mosdef-hub / foyer Goto Github PK

This will require meshing something like MDTraj's element singletons for InterMol

Current inputs can come from mbuild, parmed or openmm topologies + an xml file.

By default Forcefield.apply() yields a parmed structure but for users interested in running in OpenMM we can also easily break out earlier. This needs to be tidied up a bit and properly documented.

• .xml schema to catch typos in attributes/tags
• run SMARTS parsing on every atomtype to catch syntactic errors
• ensure ampersands only exist in SMARTS strings, then run something like this

An (optional) additional field could be added in the AtomTypes (or NonbondedForce) section of the force field XML file that includes a DOI reference for where the parameters for that atom came from.

This would look something like:
<Type name="opls_140" class="HC" element="H" mass="1.00800" def="[H][C]" desc="alkane H" doi=“10.1021/ja9621760”/>

During atomtyping, foyer could store the DOI's for each atom in the molecule and write out a text file that includes all papers that should be referenced.

This could also be extended to bonded force parameters as well.

To make it easier and reduce the chance of mistakes, we should make sure we explicitly spell out the format/functional form and units of the parameters.

Specific things to note:
distance is in nm
energy is in kJ/mol
angles are radians
harmonic bonds are defined with a factor of 0.5 in front.
Include the simple conversions from opls-style fourier series to RB.

openmm does have a good page documenting some of this; we should link to it, but also make a summary in the docs folder.

http://docs.openmm.org/7.1.0/userguide/theory.html

Within a force field XML file, torsions (and I'm pretty sure bonds and angles as well) can be defined using either "class" notation, e.g.
<Proper class1="CT" class2="CT" class3="CT" class4="CT" c0="2.9288" c1="-1.4644" c2="0.2092" c3="-1.6736" c4="0.0" c5="0.0"/>
or "type" notation, e.g.
<Proper type1="opls_961" type2="opls_136" type3="opls_136" type4="opls_136" c0="-0.987424" c1="0.08363" c2="-0.08368" c3="-0.401664" c4="1.389088" c5="0.0"/>

The practical use of having both types of torsions defined would be to have the more general case be covered using the "class" notation and defining more specific dihedrals using the "type" notation. It seems to me then that dihedrals using "type" notation should essentially override those defined using the "class" notation; however, the current functionality simply assigns to a dihedral whichever parameters are listed first in the XML. The dihedral parameters are actually assigned using OpenMM, so I'm not sure if this precedence is something they would want or not, but is something I think we should at least implement in Foyer.

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

We currently do not support all of SMARTS' features. Here we're keeping track of which portions are supported.

Left column denotes current foyer support. See here for implemented SMARTS grammar.

SMARTS Atomic Primitives

SMARTS Bond Primitives

SMARTS Logical Operators

Atomtyping a (neutral) molecule to have non-zero total charge should fail atomtyping tests.

E.g. I have a working gromacs system, how do I get into using this

As far as I can tell, there's no simple way to set cross-interactions that differ from what a combining rule would specify. In other words, I can't manually set the nonbonded interaction between atomtype i and atomtype j, which I would like to do for certain custom forcefields. I'm not sure exactly how openmm handles it (if it does) but parmed appears to

Rebuild the noarch conda package with conda 4.3 compliant settings:

build:
noarch: python

Documentation: https://www.continuum.io/blog/developer-blog/condas-new-noarch-packages

First I want to mention an issue that we found regarding generating psf file, where atom mass was not assign properly.

In addition, I would appreciate if you could add the capability to foyer forcefield reader to read
an extra parameter in nonbonded section (repulsion exponent of LJ potential) from .xml file and generate charmm style parameter file with repulsion exponent. The difference between GOMC style and standard charmm is that all energy units are in (K), and in nonbonded section, there is additional variable for n (repulsion exponent).

I attached a zip file including:
1- Mie potential implemented in xml
2- Example of GOMC charmm style parameter file
3- my python script to generate pdb and psf files for leaner alkane.
UA_atom_FF.zip

Are we able to define SMARTS for atoms attached to a hexa-coordinated atom with the SMARTS functionality we currently have? The compound of interest in this case is PF6. If all angles within the molecule are to be considered, then three fluorine atom types will be necessary. All Fx-P-Fx angles should be 180 degrees and all Fx-P-Fy angles should be 90 degrees. However, because the atom types are dependent on orientation rather than bonding, I'm not sure we're able to assign atom types with our current SMARTS functionality.

One possible solution to this involves adding the NOT operator "!" into our SMARTS logic. I'm then envisioning an XML that would look something like this:

<Type name="P" class="P" element="P"  mass="30.9738"  def="P" desc="phosphorous"/>
<Type name="F1" class="F1" element="F"  mass="18.9984"  def="FP([!%F1])([!%F1])([!%F1])([!%F1])F" desc="fluorine type 1"/>
<Type name="F2" class="F2" element="F"  mass="18.9984"  def="FP([!%F1;!%F2])([!%F1;!%F2])(F1)(F1)F" desc="fluorine type 2" overrides="F1"/>
<Type name="F3" class="F3" element="F"  mass="18.9984"  def="FP(F1)(F1)(F2)(F2)F" desc="fluorine type 3" overrides="F1,F2"/>

This should atom type two fluorines each as F1, F2, and F3. The geometry wouldn't be guaranteed, but a simple energy minimization would fix that.

I'm working on using Foyer to generate a forcefield file for PCBM. Since I'm using an all atom model, I want to use the OPLSAA forcefield file built into foyer:

pdb_path = get_fn("/Users/mikehenry/Downloads/PCBM.pdb")
untyped_pcbm = pmd.load_file(pdb_path, structure=True)
oplsaa = Forcefield(name='oplsaa')
pcbm = oplsaa.apply(untyped_pcbm)

I get an error when I apply the forcefield:
FoyerError: Found multiple types for atom 2 (carbon): ['opls_139', 'opls_713']

In order to fix this I need to do some overrides, but I'm not sure if there is a way to do that on oplsaa just with code, or if I need to modify the oplsaa.xml file.

Currently it's only supported for atomtypes

Some users may already have an atomtyped system that just needs a force field applied. Others may just want to atom type something without applying the whole force field.

It should be possible to say NeighborsExactly([141, 142, 143], 1) which is evaluated as "neighbors exactly one atom of type 141, 142 or 143."

Similarly for NeighborsAtLeast and NeighborsAtMost.

There may also be a cleaner, more natural syntax to express this but some function like these should exist.

In the ethane.py example it is possible to save .top and .gro files. It would be nice to save hoomd_xml, gsd and lammpsdata files as it is done in mbuild.

When using foyer, is there a paper we should site? Or just a URL to the github repo?

• generate blacklist graph
• check for unconnected graphs
• check for graphs without any non-recursive defs

I'm using foyer '0.2.3' and mbuld 3072c6cf8c1.

import mbuild as mb
p3ht = mb.load('p3ht.pdb')
p3ht.save("p3ht.hoomdxml", overwrite=True, forcefield_name="oplsaa")

This works (stills fails on cabon 7 since I haven't added it yet to the force field) but in order to add the parameters I need to use my own force field file, so I copied the oplsaa.xml file and renamed it myoplsaa.xml.

import mbuild as mb
p3ht = mb.load('p3ht.pdb')
p3ht.save("p3ht.hoomdxml", overwrite=True, forcefield_files="oplsaa.xml")

This now gives me a different error:

/home/mike/miniconda3/envs/p3ht-example/lib/python3.5/site-packages/plyplus/engine_ply.py in t_error(t)
     20     @staticmethod
     21     def t_error(t):
---> 22         raise TokenizeError("Illegal character in input: '%s', line: %s, %s" % (t.value[:32], t.lineno, t.type))
     23 
     24     p_error = NotImplemented

TokenizeError: Illegal character in input: 'r5]1[C;X3;r5][C;X3;r5][C;X3;r5][', line: 1, error

I deleted the offending lines in myoplsaa.xml until I got the same "FoyerError: Found no types for atom 7 (carbon)." error.

Here is a link to the diff of the files:
http://pastebin.com/zUfLRpwN

Which is a list of every line it threw up on. I'm not sure what they all have in common/why oplsaa built in force field is parsed differently than a user supplied one.

p3ht.pdb.txt

As of now it appears that all the improper dihedraltype datat is only saved into the defines attribute dictionary.

The current implementation of foyer requires an element name for all particles. This is inferred from the particle's name, such that if a particle has a name that does not match an existing element, an error occurs when element = elem.Element.getBySymbol(mb_particle.name) is executed. Would there be a way to generalize foyer such that any name could be provided for the "element" attribute? For example, if I had a box of Lennard-Jones particles, is there a way I could set their "element" to "LJ", or set the element of united-atom alkanes to "CH2" and "CH3"?

Maps rules to element types so we only run rules for the given element.

If parameters are not found for a bond/angle/torsion in a molecular system, then a warning should be issued to the user. This currently does not happen and can lead to either ambiguous errors when attempting to save atom-typed structures or (worse) data files that do not contain the expected number of bonds/angles/torsions.

The workflow we've been using with OPLS doesn't work for most forcefields. In Amber, for example, non-bonded VdW and bonded interactions are specified through atomtyping but partial charges are assigned after atomtyping using antechamber, which uses AM1-BCC. In other words, a given atomtype may have different partial charges if it's used in different molecules. (One such example is n-chloroalkanes, in which the chlorines have different charges than hydrogens but only one atomtype is used for the carbon, which can result in non-zero total charges on the molecule).

The major limitation that comes to mind is that AM1-BCC is only to be used on small molecules. A large system would need to be broken down into residues, molecules, or characteristic moieties, passed to antechamber or an equivalent, and then replicated back to the whole system.

Provide an option (possibly as an optional argument to forcefield.apply()) to write an XML file containing only the information that was used in atom typing a particular system. This could also maybe be handled via a structure.write_xml() function. This would provide an easy-to-read format for checking to make sure the atomtyping was performed correctly.

Right now we are limited in the functional forms that ParmEd supports, which are very few. Since we use the OpenMM format for our force fields, that should provide an avenue for supporting arbitrary functional forms, but it will be tricky to get this information stored into a ParmEd structure...

To better understand how atomtyping can be sped up, we should create some benchmarks. Some variables that would be worth exploring

• System size: scaling from a few hundred or thousand to hundreds of thousands or millions of atoms

• Force field size, i.e. an xml file that only contains tip3p water vs. the entire OPLS force field

• Sparse vs. dense bond networks, i.e. water vs a slab of hexa-coordinated atoms in a crystal

Might be better to store this code elsewhere & I can contribute some systems

If I have three coarse-grained atom types with "elements" _CH3, _CH2, and _CH, and a SMARTS string that looks like [_CH3;X1], the SMARTS parser fails. The parser will recognize _CH as a valid element and then yield an error because it doesn't know what to do with the 3. Since _CH3 is not a valid SMARTS if _CH is assigned as the element, there may be some way we can alter the parser such that if an invalid character is detected after an element, the element is reassigned.

Fixed by updating the writing of the files

Different levels of verbosity during atomtyping will help during debugging

We need to start making some tutorials.

I'll start with something where we build from Methane -> Ethane -> longer chain alkanes.

I'm working on a different way to apply a force field to a mBuild compound.

from foyer import Forcefield
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[$]')

oplsaa = Forcefield(name='oplsaa')
ch2 = oplsaa.apply(CH2(), output_refs=True)

  File "/home/mike/miniconda3/envs/p3ht-example/lib/python3.5/site-packages/simtk/openmm/app/topology.py", line 241, in setPeriodicBoxVectors
    raise ValueError("Periodic box vectors must be in reduced form.");

If you adjust the the particle locations to something like
self.add(mb.Particle(name='C', pos=[0,0.1,0.1]), label='C[$]')

We get back to the error in issue #104
So I think we are hitting an edge case since we have a 1 dimensional bounding box. I'm not sure the best way to fix it.

Intermol system's combination rule attribute doesn't seem to ever be given a value.

OpenMM XML files can have a section for assigning residues to particles. Right now I think foyer is ignoring this section, but it would be useful to add support to assign residues based on what's specified in the section if one is specified

For example this is making a comparison with all_patterns which only contains combinations of patterns based on element names.

Modify Foyer to accept a specific .itp file and specific rules.py file (also make it so you can easily define rules in a notebook and pass them to foyer).

In the vein of supporting custom forcefields and modifications to existing forcefields, it would be useful to be able to override bond/angle/dihedral definitions in ITP files where bond/angle/dihedral types can be determined from atom names as well as bond_types.

For example, a generic dihedral between four carbon atoms is currently identified by:
CT-CT-CT-CT
where all carbons feature a bond_type of CT. However, if one wanted to specify a custom dihedral in the case where the first carbon was of type opls_792, it would be useful to be able to identify the dihedral by:
opls_792-CT-CT-CT
rather than the currently supported method where the bond_type of opls_792 would need to be changed.