Linear Polymer (examples/pi)
=======================================================

This tutorial demonstrates the advanced **"Coarse-Graining $\rightarrow$ Pre-equilibration $\rightarrow$ Backmapping $\rightarrow$ Force Field Assignment $\rightarrow$ Write Output"** workflow provided by DoMD.

Generating equilibrated structures for long-chain polymers in atomistic detail is computationally expensive. DoMD solves this by:

1.  **CG Initialization**: Generating a specific Coarse-Grained (CG) system from SMILES.
2.  **CG Pre-equilibration**: Relaxing the system using efficient CG dynamics.
3.  **Backmapping**: Reconstructing the all-atom details from the relaxed CG configuration.
4.  **Force Field Assignment**: Assigning OPLS-AA parameters using DoMD's database and ML models.
5.  **Output Generation**: Exporting the final equilibrated structure and topology for GROMACS simulations.

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

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

First, we define the chemistry of the system. We use **ODPA** (4,4′-oxydiphthalic anhydride) and **ODA** (4,4'-Oxydianiline) as monomers.

The ``reaction_template`` is crucial. It defines how the monomers connect (using SMARTS) for both the CG topology construction and the AA reconstruction.

.. 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 Monomers
    smiA = 'Nc1ccc(Oc2ccc(N)cc2)cc1'  # ODA (Diamine)
    smiB = 'O1C(=O)c2cc(Oc3cc4C(=O)OC(=O)c4cc3)ccc2C1=O' # ODPA (Dianhydride)

    # 2. Define Reaction Rules (Imidization)
    reaction_template = {
        'B-A': {
            'cg_reactant_list': [('A', 'B')],
            # SMARTS for: Anhydride + Amine -> Imide Linkage
            'smarts': '[#7H2:1].[#6:3](=[#8:4])[#8:2][#6:5]=[#8:6]>>[#6:3](=[#8:4])[#7:1][#6:5]=[#8:6].[#8:2]',
            'prod_idx': [0]
        }
    }

    mols = {
        'A': {'smiles': smiA, 'file': None},
        'B': {'smiles': smiB, 'file': None},
    }

Step 2: Generate Initial CG System
----------------------------------

We describe the system composition using metadata. Here we define two types of polymer chains:
* **Chain Type 1**: 20 repeating units (A-B...), 20 chains.
* **Chain Type 2**: 10 repeating units, 20 chains.

Calling ``create_cg_system`` will calculate the force field parameters (using Hybrid ML/HSP) and generate the simulation files.

.. code-block:: python

    # Define chain composition
    meta1 = {
        'type': ['A','B'] * 10, # Sequence: A-B-A-B... (20 beads)
        'bond': '0-1,1-2,2-3,3-4,4-5,5-6,6-7,7-8,8-9,9-10,10-11,11-12,12-13,13-14,14-15,15-16,16-17,17-18,18-19',
        'N': 20, # Number of chains
        'is_rigid': False,
        'is_angle': True
    }

    meta2 = {
        'type': ['A','B'] * 5, # Sequence: A-B... (10 beads)
        'bond': '0-1,1-2,2-3,3-4,4-5,5-6,6-7,7-8,8-9',
        'N': 20, # Number of chains
        'is_rigid': False,
        'is_angle': True
    }

    # Generate CG input files
    create_cg_system(mols, reaction_template, [meta1, meta2])

**Outputs from Step 2:**

* ``out_chemfast_cg.xml``: Initial coordinates and topology for the CG simulation.
* ``cg_params.txt``: The derived force field parameters.

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

At this stage, you need to perform a coarse-grained simulation to relax the system.

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

   1.  Take the ``out_chemfast_cg.xml`` generated in Step 2.
   2.  Run a molecular dynamics simulation using **HOOMD-blue v1.3.3** or **GALAMOST**.
   3.  Perform energy minimization, followed by NVT/NPT equilibration to reach the desired density and chain conformation.
   4.  Save the final frame (snapshot) of the simulation as ``CG_pre_eq/final.xml``.

   This step creates a realistic, entangled polymer melt structure at a fraction of the cost of all-atom equilibration.

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

Once the CG simulation is complete and you have the ``final.xml``, DoMD can reconstruct the atomistic details. It maps the atoms back onto the bead positions and rebuilds the chemical topology based on the reaction templates.

.. code-block:: python

    # Path to the relaxed CG configuration
    xmlfile = 'CG_pre_eq/final.xml'

    # Reconstruct Topology and Coordinates
    confs, mol_graphs, box = build_aa_topology(
        mols,
        reaction_template,
        xmlfile,
        write_molecule=True
    )

Step 5: Force Field Assignment & Output
---------------------------------------

Finally, we assign OPLS-AA force field parameters to the reconstructed all-atom topology and export the files for GROMACS. DoMD handles the atom typing using its database and Graph Attention Networks (GATs).

.. code-block:: python

    # 1. Assign AA Force Field
    ffs = assign_ff_parameters(mol_graphs, strategies=['ml']*len(mol_graphs))

    # 2. Write GROMACS files
    get_gmx(
        mol_graphs,
        confs,
        ffs,
        box,
    )

**Final Outputs:**

* ``chemfast.gro``: Equilibrated all-atom coordinates.
* ``chemfast.top``: Full system topology with OPLS-AA parameters.

You can now use these files to start your simulation in GROMACS.