Skip to content

Commit

Permalink
Merge pull request #165 from Valdes-Tresanco-MS/dev
Browse files Browse the repository at this point in the history 
Improve topology cleanup
  • Loading branch information
@Valdes-Tresanco-MS
Valdes-Tresanco-MS authored Feb 25, 2022
2 parents fca8ed9 + 1865539 commit 6679503
Showing 1 changed file with 38 additions and 45 deletions.
83 changes: 38 additions & 45 deletions GMXMMPBSA/make_top.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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 '
Expand Down Expand Up @@ -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 '
Expand Down Expand Up @@ -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...')
Expand Down Expand Up @@ -456,7 +452,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':
Expand Down Expand Up @@ -498,7 +494,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':
Expand Down Expand Up @@ -582,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)
Expand Down Expand Up @@ -648,7 +643,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
Expand All @@ -661,25 +656,18 @@ def cleantop(top_file, ndx):
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 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
sol_ion = [
# 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',
Expand All @@ -697,9 +685,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'])}")
Expand Down Expand Up @@ -758,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():
Expand All @@ -772,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):

Expand Down Expand Up @@ -986,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:
Expand Down Expand Up @@ -1434,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

0 comments on commit 6679503

Please sign in to comment.