Formalizing Atom-typing and the Dissemination of Force Fields with Foyer

TL;DR

Foyer is a Python tool that formalizes atom-typing rules in a machine-readable XML format, reducing ambiguity and errors in force field parameterization for molecular simulations, thereby enhancing reproducibility.

Contribution

The paper introduces Foyer, a novel Python-based framework that standardizes force field atom-typing rules using XML and SMARTS, improving reproducibility and dissemination in molecular modeling.

Findings

Foyer successfully automates atom-typing with reduced errors.

It enhances reproducibility in force field application.

Demonstrated effectiveness in multiple use-cases.

Abstract

A key component to enhancing reproducibility in the molecular simulation community is reducing ambiguity in the parameterization of molecular models. Ambiguity in molecular models often stems from the dissemination of molecular force fields in a format that is not directly usable or is ambiguously documented via a non-machine readable mechanism. Specifically, the lack of a general tool for performing automated atom-typing under the rules of a particular force field facilitates errors in model parameterization that may go unnoticed if other researchers are unable reproduce this process. Here, we present Foyer, a Python tool that enables users to define force field atom-typing rules in a format that is both machine- and human-readable thus eliminating ambiguity in atom-typing and additionally providing a framework for force field dissemination. Foyer defines force fields in an XML format, where SMARTS strings are used to define the chemical context of a particular atom type. Herein we describe the underlying methodology of the Foyer package, highlighting its advantages over typical atom-typing approaches and demonstrate is application in several use-cases.

Figures

Tables

Table 1. Table 1: Extensions to the atom type definitions in the OpenMM XML format.

Attribute	Description	Example
def	Defines the chemical context of an atom type via SMARTS	[C;X4](H)(H)(H)C
desc	Textual description of the atom type	Alkane CH3
doi	Digital object identifier to the atom type source	10.1021/ja9621760
overrides	Atom type(s) the current rule is given precedence over	opls_136

Table 2. Table 2: 2D depictions of molecular fragments referred to in the text

Alkane			Alkene			Benzene
C, CH3	C, CH2	H	C, (R2-C=)	C, (RH-C=)	H	C	H
opls_135	opls_136	opls_140	opls_141	opls_142	opls_144	opls_145	opls_146

Table 3. Table 3: Currently implemented SMARTS atomic primitives a

Symbol	Symbol name	Atomic property requirements	Default
*	wildcard	any atom	(no default)
A	aliphatic	aliphatic	(no default)
r<n>	ring size	in smallest SSSRb ring of size <n>	any ring atom
X<n>	connectivity	<n>total connections	exactly one
#n	atomic number	atomic number <n>	(no default)

Table 4. Table 4: Extensions to SMARTS atomic primitives developed in this work

Symbol	Symbol name	Atomic property requirements	Default
_A	non-element	non-atomistic element	(no default)
%<type>	atomtype	of atomtype <type>	(no default)

Table 5. Table 5: SMARTS Logical Operators a

Symbol	Expression	Meaning
exclamation	!e1	not e1
ampersand	e1&e2	e1 and e2 (high precedence)
comma	e1,e2	e1 or e2
semicolon	e1;e2	e1 and e2 (low precedence)

Full text

\alsoaffiliation

Vanderbilt Multiscale Modeling and Simulation (MuMS) Center, Vanderbilt University, Nashville, Tennessee 37235, USA

\alsoaffiliationVanderbilt Multiscale Modeling and Simulation (MuMS) Center, Vanderbilt University, Nashville, Tennessee 37235, USA

\alsoaffiliationVanderbilt Multiscale Modeling and Simulation (MuMS) Center, Vanderbilt University, Nashville, Tennessee 37235, USA

\alsoaffiliationVanderbilt Multiscale Modeling and Simulation (MuMS) Center, Vanderbilt University, Nashville, Tennessee 37235, USA

\alsoaffiliationVanderbilt Multiscale Modeling and Simulation (MuMS) Center, Vanderbilt University, Nashville, Tennessee 37235, USA \alsoaffiliationDepartment of Chemistry, Vanderbilt University, Nashville, Tennessee 37235, United States

\alsoaffiliationVanderbilt Multiscale Modeling and Simulation (MuMS) Center, Vanderbilt University, Nashville, Tennessee 37235, USA

\alsoaffiliationVanderbilt Multiscale Modeling and Simulation (MuMS) Center, Vanderbilt University, Nashville, Tennessee 37235, USA

1 Introduction

Considerable efforts have been undertaken by many research groups to develop accurate classical force fields for a wide range of systems.1, 2, 3, 4, 5, 6, 7 Force fields are often expressed as a set of analytical functions with adjustable fitting parameters that describe the interactions between constituents of a system (often discrete atoms but, more generally, interaction sites). Classical force fields are able to achieve high accuracy by creating sets of highly specific fitting parameters (i.e., atom types), in which each atom type describes an interaction site within a different chemical context. The chemical context is typically defined by the bonded environment of an interaction site (e.g., the number of bonds and the identity of the bonded neighbors) and may also consider, among other factors, the bonded environment of the neighbors, and/or the specific molecule/structure within which the interaction site is included. Consequently, a force field may include tens or even hundreds of different atom types for a given element. For example, there are 347 different types of carbon in the OPLS forcefield parameter set distributed with GROMACS8 where each type corresponds to a carbon atom within a different chemical context (i.e., a different atom type). Thus, while force field development efforts have reduced – or in some cases completely eliminated – the need for researchers to generate their own fitting parameters, determining which parameters (i.e., atom types) to use can still be a tedious and error prone task. Failure to properly identify the chemical context and atom type of an interaction site will inevitably lead to the incorrect implementation of the force field and thus inconsistent results.

Part of the difficulty in performing atom-typing (i.e., determining which atom type applies to an interaction site) stems from the fact that there is not yet a standardized way of unambiguously expressing chemical context and parameter usage. As such, journal articles that report novel force field parameters may vary significantly in terms of their clarity. In many cases, parameters are reported in a tabular format with minimal annotations and few (if any) examples of how to appropriately assign the atom types. Since this approach does not allow for automated evaluation, different users of the force field may apply the atom types differently based on their own interpretation of the information provided. Journal articles that utilize existing force fields often do not report the specific parameters used and typically do not specify which atom types were chosen for the interaction sites, instead providing citation(s) to the source of the force field parameters. Even if the source of the parameters is clearly and fully specified, usage may again depend on the clarity of the original source(s) and the interpretation by the end user, hampering reproducibility. Force field parameter files that aggregate a large number of atom types (often thousands) into a single source suffer from some of the same issues. Often, they include only brief, unstructured – and sometimes ambiguous – annotations as to parameter usage, and may, or may not, provide clear citations of the original source of the parameters.

To apply force fields, users can perform atom-typing manually (e.g., through the creation of an atom-typed template of a molecule), although manual assignment of parameters becomes tedious and error prone for large molecules and/or complex systems. Furthermore, manual manipulation of files is not considered a good practice in terms of reproducibility 9 and manual assignment of parameters does not lend itself well to workflows such as screening, where thousands of unique systems with different chemical constituents and structures may need to be atom-typed. To avoid manual assignment, end-users often develop in-house software to apply force fields in an automated fashion, however, such software is not typically made freely available to the broader community. Without access to the same software, the exact atom-typing cannot be reproduced by others and if the source code is not made freely available, the logic used to interpret and apply the force field is unknown. Similarly, if there are errors in the software/logic, these cannot be identified. There exist a number of freely available atom-typing tools that read in a force field parameter file and execute a set of rules to apply the force field to a chemical topology 10, 11, 12, 13, 14, 15, 16, enabling the exact atom-typing process to be reproduced. However, many of these atom-typing tools are either closed-source11, 14, simulation engine-specific13, 16, and/or force field-specific12, 14, 15, which limits their utility. Furthermore, these tools almost universally rely on a rigid hierarchy of rules10, where rules must be called in a precise order such that more general atom types are only chosen when more specialized matches do not exist ($i.e.$, the order of rules defines the precedence). Maintaining, let alone constructing, these hierarchies is challenging, especially for a large number of atom types. In order to add a new atom type or correct an error in hierarchical schemes, a developer must have a complete picture of the hierarchy and know exactly where the relevant rule should be placed such that it does not inadvertently override other rules. This may impose practical limits on functionality, where, for example, a user is not able to easily extend the rules to include new atom types, or that such attempts to extend the rules result in incorrect atom-typing for other systems. For many tools, this approach is further complicated by the encoding of the hierarchy as a set of heavily nested if/else statements within the source code of the software. These hierarchies may be difficult to validate and debug, and any changes or extensions to the rules, no matter how trivial, require modification of the source code itself. Reproducibility issues may therefore arise if users make modifications or extensions to a piece of software and these changes are not made freely available to the larger community and/or incorporated into the main software distribution. This also creates a situation where there are effectively two sets of rules since there is no guarantee that the logic statements in the source code (i.e., the machine readable rules) agree with the textual annotations in the force field parameter file (i.e., the human readable rules).

Several atom-typing tools have been developed that remove the need to encode atom type usage rules within the source code itself. A unifying feature of these tools is the use of the simplified molecular-input line-entry system (SMILES 17), a language for describing chemical structures, or variants thereof, to define the chemical context of an atom type. For example, Yesselman, et al. 16 developed an atom-typing toolset for the CHARMM simulation engine, termed MATCH, that relies on assigning parameters by representing a molecule of interest as a graph and performing subgraph matching against a library of fragments with known parameters. These fragments are represented as “super smiles”, an extension of the SMILES language. By using super smiles and storing these fragments in text files separate from the software, chemical context is expressed without the need to define a rigid if/else hierarchy within the software and thus new rules (i.e., fragments encoded as super smiles) can be added without modifying the code used to evaluate them, although, the order of these fragments will still determine the final atom-type. In other work, the Enhanced Monte Carlo (EMC) software developed by in’t Veld 18 encodes chemical context of an atom type using SMILES, allowing usage rules to be defined outside of the source code, but EMC also relies on a hierarchical approach to set precedence. In recent work, Mobley and coworkers19 have developed the SMIRNOFF specification that effectively eliminates explicit atom types altogether, instead using SMIRKS (another language related to SMILES 20) to identify chemical fragments that are associated with a set of force field parameters, including bonded interactions. Similar to Yesselman, et al. 16, SMIRNOFF relies on representing the system as a graph and rules as subgraphs; in this approach rules must also be defined in a specific order to appropriately define precedence. In all cases, the use of the SMILES-based approaches not only removes the need to encode usage within the source code, but associates parameters with a human and machine readable definition of their chemical context. However, these approaches still require rules be specified in the particular order to enable correct atom-typing and since the SMILES language is only specified for atomistic system, such approaches may have limited compatibility with, e.g., coarse-grained force fields.

In this work, we present Foyer21, a Python library for performing atom-typing based upon first-order logic over graph structures, designed to address many of the aforementioned issues associated with encoding chemical context, performing atom-typing of both atomistic and non-atomistic systems, and disseminating force fields to the community. Foyer relies upon a force-field-agnostic formalism to express atom-typing and parameterization rules in a way that is expressive enough to be human readable while simultaneously being machine readable, allowing a single, unambiguous format to be constructed, suitable for use both the software and for dissemination. This logic is implemented via SMARTS20 to encode chemical context and “overrides” statements to define rule precedence. SMARTS is an extension of the SMILES language that supports substructure definitions and allows expression of greater chemical detail and logic operations within the chemical patterns. By using SMARTS to define chemical context, atom type definitions are encoded outside of the source code, and thus force fields can be created and evolved without modification to the code used to evaluate them. Foyer supports both fully atomistic and non-atomistic force fields, such as coase-grained and united-atom force fields, via an extension of SMARTS that enables user-defined “elements” (not in the periodic table) to be leveraged within the chemical context definitions. Similar to Refs. 16, 19, Foyer treats the system to be atom-typed as a graph and atom types are determined by identifying which subgraphs (specified via SMARTS) match the environment of a node in the graph (i.e., the chemical context of an interaction site). In Foyer, the method of resolving atom types differs from most tools in that it is iterative and does not rely on rigid rule hierarchies. Rule precedence is explicitly defined by the aforementioned overrides statements, thus atom-typing rules can appear in any order in the file, providing increased flexibility and eliminating physical placement in the file as a source of error. Since this iterative approach evaluates all rules (not just the first to evaluate to ”True”), automated evaluation can be used to help ensure that Foyer force field definitions (1) encompass all atom types in the force field and (2) are sufficiently descriptive without conflicting rules, both necessary conditions for publishing force fields in a way that is unambiguous and reproducible. Foyer uses an XML-formatted force file format to define force field parameters and how to apply them to chemical systems. Rather than creating a new file format, Foyer builds upon the OpenMM force field file format, extended to support the definition of the associated SMARTS definitions, overrides statements, and textual descriptions of the parameters; this format is also extended to support the definition of digital objective identifiers (DOIs) for each atom type, such that the source of the parameters can be clearly identified. By using XML, this force field file format is flexible, extensible, and by the nature of XML itself, provides clear labels/metadata of the contents in the file, improving human readability and therefore reducing ambiguity. The Foyer software provides routines that create syntactically correct input files for a variety of common simulation engines and is designed to take, as input, chemical topologies from several of other community developed tools (e.g., ParmEd22, OpenMM23, and mBuild24, 25, 26). We demonstrate that, by combining this software and annotation scheme with version control, force fields can be created and evolved in a transparent manner, allowing force field parameters and their usage to be more easily disseminated to the community. Furthermore, since the force field format is both human and machine readable, force fields can be automatically validated to ensure completeness, allowing force field parameters and their usage to be published with reduced ambiguity, helping to improve reproducibility of molecular simulations.

2 Defining chemical context and rule precedence

2.1 XML File Format

Foyer utilizes the OpenMM23 force field XML format to encode parameters, where this format is extended to allow for the definitions of chemical context and rule precedence (discussed below). To briefly summarize the OpenMM file format, atom types and forces are encoded as XML tags with various attributes defining the types of elements that they apply to (by name only), as well as the associated parameters for that interaction (e.g., the equilibrium bond length and spring constant for a harmonic bond). Listing 1 provides an example of encoding the OPLS force field parameters for linear alkanes in the OpenMM XML format (note, this Listing does not include our extensions). As shown in Listing 1, the XML format provides clear descriptions of each of the parameters/properties defined in the file (e.g., element="C" indicates the entry is defining a carbon atom), along with additional tags that provide unambiguous descriptions of the types of interactions being used (e.g., the <HarmonicBondForce> tag is used to define the use of a harmonic force to define bonds). As such, this file format includes a wealth of metadata that is both human and machine readable. For more detailed information, we refer the reader to the OpenMM manual where this force field file format is extensively documented (\urlhttp://docs.openmm.org/7.0.0/userguide/application.html#writing-the-xml-file).

The flexible nature of XML allows it to be readily extended via the addition of new tags/attributes without fundamentally changing the original format, as new tags/attributes can simply be ignored by software that does not require them. As shown in Table 1, and discussed in detail later, four new attributes have been added to the atom type entries in the existing OpenMM XML file format to enable the functionality needed to encode usage rules in Foyer: def, desc, doi, and overrides. The use of XML additionally allows sanity checks to be performed by using XML schemas to ensure the expected attributes have been provided in the file.

2.2 Using SMARTS to define chemical context

The chemical context of an interaction site is typically defined by its bonded environment, notably the number of bonds and the identities of bonded neighboring interaction sites, but may also include longer range information, such as the bonded environment of neighbors. To encode this information, Foyer utilizes SMARTS27, a language for defining chemical patterns. SMARTS is an extension of the more commonly used SMILES28 notation, providing additional tokens that enable users to express greater chemical detail and logic operations. SMARTS notation is expressed as strings that simultaneously include arbitrary chemical complexity but are concise and clear enough for human consumption, in addition to being machine readable. As an example, consider defining the chemical context of OPLS-AA atom types for carbon and hydrogen atoms in a linear alkane, as shown in Listing 2 (note, only the <AtomTypes> section of the file is shown, as this is the only section that differs from Listing 1). The reader is referred to Table 2 for a visual depiction of these atom types. To encode the chemical context, the def attribute is added to the OpenMM XML format to encode the corresponding SMARTS string. Here, the atom type that specifies the terminal methyl group, opls_135 (“-CH3”) can be expressed as C;X4(H)(H)H in the SMARTS notation. In this SMARTS notation, [C;X4] indicates that the element of interest – always the first token in the SMARTS string – is a carbon atom (i.e., C) and this carbon atom has 4 total bonds (i.e., ;X4, where ; indicates the logical operator AND). The identities of the 4 bonded neighbors are 1 carbon atom and 3 hydrogen atoms, expressed as (C)(H)(H)H. Similarly, the opls_136 atom type, which describes a methylene group in an alkane, is expressed in SMARTS notation as C;X4(C)(H)H. Here, the only change from the opls_135 SMARTS definition lies in the identity of the 4 bonded neighbors (2 carbon atoms and 2 hydrogen atoms). Increased chemical complexity can be described by adding details about each of neighboring interaction sites within SMARTS. For example, the opls_140 atom type, which describes a generic alkane hydrogen, is defined as H[C;X4] - a hydrogen atom bonded to a carbon atom with 4 bonds. Multiple valid SMARTS can be defined for each atom type, where, e.g., opls_140 could be defined simply as def="H" since there is only a single hydrogen atom type defined in Listing 2, although we note that such a simplified definition would not necessarily provide a user of the force field with a clear understanding of the chemical context for which this atom type applies. The Foyer extensions to the OpenMM XML file format also include the desc attribute (e.g., shown in Listing 2) that allows for unstructured comments to be provided in addition to the SMARTS definitions.

We note that the parser in the Foyer libraries does not currently support the full SMARTS language, instead providing support for the subset that was found to be relevant to the definition of chemical context for atom types. Table 3 lists the currently supported primitives, Table 4 shows our extensions to the language, and Table 5 outlines the logical operators supported. Increased support for the SMARTS language will be provided in future software releases, as needed.

2.3 Establishing rule precedence

Rule precedence must be established when multiple atom type definitions can apply to a given interaction site. In typical hierarchical schemes, this is determined implicitly by the order in which rules are evaluated; in general, more specific rules are evaluated first and when a match is found, the code stops evaluating rules altogether. While this approach works, it becomes more challenging to maintain the correct ordering of rules as the number of atom types grows and as chemistries become more complex and specific. Users may find it difficult, if not impossible, to make even small additions to a larger force field without breaking existing behavior. Foyer allows rule precedence to be explicitly stated via the use of the overrides attribute added to the XML file format. This allows atom type usage rules to be encoded in any order within the file, eliminating incorrectly placed rule order as a source of error. Foyer iteratively evaluates all rules on all interaction sites in the system, maintaining for each interaction site a “whitelist” consisting of rules that evaluate to True and a “blacklist” consisting of rules that have been superseded by another rule (i.e., those that appear in the overrides attribute). The set difference between the white- and blacklists of an interaction site yields the correct atom type if the force field is implemented correctly (incorrect/incomplete definition of force fields is discussed later). As an example of a system where overrides need to be defined, consider describing alkenes and benzene in a single force field file, as shown in Listing 3 (note, again only the <AtomTypes> section of the force field file is shown). The reader is again referred to Table 2 for visual depictions of the relevant atom types.

When atom-typing a benzene molecule, the carbon atoms in the ring will match the SMARTS patterns for both opls_142 (an alkene carbon) and opls_145 (a benzene carbon). Without the overrides attribute, Foyer will find that multiple atom types apply to each carbon atom. Providing the overrides attribute indicates that if the opls_145 pattern matches, it will supersede opls_142. Thus, the difference between the whitelist (containing opls_142 and opls_145) and blacklist (containing only opls_142) would be opls_145.

Note that multiple atom types can be listed in a single overrides attribute, however, the approach taken here also allows atom types to inherit overrides from the atom types they override. For example, consider a case in which atom types 1, 2, and 3 each evaluate to True for an interaction site. If atom type 3 overrides atom type 2 (i.e., adds atom type 2 to the blacklist) and atom type 2 overrides atom type 1 (i.e., adds atom type 1 to the blacklist), then atom type 3 will implicitly override atom type 1.

Additionally, in Foyer, the SMARTS grammar has been modified such that specific atom type names can also be included within the definition (see Table 4). For example, opls_146, the hydrogen atom attached to carbon atoms in a benzene ring, has the SMARTS definition [H][C;%opls_145], as shown in Listing 3. This states that the interaction site of interest is a hydrogen atom (H) and is bonded to a carbon atom that has atom type opls_145 (C;%opls_145). Because Foyer evaluates rules iteratively for each interaction site, such recursive definitions can be utilized without the need to explicitly define atom types in a chemical topology input file. For example, in this case, when Foyer identifies the interaction site of a carbon atom to be opls_145, the next iteration to evaluate the hydrogen atom will find that opls_146 now evaluates to True. Similar to how an overrides statement clearly defines precedence, this recursive definition provides a clear way to identify chemical context and the relationship between different atom types for highly specific parameters. We note, that one could also replace the recursive reference to opls_145 with its SMARTS string, although in this case it would result in a more complex, less human readable definition.

Because the logic used to define chemical context is separated from the source code used to evaluate it, one can construct a force field file that contains only the relevant subset of atom types need for a given application area. Using the above example of benzene and alkenes, if a system only contained benzene molecules, one could avoid specifying the overrides attributes altogether by simply creating a force field file containing only atom types relevant to benzene and eliminating those associated with alkenes. In many cases, considering smaller subsets is beneficial as the amount of effort required to differentiate and set rule precedence between atom types is reduced. Additionally, using smaller files reduces the likelihood of errors related to defining chemical context and rule precedence, reduces the number of test molecules with known atom types required to fully validate the rules, and increases the readability of the force field files by limiting the number of entries.

2.4 Extension of SMARTS for non-atomistic systems

Foyer is able to atom-type systems in which an interaction site does not represent a single atom with a standard element, but instead may represent a group of atoms (relevant to united-atom and coarse-grained force fields) or a generic site (relevant to simplified models). Standard SMARTS notation does not support non-atomic species due to its reliance on the presence of an element specification for each interaction site. To circumvent this limitation, the Foyer SMARTS parser allows users to define custom “elements” by prefixing their string representation with an underscore (see Table 4). For example, _CCC could represent a coarse-grained interaction site intended to model three carbon atoms. In its current implementation, Foyer makes a first pass through force field files to detect any custom element definitions. These are injected into the grammar that parses SMARTS strings and are given priority over standard elements. This allows non-atomistic and atomistic atom types to be used, either separately or together.

In practice, united-atom and coarse-grained force fields can be defined in an almost identical fashion to all-atom force fields, where the only difference is that “elements” are user-defined strings prepended with an underscore. As an example, consider an alkane modeled with the united-atom TraPPE force field. An interaction site in this force field represents both carbon and the hydrogen atoms bonded to it. Thus, this force field contains two distinct atom types, one that represents $\mathrm{CH_{3}}$ (_CH3) and one that represents $\mathrm{CH_{2}}$ (_CH2). These can be encoded as shown in Listing 4. Focusing on atom type CH3_sp3, usage is encoded with the definition [_CH3;X1][_CH3,_CH2] which states that the base “element” is _CH3 with one bond (i.e., ;X1) to either a _CH3 or a _CH2 group. In SMARTS, a comma indicates an “OR” logic statement while a semicolon is used to denote an “AND” logical statement (see Table 5 for a complete list of SMARTS logical operators). In this example, [_CH3;X1] states the element must be _CH3 “AND” have only a single bond. Atom-type _CH2_sp3, which represents a “middle” alkane carbon and its two associated hydrogen atoms, is defined similarly as _CH2;X2[_CH3,_CH2]. Here, the base “element” is a _CH2 with two bonds (i.e., ;X2), each of which may be either “element” _CH3 or _CH2.

2.5 Determining bonded parameters

Once a chemical topology is atom-typed, bonded interactions can be determined by simply searching for the matching pairs, triplets, and quartets (bonds, angles, and torsions, respectively). In many force fields, the bonded parameters are not as specific as the non-bonded interactions, and thus are not defined directly based on atom types. Thus, rather than atom types, a more general class identifier (sometimes referred to as the “bond family”) is used to identify these interactions. In Listing 5, both opls_136 and opls_962 are part of the same class “CT”. Thus a bond between opls_136 and opls_962 would have the same parameters (defined as class1="CT" class2="CT" in Listing 5) as a bond between opls_136 and opls_136 (also defined as class1="CT" class2="CT"). However, this general approach breaks down for certain chemical topologies. For example, while the atom types for carbon atoms in alkanes and perfluoralkanes are both of class “CT” and share the same bond and angle parameters for carbon atoms, they differ in terms of torsional parameters. In order to handle this conflict, many codes require users to comment out the more general set of parameters or include statements within the code that accomplish the same task. However, in this approach, one would not be able to atom-type a system composed of a mixture alkane and perfluoroalkane molecules29, since only one set of parameters can be included simultaneously. Note that while one could define a new class to, for example, differentiate between alkanes and perfluoroalkanes, this would also require defining all the interactions with existing relevant class entries and result in a force field file with many duplicate parameters sets that simply have different labels. The OpenMM format allows bonded parameters to be defined using the type attribute, referring directly to the name attribute that stores the atom type, instead of the class attribute, which allows bonded interactions to be defined with increased specificity. Additionally, mixed use of type and class in the definition of these bonded interactions is supported. Referring to Listing 5, to provide the necessary distinction between torsional parameters for perfluoroalkanes and alkanes, one could define perfluoroalkane torsions using type attributes (i.e., type1="opls_962" type2="opls_962" type3="opls_962" type4="opls_962", where opls_962 is defined in the <AtomTypes> XML section), and alkane torsions with the more general quartet for alkanes of class1="CT" class2="CT" class3="CT" class4="CT" that uses class attributes.

When iterating through bonded parameter definitions, OpenMM assigns parameters based on the first match found. In the example described above, perfluoroalkane torsional parameters would therefore need to be defined before alkane parameters in the torsional section, and thus the ordering shown in Listing 5 would result in the incorrect assignment of torsional parameters for perfluoroalkanes. While overrides statements could be used to set rule precedence for bonded topologies and thus eliminate the need to specify order in the file, additional modification to the force field file format would be required because bonded parameters do not have a “name” attribute like atom types. A more simple approach taken by Foyer is to perform a preprocessing step on the bonded parameters. This step orders bonded parameters such that the most specific cases are sorted to the top of the list to set precedence. This is accomplished by assigning a weight to each entry proportional to the number of type attributes included. For example, a torsion that defines the atom types for which it applies (i.e., has 4 type attributes) would be given the highest weight, whereas an entry that specifies only class attributes would be given the lowest. Thus, for the force field XML shown in Listing 5 Foyer would reverse the order of the two defined dihedrals during preprocessing.

3 Foyer software

The Foyer software has been developed as a Python library that can read the force field XML specification discussed above and perform atom-typing. Python allows for portability between platforms and provides a wealth of freely available modules (e.g., NumPy 30, SciPy 31, NetworkX 32) to facilitate many of the underlying operations. The source, documentation, and examples of Foyer usage are all freely available and can be found on the project’s GitHub repository \urlhttps://github.com/mosdef-hub/foyer and associated website \urlhttp://mosdef-hub.github.io/foyer/. Foyer is a core package within the Molecular Simulation Design Framework (MoSDeF), a framework/protocol aimed towards fully reproducible molecular simulations, and is hosted under the mosdef-hub GitHub organization. Additionally, tutorials to Foyer’s usage can be accessed via a separate tutorial repository \urlhttps://github.com/mosdef-hub/foyer_tutorials. Figure 1 provides an overview of the general software workflow, which will be discussed here.

3.1 Inputs and Preprocessing

Foyer accepts, as input, the XML force field file and a chemical topology for which to apply the force field. In addition to sorting bonded parameters by specificity as described in the previous section, the XML force field file undergoes a preprocessing and validation step via application of an XML schema definition. This schema definition in Foyer enforces which elements (e.g. HarmonicBondForce) are valid and how their attributes should be formatted. While this does not test the accuracy of the parameters, it does ensure that all of the expected parameters are defined. Additionally, the schema ensures that atom types are not defined more than once and that atom types referenced in other sections (e.g., <HarmonicBondForce>) are actually defined in the <AtomTypes> section. Next, the SMARTS strings defined by the def attribute for all atom types are parsed and checked for validity. This does not validate whether a SMARTS string is correctly defined for a given interaction site but simply ensures that the SMARTS string can be interpreted by Foyer and does not contain any erroneous characters. Parsing errors are captured and re-raised with error messages that allow a user to pin point the location of the problem in the XML file and within the SMARTS string. Wherever possible, Foyer attempts to provide helpful suggestions for fixing detected errors.

Input chemical topologies can be passed to Foyer through various data structures; the current version supports the OpenMM Topology object, the ParmEd Structure object22, 33, and the mBuild Compound object24, 34. Each of the OpenMM, ParmEd, and mBuild topologies support inputs from a variety of common molecular file formats, such as PDB and MOL2, and thus it is typically straightforward to convert a given system into a data structure that Foyer can accept. Regardless of the input format, once read into Foyer the chemical topology is converted to an OpenMM Topology object, in the current release of the software. The OpenMM Topology object provides a standardized data container to store the necessary system information and allows for leveraging of routines already defined within OpenMM’s library.

3.2 Atom-typing

A flowchart of Foyer’s atom-typing procedure is shown in Figure. 3. To perform atom-typing, Foyer constructs a graph of the complete system defined by the chemical topology (or alternatively a graph of each unique residue, see the Residue-based Atom-typing section below) and iteratively searches for SMARTS matches via subgraph isomorphism (where subgraphs are generated for each SMARTS definition). Graph construction and matching are performed using the NetworkX package, an open-source Python project that provides an intuitive interface for a multitude of graph-based algorithms and is the de facto standard network analysis library in Python. During this step, the iterative process of determining the atom type is undertaken, adding rules to the white and back lists for each interaction site in the system.

The implementation of the SMARTS based atom-typing scheme is comprised of several steps and internally relies on a subgraph isomorphism to detect matches as highlighted in Figure 2. First, a SMARTS string is parsed into an abstract syntax tree (AST) from which we populate a SMARTSGraph object. This class inherits from the Graph class in the NetworkX package32. Elements in this SMARTSGraph are represented as nodes and chemical bonds as edges. Inheriting from NetworkX is convenient in that it allows us to leverage most of the algorithms and visualization methods already implemented there. The primary distinguishing feature of the SMARTSGraph is the set of methods that encode the logic for matching the more complex SMARTS tokens. These methods can be directly used by NetworkX’s implementation of the VF2 subgraph isomorphism algorithm35. A thin wrapper provided by the find_matches method allows a SMARTSGraph instance to search for all subgraph isomorphisms within a bare chemical topology (an non-atom-typed graph of just elements and bonds). This method returns the indices of all elements that match the first token in the SMARTS string, which defines the atom type that we are looking for. Successfully matching elements have the atom type definition added to their whitelist and any overridden types added to their blacklist. The appropriate atom type for an interaction site is determined by examining the difference between white- and blacklists, where a sufficiently descriptive force field should yield only a single atom type as the difference between the two lists.

The use of white- and blacklists provides users with a means to validate the completeness of the chemical contexts defined by the set of SMARTS strings and overrides. For example, when considering a given atom, if multiple valid atom types are found as the difference between that atom’s white- and blacklists, this indicates the rule definitions within the force field are not sufficiently unique and likely have incomplete information provided to the override attributes. Foyer provides the list of conflicting types to aid in resolving such issues. If no atom types exist as the difference, the interaction site of interest cannot be described by the force field rules as implemented. Typically, this will require adding a new atom type or amending an existing atom type’s definition. This may also indicate, that there is an error in how rule precedence has been defined, such as, all the rules on the whitelist “overriding” each other. Note, the efficacy of this type of validation in Foyer will depend on providing a sufficient range of systems to fully explore the combinations of atom types that can be applied. As a general rule, the set of systems chosen to perform validation tests should collectively utilize all atom types defined in the force field.

3.2.1 Residue-based Atom-typing

Many systems of interest to molecular simulation contain topologies that consist of duplicates of smaller molecules or repeat units, each with identical topologies. A brute-force implementation of the atom-typing process wastes time by repeating subgraph isomorphism computation on each repeat unit. Also, this approach does not scale well with system size. To eliminate unnecessary calculations, a map of atom-typed residues is saved after each unique residue is atom-typed the first time. Then, when an identical residue is found, it copies the atom-typed information from the residue map instead of repeating the subgraph isomorphism. This feature is enabled by default but can optionally be turned off.

As a example, consider a box of $N$ hexane molecules. After the subgraph isomorphism is called on the first molecule, the result is copied and saved into a map. Then, when molecules 2 to $N$ are encountered, those results are copied into the running topology. The time it takes for the apply function to finish is timed for each case and plotted in Figure 4 relative to the brute force approach; as can be seen, the use of the residue map speeds up the overall atom-typing step by more than an order of magnitude for common system sizes.

3.3 Force field assignment and output

Once atom-typed, the OpenMM Topology now contains atom types for all particles in the system and is used to create an OpenMM System object upon which bonded parameters of the system are determined. This step can be accomplished by searching the list of bonded parameters for the appropriate pair, triplet, and quartet of atom types for bonds, angles, and dihedrals, respectively; such routines exist within OpenMM and are utilized in this context, where again we note these interactions undergo a sorting to ensure more specific definitions appear first in the file. Additionally, validation checks are performed at this time to ensure all triplets and quartets of interaction sites have had angle and dihedral (proper and improper) parameters assigned (checks for bond parameterization of interaction site pairs are performed by OpenMM in the prior step). These validation checks provide the user with an error (that can optionally be overridden) to help prevent the return of incorrectly parameterized Structures. To output the atom-typed system into a usable format for a simulation engine, the fully atom-typed and parameterized system is returned as a ParmEd Structure object. Through the use of the ParmEd Structure, Foyer has access to various additional functionality, such as I/O routines that properly parse the ParmEd Structure into common chemical file formats (MOL2, PDB, NAMD and GROMACS formats, among others22). For file formats not natively supported by ParmEd, custom I/O routines for outputting to these formats (e.g., the LAMMPS data file format) have been developed within the mBuild package.

It should be noted that by utilizing ParmEd, the force fields that Foyer currently supports must match functional forms supported by the internals of the ParmEd Structure object. For example, non-bonded interactions are currently limited to a 12-6 Lennard-Jones functional form. However, due to the large amount of force fields that utilize this functional form (e.g., Amber1, GAFF36, OPLS3, TraPPE37, 38, CHARMM2), the current version of Foyer is still widely applicable; future development will include support for additional functional forms to better accommodate the diverse force field landscape that exists.

3.4 Validating output

Foyer provides scripts to validate its output files by comparing against systems with known atom types (e.g., those determined by hand or reference molecules provided by a force field developer). Output validation requires (1) system(s) with known, validated atom types and (2) the force field XML file. The known systems are read into Foyer and atom types are determined using the rules in the XML file. The atom types generated by Foyer are then compared against the known atom-typed system(s). The pytest39 library is used to provide a clear, descriptive output of the results of these validation tests. Implementing output validation tests is particularly useful to force field developers as they ensure that the desired output is retained if a force field file is evolved through the addition of new atom type definitions or merged with a separate force field file. The utility of these validation checks relies not only upon providing accurate reference systems, but also a sufficient variety of test systems that encompass all defined atom types, as discussed previously. An example of such a validation test suite is provided as part of the Foyer template repository 40.

4 Usage Examples

Here, a basic overview of the usage of Foyer is provided. Readers are also directed to the GitHub project repository 21 and tutorial repository 41 for more thorough and current usage examples. Consider constructing a bulk system of ethane molecules and applying the OPLS force field. Listing 6 shows a simple mBuild script to load an ethane molecule and fill a 2nm x 2nm x 2nm box with 100 molecules. This defines the system’s chemical topology to which the force field will be applied. As input, the force field file is identical to Listing 1 but with the <AtomType> information from Listing 2, as Listing 2 includes the usage rule definitions. Listing 6 demonstrates two different syntaxes for applying a force field using Foyer and saving the output, in this case to the file format required by GROMACS. The second option allows different forcefields to be applied to different topologies in the system. Listing 7 shows an example of creating two separate chemical topologies in the system, and applying two different force field files to each. The two atom-typed structures that result (ethane_fluid and silica_substrate) are then combined using a simple $+$ operator and saved to any format supported by ParmEd. Note that if the surface and polymers were bonded together (e.g., to create a surface-bound monolayer), the force field files would need to be combined into a single XML document.

5 Promoting Reproducible Force Field Dissemination

Foyer force field files and associated documentation, examples, and validation tests, can be readily developed and distributed using standard software development approaches to improve quality and reproducibility. For example, the common git + GitHub/Bitbucket based distribution process allows force field creators to disseminate their force field files and associated content to the public via a version controlled repository that can be referenced from relevant publications. In this approach, a specific version of the force field used in a publication can be tagged in the git repository and a reference to this tagged version provided in the manuscript, allowing for a clear reference to the exact parameters and usage rules employed in the work. Other services, such as Zenodo 42, can additionally provide a digital object identifier (DOI) for the tagged record and a snapshot of the content of the archive. A variety of other features of this standard software development process translate well to force field development. Version control systems like git are designed to facilitate distributed, collaborative software development and allow for changes to the files in the repository to be easily tracked in a transparent manner. For example, as a force field is evolved or corrected, revisions can be easily tracked, including the author(s) responsible for the changes, and the specific differences between force field versions clearly identified using standard tools such as DIFF and through the use of descriptive “commit” statements as the content of the repository is changed. Whenever the developers wish to, they can create a new release of the force field that, as noted above, can be tagged or provided with a citable DOI. Verification and validation of a force field can also be simplified by using this software design approach, i.e., by implementing automated testing tools that can perform checks on every new iteration (i.e., commit) of the force field content to ensure errors are not introduced as the force field is changed.

To promote these practices, we have created a template git repository on GitHub which contains the basic framework needed to create, test, and publish a new force field as well as a guided tutorial that introduces users to the SMARTS based atom-typing scheme 40. This process was successfully used in recent work that derived force field parameters for perfluoropolyethers43, a novel lubricant class. The force field was published in conjunction with the manuscript and made available on GitHub (\urlhttps://github.com/mosdef-hub/forcefield_perfluoroethers). The specific version of the force field at time of publication is citable via a separate DOI44. Any adjustments or improvements to the force field can be released under a new DOI, while the old one would still exist and point to the originally published force field in order to maintain provenance.

5.1 Atom type DOI labels

While automated atom-typing and the containment of atomtypes and force field parameters within a single file helps reduce user error and promotes reproducibility, users also require knowledge of the original source of parameters, in order to ensure proper citation and validation that the parameters are appropriate for their system of interest. Foyer achieves this goal by adding a doi attribute to each Type definition within the AtomTypes block of a force field XML. Listing 8 shows the same atom-type definition for the OPLS-AA methyl carbon as in Listing 2 with the additional doi attribute providing the DOI to the original source where parameters for this atom type were derived.

This feature eliminates ambiguity concerning the origin of parameters for a particular atom type. Furthermore, Foyer automatically logs associations between DOIs and atomtypes during the atom-typing process, providing a BibTeX file featuring the full citation for the sources of all parameters applied to a particular system, along with additional notes detailing precisely which atomtypes are contained within each source. For example, Listing 9 shows the BibTeX file generated for a nitropropane molecule using the OPLS-AA force field, which includes all reference information as well as notes describing which atom type parameters were obtained from each reference.

6 Conclusion

Foyer is a Python package for atom-typing and parameterization of molecular models, yielding an output directly usable by common molecular simulation engines. Foyer extends the XML force field format introduced by OpenMM with an additional field for SMARTS definitions that define the chemical context of atom types in a single file that is both human and machine readable. Foyer uses these SMARTS definitions to assign the appropriate atom types and provides an overrides syntax for defining rule precedence that is independent of the order of the rule definitions within the file. Foyer force fields can be used along with standard software development practices for version control and DOI tagging to provide a complete record of force field development and links that researchers can include within a publication leading to the exact force field parameters used for a particular study. The extension of the XML file format to additional include DOI tags for the source of each parameter, further increases the clarity of force field usage. By promoting best practices in terms of encoding parameter usage and disseminating force fields, we believe that Foyer can play a critical role in improving reproducibility in the molecular simulation community.

7 Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant No. OCI-1047828. We also acknowledge the National Energy Research Supercomputing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.