Merge pull request #38 from open-forcefield-group/smirff_energy_check

Many (but not all) fixes relating to cross-checking force terms against AMBER-originated systems

README.md

@@ -29,9 +29,10 @@ pip install -i https://pypi.anaconda.org/OpenEye/simple OpenEye-toolkits
 Install other conda dependencies:
 ```
 conda install --yes numpy networkx
+conda install --yes -c omnia openmoltools
 ```
 
-NOTE: We'll find a better way to install these dependencies via `conda` soon.
+NOTE: We'll add a better way to install these dependencies via `conda` soon.
 
 ## Installation
 

devtools/travis-ci/install.sh

@@ -18,6 +18,7 @@ export PIP_ARGS="-U"
 export PATH=$MINICONDA_HOME/bin:$PATH
 conda update --yes conda
 conda install --yes conda-build jinja2 anaconda-client pip
+conda install --yes -c omnia openmoltools
 
 # Restore original directory
 popd

examples/README.md

@@ -3,3 +3,4 @@
 ## Manifest
 * `parm@frosst/` - example illustrating attempt to recover parm@frosst atom types
 * `AlkEtOH/` - example illustrating attempt to recover parm99 atom types for a small set of alkanes, ethers, and alcohols
+* `SMIRFF_comparison/` - temporary example (waiting for permanent home) of cross-comparison of molecule energies from SMIRFF with same molecule energies from .prmtop and .crd.

examples/SMIRFF_comparison/compare_set_energies.py

@@ -0,0 +1,33 @@
+#!/bin/env python
+
+from smarty.forcefield_utils import *
+import os
+
+# Cross-check energies of molecules from AlkEthOH set using SMIRFF xml file 
+# versus energies from AMBER .prmtop and .crd files (parm@frosst params)
+
+datapath = './AlkEthOH_inputfiles/AlkEthOH_chain_filt1' 
+molname = 'AlkEthOH_c100'
+mol_filename = os.path.join( datapath, molname+'.mol2')
+
+# Check if we have this data file; if not we have to extract the archive
+if not os.path.isfile( mol_filename):  
+    print "Extracting archived molecule files."
+    tarfile = 'AlkEthOH_inputfiles.tar.gz' 
+    os.system('tar -xf AlkEthOH_inputfiles.tar.gz')
+
+# Load OEMol
+mol = oechem.OEGraphMol()
+ifs = oechem.oemolistream(mol_filename)
+flavor = oechem.OEIFlavor_Generic_Default | oechem.OEIFlavor_MOL2_Default | oechem.OEIFlavor_MOL2_Forcefield
+ifs.SetFlavor( oechem.OEFormat_MOL2, flavor)
+oechem.OEReadMolecule(ifs, mol )
+oechem.OETriposAtomNames(mol)
+
+# Load forcefield
+forcefield = ForceField(get_data_filename('forcefield/Frosst_AlkEtOH.ffxml'))
+
+# Compare energies
+prmtop = os.path.join( datapath, molname+'.top')
+crd = os.path.join( datapath, molname+'.crd')
+results = compare_molecule_energies( prmtop, crd, forcefield, mol)

setup.py

@@ -14,7 +14,7 @@ def find_package_data(data_root, package_root):
 
 setup(
     name = "smarty",
-    version = "0.1.0",
+    version = "0.1.1",
     author = "John Chodera",
     author_email = "[email protected]",
     description = ("Automated Bayesian atomtype sampling"),

smarty/__init__.py

@@ -1,4 +1,5 @@
 from .atomtyper import *
 from .forcefield import *
+from .forcefield_utils import *
 from .sampler import *
 from .utils import *

smarty/data/forcefield/Frosst_AlkEtOH.ffxml

@@ -8,46 +8,52 @@
 <!-- This mockup is based on frcmod.Frosst_AlkEthOH, the minimal subset of parm99/parm@Frosst -->
 <!-- required for alkanes, ethers, and alcohols (excluding ketals and 3-memb-rings) -->
 
-<HarmonicBondForce length_unit="angstroms" k_unit="kilocalories_per_mole/angstrom">
-   <Bond smirks="[#6X4:1]-[#6X4:2]" length="1.526" k="310.0"/> <!-- CT-CT from frcmod.Frosst_AlkEthOH -->
-   <Bond smirks="[#6X4:1]-[#1:2]" length="1.090" k="340.0"/> <!-- CT-H_ from frcmod.Frosst_AlkEthOH -->
-   <Bond smirks="[#8:1]~[#1:2]" length="1.410" k="320.0"/> <!-- DEBUG O-H -->
-   <Bond smirks="[#6X4:1]-[O&amp;X2&amp;H1:2]" length="1.410" k="320.0"/> <!-- CT-OH from frcmod.Frosst_AlkEthOH -->
-   <Bond smirks="[#6X4:1]-[O&amp;X2&amp;H0:2]" length="1.370" k="320.0"/> <!--CT-OS from frcmod.Frosst_AlkEthOH -->
-   <Bond smirks="[#8X2:1]-[#1:2]" length="0.960" k="553.0"/> <!-- OH-HO from frcmod.Frosst_AlkEthOH -->
+<!-- WARNING: AMBER functional forms drop the factor of 2 in the bond energy term, so cross-comparing this file with a corresponding .frcmod file, it will appear that the values here are twice as large as they should be. -->
+<HarmonicBondForce length_unit="angstroms" k_unit="kilocalories_per_mole/angstrom**2">
+   <Bond smirks="[#6X4:1]-[#6X4:2]" length="1.526" k="620.0"/> <!-- CT-CT from frcmod.Frosst_AlkEthOH -->
+   <Bond smirks="[#6X4:1]-[#1:2]" length="1.090" k="680.0"/> <!-- CT-H_ from frcmod.Frosst_AlkEthOH -->
+   <Bond smirks="[#8:1]~[#1:2]" length="1.410" k="640.0"/> <!-- DEBUG O-H -->
+   <Bond smirks="[#6X4:1]-[O&amp;X2&amp;H1:2]" length="1.410" k="640.0"/> <!-- CT-OH from frcmod.Frosst_AlkEthOH -->
+   <Bond smirks="[#6X4:1]-[O&amp;X2&amp;H0:2]" length="1.370" k="640.0"/> <!--CT-OS from frcmod.Frosst_AlkEthOH -->
+   <Bond smirks="[#8X2:1]-[#1:2]" length="0.960" k="1106.0"/> <!-- OH-HO from frcmod.Frosst_AlkEthOH -->
 </HarmonicBondForce>
-<HarmonicAngleForce angle_unit="degrees" k_unit="kilocalories_per_mole/radian">
-   <Angle smirks="[a,A:1]-[#6X4:2]-[a,A:3]" angle="109.50" k="50.0"/> <!-- consensus matches all X-Csp3-X -->
-   <Angle smirks="[#1:1]-[#6X4:2]-[#1:3]" angle="109.50" k="35.0"/> <!-- H1-CT-H1 from frcmod.Frosst_AlkEthOH -->
-   <Angle smirks="[#6X4:1]-[#6X4:2]-[#6X4:3]" angle="109.50" k="40.0"/> <!-- CT-CT-CT from frcmod.Frosst_AlkEthOH -->
-   <Angle smirks="[#8X2:1]-[#6X4:2]-[#8X2:3]" angle="109.50" k="70.0"/> <!-- O_-CT-O_ from frcmod.Frosst_AlkEthOH -->
-   <Angle smirks="[#6X4:1]-[#8X2:2]-[#1:3]" angle="108.50" k="55.0"/> <!-- CT-OH-HO from frcmod.Frosst_AlkEthOH -->
-   <Angle smirks="[#6X4:1]-[#8X2:2]-[#6X4:3]" angle="109.50" k="60.0"/> <!-- CT-OS-CT from frcmod.Frosst_AlkEthOH -->
+<!-- WARNING: AMBER functional forms drop the factor of 2 in the angle energy term, so cross-comparing this file with a corresponding .frcmod file, it will appear that the values here are twice as large as they should be. -->
+<HarmonicAngleForce angle_unit="degrees" k_unit="kilocalories_per_mole/radian**2">
+   <Angle smirks="[a,A:1]-[#6X4:2]-[a,A:3]" angle="109.50" k="100.0"/> <!-- consensus matches all X-Csp3-X -->
+   <Angle smirks="[#1:1]-[#6X4:2]-[#1:3]" angle="109.50" k="70.0"/> <!-- H1-CT-H1 from frcmod.Frosst_AlkEthOH -->
+   <Angle smirks="[#6X4:1]-[#6X4:2]-[#6X4:3]" angle="109.50" k="80.0"/> <!-- CT-CT-CT from frcmod.Frosst_AlkEthOH -->
+   <Angle smirks="[#8X2:1]-[#6X4:2]-[#8X2:3]" angle="109.50" k="140.0"/> <!-- O_-CT-O_ from frcmod.Frosst_AlkEthOH -->
+   <Angle smirks="[#6X4:1]-[#8X2:2]-[#1:3]" angle="108.50" k="110.0"/> <!-- CT-OH-HO from frcmod.Frosst_AlkEthOH -->
+   <Angle smirks="[#6X4:1]-[#8X2:2]-[#6X4:3]" angle="109.50" k="120.0"/> <!-- CT-OS-CT from frcmod.Frosst_AlkEthOH -->
 </HarmonicAngleForce>
 <PeriodicTorsionForce phase_unit="degrees" k_unit="kilocalories_per_mole">
    <Proper smirks="[a,A:1]-[#6X4:2]-[#6X4:3]-[a,A:4]" periodicity1="3" phase1="0.0" k1="1.40"/> <!-- X -CT-CT-X from frcmod.Frosst_AlkEthOH -->
    <Proper smirks="[a,A:1]-[#6X4:2]-[#8X2:3]-[#1:4]" periodicity1="3" phase1="0.0" k1="0.50"/> <!--X -CT-OH-X from frcmod.Frosst_AlkEthOH -->
-   <Proper smirks="[a,A:1]-[#6X4:2]-[#8X2:3]-[!#1:4]" periodicity1="3" phase1="00" k1="1.15"/> <!-- X -CT-OS-X from frcmod.Frosst_AlkEthOH -->
+   <Proper smirks="[a,A:1]-[#6X4:2]-[#8X2:3]-[!#1:4]" periodicity1="3" phase1="0.0" k1="1.15"/> <!-- X -CT-OS-X from frcmod.Frosst_AlkEthOH -->
    <Proper smirks="[#1:1]-[#6X4:2]-[#6X4:3]-[#1:4]" periodicity1="3" phase1="0.0" k1="0.15"/> <!-- HC-CT-CT-HC from frcmod.Frosst_AlkEthOH -->
    <Proper smirks="[#1:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]" periodicity1="3" phase1="0.0" k1="0.16"/> <!-- HC-CT-CT-CT from frcmod.Frosst_AlkEthOH -->
    <Proper smirks="[#6X4:1]-[#6X4:2]-[#8X2:3]-[#1:4]" periodicity1="3" phase1="0.0" k1="0.16" periodicity2="1" phase2="0.0" k2="0.25"/> <!-- HO-OH-CT-CT from frcmod.Frosst_AlkEthOH -->
    <Proper smirks="[#6X4:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]" periodicity1="3" phase1="0.0" k1="0.18" periodicity2="2" phase2="180.0" k2="0.25" periodicity3="1" phase3="180.0" k3="0.20"/> <!-- CT-CT-CT-CT from frcmod.Frosst_AlkEthOH -->
    <Proper smirks="[#6X4:1]-[#6X4:2]-[#8X2:3]-[#6X4:4]" periodicity1="3" phase1="0.0" k1="0.383" periodicity2="2" phase2="180.0" k2="0.1"/> <!-- CT-CT-OS-CT from frcmod.Frosst_AlkEthOH -->
-   <Improper smirks="[a,A:1]~[#6X3:2]([a,A:3])~[OX1:4]" periodicity1="2" phase1="180.0" k1="10.5"/> <!-- X -X -C -O  from frcmod.Frosst_AlkEthOH -->
+   <Proper smirks="[#6X4:1]-[#8X4:2]-[#6X4:3]-[O&amp;X2&amp;H0:4]" periodicity1="3" phase1="0.0" k1="0.10" periodicity2="2" phase2="180.0" k2="0.85" periodicity3="1" phase3="180.0" k3="1.35"/> <!-- CT-OS-CT-OS from frcmod.Frosst_AlkEthOH -->
+   <Proper smirks="[#8X2:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]" periodicity1="3" phase1="0.0" k1="0.144" periodicity2="2" phase2="0.0" k2="0.175"/> <!-- O_-CT-CT-O_ from frcmod.Frosst_AlkEthOH -->
+   <Proper smirks="[#8X2:1]-[#6X4:2]-[#6X4:3]-[#1:4]" periodicity1="3" phase1="0.0" k1="0.0" periodicity2="1" phase2="0.0" k2="0.25"/> <!-- O_-CT-CT-O_ from frcmod.Frosst_AlkEthOH; discrepancy with parm@frosst with H2,H3-CT-CT-O_ per C Bayly -->
+   <Improper smirks="[a,A:1]~[#6X3:2]([a,A:3])~[OX1:4]" periodicity1="2" phase1="180.0" k1="10.5"/> <!-- X -X -C -O  from frcmod.Frosst_AlkEthOH; none in set but here as format placeholder -->
 </PeriodicTorsionForce>
+<!-- WARNING: AMBER formats typically use r_0/2=r_min/2 to describe the relevant distance parameter, where r0 = 2^(1/6)*sigma. The difference is important, and the two conventions can be used here by specifying sigma or rmin_half. -->
 <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5" sigma_unit="angstroms" epsilon_unit="kilocalories_per_mole">
    <!-- sigma is in angstroms, epsilon is in kcal/mol -->
-   <Atom smirks="[#1:1]" sigma="1.4870" epsilon="0.0157"/> <!-- making HC the generic hydrogen -->
-   <Atom smirks="[$([#1]-C):1]" sigma="1.4870" epsilon="0.0157"/> <!-- HC from frcmod.Frosst_AlkEthOH -->
-   <Atom smirks="[$([#1]-C-[#7,#8,F,#16,Cl,Br]):1]" sigma="1.3870" epsilon="0.0157"/> <!-- H1 from frcmod.Frosst_AlkEthOH -->
-   <Atom smirks="[$([#1]-C(-[#7,#8,F,#16,Cl,Br])-[#7,#8,F,#16,Cl,Br]):1]" sigma="1.2870" epsilon="0.0157"/> <!--H2 from frcmod.Frosst_AlkEthOH -->
-   <Atom smirks="[$([#1]-C(-[#7,#8,F,#16,Cl,Br])(-[#7,#8,F,#16,Cl,Br])-[#7,#8,F,#16,Cl,Br]):1]" sigma="1.1870" epsilon="0.0157"/> <!--H3 from frcmod.Frosst_AlkEthOH -->
-   <Atom smirks="[#1$(*-[#8]):1]" sigma="0.0000" epsilon="0.0000"/> <!-- HO from frcmod.Frosst_AlkEthOH -->
-   <Atom smirks="[#6:1]" sigma="1.9080" epsilon="0.1094"/> <!-- making CT the generic carbon -->
-   <Atom smirks="[#6X4:1]" sigma="1.9080" epsilon="0.1094"/> <!-- CT from frcmod.Frosst_AlkEthOH-->
-   <Atom smirks="[#8:1]" sigma="1.6837" epsilon="0.1700"/> <!-- making OS the generic oxygen -->
-   <Atom smirks="[#8X2:1]" sigma="1.6837" epsilon="0.1700"/> <!-- OS from frcmod.Frosst_AlkEthOH -->
-   <Atom smirks="[#8X2+0$(*-[#1]):1]" sigma="1.7210" epsilon="0.2104"/> <!-- OH from frcmod.Frosst_AlkEthOH -->
+   <Atom smirks="[#1:1]" rmin_half="1.4870" epsilon="0.0157"/> <!-- making HC the generic hydrogen -->
+   <Atom smirks="[$([#1]-C):1]" rmin_half="1.4870" epsilon="0.0157"/> <!-- HC from frcmod.Frosst_AlkEthOH -->
+   <Atom smirks="[$([#1]-C-[#7,#8,F,#16,Cl,Br]):1]" rmin_half="1.3870" epsilon="0.0157"/> <!-- H1 from frcmod.Frosst_AlkEthOH -->
+   <Atom smirks="[$([#1]-C(-[#7,#8,F,#16,Cl,Br])-[#7,#8,F,#16,Cl,Br]):1]" rmin_half="1.2870" epsilon="0.0157"/> <!--H2 from frcmod.Frosst_AlkEthOH -->
+   <Atom smirks="[$([#1]-C(-[#7,#8,F,#16,Cl,Br])(-[#7,#8,F,#16,Cl,Br])-[#7,#8,F,#16,Cl,Br]):1]" rmin_half="1.1870" epsilon="0.0157"/> <!--H3 from frcmod.Frosst_AlkEthOH -->
+   <Atom smirks="[#1$(*-[#8]):1]" rmin_half="0.0000" epsilon="0.0000"/> <!-- HO from frcmod.Frosst_AlkEthOH -->
+   <Atom smirks="[#6:1]" rmin_half="1.9080" epsilon="0.1094"/> <!-- making CT the generic carbon -->
+   <Atom smirks="[#6X4:1]" rmin_half="1.9080" epsilon="0.1094"/> <!-- CT from frcmod.Frosst_AlkEthOH-->
+   <Atom smirks="[#8:1]" rmin_half="1.6837" epsilon="0.1700"/> <!-- making OS the generic oxygen -->
+   <Atom smirks="[#8X2:1]" rmin_half="1.6837" epsilon="0.1700"/> <!-- OS from frcmod.Frosst_AlkEthOH -->
+   <Atom smirks="[#8X2+0$(*-[#1]):1]" rmin_half="1.7210" epsilon="0.2104"/> <!-- OH from frcmod.Frosst_AlkEthOH -->
 </NonbondedForce>
 
 </SMARFF>

smarty/forcefield.py

@@ -638,6 +638,27 @@ def createForce(self, system, topology, verbose=False, **kwargs):
 
         if verbose: print('%d bonds added' % (len(bonds)))
 
+
+        # Check that no topology bonds are missing force parameters
+        atoms = [ atom for atom in topology.atoms() ]
+        topology_bonds = ValenceDict()
+        for (atom1, atom2) in topology.bonds():
+            topology_bonds[(atom1.index,atom2.index)] = True
+        if set(bonds.keys()) != set(topology_bonds.keys()):
+            msg = 'Mismatch between bonds added and topological bonds.\n'
+            created_bondset = set(bonds.keys())
+            topology_bondset = set(topology_bonds.keys())
+            msg += 'Bonds created that are not present in Topology:\n'
+            msg += str(created_bondset.difference(topology_bondset)) + '\n'
+            msg += 'Topology bonds not assigned parameters:\n'
+            for (a1, a2) in topology_bondset.difference(created_bondset):
+                atom1 = atoms[a1]
+                atom2 = atoms[a2]
+                msg += '(%8d,%8d) : %5s %3s %3s - %5s %3s %3s' % (a1, a2, atom1.residue.index, atom1.residue.name, atom1.name, atom2.residue.index, atom2.residue.name, atom2.name)
+                msg += '\n'
+            raise Exception(msg)        
+
+
 parsers["HarmonicBondForce"] = HarmonicBondGenerator.parseElement
 
 #=============================================================================================
@@ -799,8 +820,24 @@ class NonbondedGenerator(object):
     class LennardJonesType(object):
         """A SMIRFF Lennard-Jones type."""
         def __init__(self, node, parent):
+            """Currently we support radius definition via 'sigma' or 'rmin_half'."""
             self.smirks = _validateSMIRKS(node.attrib['smirks'], node=node)
-            self.sigma = _extractQuantity(node, parent, 'sigma')
+
+            # Make sure we don't have BOTH rmin_half AND sigma
+            try:
+                a = _extractQuantity(node, parent, 'sigma')
+                a = _extractQuantity(node, parent, 'rmin_half')
+                raise Exception("Error: BOTH sigma and rmin_half cannot be specified simultaneously in the .ffxml file.")
+            except:
+                pass            
+
+            #Handle sigma
+            try: 
+                self.sigma = _extractQuantity(node, parent, 'sigma')
+            #Handle rmin_half, AMBER-style
+            except:
+                rmin_half = _extractQuantity(node, parent, 'rmin_half', unit_name='sigma_unit')
+                self.sigma = 2.*rmin_half/(2.**(1./6.))
             self.epsilon = _extractQuantity(node, parent, 'epsilon')
 
     def __init__(self, forcefield, coulomb14scale, lj14scale):