Coarse-grained

Dissipative particle dynamics (mesoscale)

Coarse beads for the mesoscale of blends, block copolymers and solutions (Groot & Warren 1997). Beads within the cut-off r_c interact by a soft conservative force a_ij(1 − r)r̂, a dissipative force and a random force related by fluctuation–dissipation (a pair thermostat that conserves momentum, so hydrodynamics survive); bonded beads by springs. Chemistry enters through a_ij: like beads 75 kT/ρ (25 at ρ = 3, water's compressibility), unlike beads a_ii + χ/0.286 at ρ = 3. A species is a bead sequence (A5B5: a symmetric diblock) with a count. Reported: kT and the pressure, the segregation order parameter ψ over the run (0 mixed, 1 separated), the structure factor of the first bead type with its peak q* and the domain spacing 2π/q*, and the frames as a structure (beads as atoms, r_c in Å for the view).

Groot & Warren 1997; Hoogerbrugge & Koelman 1992

When to use it

Morphologies over tens of nanometres and microseconds: block copolymer microdomains (lamellae, cylinders, spheres), blend phase separation, micelles, the effect of χ and composition. Map chemistry to beads (a few repeat units each) and χ from Blend phase diagram or Solvent screen.

Source

core/src/dpd.cpp (run_dpd)

Tested by

Dpd.EquationOfStateAndTemperature (ideal gas p = ρkT; water model p = 23.65 at kT = 1; (p − ρkT)/(aρ²) → 0.101 at high density); Dpd.DiblockMicrophaseSeparation (χN = 43: ψ 0.18 → 0.67, an S(q) peak)

Departure from the reference

reduced units; χ mapping for ρ = 3 only (give a_ij at other densities)

References

1. Groot, R. D., Warren, P. B., "Dissipative particle dynamics: bridging the gap between atomistic and mesoscopic simulation", J. Chem. Phys. 107, 4423–4435 (1997). doi:10.1063/1.474784
2. Hoogerbrugge, P. J., Koelman, J. M. V. A., "Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics", Europhys. Lett. 19, 155–160 (1992). doi:10.1209/0295-5075/19/3/001

Kremer–Grest bead-spring melts

Each chain is a string of beads joined by FENE springs; all beads repel through the WCA potential (Lennard-Jones cut at its minimum and shifted). The model has no chemistry, only connectivity and excluded volume, so it shows the universal physics of polymer melts — entanglements, reptation — at a fraction of the cost. CAPS builds the melt as random walks and writes a LAMMPS deck that pushes the overlaps apart before the real potential starts.

Kremer & Grest 1990

When to use it

Entanglement and dynamics studies of generic melts and rubbers, and as a fast starting point. Runs in LAMMPS (reduced units); CAPS builds, draws and writes the deck. Model resolution backmaps the beads onto an all-atom structure.

Source

core/src/kremer_grest.cpp

Tested by

KremerGrest.ChainsBoxAndDeck · the deck run in LAMMPS (P ≈ 5.0 at T = 1, ρσ³ = 0.85)

Departure from the reference

the push-off uses LAMMPS pair soft ramped to 60 ε with displacements limited, rather than the force-capped potentials of the paper

References

1. Kremer, K., Grest, G. S., "Dynamics of entangled linear polymer melts: a molecular-dynamics simulation", J. Chem. Phys. 92, 5057–5086 (1990). doi:10.1063/1.458541
2. Auhl, R., Everaers, R., Grest, G. S., Kremer, K., Plimpton, S. J., "Equilibration of long chain polymer melts in computer simulations", J. Chem. Phys. 119, 12718–12728 (2003). doi:10.1063/1.1628670

DPD chain stiffness and domains

Semiflexible DPD chains through a bending energy at every middle bead, and the morphology at the end of a run counted: the domains of the first two bead types (connected cells where one type outnumbers the other), their mean size and the largest one's share — a single domain holding nearly all of its type's cells is a continuous phase.

forces F_i = −k_θ ∂cos θ/∂r_i on the end beads, the middle bead the opposite of their sum (LAMMPS angle_style cosine); domains on cells about r_c wide joined through faces, periodic

When to use it

Rod–coil and liquid-crystalline blocks, stiff fibre-forming polymers in a soft matrix, and telling a co-continuous blend from dispersed droplets.

Source

core/src/dpd.cpp

Tested by

tests/test_dpd.cpp Dpd.StiffnessAndDomains (k_θ = 5 kT moves ⟨cos θ⟩ from −0.06 to −0.79; a strongly segregated A10/B10 blend ends as one A and one B domain)

Departure from the reference

the domain count depends on the cell size used (about r_c): very small domains merge with their neighbours' cells

DPD branched molecules and templated starts

Bead sequences with branches — a side chain in parentheses hangs off the bead before it (grafts, stars, combs) — and starts from a mesostructure of the first bead type (lamellae normal to x, cylinders along z, spheres on a cubic lattice) instead of a random mix, to test whether a morphology is stable or to reach it sooner.

P the template's period (the box over the periods), φ_A the first type's bead share; each bead placed near the bead it is bonded to, in its type's region

When to use it

Grafted rubbers (silane- or polymer-grafted), star and comb architectures, and checking a block copolymer's morphology against the one it is started in (a stable structure stays; an unstable one rearranges).

Source

core/src/dpd.cpp

Tested by

tests/test_dpd.cpp Dpd.BranchesAndTemplates (a graft's branch point has three neighbours, a star's core three; a lamellar start of a blend begins at ψ ≈ 1, a random one near 0)

Departure from the reference

the template sets the start only: the dynamics decides what survives

Iterative Boltzmann inversion

The coarse-grained bead–bead potential as a table refined until the CG model reproduces the mapped all-atom structure: from −k_BT ln g_target, each iteration runs the CG model, measures its g(r) over the same pairs as the target, and corrects the potential by k_BT ln(g/g_target). The table goes to LAMMPS as pair_style table.

a straight wall below the first resolved point, smoothed over three points, zero at the cut-off (the model's own cut-off, the table's end); the optional pressure correction ΔV = A(1 − r/r_c) with A = −0.1 k_BT sign(ΔP) min(1, 0.0003 |ΔP|/bar); pairs on other chains or more than three bonds apart

When to use it

A rubber or polymer CG model whose melt structure (and with the pressure correction, density) matches the all-atom one, for times and sizes the atoms cannot reach.

Source

core/src/ibi.cpp

Tested by

tests/test_dynamics.cpp Dynamics.TabulatedPairMatchesLennardJones (a tabulated LJ gives the analytic energy and forces) and Dynamics.IbiRecoversTheLjLiquidStructure (the residual falls from 0.0095 to 0.0025); a polystyrene CG model's energy in CAPS and in LAMMPS agree to 0.1 %

Departure from the reference

one table for every bead pair (the pooled g(r)); structure, not dynamics, is matched (CG time runs faster than real time); the model belongs to the state point it was refined at

References

1. Reith, D., Pütz, M., Müller-Plathe, F., "Deriving effective mesoscale potentials from atomistic simulations", J. Comput. Chem. 24, 1624–1636 (2003). doi:10.1002/jcc.10307

Chemistry-aware coarse-grained mapping

Beads follow the chemistry, not a count of atoms. Bonds matched by SMARTS are cut (for polyesters every ester C(=O)–O, the bond between the pattern's first and last atoms) and the connected fragments left are the beads; or fragments are given by SMARTS (an exact cover of each molecule), or listed atom by atom. Each fragment is classed by its formula, aromaticity and a Weisfeiler–Lehman key of its heavy-atom graph, and named by the first naming SMARTS that embeds in it; a class no rule names stops the mapping (or is named by its class), never merged into a neighbour. End groups stay in their bead, so mass is conserved exactly. Several systems share one type list (bead, bond, angle and dihedral types numbered alike), so PBS, PBSA and PBAT share B and its terms.

the PBS, PBSA (80/20) and PBAT (56/44) cells of DP 25: 2 beads per repeat unit, mass conserved to 10⁻⁸ g/mol

When to use it

caps cgmap (Coarse-grain page, Mapping). Trajectories are mapped frame by frame (LAMMPS dumps: xu yu zu, x y z with or without image flags, scaled); wrapped beads are made whole along each chain's bonds.

Source

core/src/cg_rules.cpp, data/cg/mapping_rules.json

Tested by

CgRules.*, bench T13

Departure from the reference

Bonded potentials by Boltzmann inversion

Bond, angle and dihedral distributions of the mapped beads, gathered over many frames of several systems (each system's histograms normalised, then weighted), are inverted per type. Only where the smoothed distribution exceeds a share of its maximum (5 %) is it inverted; outside, the potential continues with its edge slope and a stiff quadratic wall, so sparse tails — distorted contacts of an unrelaxed start — cannot leave soft spots that let beads collapse. Dihedrals use the IUPAC sign; unsampled ranges get a smooth barrier. Each table carries a split-half check (odd and even molecules inverted separately): halves more than 1 kT apart are reported as not converged. A bonded IBI step corrects the shift the non-bonded 1–3 and 1–4 terms cause.

ideal chain drawn from known potentials: recovered within 0.06 kT (bench T14: ≤ 0.1 kT)

When to use it

caps cgfit bonded / refine. Writes LAMMPS bond and angle tables (angles in degrees, the force per degree as angle_style table reads it), dihedral_style table/cut (its aat switch turns a dihedral off as an angle straightens) and GROMACS tables (which recent GROMACS versions do not run). The suggested time step is a thirtieth of the stiffest bond's period with its wall (10.5 fs for the polyester tables).

Source

core/src/cg_bonded.cpp

Tested by

CgBonded.*, bench T14

Departure from the reference

The walls and the gap barriers are CAPS's choices for unsampled ranges, not inverted data.

References

1. Tschöp, W., Kremer, K., Batoulis, J., Bürger, T., Hahn, O., "Simulation of polymer melts. I. Coarse-graining procedure for polycarbonates", Acta Polym. 49, 61–74 (1998)
2. Reith, D., Pütz, M., Müller-Plathe, F., "Deriving effective mesoscale potentials from atomistic simulations", J. Comput. Chem. 24, 1624–1636 (2003). doi:10.1002/jcc.10307

Joint per-pair iterative Boltzmann inversion

One table per bead-type pair, started from −kT ln g of the targets and refined from CG runs. Several systems are fitted at once: a pair present in several (B–B in PBS, PBSA and PBAT) gets the average of their corrections, weighted by each system's statistics at that distance and its weight, so one table set serves all of them. The pressure correction adds a linear ramp that vanishes at the cut-off; its amplitude follows the weighted mean pressure error. Convergence per pair: the relative squared deviation of g and its largest point deviation. g(r) alone hardly fixes the pressure: the ramp and the thermomechanical calibration exist for that.

LAMMPS runs the tables (pair_style table) with the bonded tables; run_ibi.sh loops LAMMPS and caps cgfit ibi-step

When to use it

caps cgfit targets, ibi-start, ibi-step (each system's dump and log), ibi-run (CAPS's engine, small systems; bonded terms as springs at the tables' wells, dihedrals left out). Pairs no system has get a zero table so LAMMPS has every coefficient (said in pair.in). Below the first resolved g the wall is linear plus quadratic (10 kcal/mol/Ų).

Source

core/src/cg_nonbonded.cpp

Tested by

CgNonbonded.*, bench T15

Departure from the reference

References

1. Reith, D., Pütz, M., Müller-Plathe, F., "Deriving effective mesoscale potentials from atomistic simulations", J. Comput. Chem. 24, 1624–1636 (2003). doi:10.1002/jcc.10307

Analytic pairs and thermomechanical calibration

IBI tables belong to one state point. For mechanics they are fitted by an analytic form (LJ 12-6, LJ 9-6 as lj/class2, Morse, or Mie n–6), over the well and the repulsion up to a few kT, by Levenberg–Marquardt with the shift at the cut-off kept. The fitted set is then scaled: every σ by one factor to reach the target density, every ε by one factor to reach the target T_g, the ratios between pairs kept. Each step is a secant in log–log from the runs so far (the first from ρ ∝ s_σ⁻³ and T_g ∝ s_ε). T_g comes from a stepwise cooling scan: a bilinear fit of the specific volume.

secant steps are limited to ×0.8–1.25 (σ) and ×0.7–1.4 (ε) per run

When to use it

caps cgfit fit, tg, calibrate (in.cg_run with ENS npt for the density, in.cg_tg for the cooling scan). Energy renormalisation (ε as a function of T) is not done.

Source

core/src/cg_nonbonded.cpp

Tested by

CgNonbonded.FitsRecoverKnownParameters, CgNonbonded.CalibrationSteps

Departure from the reference

References

1. Hsu, D. D., Xia, W., Arturo, S. G., Keten, S., "Systematic method for thermomechanically consistent coarse-graining: A universal model for methacrylate-based polymers", J. Chem. Theory Comput. 10, 2514–2527 (2014). doi:10.1021/ct500080h

Coarse-grained melts of real sequences

Chains of repeat units (each a bead sequence: BS, BA, BT) drawn with Bernoulli, first-order Markov, block, gradient, alternating or patterned sequences, monodisperse or with drawn lengths (Schulz–Zimm …), as random walks whose bonds, angles and dihedrals are drawn from the inverted distributions (their Jacobians put back), in a cubic box at the density. The walks start with the all-atom model's local structure and overlap; a soft push-off, the model's pairs with limited steps, a hot anneal and NPT follow (in.cg_equil). The report compares ⟨R²(n)⟩/n with the mapped all-atom chains and warns when the box is smaller than the chains.

the PBSA walks reproduce the mapped cell's ⟨R²(n)⟩/n within 2 % for n = 1–20

When to use it

caps cgbuild (Coarse-grain page, Build CG melt). L ≥ R_ee at least, 1.5 R_ee better: a smaller box lets chains meet their own images.

Source

core/src/cg_build.cpp

Tested by

CgBuild.*

Departure from the reference

Neighbouring angles and dihedrals are drawn independently (their correlations come back during equilibration).

References

1. Auhl, R., Everaers, R., Grest, G. S., Kremer, K., Plimpton, S. J., "Equilibration of long chain polymer melts in computer simulations", J. Chem. Phys. 119, 12718–12728 (2003). doi:10.1063/1.1628670

Entanglement length estimators

Primitive paths from CAPS's PPA (ends held, intra-chain pairs off, FENE bonds with zero rest length, K = 30 ε/σ²), a LAMMPS PPA deck with the same potentials (the two agree within 1 % on a test melt), or Z1+'s shortest paths (Z = interior kinks). N_e follows from the classical and modified single-chain-length estimators, and from M-kink and M-coil over several chain lengths, which converge fastest. M_e = N_e × the mean bead mass.

Hoy, Foteinopoulou & Kröger 2009, Eqs. (4)–(7), (13), (15); N beads per chain

When to use it

caps ppa (Coarse-grain page, Analyse). Z1+ is run when it is installed ($CAPS_Z1); otherwise config.Z1 is written for it. Kink estimators need Z1+.

Source

core/src/cg_analysis.cpp

Tested by

CgAnalysis.EstimatorFormulas, CgAnalysis.Z1ConfigAndPaths

Departure from the reference

References

1. Everaers, R., Sukumaran, S. K., Grest, G. S., Svaneborg, C., Sivasubramanian, A., Kremer, K., "Rheology and microscopic topology of entangled polymeric liquids", Science 303, 823–826 (2004). doi:10.1126/science.1091215
2. Hoy, R. S., Foteinopoulou, K., Kröger, M., "Topological analysis of polymeric melts: Chain-length effects and fast-converging estimators for entanglement length", Phys. Rev. E 80, 031803 (2009). doi:10.1103/PhysRevE.80.031803
3. Kröger, M., "Shortest multiple disconnected path for the analysis of entanglements in two- and three-dimensional polymeric systems", Comput. Phys. Commun. 168, 209–232 (2005). doi:10.1016/j.cpc.2005.01.020
4. Kröger, M., Dietz, J. D., Hoy, R. S., Luap, C., "The Z1+ package: Shortest multiple disconnected path for the analysis of entanglements in macromolecular systems", Mendeley Data (2022). doi:10.17632/m425t6xtwr.1

Uniaxial tension of coarse-grained melts

LAMMPS decks for each mode and rate: stress mode (the lateral axes at P with an uncoupled NPT, so voids and necking can form) or constant volume, after an NPT relaxation. The axial stress and its bond, angle, dihedral, pair and kinetic parts (stress/atom per term, summed; they add up to the total) are printed at even strain steps, with frames for orientation and entanglement analyses. The analysis gives the modulus, the yield (where the smoothed stress first falls 2 % below its running maximum), the softening and the strain-hardening modulus.

⟨P₂⟩ = ⟨(3 cos²β − 1)/2⟩ of the bonds against the pull; R_ee anisotropy ⟨R_z²⟩/(½(⟨R_x²⟩ + ⟨R_y²⟩))

When to use it

caps mech decks / analyze. Rates are CG time: map them to all-atom time with the factor of caps cgdyn timemap before comparing with experiment.

Source

core/src/cg_analysis.cpp

Tested by

CgAnalysis.TensionCurve

Departure from the reference

References

1. Hoy, R. S., Robbins, M. O., "Strain hardening of polymer glasses: Effect of entanglement density, temperature, and rate", J. Polym. Sci. B: Polym. Phys. 44, 3487–3500 (2006). doi:10.1002/polb.21012

Chain dynamics and the coarse-grained time scale

Mean-square displacements over all time origins at logarithmic lags: g₁ of the inner half of each chain, g₂ about the chain's centre of mass, g₃ of the centres of mass; the bond and end-to-end autocorrelations. D from g₃'s last decade, τ_R where the end-to-end correlation falls to 1/e, τ_e where g₁'s local exponent first falls below 3/8 after the Rouse ½. CG dynamics run faster than all-atom (less friction): the time factor comes from the overlap of the two g₁ curves of the same system.

s is the geometric mean of the time ratios at equal g₁; its spread says whether one factor fits all time scales

When to use it

caps cgdyn (a CG run, or a mapped all-atom trajectory), caps cgdyn timemap AA.csv CG.csv.

Source

core/src/cg_analysis.cpp

Tested by

CgAnalysis.OrientationAndDynamics

Departure from the reference

References

1. Salerno, K. M., Agrawal, A., Perahia, D., Grest, G. S., "Resolving dynamic properties of polymers through coarse-grained computational studies", Phys. Rev. Lett. 116, 058302 (2016). doi:10.1103/PhysRevLett.116.058302

Backmapping chemistry-specific beads

A fragment library from a reference all-atom cell and its mapping: each bead class (kind and place in the chain) with conformers, and every bond, angle, dihedral and improper as a template over the classes of the consecutive beads it spans. Each target bead gets a fragment of its class, its centre of mass on the bead, turned by Horn's fit so that the directions to its two chain neighbours match the reference's; the cut bonds and the terms across them come back with their types and charges. A staged relaxation follows: bonded terms only, a soft-core push-off, the full force field (minimised with Ewald when the reference uses PPPM), a short NPT.

PBSA DP-50 melt: 4000 beads → 50 580 atoms, neutral; restored ester bonds 0.71 Å off r₀ → 0.017 Å after stages 1–3

When to use it

caps backmap REF_MAP REF_DATA --cg MAP --frame FRAME --input REF.in; caps backmap check relaxed.data (bonds and angles against harmonic r₀, θ₀).

Source

core/src/cg_backmap.cpp

Tested by

python: backmap gives back every term of the cell it was mapped from

Departure from the reference

Fragments are placed rigidly: the cut bonds start stretched or compressed until the relaxation.

References

1. Horn, B. K. P., "Closed-form solution of absolute orientation using unit quaternions", Journal of the Optical Society of America A 4, 629–642 (1987). doi:10.1364/JOSAA.4.000629