millores

Author: JordiVillaFreixa

docs/episodes/04-muestreo-avanzado.md

@@ -12,19 +12,33 @@ permalink: /episodes/04-muestreo-avanzado/
 
 <!-- toc:start -->
 ## Table of contents
+- [Table of contents](#table-of-contents)
 - [Duration](#duration)
 - [Objectives](#objectives)
 - [Content](#content)
 - [Thermodynamic foundations](#thermodynamic-foundations)
 - [Official OpenMM tutorials](#official-openmm-tutorials)
 - [Simulated annealing](#simulated-annealing)
+  - [Guided demo](#guided-demo)
+  - [Exercise](#exercise)
+  - [Notebooks](#notebooks)
+- [Replica exchange (REMD)](#replica-exchange-remd)
+  - [REMD for alanine dipeptide](#remd-for-alanine-dipeptide)
+  - [Exercise](#exercise-1)
+  - [Notebooks](#notebooks-1)
 - [Gaussian accelerated MD (GaMD) with OpenMM](#gaussian-accelerated-md-gamd-with-openmm)
+  - [Guided demo: GaMD on alanine dipeptide](#guided-demo-gamd-on-alanine-dipeptide)
+  - [Exercise](#exercise-2)
+  - [Notebooks](#notebooks-2)
 - [Martini coarse-graining with OpenMM](#martini-coarse-graining-with-openmm)
+  - [Guided demo: Martini coarse-graining for the complex system](#guided-demo-martini-coarse-graining-for-the-complex-system)
+  - [Exercise](#exercise-3)
+  - [Notebooks](#notebooks-3)
 - [Umbrella sampling](#umbrella-sampling)
 - [SBMOpenMM resources](#sbmopenmm-resources)
 - [Guide scripts (OpenMM Application Layer)](#guide-scripts-openmm-application-layer)
-- [Alanine dipeptide](#alanine-dipeptide)
-- [Protein-ligand complex](#protein-ligand-complex)
+  - [Alanine dipeptide](#alanine-dipeptide)
+  - [Protein-ligand complex](#protein-ligand-complex)
 - [References](#references)
 <!-- toc:end -->
 
@@ -103,7 +117,7 @@ Script source: <a href="{{ site.baseurl }}/episodes/scripts/openmm_advanced_anne
 
 The accompanying explanation notes that the loop simply updates the `LangevinMiddleIntegrator` temperature before each batch of 1,000 steps, which mirrors the ramping schedules we describe here. Keep that block as a quick reference when tuning integrator parameters or defining temperature sequences for your exercises so the practical code stays aligned with the theory in this episode.
 
-### Replica exchange (REMD) context
+## Replica exchange (REMD)
 
 The same temperature-ramping intuition underpins replica-exchange (REMD) methods. Instead of moving a single trajectory between high and low temperatures, REMD maintains multiple replicas at different thermodynamic states and swaps them through Metropolis trials to jump across barriers. `OpenMMTools` exposes `ReplicaExchangeSampler` (see https://openmmtools.readthedocs.io/en/stable/multistate.html#replicaexchangesampler-replica-exchange-among-thermodynamic-states)
 to manage state definitions, collect swap statistics, and maintain detailed balance for arbitrary Hamiltonians.
@@ -211,256 +225,24 @@ To feed the analysis, we rely on the official scripts described in the OpenMM gu
 
 In each case, the outputs (DCD, CSV, and energy reports) are used in this episode's exercises.
 
-## Alanine dipeptide
+### Alanine dipeptide
 
-### Guided demo
-
-<!-- sync-from: docs/episodes/scripts/05-muestreo-avanzado_simple.py -->
-```python
-#!/usr/bin/env python3
-import os
-from pathlib import Path
-
-from openmm import unit, app
-import openmm as mm
-from openmm.app import PDBFile, ForceField, Simulation
-
-COURSE_DIR = Path(os.environ.get("COURSE_DIR", str(Path.home() / "Concepcion26"))).expanduser()
-DATA_DIR = COURSE_DIR / "data"
-PDB_IN = DATA_DIR / "alanine-dipeptide.pdb"
-OUT_DIR = COURSE_DIR / "results" / "05-muestreo-avanzado" / "simple"
-OUT_DIR.mkdir(parents=True, exist_ok=True)
-
-pdb = PDBFile(str(PDB_IN))
-forcefield = ForceField("amber14-all.xml", "amber14/tip3pfb.xml")
-
-system = forcefield.createSystem(
-    pdb.topology,
-    nonbondedMethod=app.NoCutoff,
-    constraints=app.HBonds,
-)
-
-# Constrain the distance between two atoms as a simple example.
-force = mm.CustomBondForce("0.5*k*(r-r0)^2")
-force.addPerBondParameter("k")
-force.addPerBondParameter("r0")
-force.addBond(0, 1, [500.0 * unit.kilojoule_per_mole / unit.nanometer**2, 0.25 * unit.nanometer])
-system.addForce(force)
-
-integrator = mm.LangevinIntegrator(300 * unit.kelvin, 1 / unit.picosecond, 2 * unit.femtoseconds)
-
-simulation = Simulation(pdb.topology, system, integrator)
-simulation.context.setPositions(pdb.positions)
-
-simulation.minimizeEnergy(maxIterations=200)
-simulation.step(2000)
-
-print("Simulation with restraint finished. Output dir:", OUT_DIR)
-```
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py">05-muestreo-avanzado_simple.py</a>
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py">05-muestreo-avanzado_simple.py</a>
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py">05-muestreo-avanzado_simple.py</a>
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py">05-muestreo-avanzado_simple.py</a>
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py">05-muestreo-avanzado_simple.py</a>
 
 Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py">05-muestreo-avanzado_simple.py</a>
 
-### Exercise
-
-- Run the simple workflow with spherical restraints and record the value of $\langle \|\mathbf{g}\|^2 \rangle$ before and after annealing.
-- Compare the average potential energy with and without CustomExternalForce.
-
-### Key points
-
-- Keep reports every 10 steps to monitor integrator drift.
-- Adjust $\beta$ to preserve numerical stability in small systems.
-
-### Notebooks and scripts
+<!-- sync-from: docs/episodes/notebooks/04-muestreo-avanzado_simple.ipynb -->
+<div class="notebook-embed"><iframe src="{{ site.baseurl }}/episodes/notebooks/rendered/04-muestreo-avanzado_simple.html" loading="lazy"></iframe><div class="notebook-links"><a href="{{ site.baseurl }}/episodes/notebooks/04-muestreo-avanzado_simple.ipynb" download>Download notebook</a> | <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py" download>Download script (.py)</a></div></div>
 
 - This notebook executes the advanced sampling routine for alanine (restraints, custom forces, reporting energies) and tracks gradient norms. (<a href="{{ site.baseurl }}/episodes/notebooks/04-muestreo-avanzado_simple.ipynb">notebook</a> | <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado_simple.py">script</a>)
 
-## Protein-ligand complex
-
-### Guided demo
-
-<!-- sync-from: docs/episodes/scripts/05-muestreo-avanzado.py -->
-```python
-#!/usr/bin/env python3
-import argparse
-import os
-import sys
-import time
-from pathlib import Path
-
-from openff.toolkit import Molecule
-from openmmforcefields.generators import SystemGenerator
-import openmm
-from openmm import app, unit, LangevinIntegrator
-from openmm.app import PDBFile, Simulation, Modeller, StateDataReporter, DCDReporter
-
-COURSE_DIR = Path(os.environ.get("COURSE_DIR", str(Path.home() / "Concepcion26"))).expanduser()
-DATA_DIR = COURSE_DIR / "data" / "complex"
-DEFAULT_PROTEIN = DATA_DIR / "protein.pdb"
-DEFAULT_LIGAND = DATA_DIR / "ligand1.mol"
-
-
-def get_platform():
-    speed = 0
-    platform = None
-    for i in range(openmm.Platform.getNumPlatforms()):
-        candidate = openmm.Platform.getPlatform(i)
-        if candidate.getSpeed() > speed:
-            platform = candidate
-            speed = candidate.getSpeed()
-    print("Using platform", platform.getName())
-    if platform.getName() in {"CUDA", "OpenCL"}:
-        platform.setPropertyDefaultValue("Precision", "mixed")
-        print("Set precision for platform", platform.getName(), "to mixed")
-    return platform
-
-
-def main() -> None:
-    parser = argparse.ArgumentParser(description="Simulate protein-ligand complex with optional solvation")
-    parser.add_argument("-p", "--protein", default=str(DEFAULT_PROTEIN), help="Protein PDB file")
-    parser.add_argument("-l", "--ligand", default=str(DEFAULT_LIGAND), help="Ligand MOL file")
-    parser.add_argument("-o", "--output", default="solvated", help="Base name for output files")
-    parser.add_argument("-s", "--steps", type=int, default=5000, help="Number of steps")
-    parser.add_argument("-z", "--step-size", type=float, default=0.002, help="Step size (ps)")
-    parser.add_argument("-f", "--friction-coeff", type=float, default=1.0, help="Friction coefficient (ps)")
-    parser.add_argument("-i", "--interval", type=int, default=1000, help="Reporting interval")
-    parser.add_argument("-t", "--temperature", type=int, default=300, help="Temperature (K)")
-    parser.add_argument("--solvate", action="store_true", help="Add solvent box")
-    parser.add_argument("--padding", type=float, default=10.0, help="Padding for solvent box (A)")
-    parser.add_argument("--water-model", default="tip3p", choices=["tip3p", "spce", "tip4pew", "tip5p", "swm4ndp"], help="Water model")
-    parser.add_argument("--positive-ion", default="Na+", help="Positive ion for solvation")
-    parser.add_argument("--negative-ion", default="Cl-", help="Negative ion for solvation")
-    parser.add_argument("--ionic-strength", type=float, default=0.0, help="Ionic strength (M)")
-    parser.add_argument("--no-neutralize", action="store_true", help="Don't neutralize")
-    parser.add_argument("-e", "--equilibration-steps", type=int, default=200, help="Equilibration steps")
-    args = parser.parse_args()
-
-    t0 = time.time()
-    out_dir = COURSE_DIR / "results" / "05-muestreo-avanzado" / "complex"
-    out_dir.mkdir(parents=True, exist_ok=True)
-    output_base = str(out_dir / args.output)
-    output_complex = output_base + "_complex.pdb"
-    output_min = output_base + "_minimised.pdb"
-    output_traj = output_base + "_traj.dcd"
-
-    print("Reading ligand")
-    ligand_mol = Molecule.from_file(args.ligand)
-
-    print("Preparing system")
-    forcefield_kwargs = {
-        "constraints": app.HBonds,
-        "rigidWater": True,
-        "removeCMMotion": False,
-        "hydrogenMass": 4 * unit.amu,
-    }
-    system_generator = SystemGenerator(
-        forcefields=["amber/ff14SB.xml", "amber/tip3p_standard.xml"],
-        small_molecule_forcefield="gaff-2.11",
-        molecules=[ligand_mol],
-        forcefield_kwargs=forcefield_kwargs,
-    )
-
-    print("Reading protein")
-    protein_pdb = PDBFile(args.protein)
-
-    print("Preparing complex")
-    modeller = Modeller(protein_pdb.topology, protein_pdb.positions)
-    lig_top = ligand_mol.to_topology()
-    modeller.add(lig_top.to_openmm(), lig_top.get_positions().to_openmm())
-
-    if args.solvate:
-        print("Adding solvent")
-        modeller.addSolvent(
-            system_generator.forcefield,
-            model=args.water_model,
-            padding=args.padding * unit.angstroms,
-            positiveIon=args.positive_ion,
-            negativeIon=args.negative_ion,
-            ionicStrength=args.ionic_strength * unit.molar,
-            neutralize=not args.no_neutralize,
-        )
-
-    with open(output_complex, "w") as outfile:
-        PDBFile.writeFile(modeller.topology, modeller.positions, outfile)
-
-    system = system_generator.create_system(modeller.topology, molecules=ligand_mol)
-    step_size = args.step_size * unit.picoseconds
-    friction = args.friction_coeff / unit.picosecond
-    temperature = args.temperature * unit.kelvin
-    duration = (step_size * args.steps).value_in_unit(unit.nanoseconds)
-
-    if system.usesPeriodicBoundaryConditions():
-        system.addForce(openmm.MonteCarloBarostat(1 * unit.atmospheres, temperature, 25))
-
-    integrator = LangevinIntegrator(temperature, friction, step_size)
-    platform = get_platform()
-
-    simulation = Simulation(modeller.topology, system, integrator, platform=platform)
-    simulation.context.setPositions(modeller.positions)
-
-    print("Minimising ...")
-    simulation.minimizeEnergy()
-
-    with open(output_min, "w") as outfile:
-        PDBFile.writeFile(
-            modeller.topology,
-            simulation.context.getState(getPositions=True, enforcePeriodicBox=True).getPositions(),
-            file=outfile,
-            keepIds=True,
-        )
-
-    simulation.context.setVelocitiesToTemperature(temperature)
-    print("Equilibrating ...")
-    simulation.step(args.equilibration_steps)
-
-    simulation.reporters.append(DCDReporter(output_traj, args.interval, enforcePeriodicBox=True))
-    simulation.reporters.append(StateDataReporter(sys.stdout, args.interval * 5, step=True, potentialEnergy=True, temperature=True))
-
-    print("Starting simulation with", args.steps, "steps ...")
-    t1 = time.time()
-    simulation.step(args.steps)
-    t2 = time.time()
-    print("Simulation complete in", round((t2 - t1) / 60, 3), "mins")
-    print("Simulation time was", round(duration, 3), "ns")
-    print("Total wall clock time was", round((t2 - t0) / 60, 3), "mins")
-
-
-if __name__ == "__main__":
-    main()
-```
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py">05-muestreo-avanzado.py</a>
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py">05-muestreo-avanzado.py</a>
-
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py">05-muestreo-avanzado.py</a>
+### Protein-ligand complex
 
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py">05-muestreo-avanzado.py</a>
 
 Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py">05-muestreo-avanzado.py</a>
 
-Script source: <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py">05-muestreo-avanzado.py</a>
-
-### Exercise
-
-- Run with `--solvate` and `--padding 12`, recording solvent density and Coulomb energy.
-- Try another water model and analyze the variation in $\Delta G$ estimated from partition differences.
-
-### Key points
-
-- Periodic solvation stabilizes $\Delta G$ on long time scales and improves annealing convergence.
-- Compare energy distributions to validate extended sampling.
+<!-- sync-from: docs/episodes/notebooks/04-muestreo-avanzado.ipynb -->
+<div class="notebook-embed"><iframe src="{{ site.baseurl }}/episodes/notebooks/rendered/04-muestreo-avanzado.html" loading="lazy"></iframe><div class="notebook-links"><a href="{{ site.baseurl }}/episodes/notebooks/04-muestreo-avanzado.ipynb" download>Download notebook</a> | <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py" download>Download script (.py)</a></div></div>
 
-### Notebooks and scripts
 
 - This notebook covers the complex system with solvation options, barostat control, and energy/force logging to validate enhanced sampling workflows. (<a href="{{ site.baseurl }}/episodes/notebooks/04-muestreo-avanzado.ipynb">notebook</a> | <a href="{{ site.baseurl }}/episodes/scripts/05-muestreo-avanzado.py">script</a>)