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 DoMD (ChemFAST) Documentation.

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.

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 = {
    'P-F': {
        'cg_reactant_list': [('P', 'F')],
        'smarts': '[CH3:1].[CH1:2]>>[C:1][C:2]',
        'prod_idx': [0]
    },
    'F-P': {
        'cg_reactant_list': [('F', 'P')],
        'smarts': '[CH3:1].[CH2:2]>>[C:1][C:2]',
        'prod_idx': [0]
    },
    'P-O': {
        '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': '[O:1].[CH3:2]>>[C:2][O:1]',
        '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…).

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

# 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.

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).

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.