polymer-build

AutoPoly: Polymer Structure Generation for LAMMPS

AutoPoly generates atomistic and coarse-grained polymer models from SMILES notation, producing ready-to-run LAMMPS simulation files. It bridges polymer chemistry (SMILES) and MD simulation (LAMMPS data files).

Always use the omnischolar conda environment, which has AutoPoly pre-installed.

Core Workflow

Every AutoPoly job follows three steps:

from AutoPoly import System, Polymer, Polymerization

# 1. Define output location
system = System(out="my_simulation")

# 2. Define polymer chains (using complement SMILES)
polymer = Polymer(
    chain_num=10,
    sequence=["CC[*]"] + ["[*]CC[*]"] * 48 + ["[*]CC"],  # DOP=50
    topology="linear",
    tacticity="atactic"
)

# 3. REQUIRED workaround for AutoPoly v1.0.0 bug (sequenceSet attribute)
polymer.sequenceSet = polymer.sequence_set

# 4. Generate LAMMPS files
Polymerization(
    name="polyethylene",
    system=system,
    model=[polymer],
    force_field="oplsaa"
)
# Output: my_simulation/polyethylene/ with system.data, system.in, system.in.settings

AutoPoly v1.0.0 bug: Atomistic Polymerization crashes with AttributeError: 'Polymer' object has no attribute 'sequenceSet'. The workaround above (polymer.sequenceSet = polymer.sequence_set) is required for all atomistic Polymer + Polymerization workflows. BeadSpringPolymer is not affected.

Complement SMILES Format

This is AutoPoly's unique notation for specifying monomer position in a chain. Understanding it is essential — most errors stem from getting this wrong.

Core principle: number of [*] wildcards = number of connections

Position	Wildcards	Example (Ethylene)	When to use
First	1 (right)	CC[*]	Chain start (linear only)
Middle	2 (both)	[*]CC[*]	Internal units + all ring positions
Last	1 (left)	[*]CC	Chain end (linear only)

Building a sequence:

# Step 1: Define middle variant (most common)
middle = "[*]CC[*]"

# Step 2: Derive first/last by removing one wildcard
first = "CC[*]"      # Remove left [*]
last = "[*]CC"       # Remove right [*]

# Step 3: Assemble
dop = 100
sequence = [first] + [middle] * (dop - 2) + [last]

Ring polymers use only middle variants (every position connects to two neighbors):

ring_sequence = ["[*]CC[*]"] * 50
polymer = Polymer(chain_num=5, sequence=ring_sequence, topology="ring")

For the full complement SMILES reference with all common monomers, read references/complement_smiles.md.

Common Monomer Quick Reference

Polymer	Abbr	First	Middle	Last
Polyethylene	PE	CC[*]	[*]CC[*]	[*]CC
Polypropylene	PP	CC(C)[*]	[*]CC([*])C	[*]CC(C)
Polystyrene	PS	CC(c1ccccc1)[*]	[*]CC([*])c1ccccc1	[*]CC(c1ccccc1)
PVC	PVC	CC(Cl)[*]	[*]CC([*])Cl	[*]CC(Cl)
PMMA	PMMA	CC(C)(C(=O)OC)[*]	[*]CC([*])(C)C(=O)OC	[*]CC(C)(C(=O)OC)
PEO	PEO	COC[*]	[*]COC[*]	[*]COC
PAN	PAN	CC(C#N)[*]	[*]CC([*])C#N	[*]CC(C#N)

Force Field Selection

Force Field	String	Best For
OPLS-AA	"oplsaa"	General organic polymers, safe default
LOPLS	"lopls"	Polymer melts, accurate densities
GAFF	"gaff"	Small molecules, functionalized monomers, mixed systems
GAFF2	"gaff2"	Updated GAFF, new projects
DREIDING	"dreiding"	Exploratory work, unusual chemistry, metals
COMPASS	"compass"	Commercial polymers, accurate mechanical properties

Decision shortcut:

• Simple hydrocarbons (PE, PP, PS) → "oplsaa"
• Functionalized monomers (PMMA, PVAc) or polymer+solvent → "gaff"
• Accurate density needed → "lopls" or "compass"
• Novel/exotic chemistry → "dreiding" (then switch for production)

GAFF/GAFF2 require separate charge calculation (AM1-BCC via Antechamber). OPLS-AA includes charges automatically.

For detailed force field comparison, read references/force_fields.md.

Common Patterns

Homopolymer

dop = 100
sequence = ["CC[*]"] + ["[*]CC[*]"] * (dop - 2) + ["[*]CC"]
polymer = Polymer(chain_num=10, sequence=sequence)
polymer.sequenceSet = polymer.sequence_set  # v1.0.0 workaround

Block Copolymer (ABA Triblock)

sequence = (
    ["CC[*]"] + ["[*]CC[*]"] * 9 +               # PE block (10 units)
    ["[*]CC([*])c1ccccc1"] * 20 +                 # PS block (20 units)
    ["[*]CC[*]"] * 9 + ["[*]CC"]                  # PE block (10 units)
)
polymer = Polymer(chain_num=5, sequence=sequence)  # DOP=40
polymer.sequenceSet = polymer.sequence_set  # v1.0.0 workaround

Alternating Copolymer

pair = ["[*]CC[*]", "[*]CC([*])c1ccccc1"]  # PE-PS alternating
sequence = ["CC[*]"] + pair * 24 + ["[*]CC(c1ccccc1)"]
polymer = Polymer(chain_num=10, sequence=sequence)
polymer.sequenceSet = polymer.sequence_set  # v1.0.0 workaround

Polymer + Solvent

from AutoPoly import Molecule

polymer = Polymer(chain_num=5, sequence=["CC[*]"] + ["[*]CC[*]"] * 48 + ["[*]CC"])
polymer.sequenceSet = polymer.sequence_set  # v1.0.0 workaround
water = Molecule(Count=100, Smiles="O", Name="water")
ethanol = Molecule(Count=20, Smiles="CCO", Name="ethanol")

Polymerization(
    name="polymer_solution",
    system=system,
    model=[polymer, water, ethanol],
    force_field="gaff"  # Use GAFF for mixed polymer+solvent
)

Ring Polymer

polymer = Polymer(
    chain_num=5,
    sequence=["[*]CC[*]"] * 30,  # ALL middle variants for rings
    topology="ring"
)
polymer.sequenceSet = polymer.sequence_set  # v1.0.0 workaround

Bead-Spring (Coarse-Grained) — Homopolymer

from AutoPoly import BeadSpringPolymer, BeadType, System

bead_A = BeadType(name="A", mass=1.0, epsilon=1.0, sigma=1.0)
system = System(out="cg_output")

polymer = BeadSpringPolymer(
    name="cg_homopolymer",
    system=system,
    n_chains=100,
    bead_types=[bead_A],
    sequence="A" * 50,       # 50 beads per chain (string format)
    topology="linear",
    pair_style="lj",         # "lj" (full LJ, cutoff 2.5σ) or "wca" (purely repulsive, cutoff 2^(1/6)σ)
    density=0.85,            # Beads per sigma^3
)
polymer.saw_generate()       # Generate with self-avoiding walk
polymer.generate_data_file() # Write in.polymer + polymer.data

Bead-Spring — Kremer-Grest (Standard Model)

polymer = BeadSpringPolymer.kremer_grest(
    name="kg_melt", system=system, n_chains=100, n_beads=50,
    topology="linear",  # density defaults to 0.74 (compress to ~0.85 during NPT)
)
polymer.saw_generate()
polymer.generate_data_file()

The kremer_grest() factory sets: FENE bonds (K=30, R0=1.5), WCA pair style (purely repulsive, cutoff 2^(1/6)σ ≈ 1.12), ε=1, σ=1, mass=1. The equilibrium bond length ~0.97σ emerges from the FENE+WCA balance — it is not an explicit parameter. Default density=0.74 for initial placement; compress to melt density ~0.85 during NPT equilibration.

Bead-Spring — AB Diblock Copolymer

bead_A = BeadType(name="A", mass=1.0, epsilon=1.0, sigma=1.0)
bead_B = BeadType(name="B", mass=1.0, epsilon=1.5, sigma=1.0)

polymer = BeadSpringPolymer(
    name="ab_diblock",
    system=system,
    n_chains=50,
    bead_types=[bead_A, bead_B],
    sequence=[("A", 20), ("B", 20)],  # Block tuple format: 20A + 20B
    density=0.85,
)
polymer.saw_generate()
polymer.generate_data_file()
# A-B cross interactions use Lorentz-Berthelot mixing automatically.
# To override: edit pair_coeff in generated in.polymer

Sequence formats for multi-bead-type CG:

• String: "AAABBB" → 3 A then 3 B
• Block tuples: [("A", 20), ("B", 30)] → 20 A then 30 B
• Explicit list: ["A", "B", "A", "B"] → alternating

Pair styles: "lj" (full LJ, cutoff 2.5σ) or "wca" (purely repulsive WCA, cutoff 2^(1/6)σ ≈ 1.12σ). Use "wca" for Kremer-Grest-style athermal melts.

Bead-Spring — Angle Potentials

from AutoPoly import AngleType

# Custom angle types for specific bead triplets
angle_AB = AngleType(triplet=("A", "B", "A"), k=25.0, theta0=120.0)

polymer = BeadSpringPolymer(
    ...,
    use_angles=True,            # Enable angle potentials
    default_k_angle=10.0,       # Default for unspecified triplets
    default_theta0=180.0,       # Default equilibrium angle
    angle_types=[angle_AB],     # Override specific triplet angles
)

Monte Carlo Placement Options (Atomistic)

For realistic initial configurations (avoiding overlaps):

Polymerization(
    name="realistic_config",
    system=system,
    model=[polymer],
    force_field="oplsaa",
    placement_method="mc_random",       # "grid" or "mc_random"
    use_mc_chain_growth=True,           # Self-avoiding walk chain growth
    mc_max_attempts=10000,              # Max placement attempts
    mc_monomer_density=0.085,           # Target density (monomers/ų)
    mc_bond_angle_min=50.0,             # Min bond deflection angle
    mc_bond_angle_max=90.0              # Max bond deflection angle
)

mc_monomer_density is in monomers/ų (not g/cm³). To convert from g/cm³:

monomers/ų = (ρ_g/cm³ × Nₐ) / (M_monomer × 10²⁴)

Example: PE at 0.85 g/cm³ with M=28.05 g/mol → 0.85 × 6.022e23 / (28.05 × 1e24) ≈ 0.0182 monomers/ų

CLI Usage

AutoPoly also provides a command-line interface:

# Show capabilities and limits
autopoly info

# Validate config without generating
autopoly validate config.json

# Generate polymer system from config
autopoly generate config.json

# Describe a SMILES structure
autopoly describe "CC[*]"

Config file format (JSON):

{
  "type": "atomistic",
  "name": "my_polymer",
  "force_field": "oplsaa",
  "polymers": [{
    "chain_num": 10,
    "sequence": ["CC[*]", "[*]CC[*]", "[*]CC[*]", "[*]CC"],
    "topology": "linear",
    "tacticity": "atactic"
  }],
  "placement_method": "mc_random",
  "mc_max_attempts": 10000
}

Agent API (Config-Driven)

For programmatic/automated use:

from AutoPoly import agent

info = agent.info()                                    # Capabilities and limits
result = agent.validate(config)                        # Validate without generating
result = agent.generate(config)                        # Generate system
desc = agent.describe_smiles("CC[*]")                  # Describe SMILES
suggestions = agent.suggest_force_field(["CC[*]"])     # Force field suggestions

Output Files

After successful generation:

output_dir/project_name/
├── moltemplate/              # Intermediate files
│   ├── monomer_*.lt          # Monomer templates
│   ├── poly_*.lt             # Chain templates
│   └── system.lt             # System definition
├── system.data               # LAMMPS data file (atoms, bonds, coordinates)
├── system.in                 # LAMMPS input script
├── system.in.settings        # Force field parameters
└── system.in.charges         # Partial charges (if applicable)

The generated system.data can be used directly with LAMMPS. The system.in provides a starting simulation script that you can customize.

Bead-Spring (CG) output is different: in.polymer + polymer.data in the output directory.

Key Constraints

• MAX_SEQUENCE_LENGTH = 10000 (max monomers per chain)
• MAX_UNIQUE_MONOMERS = 100 (max unique monomer types)
• DOP is automatically len(sequence) — no separate DOP parameter
• Wildcards must be [*] (with brackets), not bare *
• Molecule uses regular SMILES (no wildcards); Polymer uses complement SMILES (with wildcards)
• v1.0+ uses snake_case (e.g., chain_num), not CamelCase (old ChainNum)

Common Mistakes

1. Using first/last variants in ring polymers — rings need ALL middle variants
2. Wrong wildcard count — first=1, middle=2, last=1
3. Bare * instead of [*] — always use brackets
4. Using Molecule for polymers — use Polymer with complement SMILES
5. Forgetting GAFF charge calculation — GAFF needs AM1-BCC charges via Antechamber
6. Using old API — v1.0+ uses chain_num/sequence, not ChainNum/Sequence/DOP
7. BeadSpringPolymer: using chain_length — doesn't exist; chain length comes from sequence
8. BeadSpringPolymer: using write_lammps_files() — doesn't exist; use saw_generate() then generate_data_file()
9. BeadSpringPolymer: missing name and system — both are required params
10. Putting mc_monomer_density in g/cm³ — units are monomers/ų (see conversion above)
11. BeadSpringPolymer: accessing _pair_coeffs to override cross-interactions — private API; edit pair_coeff lines in the generated in.polymer file instead
12. BeadSpringPolymer: using pair_style="lj" for Kremer-Grest — KG model requires "wca" (purely repulsive); kremer_grest() factory sets this automatically
13. Atomistic: missing sequenceSet workaround — AutoPoly v1.0.0 requires polymer.sequenceSet = polymer.sequence_set before Polymerization(), otherwise crashes with AttributeError

Integration with LAMMPS Skill

AutoPoly generates the initial structure. For the simulation workflow:

1. AutoPoly → generates system.data, system.in, system.in.settings
2. LAMMPS skill → validates/customizes the simulation script, adds analysis, etc.

References

For detailed documentation on specific topics:

• references/complement_smiles.md — Full complement SMILES guide with all monomers
• references/force_fields.md — Detailed force field comparison and selection
• references/api.md — Complete API reference for all classes
• references/troubleshooting.md — Common errors and solutions