correct merging of [ atomtypes ] when bondtypes are used within.

src/pmx/ligand_alchemy.py

@@ -2870,12 +2870,43 @@ def _write_split_itp( self, fname ):
 #**********************************************#
 class _FFatom:
     def __init__(self,list):
+        #logic based on gromacs' toppush implementation: https://github.com/gromacs/gromacs/blob/31146476c633ccd15d1794ef0c4de77db9fda4ee/src/gromacs/gmxpreprocess/toppush.cpp#L340
+
+        if(len(list[3])==1 and list[3].isalpha()):
+            self.have_bonded_type=False
+            self.have_atomic_number=False
+        elif(len(list[4])==1 and list[4].isalpha()):
+            try:
+                int(list[1]) #if this succeeds, we have an atomic number in field 1
+                self.have_bonded_type=False
+                self.have_atomic_number=True
+            except ValueError: # otherwize, it's a bond type
+                self.have_bonded_type=True
+                self.have_atomic_number=False
+        elif(len(list[5])==1 and list[5].isalpha()):
+            self.have_bonded_type=True
+            self.have_atomic_number=True
+        else:
+            raise(Exception("unsupported [ atomtypes ] format"))
+
         self.type = list[0]
-        self.sigmaA = list[1]
-        self.epsA = list[2]
-        self.A = list[3]
-        self.sigmaB = list[4]
-        self.epsB = list[5]
+        shift=0
+        if(self.have_bonded_type and self.have_atomic_number):
+            self.bondtype = list[1]
+            self.anum = list[2]
+            shift=2
+        elif(self.have_bonded_type ):
+            self.bondtype = list[1]
+            shift=1
+        elif(self.have_atomic_number ):
+            self.anum = list[1]
+            shift=1
+        self.sigmaA = list[1+shift] # this should be mass
+        self.epsA = list[2+shift]   # this should be charge
+        self.A = list[3+shift]
+        self.sigmaB = list[4+shift]
+        self.epsB = list[5+shift]
+
 
 class _FFfile:
     def __init__(self, fname=None):
@@ -2890,9 +2921,10 @@ def _read_ffitp(self,file):
         l = open(file).readlines()
         toSkip = "atomtypes"
         for line in l:
-            if toSkip not in line:
-                if line.split():
-                    self.atoms.append(_FFatom(line.split()))
+            stripped_line = line.split(';')[0] #strip comments
+            if toSkip not in stripped_line:
+                if stripped_line.split():
+                    self.atoms.append(_FFatom(stripped_line.split()))
 
 def _get_FF_atoms( ffs ):
     atoms = {}
@@ -2920,5 +2952,5 @@ def _write_FF_file(atoms,file):
 def _merge_FF_files( fnameOut, ffsIn=[] ):
     atoms = _get_FF_atoms( ffsIn )
     _write_FF_file( atoms, fnameOut )
-    
+
 

tests/data/alchemy/atomtypes/ffmol1.itp

@@ -0,0 +1,2 @@
+[ atomtypes ]
+      c3      0.00000      0.00000      A      3.39967e-01      4.57730e-01  ; from gaff

tests/data/alchemy/atomtypes/ffmol2.itp

@@ -0,0 +1,3 @@
+[ atomtypes ]
+      c3      c3      0.00000      0.00000      A      3.39967e-01      4.57730e-01 ; from gaff repeated, now with bond type
+      hc      hc      0.00000      0.00000      A      2.64953e-01      6.56888e-02 ; from gaff with bond type

tests/data/alchemy/atomtypes/ffmol3.itp

@@ -0,0 +1,5 @@
+[ atomtypes ]
+      h1      h1      1      0.00000      0.00000      A      2.47135e-01      6.56888e-02 ; from gaff with bond type and atom number
+      opls_807  C807    12.0110     0.000    A    3.55000E-01   2.92880E-01                ; from opls with bond type
+; cgenff converter to gromacs does not output atomtypes, instead refering to the charmm forcefield for them
+

tests/data/alchemy/atomtypes/ffmol_ref.itp

@@ -0,0 +1,5 @@
+[ atomtypes ]
+      c3      0.00000      0.00000      A      3.39967e-01      4.57730e-01
+      hc      0.00000      0.00000      A      2.64953e-01      6.56888e-02
+      h1      0.00000      0.00000      A      2.47135e-01      6.56888e-02
+      opls_807      12.0110      0.000      A      3.55000E-01      2.92880E-01

tests/test_alchemy.py

@@ -3,6 +3,7 @@
 from pmx.forcefield import Topology
 from pmx.alchemy import mutate, gen_hybrid_top
 from pmx.gmx import set_gmxlib
+from pmx.ligand_alchemy import _merge_FF_files
 from filecmp import cmp
 
 # set GMXLIB variable
@@ -47,3 +48,18 @@ def test_gen_hybrid_top(gf, tmpdir):
     # compare
     cmp(ref_top, out_top)
     cmp(ref_itp, out_itp)
+
+def test_merging_itps(gf, tmpdir):
+    #input files
+    itp1=gf('alchemy/atomtypes/ffmol1.itp')
+    itp2=gf('alchemy/atomtypes/ffmol2.itp')
+    itp3=gf('alchemy/atomtypes/ffmol3.itp')
+    #output file
+    out_itp = str(tmpdir.join("ffmol_out.itp"))
+    # reference output file
+    ref_itp = gf('alchemy/atomtypes/ffmol_ref.itp')
+    #run merge
+    _merge_FF_files(out_itp, [itp1, itp2, itp3])
+    #raise(Exception(out_itp))
+    # compare
+    assert cmp(ref_itp, out_itp), "does not match extected output."