PEO-LiTFSI with AmberTools¶

Parameterize the ions with antechamber, grow PEO chains with tleap, and assemble a PEO-LiTFSI electrolyte at target density — a complete AmberTools workflow driven from MolPy.

External dependencies

This guide requires AmberTools (via conda), RDKit, and molcrafts-molpack. Without AmberTools, no code on this page will run.

Setting up AmberTools

Install AmberTools in a dedicated conda environment:

conda create -n AmberTools25 -c conda-forge ambertools=25
conda activate AmberTools25
# Verify installation
which antechamber # should print a path
which tleap # should print a path

MolPy's wrapper classes activate the conda environment automatically when running commands, so you do not need to keep it active in your shell. The env="AmberTools25" parameter in the code below tells the wrapper which environment to activate.

If you use a different environment name, replace "AmberTools25" throughout this guide.

Workflow overview¶

The workflow begins with parameterization of TFSI, the anion, using the standard Amber small-molecule sequence of antechamber, parmchk2, and tleap. Li⁺ is treated separately: its nonbonded parameters are taken from Åqvist (1990) and written to an frcmod file. With both ions parameterized, PEO chains are built using AmberPolymerBuilder, which wraps prepgen and tleap internally. The resulting component force fields are then merged, molpack places the molecules at the target density, and the final system is exported to LAMMPS.

Antechamber assigns GAFF types and BCC charges to TFSI¶

The Amber workflow for small molecules is: antechamber (assign types + charges) → parmchk2 (missing parameters) → tleap (topology + coordinates).

from pathlib import Path
import molpy as mp
from molpy.conformer import Conformer
from molpy.io.writers import write_pdb
from molpy.wrapper import AntechamberWrapper, Parmchk2Wrapper, TLeapWrapper

output_dir = Path("07_output")
ions_dir = output_dir / "ions"
ions_dir.mkdir(parents=True, exist_ok=True)

# Create TFSI from SMILES and generate 3D coordinates
tfsi = mp.io.read_smiles("O=S(=O)(C(F)(F)F)[N-]S(=O)(=O)C(F)(F)F")
tfsi = Conformer(add_hydrogens=False, seed=42).generate(tfsi)[0]

# Write PDB for antechamber input
write_pdb(ions_dir / "tfsi.pdb", tfsi.to_frame())

# docs: skip — AmberTools offline electrolyte workflow; not unit-tested
conda_env = "AmberTools25"

# Step 1: antechamber — assign GAFF types and BCC charges
ac = AntechamberWrapper(
 name="antechamber", workdir=ions_dir, env=conda_env, env_manager="conda"
)
ac.atomtype_assign(
 input_file=(ions_dir / "tfsi.pdb").absolute(),
 output_file=(ions_dir / "tfsi.mol2").absolute(),
 input_format="pdb",
 output_format="mol2",
 charge_method="bcc",
 atom_type="gaff2",
 net_charge=-1,
)

# Step 2: parmchk2 — generate missing parameters
parmchk2 = Parmchk2Wrapper(
 name="parmchk2", workdir=ions_dir, env=conda_env, env_manager="conda"
)
parmchk2.run(args=["-i", "tfsi.mol2", "-o", "tfsi.frcmod", "-f", "mol2", "-s", "gaff2"])

# Step 3: tleap — generate prmtop and inpcrd
leap_script = """source leaprc.gaff2
TFSI = loadmol2 tfsi.mol2
loadamberparams tfsi.frcmod
saveamberparm TFSI tfsi.prmtop tfsi.inpcrd
quit
"""
(ions_dir / "tfsi_leap.in").write_text(leap_script)

tleap = TLeapWrapper(name="tleap", workdir=ions_dir, env=conda_env, env_manager="conda")
tleap.run(args=["-f", "tfsi_leap.in"])

Li⁺ needs no charge calculation — literature parameters go directly into an frcmod file¶

Li⁺ has no bonded terms and no partial charges to compute, so antechamber is not needed. Instead, write the nonbond parameters from Åqvist (1990) directly into an frcmod file and create the prmtop with tleap.

Li⁺ nonbond parameters — Åqvist (1990), J. Phys. Chem. 94, 8021–8024, DOI: 10.1021/j100384a009. These were fitted to hydration free energies and are the standard choice for polymer electrolyte simulations with GAFF.

Parameter	Value
Rmin/2	1.137 Å
ε	0.0183 kcal/mol

# docs: skip — AmberTools offline electrolyte workflow; not unit-tested
from molpy.io import read_amber

li_dir = output_dir / "li"
li_dir.mkdir(parents=True, exist_ok=True)

# Write Åqvist (1990) frcmod — NONBON uses Rmin/2 and epsilon
li_frcmod = """Li+ Aqvist 1990 parameters
MASS
LI 6.941 0.0000000

BOND

ANGLE

DIHE

IMPROPER

NONBON
 LI 1.137 0.0183

"""
(li_dir / "li.frcmod").write_text(li_frcmod)

# Minimal mol2 for a single Li+ atom (net charge = +1)
li_mol2 = """@<TRIPOS>MOLECULE
LIT
 1 0 0 0 0
SMALL
USER_CHARGES

@<TRIPOS>ATOM
 1 LI 0.0000 0.0000 0.0000 LI 1 LIT 1.000000
@<TRIPOS>BOND
"""
(li_dir / "li.mol2").write_text(li_mol2)

# tleap: generate prmtop for Li+
li_leap = """source leaprc.gaff2
loadamberparams li.frcmod
LIT = loadmol2 li.mol2
saveamberparm LIT li.prmtop li.inpcrd
quit
"""
(li_dir / "li_leap.in").write_text(li_leap)

tleap_li = TLeapWrapper(
 name="tleap", workdir=li_dir, env=conda_env, env_manager="conda"
)
tleap_li.run(args=["-f", "li_leap.in"])

li_frame, li_ff = read_amber(li_dir / "li.prmtop", li_dir / "li.inpcrd")
print(
 f"Li+: {li_frame['atoms'].nrows} atom, charge={li_frame['atoms']['charge'][0]:.1f}"
)

MolPy chemistry defines the Amber residue variants¶

AmberTools does not have its own port or leaving-group semantics. Mark the monomer with fields.SITE and define the same Reaction used by PolymerBuilder; the Amber backend first compiles that molecular product, then translates its connection atoms and deleted atoms into prepgen HEAD, CHAIN, and TAIL inputs.

from molpy.builder.assembly import SiteMap
from molpy.conformer import Conformer

eo, _ = Conformer(add_hydrogens=True, seed=42).generate(mp.io.read_smiles("COC"))
SiteMap(eo).label_elements("C", "a", "b")

STITCH = mp.Reaction("[C;%a:1][H].[C;%b:2][H]>>[C:1][C:2]")
library = {"EO": eo}

AmberPolymerBuilder runs the full Amber pipeline internally¶

AmberPolymerBuilder wraps the monomer library, connector rules, and Amber tool chain (prepgen + tleap) into one builder that produces fully parameterized chains. Each unique chain length writes its Amber intermediate files into its own subdirectory under work_dir to prevent file conflicts.

# docs: skip — AmberTools offline electrolyte workflow; not unit-tested
from molpy.builder.polymer.ambertools import AmberPolymerBuilder

polymer_dir = output_dir / "polymer"
polymer_dir.mkdir(exist_ok=True)

builder = AmberPolymerBuilder(
 library=library,
 reaction=STITCH,
 force_field="gaff2",
 charge_method="bcc",
 env="AmberTools25",
 env_manager="conda",
 work_dir=polymer_dir,
)

result = builder.build("{[#EO]|10}")

AmberPolymerBuilder.build() internally runs antechamber, parmchk2, prepgen, and tleap. The result carries the polymer Frame, ForceField, and paths to the intermediate Amber files.

# docs: skip — AmberTools offline electrolyte workflow; not unit-tested
peo_frame = result.frame
peo_ff = result.forcefield
print(f"PEO 10-mer: {peo_frame['atoms'].nrows} atoms")

Merging three force fields before packing prevents type conflicts¶

Merging is done before packing rather than after because packing operates on coordinates only — it has no awareness of force field types. If two components share an atom type name with different parameters, a post-packing merge would silently overwrite one of them. Merging first makes any type name collision an error before coordinates are generated.

# docs: skip — AmberTools offline electrolyte workflow; packing via molpack
import numpy as np
from molpy.io import read_amber
from molpack import InsideBoxRestraint, Molpack, Target

# Read TFSI from Amber files generated in Stage 1
tfsi_frame, tfsi_ff = read_amber(
 ions_dir / "tfsi.prmtop",
 ions_dir / "tfsi.inpcrd",
)

# Merge all three force fields: PEO + TFSI + Li+
combined_ff = peo_ff.merge(tfsi_ff).merge(li_ff)

# Pack system with molpack
box_size = 60.0
box = InsideBoxRestraint([0.0, 0.0, 0.0], [box_size] * 3)
targets = [
 Target(peo_frame, count=3).with_name("peo").with_restraint(box),
 Target(li_frame, count=10).with_name("li").with_restraint(box),
 Target(tfsi_frame, count=10).with_name("tfsi").with_restraint(box),
]
system = Molpack().with_seed(12345).pack(targets, max_loops=200)
system.box = mp.Box.cubic(box_size)

Exporting skips pair_style because long-range electrostatics need it in the script¶

# docs: skip — AmberTools offline electrolyte workflow; not unit-tested
from molpy.io.writers import write_lammps_data, write_lammps_forcefield

lammps_dir = output_dir / "lammps"
lammps_dir.mkdir(exist_ok=True)
write_lammps_data(lammps_dir / "system.data", system, atom_style="full")
write_lammps_forcefield(lammps_dir / "system.ff", combined_ff, skip_pair_style=True)

skip_pair_style=True omits the pair_style line from the force-field file. This is required when using kspace (long-range electrostatics), because the pair_style must be set by the simulation input script rather than the force-field file.

Troubleshooting¶

Symptom	Check
Antechamber fails	Verify PDB has correct atom names and no duplicate IDs
TFSI charge wrong	Use `charge_method="bcc"` and verify `-nc -1`
tleap fails for Li⁺	Confirm the mol2 atom type (`LI`) matches the frcmod NONBON entry
Polymer build fails	Check port markers in monomer SMILES
Force field merge conflict	Inspect atom type names for collisions between PEO and TFSI
Packing fails	Increase box size or reduce molecule count