Brush Polymer (examples/spes_brush)
==========================================================

This tutorial demonstrates how to model complex **Solid Polymer Electrolytes (SPEs)** using DoMD. We will construct a polymer brush system consisting of a fluorinated backbone with PEO side chains and Lithium ions.

This case study highlights advanced features of DoMD:

1.  **Programmatic Topology**: Using Python loops to define complex brush/comb architectures.
2.  **Specific FF Strategies**: Applying different parameterization strategies (ML vs. Template) for polymers and ions.
3.  **Large Molecule Backmapping**: Handling spatial decomposition for efficient embedding.

.. note::
   This tutorial assumes you have configured the environment and database as described in the :doc:`../index`.

Step 1: Chemical Definitions
----------------------------

We define a system with four components:

* **F**: A fluorinated monomer (acting as part of the backbone/linker).
* **P**: A connector unit (backbone).
* **O**: PEO repeating units (side chains).
* **L**: Lithium ions.

The reaction template defines the connectivity rules for these diverse components.

.. code-block:: python

    import sys
    from domd_tools import create_cg_system, build_aa_topology, assign_ff_parameters, get_gmx
    from misc.logger import logger
    logger.setLevel('INFO')

    # 1. Define Components
    smiF = 'CC(C(F)(F)(F))C(=O)[N-]S(=O)(=O)C(F)(F)(F)' # TFSI-like / Fluorinated unit
    smiP = 'CCO' # Backbone connector
    smiO = 'CCO' # PEO unit
    smiL = '[Li+]' # Lithium Ion

    # 2. Define Reaction Rules
    reaction_template = {
        'F-P': {
            'cg_reactant_list': [('P', 'F')],
            'smarts': '[CH3,CH2:1].[CH1,CH3:2]>>[C:1][C:2]',
            'prod_idx': [0]
        },
        'O-P': {
            'cg_reactant_list': [('P', 'O')],
            'smarts': '[OH1:1].[CH3:2][CH2:3]>>[O:1][C:2][C:3]',
            'prod_idx': [0]
        },
        'O-O': {
            'cg_reactant_list': [('O', 'O')],
            'smarts': '[CH3:1].[O:2]>>[C:1][O:2]',
            'prod_idx': [0]
        },
        'S': {
            'cg_reactant_list': [('L',)],
            'smarts': '*>>*',
            'prod_idx': [0],
        }
    }

    mols = {
        'F': {'smiles': smiF, 'file': None},
        'P': {'smiles': smiP, 'file': None},
        'O': {'smiles': smiO, 'file': None},
        'L': {'smiles': smiL, 'file': None},
    }

Step 2: Generate Initial CG System (Brush)
---------------------------------------------

Unlike simple linear chains, a polymer brush has a comb-like structure. DoMD allows you to define arbitrary connectivity using Python logic.

**Target Structure:**
The code below constructs a backbone of alternating P and F units, where each P unit holds a PEO side chain (O-O-O...).

.. code-block:: text

       |
       P-O-O-O-O-O-O-O... (Side Chain)
       |
       F
       |
       P-O-O-O-O-O-O-O...
       |

.. code-block:: python

    # System Parameters
    n_chain = 500   # Number of polymer chains
    cl = 20         # Chain length (number of repeating blocks)
    n_L = 10000     # Number of Lithium ions

    # Programmatic Topology Generation
    bonds = []
    # Loop over chain length to define connectivity
    for ic in range(cl):
        # 1. Connect Backbone P (0) to first Side Chain O (1)
        bonds.append(f'{ic*15+0}-{ic*15+1}')

        # 2. Build Side Chain (O-O bonds)
        for i in range(1, 13):
             bonds.append(f'{ic*15+i}-{ic*15+i+1}')

        # 3. Connect Backbone P (0) to Backbone F (14)
        bonds.append(f'{ic*15+0}-{ic*15+14}')

    # 4. Connect Backbone F (14) of current block to P (0) of next block
    for ic in range(cl-1):
        bonds.append(f'{ic*15+14}-{ic*15+15}')

    # Create the bond string
    top_str = ','.join(bonds)

    # Define Metadata
    # Sequence: P, 13*O, F (Total 15 beads per block)
    meta1 = {
        'type': ['P','O','O','O','O','O','O','O','O','O','O','O','O','O','F'] * cl,
        'bond': top_str,
        'N': n_chain,
        'is_rigid': False,
        'is_angle': True
    }

    # Define Ions (Free particles)
    meta2 = {
        'type': ['L'],
        'N': n_L,
        'is_rigid': False,
        'is_angle': True
    }

    # Generate CG System
    create_cg_system(mols, reaction_template, [meta1, meta2], rho=0.3)

**Outputs from Step 2:**

* ``out_chemfast_cg.xml``: Initial CG configuration including the complex brush topology.
* ``cg_params.txt``: Derived force field parameters.

Step 3: CG Pre-equilibration (Simulation)
----------------------------

Run the coarse-grained simulation to relax this dense brush system.

.. important::
   **External Simulation Step:**

   1.  Take the ``out_chemfast_cg.xml`` generated in Step 2.
   2.  Run a molecular dynamics simulation (HOOMD-blue/GALAMOST).
   3.  Ensure the ions diffuse and the brushes relax from their initial generation state.
   4.  Save the final frame as ``CG_pre_eq/final.xml``.

Step 4: All-Atom Topology Building
-------------------------------

We reconstruct the all-atom coordinates. Note the use of ``chunks_per_d=5``.

.. tip::
   **Handling Large Molecules:**
   For high molecular weight polymers (like these brushes), standard embedding can be slow. Setting ``chunks_per_d=5`` splits the molecule into spatial chunks for efficient parallel reconstruction.

.. code-block:: python

    xmlfile = 'CG_pre_eq/final.xml'

    # Backmapping with spatial decomposition
    confs, mol_graphs, box = build_aa_topology(
        mols,
        reaction_template,
        xmlfile,
        chunks_per_d=5
    )

Step 5: Force Field Assignment
---------------------------------------------

For this system, we use a mixed parameterization strategy to optimize accuracy and speed:

* **Polymer Brushes (>20 atoms)**: Use **ML (Machine Learning)** for accurate charge and parameter prediction of the complex polymeric structure.
* **Lithium Ions (1 atom)**: Use **TPL (Template)** rules to strictly follow standard library parameters (e.g., OPLS-AA ion parameters).

.. code-block:: python

    strategies = []

    # Iterate through all molecules to assign strategies
    for mol in mol_graphs:
        if mol.number_of_nodes() > 20:
            strategies.append('ml')  # Use Graph Attention Network for brushes
        elif mol.number_of_nodes() == 1:
            strategies.append('tpl') # Use standard library rules for ions

    # Assign parameters
    ffs = assign_ff_parameters(mol_graphs, strategies=strategies)

    # Export to GROMACS
    get_gmx(mol_graphs, confs, ffs, box)

**Final Outputs:**
* ``spes_brush_final.gro``
* ``spes_brush_final.top``

The system is now ready for ion transport simulations in GROMACS.