From b6a49a87530c8030d66ace8300a12b0035241d9b Mon Sep 17 00:00:00 2001 From: Mario Date: Thu, 24 Feb 2022 18:01:06 -0500 Subject: [PATCH 1/5] fixed bug when include folder is not toppar or amber in topology cleanup --- GMXMMPBSA/make_top.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/GMXMMPBSA/make_top.py b/GMXMMPBSA/make_top.py index a18fb65d..6ffc56ac 100644 --- a/GMXMMPBSA/make_top.py +++ b/GMXMMPBSA/make_top.py @@ -664,14 +664,7 @@ def cleantop(top_file, ndx): # FIXME: remove minimal components, and compare with the mol_str with open(top_file) as topf: for line in topf: - if line.startswith('#include') and 'forcefield.itp' in line: - if ( - 'charmm' not in line.lower() - and 'toppar' not in line.lower() - and 'amber' not in line.lower() - ): - GMXMMPBSA_ERROR(f'Unknown force field in GROMACS topology in line:\n {line}') - elif '[ molecules ]' in line: + if '[ molecules ]' in line: molsect = True if molsect: # not copy ions and solvent From 77121d136b694cbb8319fe8d38ed0f2895b31263 Mon Sep 17 00:00:00 2001 From: Mario Date: Thu, 24 Feb 2022 18:14:33 -0500 Subject: [PATCH 2/5] fixed error when topology and index do not match --- GMXMMPBSA/make_top.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/GMXMMPBSA/make_top.py b/GMXMMPBSA/make_top.py index 6ffc56ac..23d3600a 100644 --- a/GMXMMPBSA/make_top.py +++ b/GMXMMPBSA/make_top.py @@ -456,7 +456,7 @@ def gmxtop2prmtop(self): if self.FILES.receptor_top: logging.info('A Receptor topology file was defined. Using MT approach...') logging.info('Building AMBER Receptor Topology from GROMACS Receptor Topology...') - rec_top = self.cleantop(self.FILES.receptor_top, self.indexes['REC']) + rec_top = self.cleantop(self.FILES.receptor_top, self.indexes['REC'], 'receptor') if error_info := eq_strs(rec_top, self.receptor_str): if error_info[0] == 'atoms': @@ -498,7 +498,7 @@ def gmxtop2prmtop(self): if self.FILES.ligand_top: logging.info('A Ligand Topology file was defined. Using MT approach...') logging.info('Building AMBER Ligand Topology from GROMACS Ligand Topology...') - lig_top = self.cleantop(self.FILES.ligand_top, self.indexes['LIG']) + lig_top = self.cleantop(self.FILES.ligand_top, self.indexes['LIG'], 'ligand') if error_info := eq_strs(lig_top, self.ligand_str): if error_info[0] == 'atoms': @@ -648,7 +648,7 @@ def pdb2prmtop(self): start += end @staticmethod - def cleantop(top_file, ndx): + def cleantop(top_file, ndx, id='complex'): """ Create a new top file with selected groups and without SOL and IONS :param top_file: User-defined topology file @@ -690,9 +690,17 @@ def cleantop(top_file, ndx): rtemp_top = parmed.gromacs.GromacsTopologyFile(ttp_file.as_posix()) # get the residues in the top from the com_ndx res_list = [] + if len(ndx) < len(rtemp_top.atoms): + GMXMMPBSA_ERROR(f"The {id} index has fewer atoms than the topology. Please check that the files are " + "consistent.") for i in ndx: - if rtemp_top.atoms[i - 1].residue.number + 1 not in res_list: - res_list.append(rtemp_top.atoms[i - 1].residue.number + 1) + try: + idx = rtemp_top.atoms[i - 1].residue.idx + 1 + if idx not in res_list: + res_list.append(rtemp_top.atoms[i - 1].residue.number + 1) + except IndexError: + GMXMMPBSA_ERROR(f'The atom {i} in the {id} index is not found in the topology file. Please check that ' + 'the files are consistent.') ranges = list2range(res_list) rtemp_top.strip(f"!:{','.join(ranges['string'])}") From 53455b2c97cb2a0e0496de37b9a57e3a41a0fbd8 Mon Sep 17 00:00:00 2001 From: Mario Date: Thu, 24 Feb 2022 18:15:04 -0500 Subject: [PATCH 3/5] removed CA and MG from the ion list --- GMXMMPBSA/make_top.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GMXMMPBSA/make_top.py b/GMXMMPBSA/make_top.py index 23d3600a..f59b0bf5 100644 --- a/GMXMMPBSA/make_top.py +++ b/GMXMMPBSA/make_top.py @@ -661,7 +661,7 @@ def cleantop(top_file, ndx, id='complex'): ttp_file = top_file.parent.joinpath('_temp_top.top') temp_top = ttp_file.open(mode='w') # temp_top.write('; Modified by gmx_MMPBSA\n') - # FIXME: remove minimal components, and compare with the mol_str + # TODO: keep solvent when n-wat is implemented with open(top_file) as topf: for line in topf: if '[ molecules ]' in line: @@ -672,7 +672,7 @@ def cleantop(top_file, ndx, id='complex'): # standard gmx form 'NA', 'CL', 'SOL', # charmm-GUI form ?? - 'SOD', 'Na+', 'CLA', 'Cl-', 'CAL', 'CA', 'MG', 'POT', 'K+', + 'SOD', 'Na+', 'CLA', 'Cl-', 'POT', 'K+', 'TIP3P', 'TIP3', 'TP3', 'TIPS3P', 'TIP3o', 'TIP4P', 'TIP4PEW', 'T4E', 'TIP4PD', 'TIP5P', From 9f83cba623611981c80393e7fa730606aee831a3 Mon Sep 17 00:00:00 2001 From: Mario Date: Thu, 24 Feb 2022 18:33:02 -0500 Subject: [PATCH 4/5] improved solvent and ions checking --- GMXMMPBSA/make_top.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/GMXMMPBSA/make_top.py b/GMXMMPBSA/make_top.py index f59b0bf5..097c81ae 100644 --- a/GMXMMPBSA/make_top.py +++ b/GMXMMPBSA/make_top.py @@ -394,23 +394,19 @@ def check4water(self): if counter := sum( res.name in [ - 'SOL', - 'SOD', - 'CLA', - 'TIP3P', - 'TIP4P', - 'TIPS3P', + 'SOD', 'Na+', 'CLA', 'Cl-', 'POT', 'K+', + 'SOL', 'WAT', + 'TIP3P', 'TIP3', 'TP3', 'TIPS3P', 'TIP3o', + 'TIP3P', 'TIP3', 'TP3', 'TIPS3P', 'TIP3o', + 'TIP4P', 'TIP4PEW', 'T4E', 'TIP4PD', 'TIP5P', - 'SPC', - 'SPC/E', - 'SPCE', - 'TIP3o', - 'WAT', + 'SPC', 'SPCE', + 'OPC' ] for res in self.complex_str ): - GMXMMPBSA_ERROR(f'gmx_MMPBSA does not support water molecules in any structure, but we found {counter} ' - f'molecules in the complex.') + GMXMMPBSA_ERROR(f'gmx_MMPBSA does not support water/ions molecules in any structure, but we found' + f' {counter} molecules in the complex.') def gmxtop2prmtop(self): logging.info('Using topology conversion. Setting radiopt = 0...') From 186553978437ab81b536b07d7df6440e896dd365 Mon Sep 17 00:00:00 2001 From: Mario Date: Thu, 24 Feb 2022 18:57:02 -0500 Subject: [PATCH 5/5] refactoring --- GMXMMPBSA/make_top.py | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/GMXMMPBSA/make_top.py b/GMXMMPBSA/make_top.py index 097c81ae..65d788fd 100644 --- a/GMXMMPBSA/make_top.py +++ b/GMXMMPBSA/make_top.py @@ -105,8 +105,8 @@ def __init__(self, FILES, INPUT, external_programs): self.checkFiles() def checkFiles(self): - if (not self.FILES.complex_tpr or not self.FILES.complex_index or not self.FILES.complex_trajs or not - self.FILES.complex_groups): + if (not self.FILES.complex_tpr or not self.FILES.complex_index or + not self.FILES.complex_trajs or not self.FILES.complex_groups): GMXMMPBSA_ERROR('You must define the structure, topology and index files, as well as the groups!') def buildTopology(self): @@ -161,7 +161,7 @@ def buildTopology(self): textwraped = textwrap.wrap('\t'.join(x.string for x in qm_residues), tabsize=4, width=120) logging.info(f'Selected {len(qm_residues)} residues:\n' + '\n'.join(textwraped) + '\n') - self.INPUT['qm_residues'] = ','.join(list2range(qm_residues)['string']) + self.INPUT['qm_residues'] = ','.join(list2range(qm_residues)['string']) if self.INPUT['qmcharge_com'] != rec_charge + lig_charge: logging.warning('System specified with odd number of electrons. Most likely the charge of QM region ' @@ -255,7 +255,7 @@ def gmx2pdb(self): 'fields ignore this warning') # make a temp receptor pdb (even when stability) if decomp to get correct receptor residues from complex. This - # avoid get multiples molecules from complex.split() + # avoids get multiples molecules from complex.split() if self.INPUT['decomprun'] and self.FILES.stability: self.use_temp = True logging.warning('When &decomp is defined, we generate a receptor file in order to extract interface ' @@ -578,7 +578,6 @@ def gmxtop2prmtop(self): def _split_str(self, start, r, c, basename, struct, mut_index=0): end = start + (r[1] - r[0]) mask = f'!:{start}-{end}' - # start += end str_ = self.molstr(struct) if mut_index: str_ = self.makeMutTop(str_, mut_index, True) @@ -755,7 +754,7 @@ def get_selected_residues(self, select, qm_sele=False): if [rres.chain, rres.number, rres.insertion_code] in res_selection: sele_res.append(i) if qm_sele: - rec_charge += round(sum(atm.charge for atm in com_top.residues[i -1].atoms), 0) + rec_charge += round(sum(atm.charge for atm in com_top.residues[i - 1].atoms), 0) res_selection.remove([rres.chain, rres.number, rres.insertion_code]) for j in self.resl: if j.is_receptor(): @@ -769,10 +768,7 @@ def get_selected_residues(self, select, qm_sele=False): for res in res_selection: logging.warning("We couldn't find this residue CHAIN:{} RES_NUM:{} ICODE: {}".format(*res)) sele_res.sort() - if qm_sele: - return sele_res, (rec_charge, lig_charge) - else: - return sele_res + return (sele_res, (rec_charge, lig_charge)) if qm_sele else sele_res def fixparm2amber(self, structure, str_name=None): @@ -983,7 +979,7 @@ def makeMutTop(self, wt_top, mut_index, pdb=False): at.type = h_atoms_prop['type'] at.atom_type = h_atoms_prop['atom_type'] - # change intdiel if cas_intdiel was define before end the mutation process + # change intdiel if cas_intdiel was defined before end the mutation process if self.INPUT['cas_intdiel']: if self.INPUT['gbrun']: if self.INPUT['intdiel'] != 1.0: @@ -1431,13 +1427,13 @@ def _run_tleap(self, tleap, arg1, data_path): def _set_com_order(self, REC, LIG): result = [] - l = 0 - r = 0 + l_idx = 0 + r_idx = 0 for e in self.orderl: if e in ['R', 'REC']: - result.append(REC[r]) - r += 1 + result.append(REC[r_idx]) + r_idx += 1 else: - result.append(LIG[l]) - l += 1 + result.append(LIG[l_idx]) + l_idx += 1 return result