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

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.

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.

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

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

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