Detailed protocol for working on Linux HPC systems using SLURM as a job queue manager.
Why can't I calculate MM/PBSA and MM/GBSA from MD simulation data generated using the YASARA software?
How to leverage the results generated by YASARA during a molecular dynamics simulation to calculate MM/PBSA and MM/GBSA from 12,000 snapshots on your HPC cluster. Using 64 or 128 parallel processes, the 12,000 snapshots can be processed for MM/GBSA and MM/PBSA in less than 30–45 minutes.
In implicit-solvent calculations (GB/PB), bulk solvent and ions from the
surrounding medium are removed because their effects are accounted for
by the continuum solvent model (with a dielectric constant of
\(\epsilon = 80\)) and the ionic-strength parameter
(saltcon/istrng).
WAT,
SOL, HOH, NA, CL):
These should be automatically filtered out and discarded when the
Amber topology is generated.
tleap (e.g. GAFF2, or
leaprc.water.tip3p for structural water molecules).
&pb section
to the MMPBSA.py input file. It evaluates the polar
desolvation contributions using the Poisson–Boltzmann equation
(radiopt=1, istrng=0.150).
nmode) on all 12,000 snapshots is
computationally impractical because it requires calculation of
the Hessian matrix (with \(O(N^3)\) scaling in CPU time and memory).
If NMA is required, a reduced sampling strategy
(subsampling) should be used—for example, 100–200 snapshots
uniformly distributed along the trajectory, with
--nmode-stride 60.
The next Python script (MM-PBGB-SA_v2.py) supports Poisson–Boltzmann (PB) and
Generalized Born (GB) approaches, intelligent
filtering of HETATM records and
bulk water molecules, and automatic calculation
of Interaction Entropy (IE) and
Normal Mode Analysis (NMA).
#!/usr/bin/env python3
"""
MM/GBSA over YASARA snapshots (PDB files) using AmberTools.
What it does
1. PREPARE (once, from the first snapshot): splits receptor and ligand, builds the
ligand with correct bond orders (RDKit + your SMILES), parametrizes it (GAFF2,
AM1-BCC, fixed net charge), builds the complex/receptor/ligand topologies with
tleap (protein force field + GAFF2) and CHECKS that atom count and atom order
match the PDB.
2. For EVERY selected snapshot (in parallel, one process per snapshot):
- converts the PDB to Amber coordinates (receptor atoms first, then ligand),
- relaxes ONLY THE HYDROGENS (heavy atoms restrained) with a short GB
minimization, and
- runs MM/GBSA (single trajectory) with MMPBSA.py.
3. SUMMARY: mean, standard deviation, standard errors and a plot.
Why the hydrogen relaxation is needed
Snapshots written by YASARA contain hydrogen contacts that the AMBER force field
considers clashes. In the test snapshot, evaluating it directly gave a binding
energy of +59 kcal/mol; after relaxing only the hydrogens (heavy atoms moved
0.04 A RMSD) it was -17 kcal/mol. 100 minimization steps are enough (300 steps
changed the result by 0.06 kcal/mol).
Environment (create once on the login node):
conda create -n mmgbsa -c conda-forge ambertools rdkit pandas numpy matplotlib python=3.12
Usage:
python mmgbsa_yasara.py "/path/to/snapshots/*.pdb" --workdir run1 \
--ligand-resname UNL --ligand-charge -1 \
--ligand-smiles "O=C1C(=Cc2ccc(o2)-c2ccccc2F)C(Cc2ccccc2)=NN1c1ccc(cc1)C(=O)[O-]" \
--stride 4 --processes 32
Add --prepare-only to build and check the topologies without running frames.
Add --limit N to run only the first N selected frames (useful for a timing test).
Re-running the same command resumes: frames with a result are skipped.
"""
#!/usr/bin/env python3
"""
MM/GBSA & MM/PBSA + Per-Residue Energy Decomposition + Interaction Entropy & NMA
over YASARA MD snapshots using AmberTools.
"""
import argparse
import glob
import json
import math
import multiprocessing as mp
import os
import re
import shutil
import subprocess
import sys
import time
from pathlib import Path
import numpy as np
MIN_TEMPLATE = """relax hydrogens only
&cntrl
imin=1, maxcyc={steps}, ncyc={ncyc}, ntb=0, igb={igb}, cut=999., saltcon={salt},
ntr=1, restraint_wt=100.0, restraintmask='!@H=', ntpr=50,
/
"""
MMPBSA_TEMPLATE = """MM/GBSA and MM/PBSA calculation
&general
startframe=1, endframe=999999, interval=1, verbose=2, keep_files=0,
/
&gb
igb={igb}, saltcon={salt},
/
&pb
istrng={salt}, radiopt=1, fillratio=4.0,
/
{decomp_section}
"""
NMODE_TEMPLATE = """Entropy via Normal Mode Analysis
&general
startframe=1, endframe=999999, interval=1, verbose=2, keep_files=0,
/
&nmode
nmstartframe=1, nmendframe=1, nsteps=10000, maxcyc=10000,
/
"""
RESULT_KEYS_GB = ('VDWAALS', 'EEL', 'EGB', 'ESURF', 'DELTA G gas', 'DELTA G solv', 'DELTA TOTAL')
CFG = {}
PIPELINE_TAG = 'v4-decomp-supported'
_PERM = None
def natural_key(text):
return [int(t) if t.isdigit() else t for t in re.split(r'(\d+)', text)]
def run(cmd, cwd, timeout=None):
env = dict(os.environ, OMP_NUM_THREADS='1', MKL_NUM_THREADS='1', OPENBLAS_NUM_THREADS='1')
proc = subprocess.run(cmd, cwd=str(cwd), capture_output=True, text=True, timeout=timeout, env=env)
if proc.returncode != 0:
tail = '\n'.join((proc.stdout + '\n' + proc.stderr).strip().split('\n')[-12:])
raise RuntimeError(f'Command failed ({proc.returncode}): {" ".join(cmd)}\n{tail}')
return proc
def read_pdb_lines(path):
return Path(path).read_text(errors='replace').replace('\r\n', '\n').replace('\r', '\n').split('\n')
def split_pdb_advanced(path, lig_resname, keep_hetatms=None):
keep_hetatms = set(keep_hetatms or [])
bulk_resnames = {'WAT', 'HOH', 'SOL', 'TIP3', 'TP3', 'NA', 'CL', 'K', 'MG'}
rec, lig, conect = [], [], []
for line in read_pdb_lines(path):
if line.startswith(('ATOM ', 'HETATM')):
resname = line[17:20].strip()
if resname == lig_resname:
lig.append(line)
elif line.startswith('ATOM') or resname in keep_hetatms:
rec.append(line)
elif resname in bulk_resnames:
continue
else:
rec.append(line)
elif line.startswith('CONECT'):
conect.append(line)
return rec, lig, conect
def xyz(line):
return (float(line[30:38]), float(line[38:46]), float(line[46:54]))
def write_rst7(path, coords):
flat = [c for p in coords for c in p]
if max(abs(v) for v in flat) >= 999.0:
raise ValueError('Coordinate outside range of Amber F12.7 format')
with open(path, 'w') as f:
f.write('snapshot\n')
f.write(f'{len(coords):5d}\n')
for i in range(0, len(flat), 6):
f.write(''.join(f'{v:12.7f}' for v in flat[i:i + 6]) + '\n')
def fix_receptor(rec_lines):
residues = []
index = {}
for i, line in enumerate(rec_lines):
key = (line[21], line[22:27])
if key not in index:
index[key] = len(residues)
residues.append([key, [], set()])
residues[index[key]][1].append(i)
residues[index[key]][2].add(line[12:16].strip())
new_name = {}
counts = {}
for k, (key, idxs, names) in enumerate(residues):
resname = rec_lines[idxs[0]][17:20]
out = None
if resname == 'HIS':
d, e = 'HD1' in names, 'HE2' in names
out = 'HIP' if d and e else 'HID' if d else 'HIE'
elif resname == 'ASP' and 'HD2' in names:
out = 'ASH'
elif resname == 'GLU' and 'HE2' in names:
out = 'GLH'
if out:
new_name[k] = out
counts[out] = counts.get(out, 0) + 1
out_lines = []
prev_chain = None
for i, line in enumerate(rec_lines):
if prev_chain is not None and line[21] != prev_chain:
out_lines.append('TER')
prev_chain = line[21]
res_idx = index[(line[21], line[22:27])]
if res_idx in new_name:
line = line[:17] + new_name[res_idx] + line[20:]
out_lines.append(line)
out_lines.append('TER')
out_lines.append('END')
return out_lines, counts
def build_ligand_sdf(lig_lines, conect_lines, smiles, charge, out_sdf):
from rdkit import Chem
from rdkit.Chem import AllChem
serials = {int(l[6:11]) for l in lig_lines if l[6:11].strip().isdigit()}
conect = [l for l in conect_lines if int(l[6:11]) in serials]
block = '\n'.join(lig_lines + conect + ['END']) + '\n'
tmp = Path(out_sdf).with_suffix('.pdb')
tmp.write_text(block)
pdbmol = Chem.MolFromPDBFile(str(tmp), removeHs=False, sanitize=False, proximityBonding=not conect)
tmpl = Chem.MolFromSmiles(smiles)
if Chem.GetFormalCharge(tmpl) != charge:
raise RuntimeError(f'SMILES charge ({Chem.GetFormalCharge(tmpl)}) != --ligand-charge ({charge})')
heavy = Chem.RemoveHs(pdbmol, sanitize=False)
for b in heavy.GetBonds():
b.SetBondType(Chem.BondType.SINGLE)
for a in heavy.GetAtoms():
a.SetNoImplicit(False)
a.SetFormalCharge(0)
fixed = AllChem.AssignBondOrdersFromTemplate(tmpl, heavy)
rw = Chem.RWMol(pdbmol)
idx_heavy = [a.GetIdx() for a in pdbmol.GetAtoms() if a.GetAtomicNum() > 1]
for b in fixed.GetBonds():
a1, a2 = idx_heavy[b.GetBeginAtomIdx()], idx_heavy[b.GetEndAtomIdx()]
bond = rw.GetBondBetweenAtoms(a1, a2)
bond.SetBondType(b.GetBondType())
bond.SetIsAromatic(b.GetIsAromatic())
for i, hi in enumerate(idx_heavy):
rw.GetAtomWithIdx(hi).SetFormalCharge(fixed.GetAtomWithIdx(i).GetFormalCharge())
mol = rw.GetMol()
Chem.SanitizeMol(mol)
w = Chem.SDWriter(str(out_sdf))
w.write(mol)
w.close()
def fix_mol2_charges(mol2_in, mol2_out, total_charge):
head, rest = Path(mol2_in).read_text().split('@ATOM')
body, tail = rest.split('@BOND')
rows = [r for r in body.strip('\n').split('\n') if r.strip()]
q = [float(r.split()[8]) for r in rows]
residual = round(total_charge - sum(q), 6)
k = max(range(len(q)), key=lambda i: abs(q[i]))
q[k] = round(q[k] + residual, 6)
new = []
for r, qi in zip(rows, q):
c = r.split()
new.append(f'{int(c[0]):7d} {c[1]:<8s} {float(c[2]):9.4f} {float(c[3]):9.4f} '
f'{float(c[4]):9.4f} {c[5]:<4s} {c[6]:>3s} {c[7]:<8s} {qi:10.6f}')
Path(mol2_out).write_text(head + '@ATOM\n' + '\n'.join(new) + '\n@BOND' + tail)
return abs(residual)
def read_rst7(path):
lines = Path(path).read_text().split('\n')
natoms = int(lines[1].split()[0])
vals = []
for line in lines[2:]:
vals += [float(line[i:i + 12]) for i in range(0, len(line), 12) if line[i:i + 12].strip()]
vals = vals[:3 * natoms]
return [tuple(vals[i:i + 3]) for i in range(0, len(vals), 3)]
def read_prmtop_atoms(prmtop):
sections = {}
for block in Path(prmtop).read_text().split('%FLAG ')[1:]:
head, _, body = block.partition('\n')
sections[head.strip()] = body.split('\n')[1:]
pointers = [int(x) for line in sections['POINTERS'] for x in line.split()]
natom, nres = pointers[0], pointers[11]
names = []
for line in sections['ATOM_NAME']:
names += [line[i:i + 4].strip() for i in range(0, len(line), 4) if line[i:i + 4].strip()]
first = [int(x) for line in sections['RESIDUE_POINTER'] for x in line.split()][:nres]
res_of_atom = []
for r in range(nres):
end = first[r + 1] - 1 if r + 1 < nres else natom
res_of_atom += [r] * (end - (first[r] - 1))
return names[:natom], res_of_atom
def build_permutation(prep, rec_lines, lig_lines):
names, res_of_atom = read_prmtop_atoms(prep / 'complex.prmtop')
n_rec, n_atoms = len(rec_lines), len(names)
if n_atoms != n_rec + len(lig_lines):
raise RuntimeError(f'Topology has {n_atoms} atoms but PDB has {n_rec + len(lig_lines)}')
ordinal, pdb_index = {}, {}
for j, line in enumerate(rec_lines):
key = (line[21], line[22:27])
ordinal.setdefault(key, len(ordinal))
k = (ordinal[key], line[12:16].strip())
if k in pdb_index:
raise RuntimeError(f'Duplicated atom name {k[1]} in residue {line[17:20]}{line[22:27]}')
pdb_index[k] = j
n_res_rec = len(ordinal)
perm = []
for i in range(n_atoms):
r = res_of_atom[i]
if r < n_res_rec:
k = (r, names[i])
if k not in pdb_index:
raise RuntimeError(f'Topology atom {names[i]} (residue {r + 1}) not found in PDB')
perm.append(pdb_index[k])
else:
perm.append(i)
if sorted(perm) != list(range(n_atoms)):
raise RuntimeError('Atom mapping between PDB and topology is not a permutation')
return perm
def atom_signature(lines):
import hashlib
text = '|'.join(l[12:16] + l[17:20] + l[21:27] for l in lines)
return hashlib.md5(text.encode()).hexdigest()
def get_perm():
global _PERM
if _PERM is None:
_PERM = json.loads((Path(CFG['prep_dir']) / 'perm.json').read_text())
return _PERM
def rec_fixed_atoms(lines):
return [l for l in lines if l.startswith(('ATOM', 'HETATM'))]
def prepare(args, first_pdb):
work = Path(args.workdir).resolve()
prep = work / 'prep'
prep.mkdir(parents=True, exist_ok=True)
if (prep / 'READY').exists() and (prep / 'perm.json').exists():
print('Topologies ready; skipping preparation.', flush=True)
return prep
print(f'PREPARE: using {first_pdb}', flush=True)
rec, lig, conect = split_pdb_advanced(first_pdb, args.ligand_resname, args.keep_hetatm)
print(f' Receptor atoms: {len(rec)} | Ligand atoms: {len(lig)}', flush=True)
rec_fixed, _ = fix_receptor(rec)
(prep / 'receptor.pdb').write_text('\n'.join(rec_fixed) + '\n')
build_ligand_sdf(lig, conect, args.ligand_smiles, args.ligand_charge, prep / 'lig.sdf')
run(['antechamber', '-i', 'lig.sdf', '-fi', 'mdl', '-o', 'lig_raw.mol2', '-fo', 'mol2',
'-c', 'bcc', '-nc', str(args.ligand_charge), '-at', 'gaff2',
'-rn', args.ligand_resname, '-s', '0', '-pf', 'y'], prep)
resid = fix_mol2_charges(prep / 'lig_raw.mol2', prep / 'lig.mol2', args.ligand_charge)
print(f' ligand charge rounding residual corrected: {resid:.4f} e', flush=True)
run(['parmchk2', '-i', 'lig.mol2', '-f', 'mol2', '-o', 'lig.frcmod', '-s', 'gaff2'], prep)
radii = 'mbondi3' if args.igb == 8 else 'mbondi2'
(prep / 'tleap.in').write_text(
f'source leaprc.protein.{args.protein_ff}\n'
'source leaprc.gaff2\n'
'loadamberparams lig.frcmod\n'
'LIG = loadmol2 lig.mol2\n'
'REC = loadpdb receptor.pdb\n'
f'set default PBRadii {radii}\n'
'COM = combine {REC LIG}\n'
'check COM\n'
'saveamberparm COM complex.prmtop complex.inpcrd\n'
'quit\n')
run(['tleap', '-f', 'tleap.in'], prep)
leap_log = (prep / 'leap.log').read_text() if (prep / 'leap.log').exists() else ''
m = re.findall(r'Exiting LEaP: Errors = (\d+)', leap_log)
if m and int(m[-1]) > 0:
raise RuntimeError('tleap reported errors; see prep/leap.log')
run(['ante-MMPBSA.py', '-p', 'complex.prmtop', '-c', 'complex_dry.prmtop',
'-r', 'receptor.prmtop', '-l', 'ligand.prmtop', '-n', f':{args.ligand_resname}',
'--radii', radii], prep)
perm = build_permutation(prep, rec_fixed_atoms(rec_fixed), lig)
coords_pdb = [xyz(l) for l in rec] + [xyz(l) for l in lig]
coords_top = read_rst7(prep / 'complex.inpcrd')
worst = max(math.dist(coords_top[i], coords_pdb[perm[i]]) for i in range(len(perm)))
if worst > 0.01:
raise RuntimeError(f'Atom mapping check failed (max deviation {worst:.3f} A)')
print(f' validation OK: {len(perm)} atoms mapped between PDB and topology '
f'(max deviation {worst:.4f} A)', flush=True)
(prep / 'perm.json').write_text(json.dumps(perm))
(prep / 'READY').write_text('ok\n')
return prep
def _init_worker(cfg):
global CFG
CFG = cfg
def parse_mmpbsa_output(dat_path):
text = Path(dat_path).read_text()
blocks = text.split('Differences (Complex - Receptor - Ligand):')
gb_data, pb_data = {}, {}
if len(blocks) > 1:
gb_text = blocks[1]
for key in RESULT_KEYS_GB:
m = re.search(r'^' + re.escape(key) + r'\s+(-?[0-9.]+)', gb_text, re.M)
if m:
gb_data['GB_' + key] = float(m.group(1))
if len(blocks) > 2:
pb_text = blocks[2]
for key in RESULT_KEYS_GB:
m = re.search(r'^' + re.escape(key) + r'\s+(-?[0-9.]+)', pb_text, re.M)
if m:
pb_data['PB_' + key] = float(m.group(1))
return {**gb_data, **pb_data}
def parse_decomp_output(fdir):
"""
Parses per-residue energy decomposition from FINAL_DECOMP.csv or FINAL_DECOMP.dat.
Returns dict: { residue_name: {'vdw': float, 'elec': float, 'pol': float, 'npol': float, 'total': float} }
"""
csv_path = Path(fdir) / 'FINAL_DECOMP.csv'
dat_path = Path(fdir) / 'FINAL_DECOMP.dat'
decomp_results = {}
if csv_path.exists():
text = csv_path.read_text(errors='replace')
if 'Differences (Complex - Receptor - Ligand)' in text:
text = text.split('Differences (Complex - Receptor - Ligand)')[1]
elif 'DELTA' in text:
text = text.split('DELTA')[1]
for line in text.splitlines():
parts = [p.strip('" ').strip() for p in line.split(',')]
if len(parts) >= 7 and parts[0] not in ('Residue', 'Residue Name', 'Total Energy Decomposition', ''):
try:
res_name = parts[0]
vdw = float(parts[2])
elec = float(parts[3])
pol = float(parts[4])
npol = float(parts[5])
tot = float(parts[6])
decomp_results[res_name] = {'vdw': vdw, 'elec': elec, 'pol': pol, 'npol': npol, 'total': tot}
except ValueError:
continue
elif dat_path.exists():
text = dat_path.read_text(errors='replace')
if 'Differences (Complex - Receptor - Ligand):' in text:
text = text.split('Differences (Complex - Receptor - Ligand):')[1]
if 'Total Energy Decomposition:' in text:
text = text.split('Total Energy Decomposition:')[1]
if 'Sidechain Energy Decomposition:' in text:
text = text.split('Sidechain Energy Decomposition:')[0]
for line in text.splitlines():
if '|' in line:
parts = [p.strip() for p in line.split('|')]
if len(parts) >= 7 and not parts[0].startswith(('Residue', '-')):
try:
res_name = parts[0]
vdw = float(parts[2])
elec = float(parts[3])
pol = float(parts[4])
npol = float(parts[5])
tot = float(parts[6])
decomp_results[res_name] = {'vdw': vdw, 'elec': elec, 'pol': pol, 'npol': npol, 'total': tot}
except ValueError:
continue
return decomp_results
def process_frame(job):
idx, path = job
fdir = Path(CFG['frames_dir']) / f'frame_{idx:05d}'
result_file = fdir / 'result.json'
if result_file.exists():
cached = json.loads(result_file.read_text())
if cached.get('pipeline') == PIPELINE_TAG:
return cached
t0 = time.time()
res = {'index': idx, 'file': str(path), 'status': 'OK', 'message': '', 'pipeline': PIPELINE_TAG}
try:
fdir.mkdir(parents=True, exist_ok=True)
rec, lig, _ = split_pdb_advanced(path, CFG['ligand_resname'], CFG['keep_hetatm'])
if len(rec) != CFG['n_rec'] or len(lig) != CFG['n_lig']:
raise RuntimeError(f'Atom counts differ: receptor {len(rec)} (exp {CFG["n_rec"]}), ligand {len(lig)} (exp {CFG["n_lig"]})')
if atom_signature(rec + lig) != CFG['signature']:
raise RuntimeError('Atom signature mismatch between snapshot and template')
coords_pdb = [xyz(l) for l in rec] + [xyz(l) for l in lig]
write_rst7(fdir / 'start.rst7', [coords_pdb[j] for j in get_perm()])
prep = Path(CFG['prep_dir'])
(fdir / 'min.in').write_text(MIN_TEMPLATE.format(
steps=CFG['min_steps'], ncyc=max(1, CFG['min_steps'] // 2), igb=CFG['igb'], salt=CFG['saltcon']))
run(['sander', '-O', '-i', 'min.in', '-o', 'min.out', '-p', str(prep / 'complex.prmtop'),
'-c', 'start.rst7', '-ref', 'start.rst7', '-r', 'relaxed.rst7'], fdir, CFG['timeout'])
decomp_sec = ""
if CFG['run_decomp']:
print_res_line = f' print_res="{CFG["print_res"]}",\n' if CFG['print_res'] else ''
decomp_sec = f"""&decomp\n idecomp=2, dec_verbose=1, csv_format=1,\n{print_res_line}/"""
(fdir / 'mmpbsa.in').write_text(MMPBSA_TEMPLATE.format(igb=CFG['igb'], salt=CFG['saltcon'], decomp_section=decomp_sec))
cmd = ['MMPBSA.py', '-O', '-i', 'mmpbsa.in', '-o', 'FINAL.dat']
if CFG['run_decomp']:
cmd += ['-do', 'FINAL_DECOMP.dat', '-deo', 'FINAL_DECOMP.csv']
cmd += ['-cp', str(prep / 'complex.prmtop'), '-rp', str(prep / 'receptor.prmtop'),
'-lp', str(prep / 'ligand.prmtop'), '-y', 'relaxed.rst7']
run(cmd, fdir, CFG['timeout'])
res.update(parse_mmpbsa_output(fdir / 'FINAL.dat'))
if CFG['run_decomp']:
res['decomp'] = parse_decomp_output(fdir)
if CFG['run_nmode'] and (idx % CFG['nmode_stride'] == 0):
(fdir / 'nmode.in').write_text(NMODE_TEMPLATE)
run(['MMPBSA.py', '-O', '-i', 'nmode.in', '-o', 'NMODE.dat',
'-cp', str(prep / 'complex.prmtop'), '-rp', str(prep / 'receptor.prmtop'),
'-lp', str(prep / 'ligand.prmtop'), '-y', 'relaxed.rst7'], fdir, CFG['timeout'] * 2)
nmode_txt = (fdir / 'NMODE.dat').read_text()
m = re.search(r'DELTA S total=\s*(-?[0-9.]+)', nmode_txt)
if m:
res['NMODE_TDS'] = float(m.group(1))
except Exception as e:
res['status'] = 'ERROR'
res['message'] = str(e)[:500]
res['seconds'] = round(time.time() - t0, 1)
result_file.write_text(json.dumps(res))
if res['status'] == 'OK' and not CFG['keep_files']:
for name in ('start.rst7', 'relaxed.rst7', 'min.out', 'min.in', 'mmpbsa.in', 'nmode.in', 'FINAL_DECOMP.dat', 'FINAL_DECOMP.csv'):
(fdir / name).unlink(missing_ok=True)
return res
def compute_interaction_entropy(df, temp=298.15):
R = 0.0019872042
beta = 1.0 / (R * temp)
E_gas = df['GB_VDWAALS'] + df['GB_EEL']
E_gas_mean = E_gas.mean()
dE = E_gas - E_gas_mean
exp_dE = np.exp(beta * dE)
minus_T_dS = (1.0 / beta) * np.log(exp_dE.mean())
return minus_T_dS
def summarize(workdir, results, temp=298.15):
import pandas as pd
df = pd.DataFrame(results).sort_values('index')
csv_path = Path(workdir) / 'results_frames.csv'
df.to_csv(csv_path, index=False)
ok = df[df['status'] == 'OK']
lines = [f'Frames Processed: {len(df)} | OK: {len(ok)} | Errors: {len(df) - len(ok)}']
if len(ok) >= 2:
gb_total = ok['GB_DELTA TOTAL'].to_numpy()
pb_total = ok['PB_DELTA TOTAL'].to_numpy()
tds_ie = compute_interaction_entropy(ok, temp=temp)
abs_gb_ie = gb_total.mean() + tds_ie
abs_pb_ie = pb_total.mean() + tds_ie
lines += [
"--------------------------------------------------",
f"MM/GBSA Delta H (mean): {gb_total.mean():.2f} +/- {gb_total.std():.2f} kcal/mol",
f"MM/PBSA Delta H (mean): {pb_total.mean():.2f} +/- {pb_total.std():.2f} kcal/mol",
"--------------------------------------------------",
f"Interaction Entropy (-T*dS_IE): {tds_ie:.2f} kcal/mol (at {temp} K)",
f"Delta G estimate (MM/GBSA + IE): {abs_gb_ie:.2f} kcal/mol",
f"Delta G estimate (MM/PBSA + IE): {abs_pb_ie:.2f} kcal/mol",
"--------------------------------------------------"
]
for label, arr in (('MM/GBSA', gb_total), ('MM/PBSA', pb_total)):
nb = 5 if len(arr) >= 10 else 2
sem_block = np.std([b.mean() for b in np.array_split(arr, nb)], ddof=1) / np.sqrt(nb)
half = len(arr) // 2
lines.append(f"{label}: SEM (block averaging, {nb} blocks) {sem_block:.2f} | "
f"first half {arr[:half].mean():.2f} vs second half {arr[half:].mean():.2f}")
# Per-residue Decomposition Summary
if 'decomp' in ok.columns:
residue_data = {}
for res_dict in ok['decomp'].dropna():
if isinstance(res_dict, dict):
for r_name, vals in res_dict.items():
if r_name not in residue_data:
residue_data[r_name] = {'total': [], 'vdw': [], 'elec': [], 'solv': []}
residue_data[r_name]['total'].append(vals['total'])
residue_data[r_name]['vdw'].append(vals['vdw'])
residue_data[r_name]['elec'].append(vals['elec'])
residue_data[r_name]['solv'].append(vals['pol'] + vals['npol'])
decomp_rows = []
for r_name, data in residue_data.items():
n_f = len(data['total'])
if n_f == 0:
continue
tot_arr = np.array(data['total'])
std_tot = tot_arr.std(ddof=1) if n_f > 1 else 0.0
decomp_rows.append({
'Residue': r_name,
'Delta_G_mean': tot_arr.mean(),
'SD': std_tot,
'SEM': std_tot / math.sqrt(n_f),
'VdW_mean': np.array(data['vdw']).mean(),
'Elec_mean': np.array(data['elec']).mean(),
'Solv_mean': np.array(data['solv']).mean()
})
if decomp_rows:
df_decomp = pd.DataFrame(decomp_rows).sort_values('Delta_G_mean')
decomp_csv_path = Path(workdir) / 'results_decomp_residues.csv'
df_decomp.to_csv(decomp_csv_path, index=False)
hotspots = df_decomp[df_decomp['Delta_G_mean'] <= -0.5]
lines.append("--------------------------------------------------")
lines.append(f"PER-RESIDUE DECOMPOSITION: Top Hotspots (Delta G <= -0.5 kcal/mol)")
lines.append(f"{'Residue':<14s} | {'Delta G (mean)':<14s} | {'SD':<6s} | {'VdW':<7s} | {'Elec':<7s} | {'Solv':<7s}")
lines.append("-" * 75)
for _, row in hotspots.head(15).iterrows():
lines.append(f"{row['Residue']:<14s} | {row['Delta_G_mean']:8.2f} kcal/mol | {row['SD']:6.2f} | {row['VdW_mean']:7.2f} | {row['Elec_mean']:7.2f} | {row['Solv_mean']:7.2f}")
lines.append("--------------------------------------------------")
try:
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
top_plot = df_decomp.head(20)
fig, ax = plt.subplots(figsize=(10, 5))
ax.bar(top_plot['Residue'], top_plot['Delta_G_mean'], yerr=top_plot['SEM'], capsize=3, color='#1f77b4', edgecolor='black')
ax.axhline(0, color='grey', linewidth=0.8, linestyle='--')
ax.set_ylabel('Delta G contribution (kcal/mol)')
ax.set_xlabel('Residue')
ax.set_title('Per-Residue Energy Decomposition (Top 20 Hotspots)')
plt.xticks(rotation=45, ha='right')
fig.tight_layout()
fig.savefig(Path(workdir) / 'decomp_hotspots.png', dpi=150)
lines.append("Plot generated: decomp_hotspots.png")
except Exception as e:
lines.append(f"(Decomp plot skipped: {e})")
summary_text = '\n'.join(lines)
(Path(workdir) / 'summary.txt').write_text(summary_text + '\n')
print('=' * 60)
print(summary_text)
print(f'\nFull results saved to: {csv_path}')
print('=' * 60)
def main():
ap = argparse.ArgumentParser()
ap.add_argument('frames', help='Glob pattern for snapshot PDBs')
ap.add_argument('--workdir', required=True)
ap.add_argument('--ligand-resname', required=True)
ap.add_argument('--ligand-charge', type=int, required=True)
ap.add_argument('--ligand-smiles', required=True)
ap.add_argument('--keep-hetatm', nargs='*', help='Non-standard HETATM resnames to keep in receptor')
ap.add_argument('--protein-ff', default='ff14SB')
ap.add_argument('--igb', type=int, default=2)
ap.add_argument('--saltcon', type=float, default=0.150)
ap.add_argument('--min-steps', type=int, default=100)
ap.add_argument('--stride', type=int, default=1)
ap.add_argument('--decomp', action='store_true', help='Enable per-residue energy decomposition')
ap.add_argument('--print-res', default=None, help='Restrict decomposition printing (e.g. "within 6.0" or "1-300")')
ap.add_argument('--run-nmode', action='store_true', help='Enable NMA entropy calculation on a subset')
ap.add_argument('--nmode-stride', type=int, default=100, help='Stride for NMA frames')
ap.add_argument('--processes', type=int, default=None)
ap.add_argument('--timeout', type=int, default=3600)
ap.add_argument('--keep-files', action='store_true')
ap.add_argument('--limit', type=int, default=None, help='only the first N selected frames (timing test)')
ap.add_argument('--prepare-only', action='store_true')
args = ap.parse_args()
for tool in ('tleap', 'antechamber', 'parmchk2', 'sander', 'MMPBSA.py', 'ante-MMPBSA.py'):
if shutil.which(tool) is None:
sys.exit(f'Error: {tool} not found. Activate the environment with AmberTools first.')
files = sorted(glob.glob(args.frames), key=natural_key)
if not files:
sys.exit(f'No files match {args.frames}')
selected = files[::args.stride]
if args.limit:
selected = selected[:args.limit]
print(f'{len(files)} snapshots found; {len(selected)} selected (stride {args.stride}).', flush=True)
work = Path(args.workdir).resolve()
prep = prepare(args, files[0])
if args.prepare_only:
print('Prepare-only: done.')
return
rec, lig, _ = split_pdb_advanced(files[0], args.ligand_resname, args.keep_hetatm)
cfg = dict(frames_dir=str(work / 'frames'), prep_dir=str(prep), ligand_resname=args.ligand_resname,
keep_hetatm=args.keep_hetatm, n_rec=len(rec), n_lig=len(lig), signature=atom_signature(rec + lig),
igb=args.igb, saltcon=args.saltcon, min_steps=args.min_steps, timeout=args.timeout,
keep_files=args.keep_files, run_decomp=args.decomp, print_res=args.print_res,
run_nmode=args.run_nmode, nmode_stride=args.nmode_stride)
procs = args.processes or int(os.environ.get('SLURM_CPUS_PER_TASK', 0)) or os.cpu_count() or 1
jobs = [(i, f) for i, f in zip(range(0, len(files), args.stride), selected)]
procs = max(1, min(procs, len(jobs)))
print(f'Running MM/GBSA + MM/PBSA (Decomp={args.decomp}) on {len(jobs)} frames with {procs} processes...', flush=True)
results, start = [], time.time()
with mp.Pool(procs, initializer=_init_worker, initargs=(cfg,)) as pool:
for n, res in enumerate(pool.imap_unordered(process_frame, jobs), start=1):
results.append(res)
if res['status'] != 'OK':
print(f' frame {res["index"]}: ERROR {res["message"].splitlines()[0][:150]}', flush=True)
if n % max(1, len(jobs) // 50) == 0 or n == len(jobs):
el = time.time() - start
print(f' {n}/{len(jobs)} frames done ({el / 60:.1f} min, ~{el / n * (len(jobs) - n) / 60:.0f} min left)', flush=True)
summarize(work, results)
if __name__ == '__main__':
main()
Slurm Submission Script (submit_mmpbsa.sbatch)
#!/bin/bash
#SBATCH --job-name=MM_PBSA_12k
#SBATCH --qos=long
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=64
#SBATCH --mem=100G
#SBATCH --time=12:00:00
#SBATCH --output=mmpbsa_%j.log
# Activar el entorno de conda con AmberTools y RDKit
source ~/miniconda3/etc/profile.d/conda.sh
conda activate mmgbsa
# Ejecutar script
python3 MM-PBGB-SA_v2.py "/home/jant.encinar/runMD/MD01-bTrCp1-GS143_03md_fin/*.pdb" \
--limit 12000 \
--workdir run_12000_frames \
--ligand-resname UNL \
--ligand-charge -1 \
--ligand-smiles "O=C1C(=Cc2ccc(o2)-c2ccccc2F)C(Cc2ccccc2)=NN1c1ccc(cc1)C(=O)[O-]" \
--stride 1 \
--processes 64
| Method | Purpose | Computational Cost for 12,000 Snapshots |
|---|---|---|
| MM/GBSA (\(\Delta H\)) | Generalized Born implicit-solvent calculation. | ~15 min (64 cores) |
| MM/PBSA (\(\Delta H\)) | Poisson–Boltzmann implicit-solvent calculation. | ~35 min (64 cores) |
| Interaction Entropy (\(-T\Delta S_{\mathrm{IE}}\)) | Entropy estimated directly from interaction-energy fluctuations. | 0 seconds (immediate post-processing) |
| Normal Mode Analysis (NMA) (\(-T\Delta S_{\mathrm{NMA}}\)) | Entropy estimated from normal modes (optional subsampling). | Optional; enabled with --run-nmode --nmode-stride 100 |
The results and files generated during script execution are organized
within the working directory specified by --workdir.
They are divided into four main categories:
--workdir)summary.txt:
A plain-text file containing the final quantitative summary.
It includes:
VDWAALS, EEL, EGB, ESURF, PB).--run-nmode option was enabled.results_frames.csv:
A CSV-formatted table containing one row for each of the
12,000 snapshots. Its columns include:
index, file, status,
all individual GB and PB energy terms, computation time in
seconds (seconds), and the error message if a
frame fails.
--workdir/prep/)These files are generated only at the beginning of the calculation, using the first snapshot, and provide the fixed structural basis for the entire analysis:
receptor.pdb:
Receptor PDB structure with solvent molecules removed and
standardized residue protonation states (HIS, ASP, GLU)
for Amber.
lig.sdf:
Three-dimensional ligand structure in SDF format, with bond
ordering corrected using RDKit based on the supplied SMILES.
lig.mol2 and
lig.frcmod:
Ligand force-field parameters generated using
antechamber (GAFF2, AM1-BCC charges) and
parmchk2.
complex.prmtop and
complex.inpcrd:
Amber topology and coordinate files for the complete
solvent-free complex, generated using tleap.
complex_dry.prmtop,
receptor.prmtop, and
ligand.prmtop:
Topology files for three independent groups, generated using
ante-MMPBSA.py for the single trajectory.
tleap.in and
leap.log:
tleap input script and execution log, respectively.
READY:
A text-based flag file indicating that the topologies have
been successfully constructed. If the script is executed
again, it detects this file and skips the preparation stage.
--workdir/frames/frame_XXXXX/)
A separate directory is created for each processed frame
(e.g. frame_00000, frame_00001, etc.):
result.json:
A lightweight JSON file storing the energy-result dictionary
for the corresponding snapshot. It serves as a checkpoint:
if the job is interrupted and subsequently restarted, the
script reads these files and skips snapshots that have
already been processed.
--keep-files flag is
specified. This helps save disk space and inodes on the cluster.
The files include:
start.rst7: Amber coordinates converted from the original PDB file.relaxed.rst7: Coordinates after hydrogen minimization using sander.min.in and min.out: Input and output files for the minimization step.mmpbsa.in and FINAL.dat: MMPBSA.py command/input file and raw output, respectively.nmode.in and NMODE.dat: Normal Mode Analysis input and output files (only for selected frames when --run-nmode is enabled).mmpbsa_<JOBID>.log:
Standard output file generated by the Slurm job scheduler
(specified using #SBATCH --output).
It records real-time execution progress, including the
percentage of completed frames and the estimated remaining time.
.sim files;
instead, it uses snapshots that have been previously exported in
PDB format (*.pdb).
The .sim files correspond to YASARA's proprietary format
for storing the complete state of a simulation scene. The underlying
AmberTools utilities executed by the script (tleap,
sander, and MMPBSA.py) only support
coordinates in standard PDB format or Amber trajectory formats
(.nc, .rst7, .mdcrd).
Performing the analysis using PDB files allows all 12,000 snapshots to be processed in a massively parallel manner on the HPC cluster, without requiring active YASARA licenses during the MM/PBSA calculation stage.
saltcon=0.150).
split_pdb_advanced filter,
which automatically detects and discards water molecules
(WAT, HOH, SOL) and
free ions (NA, CL, etc.) before
constructing the Amber topologies.
If your system contains any of the following:
In such cases, retain only these specific water molecules or ions
and remove all remaining bulk solvent. To ensure that the script includes
them in the receptor topology, specify their residue names using the
--keep-hetatm flag (for example:
--keep-hetatm MG).