Martini coarse-graining for the complex system

This notebook demonstrates how to convert the protein-ligand complex into Martini and then run a short OpenMM trajectory on the coarse-grained structure. The conversion command uses martini_openmm; the ensuing cells assume the .top/.grooutputs live under yourCOURSE_DIR/results/05-muestreo-avanzado/martini/` tree, just as the script logs.

Table of contents

• Underlying Martini math
• Step 2
• Step 3
• Step 4

Underlying Martini math

Martini lumps four heavy atoms (on average) into each bead, then rescales the interaction potentials so the coarse-grained Gibbs free energy matches the atomistic reference. The nonbonded part keeps a Lennard–Jones form

$$ V_{ ext{LJ}}(r) = 4 arepsilon\left[\left( rac{\sigma}{r} ight)^{12} - \left( rac{\sigma}{r} ight)^6 ight], $$

with σ set by the bead radius (≈0.47 nm for polar beads) and ε tuned by bead type (polar, apolar, charged). Martini applies the Lorentz–Berthelot mixing rules

$$ \sigma_{ij} =  rac{\sigma_i + \sigma_j}{2}, \quad  arepsilon_{ij} = \sqrt{ arepsilon_i  arepsilon_j}, $$

so cross terms remain smooth. Coulomb interactions use scaled charges (typically ±0.5) and the relative dielectric constant ε_r=15, yielding

$$ V_{ ext{Coul}}(r) =  rac{1}{4\pi arepsilon_0} rac{q_i q_j}{ arepsilon_r r}. $$

We also keep bonded terms (bonds, angles, dihedrals) to preserve chain stiffness; the Martini setup script writes the appropriate martini_openmm topology that specifies all these parameters.

from pathlib import Path
import os

COURSE_DIR = Path(os.environ.get('COURSE_DIR', str(Path.home() / 'Concepcion26'))).expanduser()
COMPLEX_DIR = COURSE_DIR / 'data' / 'complex'
MARTINI_DIR = COURSE_DIR / 'results' / '05-muestreo-avanzado' / 'martini'
MARTINI_DIR.mkdir(parents=True, exist_ok=True)
MARTINI_TOP = MARTINI_DIR / 'complex-martini.top'
MARTINI_GRO = MARTINI_DIR / 'complex-martini.gro'

print('Course directory:', COURSE_DIR)
print('Complex directory contents:')
for path in sorted(COMPLEX_DIR.glob('*')):
    print(' ', path.name)
print('Martini outputs:', MARTINI_TOP.name, MARTINI_GRO.name)

Course directory: /home/jordivilla/Concepcion26
Complex directory contents:
  ligand1.mol
  ligand1.sdf
  protein.pdb
  protein2-ligand2.pdb
  protein_orig.pdb
Martini outputs: complex-martini.top complex-martini.gro

Step 2

command = [
    'martini_openmm',
    '--protein', str(COMPLEX_DIR / 'protein.pdb'),
    '--ligand', str(COMPLEX_DIR / 'ligand1.mol'),
    '--output-top', str(MARTINI_TOP),
    '--output-gro', str(MARTINI_GRO),
]
print('Run this command to build the Martini topology:')
print(' ', ' '.join(command))

Run this command to build the Martini topology:
  martini_openmm --protein /home/jordivilla/Concepcion26/data/complex/protein.pdb --ligand /home/jordivilla/Concepcion26/data/complex/ligand1.mol --output-top /home/jordivilla/Concepcion26/results/05-muestreo-avanzado/martini/complex-martini.top --output-gro /home/jordivilla/Concepcion26/results/05-muestreo-avanzado/martini/complex-martini.gro

Step 3

import shutil
import subprocess

if shutil.which('martini_openmm'):
    print('martini_openmm available; running conversion now...')
    subprocess.run(command, check=True)
else:
    print('martini_openmm is not installed. Follow https://github.com/maccallumlab/martini_openmm to get it.')

martini_openmm is not installed. Follow https://github.com/maccallumlab/martini_openmm to get it.

Step 4

from openmm import app, unit
from openmm.app import GromacsTopFile, GromacsGroFile, Simulation, DCDReporter
from openmm import LangevinIntegrator

if MARTINI_TOP.exists() and MARTINI_GRO.exists():
    top = GromacsTopFile(str(MARTINI_TOP))
    gro = GromacsGroFile(str(MARTINI_GRO))
    system = top.createSystem(nonbondedMethod=app.NoCutoff, constraints=app.HBonds)
    integrator = LangevinIntegrator(300 * unit.kelvin, 1 / unit.picosecond, 20 * unit.femtoseconds)
    simulation = Simulation(top.topology, system, integrator)
    simulation.context.setPositions(gro.getPositions())
    simulation.minimizeEnergy()
    simulation.reporters.append(DCDReporter(str(MARTINI_DIR / 'martini_bg.dcd'), 100, enforcePeriodicBox=True))
    simulation.step(500)
    print('Martini simulation finished; DCD saved to', MARTINI_DIR / 'martini_bg.dcd')
else:
    print('Missing Martini files. Run the conversion first, then rerun this cell.')

Missing Martini files. Run the conversion first, then rerun this cell.