how to build polymer melts?

A LAMMPS user recently emailed me to ask how to build a large all-atom
polymer melt. (In his case, his monomers are not simple
coarse-grained one-atom beads. Instead, he is trying to build an
all-atom phenolic resin.)

I'm not sure moltemplate is the tool for the job. I've received this
question several times and I feel bad because I don't know the answer.
So I thought I'd put it to the LAMMPS community. To be frank, tools
like moltemplate are more useful when the system being simulated
consists of many copies of the same (mostly rigid) molecule (although
can customize individual molecules).

I told him that once he has constructed a reasonably well-equilibrated
polymer melt, it should be easy to add cross-links using "fix
bond/create" (if needed).

To build a polymer melt, I suggested that the he could:

1) use a molecule builder to create many copies of a monomer and
simulate the process of (step-growth) polymerization using LAMMPS.
2) use a molecule builder to create many polymers, and let them pass
through each other using small "sigma" Lennard-Jones parameters (or
"pair_style table") to weaken steric repulsion. (And then add
cross-linking if needed.)
3) Create a polymer melt of simple one-bead chains (by other means).
Later he could replace each bead in the melt with an all-atom monomer
oriented in the direction of the separation between neighboring beads
(using a molecule-builder like moltemplate or some other script).

However I suspect there are other ways...

But I am wondering if anybody knows of a good paper which has a method
section containing a good protocol to build polymer melts.
I figured somebody here would have some good advice.

Hi Andrew,

If you mean creating a crosslinked polymer network in atomic scale (polymer melt) in order to simulate it in LAMMPS, the first issue is to make a lammps data file for liquid phase polymer structure.

I suggest that he could use

  1. “Jmol” to create each monomer molecule and save it as a state coordinate file (e.g. xyz file)
  2. then randomly mix them using “Packmol” and save the coordinate file like
  3. Import the system file in VMD in order to create a lammps data file using “topo tool”
  4. assign force field parameters indicated inside of the lammps data file created in step (4)
  5. then use “fix create/bond” command and “delete_bond” command (step-growth = close distance criteria) in order to create polymer network

This method could be applied to epoxy which is polymerized by step-growth, but not for free radical polymerization.

See the “Relative reactivity volume criterion for cross-linking: Application to vinyl ester resin molecular dynamics simulations” for free radical polymerization modeling done by Materials Studio.


Changwoon Jang

Thank you Changwoon for your detailed instructions!

   Links to the relevant lammps commands:
( ?)

   Incidentally, there is a way to add/remove angle/dihedral
interactions as the bonds are forming/breaking. But it is not trivial
to do. There is a long thread discussing various ways to do that
   (read the follow-up messages on that thread. Note: the bug in has since been fixed, but the URL has moved to

Hi Andrew,

   The "polymer melt" problem is too nebulous to be solved with any off-the-shelf software. Do you want to build 3D cross linked systems or a bucket of oligomers. Very different problems. The bucket of oligomers problem is very similar to the
bucket of solvent problem and can be solved by similar means (ie build a single molecule and replicate it with packmol) but
has the specific issue of how you get long polymer chains to intertwine. The 3D cross linked problem is one I think can only be solved by writing code specific to the particular polymer.

   I have worked on the bucket of oligomers problem and that solution can be used to solve the 3D cross link problem with
some extra programming to make the cross links. I will first describe the code for coordinate generation and then I will
some examples to show how it has been used.

  I have developed a linear proximity embedding method which takes a connectivity table for a collection of atoms and
generates a set of 3D coordinates using standard bond lengths, angles and dihedrals. The connectivity table is converted to
a matrix of minimum and maximum distances between each pair of atoms - 1,2 distances are bond lengths; 1,3 distances are angles; 1,4 are the 0 and 180 degree dihedrals; non-bonded pairs - min is 1/4 sum of the vdw radii, max is the length of the side of a cube - length = cube root(mass of system/density desired). The program then iterates through all pairs of atoms solving the disance constraints several times using a learning correction after each cycle. I have published the code as part of the PubChem work I did for the ChemBioGrid project at Indiana University and you can download the code by searching for smi2sdf (that code took a smiles string can generated a sdf file of 3d coordinates). The important subroutine is: iterate. I think there may also be a Java version written by Rajarshi Guha. The code is open source and available for use. Agriafotis has published several papers on a stochastic proximity embedding method which is similar.

Some examples will show how this code is used:

3D cross linked model for oil shale kerogen:

The polymer unit starts with a small ring (single or fused, anything from benzene to a steroid) which we call a core. A long chain aliphatic it attached to the core which we call a linker. The core-linker unit is extended (polymerized) as desired. Next side chains are randomly added to cores (0, 1 or 2 to each core). The side chains are aliphatic chains similar to the links but only attached to a core at one end. This is all done to generate a connectivity table (a list of all the atoms and bonds) and there are no coordinates at this point. The connectivity table is sent to the coordinate generation routine and we get back a set of 3D coordinates that fit within the prescibed volume and satisfy the rules of organic chemistry. At this point we look at the ends of the side chains (easy to find since they are the only terminal aliphatic atoms present), check to see if any are within bonding distance of a different core (not the one the side chain attached to) and if so make a bond - thus the cross link. We can next decorate the structure with oxygen, nitrogen and/or sulfur to make other functional groups, save it as a Mol2 file and then convert the Mol2 file into a lammps input file along with the parameters/data for ReaxFF, OplsUA or OplsAA simulations. A 10 core model (~800 heavy atoms) takes less than 5 seconds. A 300 core model (10,000 heavy atoms) takes and hour and the largest model I build was 1000 cores and took a day on a desk top linux workstation. The particular models I was building had link lengths of 13 +- 6 carbons and are reasonably good models for Green River Oil Shales.

Polymer melts

A much simpler problem. Build a single oligomer - a 50 mer of polyethylene oxide (C-C-O ---) using your favorite molecular editor (I use Pcmodel) and save it as a mol2 file. Read x copies of the mol2 file into the program, keep the connectivity tables and discard the coordinates - basically concatenate the connectivity tables. Run the connectivity tables through the coordinate generator and save the 3d coordinates. The resulting bucket of oligomers has random orientations of all the chains, intertwining of the chains and random orientations of all the internal dihedral angles!!!! I have used the replicate a single molecule method in Pcmodel for many years and it works great for small, spherical molecules but long, skinny molecules are a real problem = you can either align them all the same way and avoid bad collisions, or randomly orient the molecules and have failure when the chains cross. This method solves that problem, gives you a relaxed structure to start with and save enormous amounts of time if you use Towhee.


We were interested in studing the reactions of single metal atoms with hydrocarbons using ReaxFF. So using the above methods we read in - 1,6-diphenyl hexane, Ni atoms, H2 - repeat for as many times as you wish - generate the coordinates and you get an simulation cell with randomly dispearsed metal atoms. We could build systems of hydrocarbon only, hydrocarbon plus Ni and hydrocarbon plus Ni plus H2. Conversion to lammps input files and simulation. We think the current ReaxFF parameter set for CHNi under estimates the reactivity of single Ni atoms. What you can see from watching movies of the simulation is the aggolmeration of the Ni atoms as the simulation proceeds, which is a known mechanism for the loss of catalytic activity.

Please feel free to get the code and use it in Moltemplate and send me an email if you have questions.