Polymer Network Analysis

Hi LAMMPS Users,

My name is Sam Edgecombe and I work as a researcher at the department of Theoretical Chemistry at Lund University in Sweden. I am researching into the swelling and mechanical properties of polymer networks

I am fairly new to LAMMPS, but have managed to create some end-linked polymer networks using the the fix bond/create command, but now I need to analyse the network for active strands, number of loose ends, number of loops etc. I wonder if anyone using lammps knows how to/has some source code in order to do this kind of analysis?

If so then I would be grateful to know how to do this.

Yours

Sam Edgecombe

You can create a restart file at the end of your simulation and use
the "restart2data" utility to create a lammps DATA file, and then load
it with some analysis program (like networkx). Restart2data is
discussed here:
http://lammps.sandia.gov/doc/Section_tools.html#restart

Then you can use a text editor to copy the "Bonds" section of the data
file and save it to another file. If you have to do this many times,
you can use:
   extract_lammps_data.py "Bonds" < DATA_FILE > bonds.dat
(See attached python script. Alternately I bet there's a nicer way to
do this using pizza.py)
  The "Bonds" section of a LAMMPS data file is a 4-column text file
containing a list of bonds ("edges"), and bond-types
("edge-attributes"), and the atom pairs. This text should be very
easy to parse. (See example script below.) If you are comfortable
with python you can try loading this file in python in networkx:
http://networkx.lanl.gov/

  For example, to load the bonds data from a 4-column ascii file, try this:

import networkx as nx
fbonds = open('bonds.dat', 'r') #<--contains only the Bonds section
G=nx.DiGraph()

for line in fbonds:
    tokens=line.strip().split()
    if len(tokens) == 0:
        continue
    # The bond-type is in column 2.
    bondtype = int(tokens[1])
    # The atom-IDs are in columns 3 and 4.
    atomid1 = int(tokens[2])
    atomid2 = int(tokens[3])
    G.add_edge(atomid1, atomid2, type=bondtype)

    # ...or just use "G.add_edge(atomid1, atomid2)"
    # if you don't care bout bond type
    # (If multiple bonds connect same pair, try using
    # MultiDiGraph instead of DiGraph. Hopefully not.)

For more information look here:
http://networkx.lanl.gov/reference/generated/networkx.DiGraph.add_edge.html

Once loaded this way, there are many things you can do.

   ...If anyone needs it, there is a way to add angular interactions
shortly after bonds have been created. I have a script which can add
angle, dihedral, and improper interactions automatically to a data
file (according to atom&bond type). This could be run periodically
from within LAMMPS using the "shell" command. For a simulation which
uses ""fix bond/create"", this would automatically add angular
interactions after the bonds had been formed. (Warning: this is slow
and you should probably not run it more often than once every 1000
timesteps. It does not yet work with "fix bond/break", but I can add
that feature.)

Cheers!
Andrew

P.S.
My one complaint with "networkx" is that the subgraph
search/isomorphism algorithm, ("VF2"), is very slow (~O(n1+n2), per
iteration), when one graph is much smaller than the other (n2<<n1).
For this specific case, I have a faster (~O(min(n1,n2)), per
iteration) version I can share if anybody needs it.

extract_lammps_data.py (4.29 KB)