From f34a761126fe1b0d665a0d1cafe23620e7c01b1b Mon Sep 17 00:00:00 2001 From: Yuriy Khalak Date: Wed, 27 Oct 2021 14:13:26 +0200 Subject: [PATCH] allow merging [ atomtypes ] when bondtypes are defined. --- src/pmx/ligand_alchemy.py | 50 ++++++++++++++++++---- tests/data/alchemy/atomtypes/ffmol1.itp | 2 + tests/data/alchemy/atomtypes/ffmol2.itp | 3 ++ tests/data/alchemy/atomtypes/ffmol3.itp | 5 +++ tests/data/alchemy/atomtypes/ffmol_ref.itp | 5 +++ tests/test_alchemy.py | 16 +++++++ 6 files changed, 72 insertions(+), 9 deletions(-) create mode 100644 tests/data/alchemy/atomtypes/ffmol1.itp create mode 100644 tests/data/alchemy/atomtypes/ffmol2.itp create mode 100644 tests/data/alchemy/atomtypes/ffmol3.itp create mode 100644 tests/data/alchemy/atomtypes/ffmol_ref.itp diff --git a/src/pmx/ligand_alchemy.py b/src/pmx/ligand_alchemy.py index d91cebef..a5fb4aef 100644 --- a/src/pmx/ligand_alchemy.py +++ b/src/pmx/ligand_alchemy.py @@ -2792,12 +2792,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): @@ -2812,9 +2843,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 = {} @@ -2842,5 +2874,5 @@ def _write_FF_file(atoms,file): def _merge_FF_files( fnameOut, ffsIn=[] ): atoms = _get_FF_atoms( ffsIn ) _write_FF_file( atoms, fnameOut ) - + diff --git a/tests/data/alchemy/atomtypes/ffmol1.itp b/tests/data/alchemy/atomtypes/ffmol1.itp new file mode 100644 index 00000000..8a95a993 --- /dev/null +++ b/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 diff --git a/tests/data/alchemy/atomtypes/ffmol2.itp b/tests/data/alchemy/atomtypes/ffmol2.itp new file mode 100644 index 00000000..cb9c5f22 --- /dev/null +++ b/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 diff --git a/tests/data/alchemy/atomtypes/ffmol3.itp b/tests/data/alchemy/atomtypes/ffmol3.itp new file mode 100644 index 00000000..7f9ff70a --- /dev/null +++ b/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 + diff --git a/tests/data/alchemy/atomtypes/ffmol_ref.itp b/tests/data/alchemy/atomtypes/ffmol_ref.itp new file mode 100644 index 00000000..20b95b36 --- /dev/null +++ b/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 diff --git a/tests/test_alchemy.py b/tests/test_alchemy.py index a8459d07..861fb50a 100644 --- a/tests/test_alchemy.py +++ b/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."