Run simulation
==============

The ``mdrun`` program executes **molecular dynamics simulations** using the
DROPPS/CGPS OpenMM-based backend.  
It consumes a runtime file generated by ``dps grompp`` and produces:

- trajectory files (XTC)  
- energy & log files  
- checkpoint files  
- final wrapped configuration (PDB)  

This is analogous to ``gmx mdrun`` in GROMACS.

Implemented in ``mdrun.py``.


Overview
--------

``dps mdrun``:

1. Loads **all runtime data** (system, topology, integrator parameters, positions)
   from a `.tpr` pickle generated by ``dps grompp``.
2. Sets up the **OpenMM integrator** (Langevin).  
3. Automatically selects the **best available compute platform**:
   ``CUDA → OpenCL → CPU``.
4. Applies optional:
   - energy minimization  
   - velocity initialization  
   - center-of-mass motion removal  
   - pressure coupling (MC barostat)  
5. Configures **reporters** for:
   - screen log  
   - file log  
   - trajectory (XTC)  
   - checkpoint  
6. Runs warming phase (if enabled).  
7. Runs **production dynamics**.  
8. Writes the **final wrapped PDB file** to disk.


Usage
-----

Run a simulation:

.. code-block:: bash

   dps mdrun -s run.tpr -o output/run

Resume from a checkpoint:

.. code-block:: bash

   dps mdrun -s run.tpr -cpi output/run.chk -o output/run2


Arguments
---------

Required
--------

.. option:: -s, --run-input FILE

   Runtime file produced by ``dps grompp`` (Python pickle format).

.. option:: -o, --output-prefix PREFIX

   Prefix for all output files:

   - ``PREFIX.log``  
   - ``PREFIX.xtc``  
   - ``PREFIX.chk``  
   - ``PREFIX.pdb``

Optional
--------

.. option:: -cpi, --checkpoint FILE

   Checkpoint file from a previous run.  
   If provided:

   - **no** minimization is performed  
   - **no** warming is performed  
   - velocities & box come from checkpoint


Simulation Workflow
-------------------

### 1. Load runtime data

.. code-block:: python

   with open(args.run_input, "rb") as f:
       loaded = pickle.load(f)

Loaded fields include:

- ``parameters`` (MDP dictionary)  
- ``mdsystem`` (OpenMM System)  
- ``mdtopology``  
- ``positions``  
- ``pdb_raw`` (for final PDB writing)  

### 2. Set up Langevin integrator

.. code-block:: python

   temperature = parameters["production_temperature"] * kelvin
   friction    = parameters["friction"] / picosecond
   timestep    = parameters["dt"] * picosecond

### 3. Automatic platform selection

Preferred order:

.. code-block:: python

   ['CUDA', 'OpenCL', 'CPU']

Platform-specific features:

- CUDA: ``DeviceIndex=0``, ``Precision=mixed``  
- OpenCL: similar  
- CPU: no properties

Aborts if no platform is available.

### 4. Center-of-mass motion removal

If ``comm_mode = Linear``:

.. code-block:: python

   mdsystem.addForce(openmm.CMMotionRemover(...))

Otherwise a warning is printed.

### 5. Pressure coupling (optional)

If ``pcoulp = True``:

.. code-block:: python

   mdsystem.addForce(openmm.MonteCarloBarostat(pressure, temperature, tau_P))

### 6. Energy minimization

Performed only if:

- ``minimize = True`` in parameters  
- **no** checkpoint is used  

Done using:

.. code-block:: python

   simulation.minimizeEnergy(force_tol, max_steps)

Maximum force and RMS force are reported.

### 7. Reporters

Supported reporters:

- **screen log** (printed to STDOUT)  
- **file log** (written to ``PREFIX.log``)  
- **XTC trajectory** (``PREFIX.xtc``)  
- **checkpoint files** (``PREFIX.chk``)

Configuration is driven by parameters:

- ``nst_screenlog``  
- ``screenlog_grps``  
- ``nst_filelog``  
- ``filelog_grps``  
- ``nst_xout``  
- ``nst_cp``  

### 8. Velocity initialization

If ``gen_vel = True``:

.. code-block:: python

   simulation.context.setVelocitiesToTemperature(initial_T)

### 9. Warming schedule (optional)

If:

- initial temperature < production temperature  
- **no** checkpoint is used  

Then temperature is raised gradually:

.. code-block:: python

   for T in range(T_initial, T_final):
       simulation.integrator.setTemperature(T)
       simulation.step(step_per_kelvin)

Temperature update every 10 K printed.

### 10. Production dynamics

Simulation length:

.. code-block:: python

   parameters["nsteps"]

Temperature is set to the final production temperature.

If a checkpoint is supplied:

.. code-block:: python

   simulation.loadCheckpoint(args.checkpoint)

### 11. Final wrapping and PDB writing

Final coordinates:

.. code-block:: python

   positions = simulation.context.getState(getPositions=True)
                    .getPositions(asNumpy=True).value_in_unit(nanometer)

Box vectors:

.. code-block:: python

   box_vec = [Lx, Ly, Lz]

Coordinates wrapped:

.. code-block:: python

   wrapped = positions % box_vec

Final PDB written with preserved metadata from ``pdb_raw``:

.. code-block:: python

   write_pdb(PREFIX.pdb, box_vec, ...)