Merge pull request #124 from MobleyLab/master

Utilize ParmEd to simplify GROMACS workflows and parameter conversion; fix ParmEd 2.0 compatibility

devtools/conda-recipe/meta.yaml

@@ -20,6 +20,8 @@ requirements:
     - pandas    
     - openmm
     - ambermini
+    - pytables
+    - parmed
 #    - rdkit    # rdkit is an optional dependency, may want to comment this out for the release version.
   run:
     - python
@@ -31,6 +33,8 @@ requirements:
     - scipy
     - openmm
     - ambermini
+    - pytables
+    - parmed
 #    - rdkit    # rdkit is an optional dependency, may want to comment this out for the release version.
 
 test:

openmoltools/tests/test_gromacs.py

@@ -31,7 +31,7 @@ def test_gromacs_merge():
         utils.convert_via_acpype( "benzene", prmtop_filename2, crd_filename2, out_top = top_filename2, out_gro = gro_filename2 )
 
         #Merge topologies
-        gromacs.merge_topologies( [ top_filename1, top_filename2], './combined.top', 'combined', molecule_numbers = [1, 5] )
+        gromacs.merge_topologies( [ top_filename1, top_filename2], './combined.top', 'combined', molecule_numbers = [1, 5], molecule_names = ['etoh', 'benzene'] )
 
         #Test editing of molecule numbers in topology file
         gromacs.change_molecules_section( './combined.top', './edited.top', ['etoh', 'benzene'], [10, 20] )

openmoltools/tests/test_openeye.py

@@ -28,7 +28,7 @@
     HAVE_OE = False
 
 try:
-    import chemistry
+    import parmed 
     HAVE_PARMED = True
 except ImportError:
     HAVE_PARMED = False

openmoltools/tests/test_utils.py

@@ -7,7 +7,10 @@
 import simtk.unit as u
 from simtk.openmm import app
 import simtk.openmm as mm
+import simtk.openmm.openmm as mmmm
 from distutils.spawn import find_executable
+import parmed
+
 
 HAVE_RDKIT = True
 try:
@@ -50,6 +53,44 @@ def test_acpype_conversion():
         prmtop, inpcrd = utils.run_tleap(molecule_name, gaff_mol2_filename, frcmod_filename)
         out_top, out_gro = utils.convert_via_acpype( molecule_name, prmtop, inpcrd ) 
 
+def test_parmed_conversion():
+    molecule_name = 'sustiva'
+    input_filename = utils.get_data_filename("chemicals/sustiva/sustiva.mol2")
+    with utils.enter_temp_directory(): # Prevents creating tons of GAFF files everywhere.
+        #Make sure conversion runs
+        gaff_mol2_filename, frcmod_filename = utils.run_antechamber(molecule_name, input_filename, charge_method=None)
+        prmtop, inpcrd = utils.run_tleap(molecule_name, gaff_mol2_filename, frcmod_filename)
+        out_top, out_gro = utils.amber_to_gromacs( molecule_name, prmtop, inpcrd, precision = 8 ) 
+
+        #Test energies before and after conversion
+        #Set up amber system
+        a = parmed.amber.AmberParm( prmtop, inpcrd )
+        ambersys = a.createSystem()
+        ambercon = mmmm.Context( ambersys, mm.VerletIntegrator(0.001))
+        ambercon.setPositions( a.positions )
+        #Set up GROMACS system
+        g = parmed.load_file( out_top )
+        gro = parmed.gromacs.GromacsGroFile.parse( out_gro ) 
+        g.box = gro.box
+        g.positions = gro.positions
+        gromacssys = g.createSystem()
+        gromacscon = mmmm.Context( gromacssys, mm.VerletIntegrator(0.001))
+        gromacscon.setPositions( g.positions ) 
+
+        #Check energies
+        a_energies = parmed.openmm.utils.energy_decomposition( a, ambercon )    
+        g_energies = parmed.openmm.utils.energy_decomposition( g, gromacscon )
+        #Check components
+        tolerance = 1e-5
+        ok = True
+        for key in a_energies.keys():
+            diff = np.abs(a_energies[key] - g_energies[key] )
+            if diff/np.abs(a_energies[key]) > tolerance:
+                ok = False
+                print("In testing AMBER to GROMACS conversion, %s energy differs by %.5g, which is more than a fraction %.2g of the total, so conversion appears not to be working properly." % ( key, diff, tolerance) )
+        if not ok:
+            raise(ValueError("AMBER to GROMACS conversion yields energies which are too different.")) 
+
 
 @skipIf(SKIP_CHECKMOL, "Skipping testing of checkmol descriptors since checkmol is not found (under that name)." )
 def test_checkmol_descriptors():

openmoltools/utils.py

@@ -78,8 +78,14 @@ def convert_via_acpype( molecule_name, in_prmtop, in_crd, out_top = None, out_gr
         GROMACS topology file produced by acpype 
     out_gro : str
         GROMACS coordinate file produced by acpype
+
+    Notes
+    -----
+    Deprecated. Please use ParmEd (especially amber_to_gromacs) instead.
     """
 
+    print("WARNING: Use of acpype for conversion is deprecated. ParmEd is preferred; please use amber_to_gromacs instead.") 
+
     #Create output file names if needed
     if out_top is None:
         out_top = "%s.top" % molecule_name        
@@ -406,13 +412,13 @@ def randomize_mol2_residue_names(mol2_filenames):
     """Find unique residue names for a list of MOL2 files.  Then
     re-write the MOL2 files using ParmEd with the unique identifiers.
     """
-    import chemistry    
+    import parmed   
     names = get_unique_names(len(mol2_filenames))
 
     for k, filename in enumerate(mol2_filenames):
-        struct = chemistry.load_file(filename)
+        struct = parmed.load_file(filename)
         struct.name = names[k]
-        mol2file = chemistry.formats.Mol2File
+        mol2file = parmed.formats.Mol2File
         mol2file.write(struct, filename)
 
 def get_checkmol_descriptors( molecule_filename, executable_name = 'checkmol' ):
@@ -479,3 +485,63 @@ def get_checkmol_descriptors( molecule_filename, executable_name = 'checkmol' ):
 
     return descriptors
 
+def amber_to_gromacs( molecule_name, in_prmtop, in_crd, out_top = None, out_gro = None, precision = None): 
+    """Use ParmEd to convert AMBER prmtop and crd files to GROMACS format.
+
+    Requires
+    --------
+    Currently requires ParmEd v2.0 beta1 or later.
+
+
+    Parameters
+    ----------
+    molecule_name : str
+        String specifying name of molecule
+    in_prmtop : str
+        String specifying path to AMBER-format parameter/topology (parmtop) file
+    in_crd : str
+        String specifying path to AMBER-format coordinate file
+    out_top : str, optional, default = None
+        String specifying path to GROMACS-format topology file which will be written out. If none is provided, created based on molecule_name.
+    out_gro : str, optional, default = None
+        String specifying path to GROMACS-format coordinate (.gro) file which will be written out. If none is provided, created based on molecule_name.
+    precision : int, optional, default = None
+        If not none, set the precision of the coordinates in the written .gro file to the specified number of decimal places.
+
+    Returns
+    -------
+    out_top : str
+        GROMACS topology file produced by ParmEd
+    out_gro : str
+        GROMACS coordinate file produced by ParmEd
+
+    Notes
+    -----
+        molecule_name is not currently used except to generate output file names if gro/top file names are not provided. It is an argument partly for API consistency.
+    """
+    #Create output file names if needed
+    if out_top is None:
+        out_top = "%s.top" % molecule_name        
+    if out_gro is None:
+        out_gro = "%s.gro" % molecule_name
+
+    #Check precision
+    if precision is not None:
+        assert isinstance(precision, int), "Precision %s is not an integer." % precision
+
+    #Import ParmEd
+    import parmed
+
+    #Read AMBER to ParmEd object
+    structure = parmed.amber.AmberParm( in_prmtop, in_crd )
+    #Make GROMACS topology
+    gromacs_topology = parmed.gromacs.GromacsTopologyFile.from_structure( structure )
+    #Write
+    parmed.gromacs.GromacsTopologyFile.write( gromacs_topology, out_top )
+    if precision == None:    
+        parmed.gromacs.GromacsGroFile.write( gromacs_topology, out_gro )
+    else:
+        parmed.gromacs.GromacsGroFile.write( gromacs_topology, out_gro, precision = precision )
+
+    return out_top, out_gro
+

setup.py

@@ -24,8 +24,8 @@
 
 
 ##########################
-VERSION = "0.6.8.dev0"
-ISRELEASED = False
+VERSION = "0.6.8"
+ISRELEASED = True
 __version__ = VERSION
 ##########################