Crosslinked Networks (examples/spes_network)

This tutorial demonstrates how to model a crosslinked Solid Polymer Electrolyte (SPE) network. Unlike previous tutorials where chains were pre-built, here we simulate the in-situ polymerization process.

The system consists of:

  • Precursors: Mono-functional acrylates (A-B) and Di-functional crosslinkers (A-C-A).

  • Electrolyte: Lithium salt (LiTFSI) and Succinonitrile (SN) solvent.

  • Reaction: Vinyl polymerization of the acrylate groups to form a 3D network.

Key features:

  1. Reactive CG Simulation: Using DoMD to define reaction templates that are executed dynamically during the CG simulation (via HOOMD-blue/PyGAMD).

  2. Regiospecific Control: Defining strict SMARTS patterns to control reaction products (e.g., Head-to-Tail connectivity).

  3. Hybrid Force Field Strategy: Applying different parameterization methods based on molecule size and type.

Note

This tutorial assumes you have configured the environment and database as described in the DoMD (ChemFAST) Documentation.

Step 1: Chemical Definitions & Reaction Rules

We define the components and the reaction rules.

  • A: Acrylic Acid group (C=C center).

  • B: Butanol tail.

  • C: PEG crosslinker chain.

  • L/T: Li+ and TFSI-.

  • S: Succinonitrile (Solvent).

Important

Regiospecificity: Bridging Isotropic Beads and Anisotropic Chemistry

In Coarse-Grained (CG) simulations, reactions (e.g., A + A -> A-A) are typically triggered by spatial proximity, treating beads as isotropic spheres without inherent directionality. However, the resulting all-atom chemistry is strictly regiospecific.

For vinyl polymerization, the connectivity direction determines the microstructure. While Head-to-Tail (H-T) is common, Head-to-Head (H-H) or Tail-to-Tail (T-T) linkages are chemically valid structures that may occur under specific synthetic conditions:

  • Head-to-Tail (H-T): -CH2-CH(R)-CH2-CH(R)-

  • Head-to-Head (H-H): -CH2-CH(R)-CH(R)-CH2-

  • Tail-to-Tail (T-T): -CH(R)-CH2-CH2-CH(R)-

The “What You See Is What You Get” Principle:

ChemFAST does not enforce “standard” chemistry; it strictly executes the topology defined by your inputs. Therefore:

  1. Define SMARTS Precisely: You must write the SMARTS pattern to explicitly define the atom mapping for the desired topology (H-T, H-H, or T-T).

  2. Match Reaction Recording: Ensure that the reactant indices recorded by the CG engine (HOOMD-blue/PyGAMD) match the order of reactants in your SMARTS.

  • Example (H-T): [C:1]=[C:2](R).[C:3]=[C:4](R) >> [C:1][C:2](R)[C:3][C:4](R)

  • This SMARTS specifically instructs the software to connect the “Tail” ([C:2]) of the first monomer to the “Head” ([C:3]) of the second.

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
smiA = 'C=CC(=O)O'   # Acrylate monomer (reactive)
smiB = 'OCCCC'       # Butyl tail
smiC = 'OCCOCCOCCOCCO' # PEG Crosslinker
smiL = '[Li+]'       # Lithium
smiT = 'C(F)(F)(F)S(=O)(=O)[N-]S(=O)(=O)C(F)(F)(F)' # TFSI
smiS = 'N#CCCC#N'    # Succinonitrile

# 2. Define Reaction Rules
reaction_template = {
    'A-A': {
        'cg_reactant_list': [('A', 'A')],
        'smarts': '[C:1]=[C:2]C(=O)O.[C:3]=CC(=O)O>>[C:1][C:2](C(=O)O)[C:3]=CC(=O)O',
        'prod_idx': [0]
    },
    'A-B': {
        'cg_reactant_list': [('A', 'B'),],
        'smarts': 'C(=O)[#8H1:1].[#6:2][#8H1:3]>>C(=O)[#8:1][#6:2].[#8:3]',
        'prod_idx': [0]
    },
    'A-C': {
        'cg_reactant_list': [('A', 'C'),],
        'smarts': 'C(=O)[#8H1:1].[#6:2][#8H1:3]>>C(=O)[#8:1][#6:2].[#8:3]',
        'prod_idx': [0]
    },
    'S': {
        'cg_reactant_list': [('L',), ('T',), ('S', )],
        'smarts': '*>>*',
        'prod_idx':[0],
    }

}


mols = {
    'A': {'smiles': smiA, 'file': None},
    'B': {'smiles': smiB, 'file': None},
    'C': {'smiles': smiC, 'file': None},
    'L': {'smiles': smiL, 'file': None},
    'T': {'smiles': smiT, 'file': None},
    'S': {'smiles': smiS, 'file': None},
}

Step 2: Generate Initial CG System

We initialize the system with precursors rather than a fully formed network. The network will be formed during the CG simulation step.

  • Meta 1: A-B (Butyl Acrylate monomer).

  • Meta 2: A-C-A (PEG Diacrylate crosslinker).

  • Meta 3/4/5: Solvents and Ions.

# System Composition
n_ba = 30000        # Number of Butyl Acrylate precursors
n_pegda = int(n_ba * 0.01) # 1% Crosslinker
n_sn = int(n_ba * 12/7)    # Solvent
n_L = int(2 * n_ba/7)      # LiTFSI

# Define Precursor Topologies
meta1 = { 'type': ['A','B'], 'bond': '0-1', 'N': n_ba, 'is_rigid': False, 'is_angle': True }
meta2 = { 'type': ['A','C','A'], 'bond': '0-1,1-2', 'N': n_pegda, 'is_rigid': False, 'is_angle': True }
meta3 = { 'type': ['S'], 'N': n_sn, 'is_rigid': False, 'is_angle': True }
meta4 = { 'type': ['L'], 'N': n_L , 'is_rigid': False, 'is_angle': True }
meta5 = { 'type': ['T'], 'N': n_L , 'is_rigid': False, 'is_angle': True } # Corrected type to 'T'

metas = [meta1, meta2, meta3, meta4, meta5]

# Generate Initial "Soup"
create_cg_system(mols, reaction_template, metas)

Outputs: * out_chemfast_cg.xml: Contains the unconnected precursors mixed with solvents/ions.

Step 3: Reactive CG Simulation (Polymerization)

This is the most critical step. You must run a simulation that supports topology changes (bond formation).

External Protocol (HOOMD-blue / PyGAMD):

  1. Load out_chemfast_cg.xml.

  2. Define Reaction Criteria: Set the reaction radius and probability for the A-A reaction defined in Step 1.

  3. Run Dynamics: As beads diffuse and collide, ‘A’ beads will bond to form the polyacrylate backbone, crosslinked by ‘C’ chains.

  4. Equilibrate: After the reaction reaches the desired conversion (e.g., 90%), continue relaxing the network.

  5. Save the final crosslinked network as CG_pre_eq/final.xml.

Polymerization at CG scale to generate network topology

Step 4: All-Atom Topology Building

DoMD reconstructs the complex, irregular network topology from the reacted CG snapshot.

# Path to the reacted, crosslinked network snapshot
xmlfile = 'CG_pre_eq/final.xml'

# Reconstruct AA Topology
# chunks_per_d=5 is highly recommended for infinite networks to save memory
confs, mol_graphs, box = build_aa_topology(
    mols,
    reaction_template,
    xmlfile,
    chunks_per_d=5
)

Step 5: Force Field Assignment

A crosslinked network is essentially one giant molecule (Infinite Cluster). Standard template matching can fail on such irregular structures. We use a Size-Dependent Hybrid Strategy:

  • Giant Molecules (>20 atoms): The Polymer Network. Use ML (GATs) for robust parameterization of the irregular backbone.

  • Single Ions (1 atom): Li+. Use TPL (Templates) for standard OPLS parameters.

  • Small Molecules: TFSI / Solvent. Use Hybrid (DB + ML fallback).

strategies = []
for mol in mol_graphs:
    if mol.number_of_nodes()>20:
        strategies.append('ml')
    elif mol.number_of_nodes()==1:
        strategies.append('tpl')
    elif mol.number_of_nodes()==10:
        strategies.append('ml')
    else:
        strategies.append('hybrid')

Final Outputs: The generated chemfast.gro and chemfast.top represent a chemically specific, crosslinked polymer electrolyte network swollen with solvent and ions, ready for conductivity simulations.