drug-protein-ligand-md

Goal

To run a complete protein-ligand molecular dynamics simulation using OpenMM, starting from a system bundle produced by drug-complex-system-builder. The workflow includes:

Energy minimization
NVT equilibration with positional restraints on heavy atoms
NPT equilibration with restraints gradually released
NPT production run

The output is a DCD trajectory + final state checkpoint suitable for drug-trajectory-analysis.

Instructions

1. Prepare inputs

Required from drug-complex-system-builder:

system.xml: serialized OpenMM System
complex_solvated.pdb: solvated complex PDB (used as topology reference)

2. Run the simulation

# Env: drugmd-agent
python .agents/skills/drug-protein-ligand-md/scripts/run_md.py \
  --system_xml md/system/system.xml \
  --input_pdb md/system/complex_solvated.pdb \
  --temperature 300 \
  --pressure 1.0 \
  --timestep 4.0 \
  --minimize_steps 5000 \
  --equil_nvt_steps 25000 \
  --equil_npt_steps 50000 \
  --production_steps 2500000 \
  --restraint_k 50.0 \
  --reporting_interval 5000 \
  --checkpoint_interval 25000 \
  --output_dir md/run/

Key parameters:

--temperature: simulation temperature in Kelvin (default: 300).
--pressure: target pressure in atm (default: 1.0).
--timestep: integration timestep in fs (default: 4.0). 4 fs is safe with hydrogen mass repartitioning (HMR) from the system builder; use 2 fs without HMR.
--minimize_steps: max minimization steps (default: 5000). Set to 0 to skip.
--equil_nvt_steps: NVT equilibration steps with restraints on protein/ligand heavy atoms (default: 25000 = 100 ps at 4 fs).
--equil_npt_steps: NPT equilibration steps with restraints released (default: 50000 = 200 ps).
--production_steps: production NPT steps (default: 2500000 = 10 ns at 4 fs).
--restraint_k: restraint force constant for equilibration in kJ/mol/nm^2 (default: 50.0).
--reporting_interval: write trajectory frame every N steps (default: 5000 = 20 ps).
--checkpoint_interval: write checkpoint every N steps (default: 25000).

3. Output files

The script produces:

md/run/minimized.pdb: structure after energy minimization
md/run/nvt_equilibration.log: energy/temperature log during NVT equilibration
md/run/npt_equilibration.log: energy/temperature/density log during NPT equilibration
md/run/production.dcd: production trajectory (DCD format)
md/run/production.log: production energy/temperature/density log
md/run/final_state.xml: serialized simulation state for restarts
md/run/md_provenance.json: all simulation parameters and timing

4. Running replicates

For statistical confidence, run multiple independent replicates with different random seeds:

# Env: drugmd-agent
for i in 1 2 3; do
  python .agents/skills/drug-protein-ligand-md/scripts/run_md.py \
    --system_xml md/system/system.xml \
    --input_pdb md/system/complex_solvated.pdb \
    --production_steps 2500000 \
    --seed $((42 + i)) \
    --output_dir md/rep_${i}/
done

5. Quick validation checks

After the run, verify:

Temperature fluctuates around the target (check production.log)
Density is stable around 1.0 g/cm^3 for aqueous systems
Total energy does not drift monotonically
The ligand remains in the binding pocket (use drug-trajectory-analysis)

Examples

Example: 10 ns production MD of TYK2 complex

# Env: drugmd-agent
python .agents/skills/drug-protein-ligand-md/scripts/run_md.py \
  --system_xml tyk2/md/system/system.xml \
  --input_pdb tyk2/md/system/complex_solvated.pdb \
  --temperature 300 \
  --timestep 4.0 \
  --production_steps 2500000 \
  --output_dir tyk2/md/run/

Example: short 1 ns refinement for pose assessment

# Env: drugmd-agent
python .agents/skills/drug-protein-ligand-md/scripts/run_md.py \
  --system_xml md/system/system.xml \
  --input_pdb md/system/complex_solvated.pdb \
  --production_steps 250000 \
  --reporting_interval 2500 \
  --output_dir md/short_refine/

Constraints

Environment: Requires drugmd-agent.
GPU acceleration: The script auto-detects CUDA GPUs. Without a GPU, simulations will run on CPU (significantly slower; consider reducing production_steps for testing).
Timestep: 4 fs requires hydrogen mass repartitioning (HMR) in the system. The drug-complex-system-builder applies HMR by default. If using a system without HMR, set --timestep 2.0.
Trajectory size: DCD files grow ~1 MB per 1000 frames for a typical 50k-atom system. A 10 ns run at 20 ps intervals produces ~500 frames (~500 MB).
Restarts: use --restart_from with a saved state XML to continue a simulation.

References

Eastman, P.; Galvelis, R.; Pelaez, R. P.; et al. OpenMM 8: Molecular Dynamics Simulation with Machine Learning Potentials. J. Phys. Chem. B 2024, 128(1), 109-116. https://doi.org/10.1021/acs.jpcb.3c06662
Hopkins, C. W.; Le Grand, S.; Walker, R. C.; Roitberg, A. E. Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning. J. Chem. Theory Comput. 2015, 11, 1864-1874. https://doi.org/10.1021/ct5010406

Author: Matthew Cox Contact: GitHub @mcox3406