From 05bc24ba1c45a046e16b1a93b42e1f0e25ae377b Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Mon, 8 Jun 2026 17:08:38 +0200 Subject: [PATCH] Add membrane tutorial with packmol-memgen setup --- membranes/packmol_memgen_setup.ipynb | 610 +++++++++++++++++++++++++++ 1 file changed, 610 insertions(+) create mode 100644 membranes/packmol_memgen_setup.ipynb diff --git a/membranes/packmol_memgen_setup.ipynb b/membranes/packmol_memgen_setup.ipynb new file mode 100644 index 0000000..3f6b945 --- /dev/null +++ b/membranes/packmol_memgen_setup.ipynb @@ -0,0 +1,610 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "814578a2-71ad-4d3d-9004-221a0f66fe87", + "metadata": {}, + "source": [ + "# Membrane Protein RBFE with packmol-memgen and OpenFE\n", + "\n", + "This tutorial walks through setting up a Relative Binding Free Energy (RBFE) simulation for a membrane-embedded protein using:\n", + "- `packmol-memgen`: to build the lipid bilayer system\n", + "- `OpenMM`: to check if we can run a few steps of MD with the system\n", + "- `openfe`: to define and run the RBFE network\n", + "\n", + "**System**: A2A adenosine receptor (GPCR) in a pure POPC membrane\n", + "**Ligands**: Two (or more) congeneric small molecules whose relative binding affinities we wish to predict" + ] + }, + { + "cell_type": "markdown", + "id": "7bfa183b-2087-44db-8420-92e47c570cea", + "metadata": {}, + "source": [ + "## Imports" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "e4e1763e-bdb2-46ca-a2e0-c1401af9c728", + "metadata": {}, + "outputs": [], + "source": [ + "import pathlib\n", + "import subprocess\n", + "import re\n", + "from openmm import unit, Vec3\n", + "import xml.etree.ElementTree as ET\n", + "import openmmforcefields\n", + "from collections import defaultdict\n", + "import openmm\n", + "from openmm import app\n", + "from openmm.app import (\n", + " PDBFile, ForceField, Simulation, PDBReporter, \n", + " StateDataReporter, PME, HBonds, Topology, Modeller,\n", + ")\n", + "from openmm import LangevinMiddleIntegrator, MonteCarloMembraneBarostat, unit\n", + "import openmm\n", + "from sys import stdout" + ] + }, + { + "cell_type": "markdown", + "id": "97214289-05af-4a04-8e7b-3487a757af36", + "metadata": {}, + "source": [ + "## Contents\n", + "\n", + "1. [System preparation with packmol-memgen](#1-system-preparation-with-packmol-memgen)\n", + "2. [Fix common PDB issues](#2-fix-common-pdb-issues)\n", + "3. [Parameterise and solvate with OpenMM / OpenFF](#3-parameterise-and-solvate-with-openmm--openff)\n", + "4. [Equilibration protocol](#4-equilibration-protocol)\n", + "5. [Load the equilibrated system into OpenFE](#5-load-the-equilibrated-system-into-openfe)\n", + "6. [Prepare ligands and the RBFE network](#6-prepare-ligands-and-the-rbfe-network)\n", + "7. [Run the RBFE campaign](#7-run-the-rbfe-campaign)\n", + "8. [Analyse results](#8-analyse-results)" + ] + }, + { + "cell_type": "markdown", + "id": "9f1fb27a-81d1-4c31-a5fc-933d9a9e20ee", + "metadata": {}, + "source": [ + "## 1. System preparation with packmol-memgen\n", + "\n", + "We use `packmol-memgen` to embed the A2A receptor into a pure POPC bilayer. \n", + "The receptor PDB (`a2a/a2a.pdb`) should already be:\n", + "- Oriented along the membrane normal (PPM server or OPM database)\n", + "- Protonation states assigned\n", + "\n", + "> **Note:** We pass `--preoriented` because the PDB was downloaded from the OPM database, \n", + "> which already aligns the protein relative to the membrane midplane. \n", + "> Remove this flag if your structure has not been pre-oriented." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a08b8f1d-7de9-45f3-9065-684cc936837b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Processing segment 2 of 9..\n", + "Processing segment 3 of 9..\n", + "Processing segment 4 of 9\n", + "Processing segment 5 of 9\n", + "Processing segment 6 of 9\n", + "Processing segment 7 of 9\n", + "Processing segment 8 of 9\n", + "Processing segment 9 of 9\u001b[F\n", + "All-together Packing..........\n", + "########################################################################################################################################\n", + "########################################################################################################################################\n", + "DONE!\n", + "\n" + ] + } + ], + "source": [ + "cmd = [\n", + " \"packmol-memgen\",\n", + " \"--pdb\", \"a2a/a2a.pdb\",\n", + " \"--lipids\", \"POPC\",\n", + " \"--preoriented\",\n", + " \"--dist\", \"17\", # minimum distance (Å) from protein to box edge in xy\n", + " \"--dist_wat\", \"15\", # water layer thickness (Å) on each side\n", + " \"--salt\",\n", + " \"--salt_c\", \"K+\", # physiological cation\n", + " \"--saltcon\", \"0.15\", # 150 mM KCl\n", + " \"--nottrim\", # keep all atoms; don't trim lipids that overlap with protein\n", + " \"--overwrite\",\n", + " \"--notprotonate\", # protonation already handled\n", + " \"--charmm\", # write CHARMM-compatible PDB (lipid as single residue)\n", + "]\n", + "\n", + "result = subprocess.run(cmd, capture_output=True, text=True)\n", + "print(result.stdout[-3000:]) # last 3000 chars of stdout\n", + "if result.returncode != 0:\n", + " print(\"STDERR:\", result.stderr[-2000:])" + ] + }, + { + "cell_type": "markdown", + "id": "1adc9fdc-e7c1-428c-a757-cc6f1ed6fed4", + "metadata": {}, + "source": [ + "### Key flags explained\n", + "\n", + "| Flag | Purpose |\n", + "|---|---|\n", + "| `--preoriented` | Skip PPM alignment; the PDB is already membrane-oriented |\n", + "| `--dist 17` | ~17 Å padding between protein and periodic box edge in the xy-plane |\n", + "| `--dist_wat 15` | Water thickness above/below the bilayer |\n", + "| `--nottrim` | Prevents clipping of lipid tails that partially overlap the protein; important for GPCRs with deep binding pockets |\n", + "| `--charmm` | Writes a CHARMM-formatted PDB; atom naming compatible with CHARMM36 lipid FF |\n", + "| `--notprotonate` | Skip internal protonation step (we have already handled this) |\n", + "\n", + "The main output we need is **`bilayer_a2a_prepped.pdb`**. \n", + "packmol-memgen also writes a log file containing the box vectors — we will need these shortly." + ] + }, + { + "cell_type": "markdown", + "id": "49cd9de0-456e-4dc8-8c1d-331aedd0d8b0", + "metadata": {}, + "source": [ + "## 2. Fix common PDB issues and prepare the topology\n", + "\n", + "### 2a. Parse box dimensions from the log file and patch the CRYST1 record\n", + "packmol-memgen writes box dimensions to the log file but not to the CRYST1 record. Extract the dimensions from the log and patch them in before doing anything else." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "ae62cc19-c008-4660-ade5-7fb943353806", + "metadata": {}, + "outputs": [], + "source": [ + "def parse_box_from_log(log_path: str) -> tuple[float, float, float]:\n", + " \"\"\"\n", + " Extract box dimensions (Å) from a packmol-memgen log file.\n", + " \n", + " Looks for lines of the form:\n", + " x_len = 88.52550719674292\n", + " y_len = 88.52550719674292\n", + " z_len = 108.602\n", + " \"\"\"\n", + " dims = {}\n", + " pattern = re.compile(r\"(x_len|y_len|z_len)\\s*=\\s*(\\d+\\.\\d+)\")\n", + " with open(log_path) as fh:\n", + " for line in fh:\n", + " m = pattern.search(line)\n", + " if m:\n", + " dims[m.group(1)] = float(m.group(2))\n", + " if len(dims) == 3:\n", + " break\n", + " if len(dims) < 3:\n", + " raise ValueError(\n", + " f\"Could not find x_len/y_len/z_len in {log_path}. \"\n", + " f\"Found: {dims}\"\n", + " )\n", + " return dims[\"x_len\"], dims[\"y_len\"], dims[\"z_len\"]\n", + "\n", + "def patch_cryst1(pdb_in: str, pdb_out: str, a: float, b: float, c: float) -> None:\n", + " \"\"\"Write / overwrite the CRYST1 record with correct box vectors.\"\"\"\n", + " cryst1 = f\"CRYST1{a:9.3f}{b:9.3f}{c:9.3f} 90.00 90.00 90.00 P 1 1\\n\"\n", + " lines = open(pdb_in).readlines()\n", + " lines = [l for l in lines if not l.startswith(\"CRYST1\")]\n", + " lines.insert(0, cryst1)\n", + " with open(pdb_out, \"w\") as fh:\n", + " fh.writelines(lines)\n", + " print(f\"Written {pdb_out} with CRYST1: a={a}, b={b}, c={c}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "3ab57fef-f4bf-4d0f-916a-dc4382f7fc59", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Written a2a/bilayer_a2a_fixed.pdb with CRYST1: a=88.52550719674292, b=88.52550719674292, c=108.602\n" + ] + } + ], + "source": [ + "a, b, c = parse_box_from_log(\"a2a/packmol-memgen.log\")\n", + "patch_cryst1(\"a2a/bilayer_a2a.pdb\", \"a2a/bilayer_a2a_fixed.pdb\", a, b, c)" + ] + }, + { + "cell_type": "markdown", + "id": "57396e12-e567-4810-bd76-c5b9192c40fb", + "metadata": {}, + "source": [ + "### 2b. Fix the element column for ions\n", + "\n", + "The last column of the PDB ATOM/HETATM record may contain a single-character element symbol instead of the correct two-character one. For example, the sodium ion line looks like:\n", + "\n", + "`ATOM 4654 Na NA A 298 -2.406 -0.748 -2.531 1.00 0.00 N`\n", + "\n", + "The element column reads N instead of NA, causing OpenMM to misidentify the ion as nitrogen." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "06347612-7795-46eb-8476-7368a135454b", + "metadata": {}, + "outputs": [], + "source": [ + "def fix_ion_elements(pdb_in: str, pdb_out: str) -> None:\n", + " \"\"\"Correct the element column (cols 77-78) for ions that packmol-memgen\n", + " writes with a truncated element symbol.\n", + "\n", + " Args:\n", + " pdb_in: Path to the input PDB file.\n", + " pdb_out: Path to the output PDB file.\n", + " \"\"\"\n", + " # (residue_name, atom_name) → correct two-character element string\n", + " fixes = {\n", + " (\"NA\", \"Na\"): \"NA\",\n", + " (\"CL\", \"Cl\"): \"CL\",\n", + " (\"K+\", \"K\" ): \" K\",\n", + " (\"MG\", \"Mg\"): \"MG\",\n", + " }\n", + " corrected = 0\n", + " out_lines = []\n", + " with open(pdb_in) as fh:\n", + " for line in fh:\n", + " if line.startswith((\"ATOM\", \"HETATM\")) and len(line) >= 78:\n", + " res_name = line[17:20].strip()\n", + " atom_name = line[12:16].strip()\n", + " key = (res_name, atom_name)\n", + " if key in fixes and line[76:78].strip() != fixes[key].strip():\n", + " line = line[:76] + fixes[key] + line[78:]\n", + " corrected += 1\n", + " out_lines.append(line)\n", + " with open(pdb_out, \"w\") as fh:\n", + " fh.writelines(out_lines)\n", + " print(f\"Fixed element column for {corrected} ion records → {pdb_out}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "6f5c84fb-9767-4237-94c7-c8cfd6a114f4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Fixed element column for 1 ion records → a2a/bilayer_a2a_fixed.pdb\n" + ] + } + ], + "source": [ + "fix_ion_elements(\"a2a/bilayer_a2a_fixed.pdb\", \"a2a/bilayer_a2a_fixed.pdb\")" + ] + }, + { + "cell_type": "markdown", + "id": "328ebd41-8e1d-4696-a8eb-30af0c39a3b6", + "metadata": {}, + "source": [ + "### 2d. Add lipid bond information and prepare the topology\n", + "\n", + "packmol-memgen does not write CONECT records for lipid residues, which is a known issue. Lipids are heterogens and are required by the PDB spec to use CONECT records and files that omit them are nonstandard.\n", + "The fix is to load the bond definitions from the Amber lipid17 force field XML before reading the PDB, so that `createStandardBonds()` can apply them. There is one complication: lipid17_merged.xml stores bonds with atomName1/atomName2 attributes, but Topology.loadBondDefinitions() expects from/to attributes. We therefore need to convert the format first. See https://github.com/openmm/openmm/issues/4997#issuecomment-3290668970 for a discussion on that." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "1e3caf73-5ec0-4a88-a892-ba80d4139b0b", + "metadata": {}, + "outputs": [], + "source": [ + "_RESNAME_FIXES: dict[str, str] = {\n", + " \"POP\": \"POPC\",\n", + "}\n", + "\n", + "def extract_bond_info(xml_path: str, output_path: str) -> None:\n", + " \"\"\"Convert bond definitions from Amber FF XML format to the\n", + " from/to format expected by Topology.loadBondDefinitions().\n", + "\n", + " Amber FF XML uses: \n", + " loadBondDefinitions expects: \n", + "\n", + " Args:\n", + " xml_path: Path to the Amber FF XML (e.g. lipid17_merged.xml).\n", + " output_path: Path to write the converted XML.\n", + " \"\"\"\n", + " tree = ET.parse(xml_path)\n", + " root = tree.getroot()\n", + " residues = root.find(\"Residues\")\n", + " with open(output_path, \"w\") as f:\n", + " f.write(\"\\n\")\n", + " for residue in residues.findall(\"Residue\"):\n", + " res_name = residue.get(\"name\")\n", + " f.write(f' \\n')\n", + " for bond in residue.findall(\"Bond\"):\n", + " atom1 = bond.get(\"atomName1\")\n", + " atom2 = bond.get(\"atomName2\")\n", + " f.write(f' \\n')\n", + " f.write(' \\n')\n", + " f.write(\"\\n\")\n", + " print(f\"Bond definitions saved to: {output_path}\")\n", + "\n", + "\n", + "def prepare_membrane_topology(\n", + " input_pdb: str,\n", + " output_pdb: str,\n", + " verbose: bool = True,\n", + ") -> None:\n", + " \"\"\"Prepare a packmol-memgen --charmm PDB for use with OpenMM.\n", + "\n", + " Topology.loadBondDefinitions() must be called with the converted lipid\n", + " bond XML before calling this function, so that createStandardBonds()\n", + " can find and apply the lipid bond templates.\n", + "\n", + " This function also renames POP → POPC to match the force field template\n", + " names, since packmol-memgen truncates POPC to POP in CHARMM output.\n", + "\n", + " Args:\n", + " input_pdb: Path to the input PDB (with CRYST1 already patched).\n", + " output_pdb: Path to write the corrected PDB.\n", + " verbose: If True, print a per-residue-type bond count summary.\n", + " \"\"\"\n", + " pdb = PDBFile(input_pdb)\n", + "\n", + " if pdb.topology.getPeriodicBoxVectors() is None:\n", + " raise ValueError(\n", + " \"No periodic box vectors found. \"\n", + " \"Run patch_cryst1() before calling prepare_membrane_topology().\"\n", + " )\n", + "\n", + " for res in pdb.topology.residues():\n", + " if res.name in _RESNAME_FIXES:\n", + " res.name = _RESNAME_FIXES[res.name]\n", + "\n", + " modeller = Modeller(pdb.topology, pdb.positions)\n", + " modeller.topology.createStandardBonds()\n", + "\n", + " with open(output_pdb, \"w\") as fh:\n", + " PDBFile.writeFile(modeller.topology, modeller.positions, fh, keepIds=True)\n", + "\n", + " if verbose:\n", + " counts: dict[str, int] = defaultdict(int)\n", + " for bond in modeller.topology.bonds():\n", + " counts[bond[0].residue.name] += 1\n", + " print(f\"{'Residue':<8} {'Total bonds':>12}\")\n", + " print(\"-\" * 22)\n", + " for name, count in sorted(counts.items()):\n", + " print(f\"{name:<8} {count:>12}\")\n", + " print(f\"\\nSaved : {output_pdb}\")" + ] + }, + { + "cell_type": "markdown", + "id": "566c4f87-e783-4bbf-961e-334bf882ac10", + "metadata": {}, + "source": [ + "Now run the full pipeline. The order matters: `loadBondDefinitions` must be called before PDBFile:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "af484c5e-e74d-45cd-8beb-d2bbfa9c5afb", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Bond definitions saved to: a2a/lipid17_bonds.xml\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/hannahbaumann/.local/share/mamba/envs/openfe_env/lib/python3.12/site-packages/openmm/app/internal/pdbstructure.py:441: UserWarning: WARNING: two consecutive residues with same number (ATOM 3199 N NME A 210 17.633 8.568 -23.707 1.00 0.00 N , ATOM 3198 HD23 LEU A 210 13.213 5.988 -24.001 1.00 0.00 H )\n", + " warnings.warn(\"WARNING: two consecutive residues with same number (%s, %s)\" % (atom, self._current_residue.atoms[-1]))\n", + "/Users/hannahbaumann/.local/share/mamba/envs/openfe_env/lib/python3.12/site-packages/openmm/app/internal/pdbstructure.py:441: UserWarning: WARNING: two consecutive residues with same number (ATOM 3211 N GLU A 211 9.982 5.977 -38.597 1.00 0.00 N , ATOM 3210 HH33 ACE A 211 11.950 4.338 -37.817 1.00 0.00 H )\n", + " warnings.warn(\"WARNING: two consecutive residues with same number (%s, %s)\" % (atom, self._current_residue.atoms[-1]))\n", + "/Users/hannahbaumann/.local/share/mamba/envs/openfe_env/lib/python3.12/site-packages/openmm/app/internal/pdbstructure.py:441: UserWarning: WARNING: two consecutive residues with same number (ATOM 4648 N NME A 297 -24.566 0.170 -17.915 1.00 0.00 N , ATOM 4647 HG SER A 297 -23.441 -0.799 -19.073 1.00 0.00 H )\n", + " warnings.warn(\"WARNING: two consecutive residues with same number (%s, %s)\" % (atom, self._current_residue.atoms[-1]))\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Residue Total bonds\n", + "----------------------\n", + "ACE 6\n", + "ALA 390\n", + "ARG 312\n", + "ASN 182\n", + "ASP 48\n", + "CYS 150\n", + "GLN 136\n", + "GLU 105\n", + "GLY 128\n", + "HIS 110\n", + "HOH 32374\n", + "ILE 551\n", + "LEU 684\n", + "LYS 110\n", + "MET 102\n", + "NME 10\n", + "PHE 399\n", + "POPC 28462\n", + "PRO 180\n", + "SER 165\n", + "THR 154\n", + "TRP 156\n", + "TYR 220\n", + "VAL 416\n", + "\n", + "Saved : a2a/bilayer_a2a_prepped.pdb\n" + ] + } + ], + "source": [ + "# Find the lipid17 XML bundled with openmmforcefields\n", + "ff_xml = str(\n", + " pathlib.Path(openmmforcefields.__file__).parent\n", + " / \"ffxml\" / \"amber\" / \"lipid17_merged.xml\"\n", + ")\n", + "\n", + "# Step 1: convert atomName1/atomName2 → from/to format\n", + "extract_bond_info(ff_xml, \"a2a/lipid17_bonds.xml\")\n", + "\n", + "# Step 2: register with OpenMM — must happen before PDBFile() is called\n", + "Topology.loadBondDefinitions(\"a2a/lipid17_bonds.xml\")\n", + "\n", + "# Step 3: rename residues, build bonds, write corrected PDB\n", + "prepare_membrane_topology(\n", + " input_pdb=\"a2a/bilayer_a2a_fixed.pdb\",\n", + " output_pdb=\"a2a/bilayer_a2a_prepped.pdb\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "41fe7285-54e7-4989-b5be-86aef224888e", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/hannahbaumann/.local/share/mamba/envs/openfe_env/lib/python3.12/site-packages/openmm/app/internal/pdbstructure.py:441: UserWarning: WARNING: two consecutive residues with same number (HETATM 3199 N NME A 210 17.633 8.568 -23.707 1.00 0.00 N , ATOM 3198 HD23 LEU A 210 13.213 5.988 -24.001 1.00 0.00 H )\n", + " warnings.warn(\"WARNING: two consecutive residues with same number (%s, %s)\" % (atom, self._current_residue.atoms[-1]))\n", + "/Users/hannahbaumann/.local/share/mamba/envs/openfe_env/lib/python3.12/site-packages/openmm/app/internal/pdbstructure.py:441: UserWarning: WARNING: two consecutive residues with same number (ATOM 3212 N GLU A 211 9.982 5.977 -38.597 1.00 0.00 N , HETATM 3211 H3 ACE A 211 11.950 4.338 -37.817 1.00 0.00 H )\n", + " warnings.warn(\"WARNING: two consecutive residues with same number (%s, %s)\" % (atom, self._current_residue.atoms[-1]))\n", + "/Users/hannahbaumann/.local/share/mamba/envs/openfe_env/lib/python3.12/site-packages/openmm/app/internal/pdbstructure.py:441: UserWarning: WARNING: two consecutive residues with same number (HETATM 4649 N NME A 297 -24.566 0.170 -17.915 1.00 0.00 N , ATOM 4648 HG SER A 297 -23.441 -0.799 -19.073 1.00 0.00 H )\n", + " warnings.warn(\"WARNING: two consecutive residues with same number (%s, %s)\" % (atom, self._current_residue.atoms[-1]))\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Minimizing energy...\n", + "Potential energy after minimisation: -936771.6 kJ/mol\n", + "Running short NPT (500 steps)...\n" + ] + } + ], + "source": [ + "pdb = PDBFile('a2a/bilayer_a2a_prepped.pdb')\n", + "\n", + "forcefield = ForceField(\n", + " 'amber/protein.ff14SB.xml',\n", + " 'amber/lipid17_merged.xml',\n", + " 'amber/tip3p_standard.xml',\n", + ")\n", + "\n", + "system = forcefield.createSystem(\n", + " pdb.topology,\n", + " nonbondedMethod=PME,\n", + " nonbondedCutoff=1.0 * unit.nanometer,\n", + " constraints=HBonds,\n", + " rigidWater=True,\n", + ")\n", + "\n", + "# MonteCarloMembraneBarostat is required for lipid bilayers —\n", + "# it allows independent xy and z fluctuations.\n", + "# Surface tension of 0 bar·nm = tensionless ensemble.\n", + "system.addForce(MonteCarloMembraneBarostat(\n", + " 1.0 * unit.bar,\n", + " 0.0 * unit.bar * unit.nanometer,\n", + " 300 * unit.kelvin,\n", + " MonteCarloMembraneBarostat.XYIsotropic,\n", + " MonteCarloMembraneBarostat.ZFree,\n", + " 50,\n", + "))\n", + "\n", + "# 4 fs timestep enabled by HMR above\n", + "integrator = LangevinMiddleIntegrator(\n", + " 300 * unit.kelvin,\n", + " 1.0 / unit.picosecond,\n", + " 2.0 * unit.femtoseconds,\n", + ")\n", + "\n", + "simulation = Simulation(pdb.topology, system, integrator)\n", + "simulation.context.setPositions(pdb.positions)\n", + "simulation.context.setPeriodicBoxVectors(*pdb.topology.getPeriodicBoxVectors())\n", + "\n", + "# Energy minimisation — membrane systems often have clashes from packing\n", + "print(\"Minimizing energy...\")\n", + "simulation.minimizeEnergy(maxIterations=1000)\n", + "\n", + "# Check the energy is sensible after minimisation — NaN or very large\n", + "# values (> 1e6 kJ/mol) indicate something is wrong\n", + "state = simulation.context.getState(getEnergy=True)\n", + "e = state.getPotentialEnergy().value_in_unit(unit.kilojoules_per_mole)\n", + "print(f\"Potential energy after minimisation: {e:.1f} kJ/mol\")\n", + "if e > 1e7 or e != e: # NaN check: NaN != NaN is True\n", + " raise RuntimeError(f\"Energy looks wrong after minimisation: {e} kJ/mol\")\n", + "\n", + "# Short NPT run \n", + "print(\"Running short NPT (500 steps)...\")\n", + "simulation.step(500)\n", + "\n", + "# Save final state\n", + "state = simulation.context.getState(getPositions=True, enforcePeriodicBox=True)\n", + "positions = state.getPositions()\n", + "box_vectors = state.getPeriodicBoxVectors()\n", + "pdb.topology.setPeriodicBoxVectors(box_vectors)\n", + "\n", + "with open('a2a/test_output.pdb', 'w') as f:\n", + " PDBFile.writeFile(pdb.topology, positions, f)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11ee11c3-992e-42e2-b93a-38e92b09ddf4", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}