Density profile
===============

The ``density`` tool computes **mass density** or **charge density** profiles of
a molecular system along any of the three axes (x, y, z).  
It supports:

- slab-density calculations  
- dense-phase recentering  
- multi-group density reporting  
- customizable time windows  
- mass or charge density modes  
- interactive index-group selection  

This tool is implemented in ``density.py``.


Overview
--------

The program performs the following tasks:

1. Load a molecular trajectory (TPR + XTC + optional NDX).  
2. Determine time windows for analysis.  
3. Select a group of atoms that defines the **dense phase**.  
4. Select groups of atoms for which density will be computed.  
5. For each time window:  
   - compute linear-density profiles  
   - optionally shift density profiles to recenter dense phase  
6. Accumulate and average density profiles across windows.  
7. Write one or more averaged density profiles into an XVG file.

A slab LLPS simulation typically shows:

- thick dense-phase region  
- low-density dilute region  
- need for recentering due to diffusion / drift

This program automates the centering process frame-by-frame or interval-by-interval.


Usage
-----

.. code-block:: bash

   dps density -s run.tpr -f run.xtc -o density.xvg -axis z -type mass


Arguments
---------

Required
--------

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

   TPR file containing topology, box size, chain info, atom masses/charges.

.. option:: -f, --input XTC

   Trajectory file to compute density from.

.. option:: -o, --output FILE

   Name of the XVG output file.  
   If not ending with ``.xvg``, the extension will be added automatically.


Optional
--------

Axis selection
~~~~~~~~~~~~~~

.. option:: -x, --axis {x,y,z}

   Axis along which density is calculated. Default: ``z``.

Density type
~~~~~~~~~~~~

.. option:: -tp, --type {mass, charge}

   - ``mass``: mass density (mg/mL)  
   - ``charge``: charge density (e/mol/mL)

Time window
~~~~~~~~~~~

.. option:: -b, --start-time INT  
.. option:: -e, --end-time INT  

   Start and end time (ns). Default: first → last frame.

.. option:: -dt, --delta-time INT

   Interval (ns) between *recenter attempts*.  
   Default: ``1`` ns.

Selection
~~~~~~~~~

.. option:: -selfit, --selection-fit INT

   Index group ID to use for dense-phase detection.

.. option:: -sel, --selection-calculate INT INT ...

   List of index groups for which density profiles are computed.

If not provided, the program enters an **interactive group selection** mode.

Centering
~~~~~~~~~

.. option:: -nc, --no-center

   Disable centering.  
   When disabled, raw density profiles are used.

.. option:: -t, --dense-phase-threshold FLOAT

   Threshold for defining the dense phase.  
   Dense bins are those with density > ``threshold × max_density``.  
   Default: ``0.5``.


Program Workflow
----------------

(1) **Load trajectory**

.. code-block:: python

   trajectory = trajectory_class(tpr, ndx, xtc)

If loading fails, the program prints an error and exits.

(2) **Time and frame range**

The program determines the frame indices:

::

   start_frame = (start_time - time_init) / dt
   end_frame   = (end_time   - time_init) / dt
   center_interval_frame = (delta_time / dt)

(3) **Group selection**

- Fit group: used to detect dense-phase center  
- Calculation groups: used to compute density curves

Interactive examples:

- ``group 0``  
- ``group 18 19 20``  
- ``abbr ALA GLY PRO``  
- ``splitch``

(4) **Density calculation (MDAnalysis LinearDensity)**

For each time window:

.. code-block:: python

   lin.LinearDensity(selection, binsize=0.5).run(start=F1, stop=F2)

Density is taken from:

- ``mass_density`` (×1000 → mg/mL)  
- ``charge_density`` (e / mol / mL)

(5) **Recentering (default enabled)**

The dense-phase region is detected using:

- threshold × max-density  
- periodic-connected regions  
- largest continuous dense region  
- compute center-of-dense-phase  
- compute shift to move dense phase to box center  
- apply shift with ``np.roll``

Function:

.. code-block:: python

   shift_density_center(density_profiles, threshold, reference_profile)

(6) **Averaging**

All density curves are summed and divided by number of windows.

(7) **Writing output**

.. code-block:: python

   write_xvg(output_file, x_axis, density_profiles, legends=[group names])

Output labels:

- X-axis: position along selected axis (nm)  
- Y-axis: mass or charge density  


Output Format
-------------

The output XVG contains:

- one density profile per selected calculation group  
- evenly spaced bin centers  
- averaged (and recentered) values across all windows  

Y-axis units:

- mass density: ``mg/mL``  
- charge density: ``e/mol/mL``


Example
-------

Compute mass density of two groups and recenter automatically:

.. code-block:: bash

   dps density \
       -s run.tpr \
       -f run.xtc \
       -n run.ndx \
       -o density_mass.xvg \
       -axis z \
       -type mass \
       -selfit 0 \
       -sel 1 2 \
       -b 10 -e 100 -dt 2

Compute charge density along x-axis, no centering:

.. code-block:: bash

   dps density \
       -s run.tpr \
       -f run.xtc \
       -o charge.xvg \
       -type charge \
       -axis x \
       -nc


Error Messages
--------------

**“YOU MUST BE KIDDING ME.”**  
Start/end times exceed trajectory range.

**“Exception occurred when trying to open trajectory file …”**  
The TPR or XTC file failed to open.

**“Exception occurred when trying to open text file …”**  
Output file could not be created (permissions, path, etc.).

**Interactive selection errors**  
Any malformed selection (invalid group number, empty group) triggers an error.


Summary
-------

``dps density`` computes high-resolution slab density profiles with:

- dense-phase-aware recentering  
- multi-group density reporting  
- mass and charge density modes  
- interactive or scripted workflows  

It is a key analysis tool for LLPS simulations, where shifting and averaging
density profiles is critical for phase-behavior characterization.