Core-Shell Nanoparticles (examples/core_shell_nps)
=================================================================

This tutorial demonstrates the modeling of a co-assembly system involving a Block Copolymer (**PEO-b-PAsp(DET)**) and a rigid biological macromolecule (**CRISPR-Cas9**, PDB: 4OO8).

Key features of this tutorial:

1.  **Rigid Body Integration**: Importing an external PDB file (``4OO8_1_hs.pdb``) to represent a complex biological macromolecule (RNP) as a rigid bead in the CG simulation.
2.  **Implicit Solvent Modeling**: Controlling the assembly behavior (collapse vs. swelling) by tuning the attraction strength scaling factor ($\alpha$).
3.  **Specific Backmapping**: Reconstructing all-atom coordinates while preserving the internal structure of the rigid protein.

.. note::
   **Force Field Note**: While we assign OPLS-AA parameters here, standard OPLS-AA may require specific tuning for complex proteins/nucleic acids (e.g. OPLS-AA/M). In this context, the force field assignment primarily serves to construct a valid topology for further relaxation or specific force field replacement (e.g., CHARMM, AMBER) later.

Step 1: Chemical Definitions & Rigid Input
------------------------------------------

We define the block copolymer components (PEO and Asp-derivative) and the rigid protein.

* **O**: PEO unit (Shell-forming block).
* **N**: Functionalized Aspartame unit (Core-forming block).
* **R**: CRISPR-Cas9 RNP, defined by an external PDB file.

.. code-block:: python

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

    # 1. Define SMILES
    smiO = 'OCC' # Polyethylene oxide
    smiN = 'OC(=O)C(CC(=O)NCCNCCN)N' # Functionalized Asp monomer (DET)
    smiR = 'CC' # Placeholder SMILES for the rigid center

    # 2. Define Reaction (Polymerization)
    reaction_template = {
        'O-O': {
            'cg_reactant_list': [('O', 'O')],
            'smarts': '[CH3:1].[OH1:2]>>[C:1][O:2]',
            'prod_idx': [0]
        },
        'N-N': {
            'cg_reactant_list': [('N', 'N')],
            'smarts': '[C:1](=[O:2])[OH1:3].[NH2:4][CH1:5]>>[C:1](=[O:2])[N:4][C:5].[O:3]',
            'prod_idx': [0]
        },
        'O-N': {
            'cg_reactant_list': [('O', 'N')],
            'smarts': '[CH3:1].[NH2:2][CH1:3]>>[C:1][N:2][C:3]',
            'prod_idx': [0]
        },
    }

    # 3. Load Molecules with Rigid Body PDB
    mols = {
        'O': {'smiles': smiO, 'file': None},
        'N': {'smiles': smiN, 'file': None},
        # Crucial: 'file' points to the PDB, 'is_rigid' marks it as a rigid body
        'R': {'smiles': smiR, 'file': '4OO8_1_hs.pdb', 'is_rigid': True},
    }

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

We construct a system containing 100 block copolymer chains and 5 rigid RNP complexes.

* **Block Copolymer**: A linear chain of 272 'O' units and 78 'N' units.
* **Rigid Protein**: 5 isolated rigid bodies.

.. code-block:: python

    # Define Block Copolymer Topology (O-O-...-O-N-N-...-N)
    n_polymer = 100
    bonds = [f'{i}-{i+1}' for i in range(272+78-1)]
    bondstring = ','.join(bonds)

    # Flexible Chain Metadata
    meta1 = {
        'type': ['O']*272 + ['N']*78,
        'bond': bondstring,
        'N': n_polymer,
        'is_rigid': False,
        'is_angle': True
    }

    # Rigid Body Metadata
    n_R = 5
    meta2 = {
        'type': ['R'],
        'N': n_R,
        'is_rigid': True, # Must match the molecule definition
        'is_angle': True
    }

    metas = [meta1, meta2]

    # Generate System
    # Note: iterations=3 helps converge the initial parameter estimation
    create_cg_system(mols, reaction_template, metas, boxl=150, iterations=3)

**Outputs:**

* ``out_chemfast_cg.xml``: System topology.
* ``rigid_meta.pkl``: **Crucial**. Contains the mapping between the CG rigid bead and the atoms in the PDB file.

Step 3: Implicit Solvent Assembly (CG Simulation)
----------------------------------------------

This step is performed in the MD engine (HOOMD/GALAMOST). Since we are using implicit solvent, the assembly behavior is controlled by scaling the Lennard-Jones attraction strength ($\alpha$).

.. tip::
   **$\alpha$ Tuning of Lennard-Jones Potential**
   The nonbond Lennard-Jones potential is below:

   .. math::

       U_{LJ}(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \alpha \left( \frac{\sigma}{r} \right)^6 \right]

   while $\epsilon$ and $\sigma$ are determined by the HSP, $\alpha$ is a user-defined scaling factor to tune the attraction strength.
   The scaling parameter $\alpha$ directly correlates with the solvent quality and polymer scaling laws. The empirical guidelines for tuning $\alpha$ are derived from the study *"Brownian dynamics simulation study on the self-assembly of incompatible star-like block copolymers in dilute solution"* (Li, B. et al., *Phys. Chem. Chem. Phys.*, **2012**, 14, 4964-4970).

   * **$\alpha < 0.4$ (Hydrophilic/Good Solvent)**: Weak attraction. The polymer chains swell ($R_g \sim N^\nu, \nu > 0.5$).
   * **$\alpha > 0.6$ (Hydrophobic/Poor Solvent)**: Strong attraction. The polymer chains collapse ($R_g \sim N^\nu, \nu < 0.5$), driving aggregation.

**Defining Interactions for Core-Shell Assembly:**

For a successful Core-Shell assembly loaded with a rigid cargo (RNP), users are expected to define the interaction potentials based on the specific physicochemical properties of the cargo (e.g., hydrophobicity, charge).

1.  **Core Formation**: The 'N' block (PAsp) should be hydrophobic ($\alpha > 0.6$) to form the core.
2.  **Shell Stabilization**: The 'O' block (PEO) should be hydrophilic ($\alpha < 0.4$) to stabilize the nanoparticle in solution.
3.  **Cargo Loading (Rigid Body 'R')**:
    * If 'R' is hydrophobic or electrostatically attracted to the core, set **high affinity** between R and N (e.g., $\alpha_{R-N} \approx 1.0$).
    * Ensure 'R' has **low affinity** with the shell 'O' (e.g., $\alpha_{R-O} \approx 0.3$) to avoid surface adsorption and promote encapsulation.

**Protocol:**
1.  Initialize with ``out_chemfast_cg.xml``.
2.  Apply the interaction matrix defined above.
3.  Run the simulation until core-shell nanoparticles form and equilibrate.
4.  Save the snapshot as ``CG_pre_eq/final.xml``.

Step 4: All-Atom Topology Building (with Rigid Meta)
-----------------------------------

When backmapping systems with rigid bodies, we must load the ``rigid_meta.pkl`` generated in Step 2. This ensures the protein atoms are placed correctly relative to the rigid bead's center of mass and orientation.

.. code-block:: python

    # Path to equilibrated CG snapshot
    xmlfile = 'CG_pre_eq/final.xml'

    # Load rigid body metadata
    rigid_meta = pickle.load(open('rigid_meta.pkl', 'rb'))

    # Backmap
    # We pass rigid_meta to the builder so it knows how to reconstruct the protein
    confs, mol_graphs, box = build_aa_topology(
        mols,
        reaction_template,
        xmlfile,
        chunks_per_d=5,
        rigid_meta=rigid_meta
    )

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

We use a pure ML strategy here for efficiency on the large polymer systems.

.. code-block:: python

    # Assign Parameters using ML strategy for all molecules
    strategies = ['ml'] * len(confs)

    ffs = assign_ff_parameters(mol_graphs, strategies=strategies)

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

**Final Outputs:**
The resulting ``chemfast.gro`` and ``chemfast.top`` will contain the assembled core-shell nanoparticle with the atomistic RNP structure embedded within the PAsp core.