{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular Dynamics Preparation\n",
    "\n",
    "This notebook demonstrates preparing Auto3D structures for MD simulations:\n",
    "\n",
    "1. **Topology generation** - GROMACS, Amber, CHARMM formats\n",
    "2. **Force field assignment** - GAFF, CGenFF, OPLS\n",
    "3. **Solvation** - box setup and ion addition\n",
    "4. **Minimization** - pre-equilibration\n",
    "\n",
    "## Workflow\n",
    "\n",
    "```\n",
    "Auto3D (3D structure) → Topology/Parameters → Solvation → Minimization → MD\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import tempfile\n",
    "import subprocess\n",
    "from pathlib import Path\n",
    "\n",
    "import Auto3D\n",
    "from Auto3D import Auto3DOptions, main\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "\n",
    "print(f\"Auto3D version: {Auto3D.__version__}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Generate Optimized Structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Small drug molecule for MD\n",
    "molecule = {\"ibuprofen\": \"CC(C)CC1=CC=C(C=C1)C(C)C(=O)O\"}\n",
    "\n",
    "with tempfile.NamedTemporaryFile(mode='w', suffix='.smi', delete=False) as f:\n",
    "    for name, smi in molecule.items():\n",
    "        f.write(f\"{smi} {name}\\n\")\n",
    "    input_file = f.name\n",
    "\n",
    "# Generate structure\n",
    "if __name__ == \"__main__\":\n",
    "    config = Auto3DOptions(\n",
    "        path=input_file,\n",
    "        k=1,\n",
    "        optimizing_engine=\"AIMNET\",\n",
    "        use_gpu=True,\n",
    "    )\n",
    "    \n",
    "    sdf_output = main(config)\n",
    "    print(f\"Output: {sdf_output}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Convert to PDB Format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sdf_to_pdb(sdf_path, pdb_path=None):\n",
    "    \"\"\"\n",
    "    Convert SDF to PDB using RDKit.\n",
    "    \"\"\"\n",
    "    if pdb_path is None:\n",
    "        pdb_path = sdf_path.replace('.sdf', '.pdb')\n",
    "    \n",
    "    mol = next(Chem.SDMolSupplier(sdf_path, removeHs=False))\n",
    "    \n",
    "    if mol is not None:\n",
    "        Chem.MolToPDBFile(mol, pdb_path)\n",
    "        return pdb_path\n",
    "    \n",
    "    return None\n",
    "\n",
    "\n",
    "if 'sdf_output' in dir():\n",
    "    pdb_output = sdf_to_pdb(sdf_output)\n",
    "    print(f\"PDB output: {pdb_output}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. GROMACS Topology with ACPYPE\n",
    "\n",
    "ACPYPE (AnteChamber PYthon Parser interfacE) generates GROMACS topologies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_gromacs_topology(pdb_path, charge=0, multiplicity=1):\n",
    "    \"\"\"\n",
    "    Generate GROMACS topology using ACPYPE.\n",
    "    \n",
    "    Requires: acpype installed (`pip install acpype` or `conda install -c conda-forge acpype`)\n",
    "    \n",
    "    Args:\n",
    "        pdb_path: Input PDB file\n",
    "        charge: Net molecular charge\n",
    "        multiplicity: Spin multiplicity (1 for singlet)\n",
    "    \"\"\"\n",
    "    output_dir = Path(pdb_path).parent / \"gromacs\"\n",
    "    output_dir.mkdir(exist_ok=True)\n",
    "    \n",
    "    cmd = f\"acpype -i {pdb_path} -n {charge} -m {multiplicity} -o gmx -b LIG\"\n",
    "    \n",
    "    print(f\"Running: {cmd}\")\n",
    "    print(\"\\nThis generates:\")\n",
    "    print(\"  - LIG.acpype/LIG_GMX.gro  (coordinates)\")\n",
    "    print(\"  - LIG.acpype/LIG_GMX.top  (topology)\")\n",
    "    print(\"  - LIG.acpype/LIG_GMX.itp  (include topology)\")\n",
    "    print(\"  - LIG.acpype/posre_LIG.itp (position restraints)\")\n",
    "    \n",
    "    try:\n",
    "        result = subprocess.run(cmd, shell=True, capture_output=True, text=True)\n",
    "        return result.stdout\n",
    "    except Exception as e:\n",
    "        return f\"Error: {e}. Install ACPYPE: pip install acpype\"\n",
    "\n",
    "\n",
    "print(\"GROMACS topology generation:\")\n",
    "print(\"-\" * 40)\n",
    "print(generate_gromacs_topology.__doc__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. GROMACS Solvation and Box Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_gromacs_box_commands(gro_file, top_file, box_size=3.0):\n",
    "    \"\"\"\n",
    "    Generate GROMACS commands for solvation.\n",
    "    \"\"\"\n",
    "    commands = f\"\"\"\n",
    "# 1. Create box (cubic, {box_size} nm buffer)\n",
    "gmx editconf -f {gro_file} -o box.gro -c -d {box_size} -bt cubic\n",
    "\n",
    "# 2. Solvate with water (TIP3P)\n",
    "gmx solvate -cp box.gro -cs spc216.gro -o solvated.gro -p {top_file}\n",
    "\n",
    "# 3. Add ions (neutralize system)\n",
    "gmx grompp -f ions.mdp -c solvated.gro -p {top_file} -o ions.tpr\n",
    "gmx genion -s ions.tpr -o system.gro -p {top_file} -pname NA -nname CL -neutral\n",
    "\n",
    "# 4. Energy minimization\n",
    "gmx grompp -f em.mdp -c system.gro -p {top_file} -o em.tpr\n",
    "gmx mdrun -v -deffnm em\n",
    "\"\"\"\n",
    "    return commands.strip()\n",
    "\n",
    "\n",
    "print(\"GROMACS solvation workflow:\")\n",
    "print(create_gromacs_box_commands(\"LIG_GMX.gro\", \"LIG_GMX.top\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_em_mdp():\n",
    "    \"\"\"\n",
    "    Create energy minimization MDP file for GROMACS.\n",
    "    \"\"\"\n",
    "    mdp = \"\"\"\n",
    "; em.mdp - Energy minimization parameters\n",
    "integrator  = steep         ; Steepest descent\n",
    "emtol       = 1000.0        ; Stop when max force < 1000 kJ/mol/nm\n",
    "emstep      = 0.01          ; Step size\n",
    "nsteps      = 50000         ; Max steps\n",
    "\n",
    "; Neighbor searching\n",
    "cutoff-scheme   = Verlet\n",
    "nstlist         = 10\n",
    "rlist           = 1.2\n",
    "\n",
    "; Electrostatics\n",
    "coulombtype     = PME\n",
    "rcoulomb        = 1.2\n",
    "\n",
    "; Van der Waals\n",
    "vdwtype         = Cut-off\n",
    "rvdw            = 1.2\n",
    "\n",
    "; Temperature/pressure (not used in EM)\n",
    "tcoupl          = no\n",
    "pcoupl          = no\n",
    "\"\"\"\n",
    "    return mdp.strip()\n",
    "\n",
    "\n",
    "print(\"Energy minimization MDP:\")\n",
    "print(create_em_mdp())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Amber Topology with Antechamber"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_amber_topology(pdb_path, charge=0):\n",
    "    \"\"\"\n",
    "    Generate Amber topology using Antechamber.\n",
    "    \n",
    "    Requires: AmberTools installed\n",
    "    \"\"\"\n",
    "    base = Path(pdb_path).stem\n",
    "    \n",
    "    commands = f\"\"\"\n",
    "# 1. Generate mol2 with AM1-BCC charges\n",
    "antechamber -i {pdb_path} -fi pdb -o {base}.mol2 -fo mol2 -c bcc -s 2 -nc {charge}\n",
    "\n",
    "# 2. Check for missing parameters\n",
    "parmchk2 -i {base}.mol2 -f mol2 -o {base}.frcmod\n",
    "\n",
    "# 3. Create Amber topology with tleap\n",
    "tleap -f - <<EOF\n",
    "source leaprc.gaff2\n",
    "LIG = loadmol2 {base}.mol2\n",
    "loadamberparams {base}.frcmod\n",
    "saveamberparm LIG {base}.prmtop {base}.inpcrd\n",
    "savepdb LIG {base}_amber.pdb\n",
    "quit\n",
    "EOF\n",
    "\"\"\"\n",
    "    return commands.strip()\n",
    "\n",
    "\n",
    "print(\"Amber topology generation:\")\n",
    "print(\"-\" * 40)\n",
    "print(generate_amber_topology(\"ligand.pdb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. OpenMM Python Interface\n",
    "\n",
    "OpenMM provides a Python API for MD simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_openmm_script(pdb_path, output_prefix=\"md\"):\n",
    "    \"\"\"\n",
    "    Create OpenMM simulation script.\n",
    "    \n",
    "    Requires: openmm, openmmforcefields\n",
    "    \"\"\"\n",
    "    script = f'''\n",
    "\"\"\"OpenMM simulation script for small molecule.\"\"\"\n",
    "from openmm.app import *\n",
    "from openmm import *\n",
    "from openmm.unit import *\n",
    "from openmmforcefields.generators import GAFFTemplateGenerator\n",
    "\n",
    "# Load structure\n",
    "pdb = PDBFile(\"{pdb_path}\")\n",
    "\n",
    "# Create force field with GAFF for small molecules\n",
    "from openff.toolkit.topology import Molecule\n",
    "molecule = Molecule.from_file(\"{pdb_path}\")\n",
    "\n",
    "gaff = GAFFTemplateGenerator(molecules=molecule)\n",
    "forcefield = ForceField(\"amber14-all.xml\", \"amber14/tip3pfb.xml\")\n",
    "forcefield.registerTemplateGenerator(gaff.generator)\n",
    "\n",
    "# Create system\n",
    "modeller = Modeller(pdb.topology, pdb.positions)\n",
    "modeller.addSolvent(forcefield, model=\"tip3p\", padding=1.0*nanometer)\n",
    "\n",
    "system = forcefield.createSystem(\n",
    "    modeller.topology,\n",
    "    nonbondedMethod=PME,\n",
    "    nonbondedCutoff=1.0*nanometer,\n",
    "    constraints=HBonds\n",
    ")\n",
    "\n",
    "# Setup simulation\n",
    "integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)\n",
    "simulation = Simulation(modeller.topology, system, integrator)\n",
    "simulation.context.setPositions(modeller.positions)\n",
    "\n",
    "# Minimize\n",
    "print(\"Minimizing...\")\n",
    "simulation.minimizeEnergy()\n",
    "\n",
    "# Equilibrate\n",
    "print(\"Equilibrating...\")\n",
    "simulation.context.setVelocitiesToTemperature(300*kelvin)\n",
    "simulation.step(10000)  # 20 ps\n",
    "\n",
    "# Production\n",
    "print(\"Running production...\")\n",
    "simulation.reporters.append(DCDReporter(\"{output_prefix}.dcd\", 1000))\n",
    "simulation.reporters.append(StateDataReporter(\"{output_prefix}.log\", 1000, \n",
    "    step=True, potentialEnergy=True, temperature=True))\n",
    "simulation.step(50000)  # 100 ps\n",
    "\n",
    "print(\"Done!\")\n",
    "'''\n",
    "    return script.strip()\n",
    "\n",
    "\n",
    "print(\"OpenMM simulation script:\")\n",
    "print(\"-\" * 40)\n",
    "print(\"Requires: pip install openmm openmmforcefields openff-toolkit\")\n",
    "print(\"\\nScript generated for ligand.pdb:\")\n",
    "print(create_openmm_script(\"ligand.pdb\")[:500] + \"...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Best Practices\n",
    "\n",
    "### Force Field Selection\n",
    "\n",
    "| Force Field | Use Case | Notes |\n",
    "|-------------|----------|-------|\n",
    "| GAFF/GAFF2 | Drug-like molecules | Amber compatible |\n",
    "| CGenFF | Drug-like molecules | CHARMM compatible |\n",
    "| OPLS-AA | General organics | Good for liquids |\n",
    "| OpenFF | Modern drug molecules | Machine-learned |\n",
    "\n",
    "### Tips\n",
    "\n",
    "1. **Protonation**: Match pH of your system\n",
    "2. **Charges**: AM1-BCC or RESP for accuracy\n",
    "3. **Box size**: At least 1.0 nm buffer for small molecules\n",
    "4. **Equilibration**: NVT then NPT before production"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This tutorial covered:\n",
    "\n",
    "1. **Auto3D → PDB** conversion\n",
    "2. **GROMACS** topology with ACPYPE\n",
    "3. **Amber** topology with Antechamber\n",
    "4. **OpenMM** Python interface\n",
    "5. **Solvation** and box setup\n",
    "\n",
    "Key workflow:\n",
    "```\n",
    "SMILES → Auto3D (optimized 3D) → PDB → Topology → Solvation → MD\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cleanup\n",
    "if 'input_file' in dir() and os.path.exists(input_file):\n",
    "    os.unlink(input_file)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}