Molecular Dynamics Preparation

This notebook demonstrates preparing Auto3D structures for MD simulations:

  1. Topology generation - GROMACS, Amber, CHARMM formats

  2. Force field assignment - GAFF, CGenFF, OPLS

  3. Solvation - box setup and ion addition

  4. Minimization - pre-equilibration

Workflow

Auto3D (3D structure) → Topology/Parameters → Solvation → Minimization → MD

[ ]:

import os
import tempfile
import subprocess
from pathlib import Path

import Auto3D
from Auto3D import Auto3DOptions, main
from rdkit import Chem
from rdkit.Chem import AllChem

print(f"Auto3D version: {Auto3D.__version__}")

1. Generate Optimized Structure

[ ]:

# Small drug molecule for MD
molecule = {"ibuprofen": "CC(C)CC1=CC=C(C=C1)C(C)C(=O)O"}

with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:
    for name, smi in molecule.items():
        f.write(f"{smi} {name}\n")
    input_file = f.name

# Generate structure
if __name__ == "__main__":
    config = Auto3DOptions(
        path=input_file,
        k=1,
        optimizing_engine="AIMNET",
        use_gpu=True,
    )

    sdf_output = main(config)
    print(f"Output: {sdf_output}")

2. Convert to PDB Format

[ ]:

def sdf_to_pdb(sdf_path, pdb_path=None):
    """
    Convert SDF to PDB using RDKit.
    """
    if pdb_path is None:
        pdb_path = sdf_path.replace('.sdf', '.pdb')

    mol = next(Chem.SDMolSupplier(sdf_path, removeHs=False))

    if mol is not None:
        Chem.MolToPDBFile(mol, pdb_path)
        return pdb_path

    return None


if 'sdf_output' in dir():
    pdb_output = sdf_to_pdb(sdf_output)
    print(f"PDB output: {pdb_output}")

3. GROMACS Topology with ACPYPE

ACPYPE (AnteChamber PYthon Parser interfacE) generates GROMACS topologies.

[ ]:

def generate_gromacs_topology(pdb_path, charge=0, multiplicity=1):
    """
    Generate GROMACS topology using ACPYPE.

    Requires: acpype installed (`pip install acpype` or `conda install -c conda-forge acpype`)

    Args:
        pdb_path: Input PDB file
        charge: Net molecular charge
        multiplicity: Spin multiplicity (1 for singlet)
    """
    output_dir = Path(pdb_path).parent / "gromacs"
    output_dir.mkdir(exist_ok=True)

    cmd = f"acpype -i {pdb_path} -n {charge} -m {multiplicity} -o gmx -b LIG"

    print(f"Running: {cmd}")
    print("\nThis generates:")
    print("  - LIG.acpype/LIG_GMX.gro  (coordinates)")
    print("  - LIG.acpype/LIG_GMX.top  (topology)")
    print("  - LIG.acpype/LIG_GMX.itp  (include topology)")
    print("  - LIG.acpype/posre_LIG.itp (position restraints)")

    try:
        result = subprocess.run(cmd, shell=True, capture_output=True, text=True)
        return result.stdout
    except Exception as e:
        return f"Error: {e}. Install ACPYPE: pip install acpype"


print("GROMACS topology generation:")
print("-" * 40)
print(generate_gromacs_topology.__doc__)

4. GROMACS Solvation and Box Setup

[ ]:

def create_gromacs_box_commands(gro_file, top_file, box_size=3.0):
    """
    Generate GROMACS commands for solvation.
    """
    commands = f"""
# 1. Create box (cubic, {box_size} nm buffer)
gmx editconf -f {gro_file} -o box.gro -c -d {box_size} -bt cubic

# 2. Solvate with water (TIP3P)
gmx solvate -cp box.gro -cs spc216.gro -o solvated.gro -p {top_file}

# 3. Add ions (neutralize system)
gmx grompp -f ions.mdp -c solvated.gro -p {top_file} -o ions.tpr
gmx genion -s ions.tpr -o system.gro -p {top_file} -pname NA -nname CL -neutral

# 4. Energy minimization
gmx grompp -f em.mdp -c system.gro -p {top_file} -o em.tpr
gmx mdrun -v -deffnm em
"""
    return commands.strip()


print("GROMACS solvation workflow:")
print(create_gromacs_box_commands("LIG_GMX.gro", "LIG_GMX.top"))

[ ]:

def create_em_mdp():
    """
    Create energy minimization MDP file for GROMACS.
    """
    mdp = """
; em.mdp - Energy minimization parameters
integrator  = steep         ; Steepest descent
emtol       = 1000.0        ; Stop when max force < 1000 kJ/mol/nm
emstep      = 0.01          ; Step size
nsteps      = 50000         ; Max steps

; Neighbor searching
cutoff-scheme   = Verlet
nstlist         = 10
rlist           = 1.2

; Electrostatics
coulombtype     = PME
rcoulomb        = 1.2

; Van der Waals
vdwtype         = Cut-off
rvdw            = 1.2

; Temperature/pressure (not used in EM)
tcoupl          = no
pcoupl          = no
"""
    return mdp.strip()


print("Energy minimization MDP:")
print(create_em_mdp())

5. Amber Topology with Antechamber

[ ]:

def generate_amber_topology(pdb_path, charge=0):
    """
    Generate Amber topology using Antechamber.

    Requires: AmberTools installed
    """
    base = Path(pdb_path).stem

    commands = f"""
# 1. Generate mol2 with AM1-BCC charges
antechamber -i {pdb_path} -fi pdb -o {base}.mol2 -fo mol2 -c bcc -s 2 -nc {charge}

# 2. Check for missing parameters
parmchk2 -i {base}.mol2 -f mol2 -o {base}.frcmod

# 3. Create Amber topology with tleap
tleap -f - <<EOF
source leaprc.gaff2
LIG = loadmol2 {base}.mol2
loadamberparams {base}.frcmod
saveamberparm LIG {base}.prmtop {base}.inpcrd
savepdb LIG {base}_amber.pdb
quit
EOF
"""
    return commands.strip()


print("Amber topology generation:")
print("-" * 40)
print(generate_amber_topology("ligand.pdb"))

6. OpenMM Python Interface

OpenMM provides a Python API for MD simulations.

[ ]:

def create_openmm_script(pdb_path, output_prefix="md"):
    """
    Create OpenMM simulation script.

    Requires: openmm, openmmforcefields
    """
    script = f'''
"""OpenMM simulation script for small molecule."""
from openmm.app import *
from openmm import *
from openmm.unit import *
from openmmforcefields.generators import GAFFTemplateGenerator

# Load structure
pdb = PDBFile("{pdb_path}")

# Create force field with GAFF for small molecules
from openff.toolkit.topology import Molecule
molecule = Molecule.from_file("{pdb_path}")

gaff = GAFFTemplateGenerator(molecules=molecule)
forcefield = ForceField("amber14-all.xml", "amber14/tip3pfb.xml")
forcefield.registerTemplateGenerator(gaff.generator)

# Create system
modeller = Modeller(pdb.topology, pdb.positions)
modeller.addSolvent(forcefield, model="tip3p", padding=1.0*nanometer)

system = forcefield.createSystem(
    modeller.topology,
    nonbondedMethod=PME,
    nonbondedCutoff=1.0*nanometer,
    constraints=HBonds
)

# Setup simulation
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)

# Minimize
print("Minimizing...")
simulation.minimizeEnergy()

# Equilibrate
print("Equilibrating...")
simulation.context.setVelocitiesToTemperature(300*kelvin)
simulation.step(10000)  # 20 ps

# Production
print("Running production...")
simulation.reporters.append(DCDReporter("{output_prefix}.dcd", 1000))
simulation.reporters.append(StateDataReporter("{output_prefix}.log", 1000,
    step=True, potentialEnergy=True, temperature=True))
simulation.step(50000)  # 100 ps

print("Done!")
'''
    return script.strip()


print("OpenMM simulation script:")
print("-" * 40)
print("Requires: pip install openmm openmmforcefields openff-toolkit")
print("\nScript generated for ligand.pdb:")
print(create_openmm_script("ligand.pdb")[:500] + "...")

7. Best Practices

Force Field Selection

Force Field

Use Case

Notes

GAFF/GAFF2

Drug-like molecules

Amber compatible

CGenFF

Drug-like molecules

CHARMM compatible

OPLS-AA

General organics

Good for liquids

OpenFF

Modern drug molecules

Machine-learned

Tips

  1. Protonation: Match pH of your system

  2. Charges: AM1-BCC or RESP for accuracy

  3. Box size: At least 1.0 nm buffer for small molecules

  4. Equilibration: NVT then NPT before production

Summary

This tutorial covered:

  1. Auto3D → PDB conversion

  2. GROMACS topology with ACPYPE

  3. Amber topology with Antechamber

  4. OpenMM Python interface

  5. Solvation and box setup

Key workflow:

SMILES → Auto3D (optimized 3D) → PDB → Topology → Solvation → MD

[ ]:

# Cleanup
if 'input_file' in dir() and os.path.exists(input_file):
    os.unlink(input_file)