LAMMPS reactive deformation of a single polyethylene chain

From EVOCD
Jump to: navigation, search

Contents

Abstract

Figure 1. Ruptured polymer chain visualized in OVITO.

In this tutorial, molecular dynamics simulation in LAMMPS is used to show the deformation of a single polymer chain using the Modified Embedded Atom Method (MEAM) potential for hydrocarbons (C/H). This example script first performs an energy minimization and equilibration on the polymer chain. Then, the chain is deformed by fixing one end of the chain to an imaginary wall while applying a velocity to the other end of the chain. With the utilization of the reactive MEAM potential we can directly observe bond rupture caused by deformation. Using OVITO, the deformation and chain rupture can be visualized.

Author(s): Andrew Bowman, Sungkwang Mun.

Corresponding Author: Andrew Bowman



Input File

Figure 2. Undeformed polymer chain visualized in OVITO.

The polymer chain of 62 atoms (20 carbon, 42 hydrogen) was prepared in Accelrys Materials Studio. For the MEAM potential, each atom, including hydrogen atoms, must be explicitly represented in the input file. Additional open source polymer creation tools are available or a custom code can be used. For this simulation the polymer was built along the y-direction as indicated in Fig. 2.

The input data file is shown below:

# nc1cl10 

62 atoms
-2.377300 5.360300 xlo xhi
0.590000 21.363700 ylo yhi
-2.163700 3.024900 zlo zhi

2 atom types

Masses

1 12.0111
2 1.0079

Atoms

    1    1     1.000000     1.540000     0.000000
    2    1     2.456100     2.041400     0.000000
    3    1     1.538500     4.082700    -1.130600
    4    1     2.456100     3.581400     0.000000
    5    1     1.538500     5.622700    -1.130600
    6    1     0.370100     6.124100    -0.261600
    7    1     0.370100     7.664100    -0.261600
    8    1    -0.875500     8.165500     0.492500
    9    1    -0.875500     9.705500     0.492500
   10    1    -1.359100    10.206800    -0.881000
   11    1    -1.359100    11.746800    -0.881000
   12    1     0.095700    12.248200    -0.944400
   13    1     0.095700    13.788200    -0.944400
   14    1     1.468100    14.289600    -0.458000
   15    1     2.898200    16.331000    -0.183900
   16    1     1.468100    15.829600    -0.458000
   17    1     2.898200    17.871000    -0.183900
   18    1     3.201500    18.372300     1.240300
   19    1     4.452800    20.413700     1.985000
   20    1     3.201500    19.912300     1.240300
   21    2     1.000000     0.590000     0.000000
   22    2     0.480000     1.898100     0.900700
   23    2     0.480000     1.898100    -0.900700
   24    2     2.964000     1.666300    -0.900700
   25    2     2.964000     1.666300     0.900700
   26    2     1.910100     3.724600    -2.102000
   27    2     0.511400     3.724600    -0.966700
   28    2     2.076800     3.956500     0.962000
   29    2     3.475500     3.956500    -0.173300
   30    2     1.418200     5.980900    -2.163700
   31    2     2.493300     5.980900    -0.718200
   32    2     0.500100     5.749000     0.764300
   33    2    -0.575000     5.749000    -0.681200
   34    2     0.348500     8.022200    -1.301400
   35    2     1.281400     8.022200     0.239700
   36    2    -0.843600     7.790400     1.526100
   37    2    -1.776500     7.790400    -0.015000
   38    2     0.146800    10.063600     0.683900
   39    2    -1.552500    10.063600     1.282100
   40    2    -2.377300     9.831700    -1.061000
   41    2    -0.678100     9.831700    -1.659200
   42    2    -1.839400    12.105000     0.041600
   43    2    -1.917800    12.105000    -1.758200
   44    2     0.563900    11.873100    -1.866400
   45    2     0.642300    11.873100    -0.066600
   46    2    -0.695400    14.146300    -0.269100
   47    2    -0.093600    14.146300    -1.967100
   48    2     2.247700    13.914500    -1.137300
   49    2     1.646000    13.914500     0.560700
   50    2     3.239400    15.972800     0.798600
   51    2     3.578500    15.972800    -0.970600
   52    2     1.138800    16.204700    -1.438200
   53    2     0.799700    16.204700     0.331100
   54    2     1.908900    18.229100    -0.504900
   55    2     3.670800    18.229100    -0.880100
   56    2     4.188300    17.997200     1.549400
   57    2     2.426400    17.997200     1.924700
   58    2     4.452800    21.363700     1.985000
   59    2     4.439000    20.055600     3.024900
   60    2     5.360300    20.055600     1.476900
   61    2     3.225700    20.287400     0.206500
   62    2     2.304400    20.287400     1.754500

LAMMPS Script

Below is the script used for the actual simulation. The Aug 2015 version of LAMMPS was used. The MEAM potential for saturated hydrocarbons was used for this simulation. The required library.meam and CH.meam parameters can be found here. Minimization and equilibration were performed first, then atom groups were defined to allow for chain stretching at a constant velocity. The compute " stress/atom" was used to visualize the stress distribution along the chain. This input script was run using the February 2014 version of LAMMPS. Changes in some commands in more recent versions may require revision of the input script. This script runs the simulation with the previously discussed input file.

#Initialization

units metal
boundary s s s
atom_style atomic
variable filename string nc1cl10
read_data ${filename}.dat

# Meam Potential for Saturated Hydrocarbons Information

pair_style meam
pair_coeff * * library.meam C H CH.meam C H

##################################################################################
#Minimization

minimize 0.0 1.0e-8 100000 100000

# Equilibration

reset_timestep 0
velocity 		all create 0.0 1231 dist gaussian
fix 			1 all nve
fix 			2 all langevin 0 0 0.001 904297
thermo_style 	custom step etotal press vol temp density pxx pyy pzz
thermo	  	100
timestep 		0.0003
run 			20000
unfix 1
unfix 2

######################################################################
# Deformation

run 0
reset_timestep 0

#Define Variables

variable 	yl equal "ylo+1"
variable 	yh equal "yhi-1"
variable 	yh2 equal "yhi*3"
variable 	y_step equal "yhi/2"

variable 	p3 equal "ly"
variable 	L0 equal ${p3}
variable 	strain equal "(ly - v_L0)/v_L0"
variable 	p1 equal v_strain
variable 	p2 equal "-pyy/10"

# define groups

region   	1 block INF INF INF ${yl} INF INF
group		wall region 1
region	2 block INF INF ${yh} INF INF INF
group		pull region 2
group		boundary union wall pull
group		mobile subtract all wall pull

# computes

compute peratom all stress/atom NULL pair	

# initial velocities

velocity	mobile create 0.01 887723 
velocity	pull set 0.0 .5 0.0
velocity	mobile ramp vy 0 ${y_step} y 0 ${yh2} sum yes

# fixes

fix	f1 all nve
fix 	ftemp all langevin 25 25 0.05 904297 
fix	f2 boundary setforce NULL 0.0 0.0

# Dump
dump 100 all custom 1000 stress_per_atom.txt id type x y z c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]
fix def1 all print 100 "${p1} ${p2}" file Stress_strain.def1.txt screen no

thermo_style custom step temp pyy lx ly lz v_strain
thermo 1000
timestep 0.0005
run   50000

print "All Done"

LAMMPS Logfile

Here is the logfile produced by LAMMPS during the simulation.

LAMMPS (1 Feb 2014)
#Initialization

units metal
boundary s s s
atom_style atomic
variable filename string nc1cl10
read_data ${filename}.dat
read_data nc1cl10.dat
  orthogonal box = (-2.3773 0.59 -2.1637) to (5.3603 21.3637 3.0249)
  2 by 5 by 1 MPI processor grid
  reading atoms ...
  62 atoms

# Meam Potential for Saturated Hydrocarbons Information

pair_style meam
pair_coeff * * library.meam C H CH.meam C H

dump 	d1 all movie 100 movie.mp4 type type size 640 480 zoom 2.5

##################################################################################
#Minimization

minimize 0.0 1.0e-8 100000 100000
WARNING: Resetting reneighboring criteria during minimization (../min.cpp:173)
Memory usage per processor = 5.7362 Mbytes
Step Temp E_pair E_mol TotEng Press Volume 
       0            0   -256.88642            0   -256.88642    16039.246    834.50871 
    2275            0   -258.98482            0   -258.98482   -15.964463    880.60809 
Loop time of 9.09964 on 10 procs for 2275 steps with 62 atoms

Minimization stats:
  Stopping criterion = linesearch alpha is zero
  Energy initial, next-to-last, final = 
        -256.886423993     -258.984817441     -258.984817441
  Force two-norm initial, final = 17.9272 0.00267172
  Force max component initial, final = 8.69935 0.00144996
  Final line search alpha, max atom move = 0.00195312 2.83195e-06
  Iterations, force evaluations = 2275 18266

Pair  time (%) = 5.0886 (55.9209)
Neigh time (%) = 3.75032e-05 (0.00041214)
Comm  time (%) = 1.83001 (20.1108)
Outpt time (%) = 0.740981 (8.14297)
Other time (%) = 1.44001 (15.8249)

Nlocal:    6.2 ave 13 max 0 min
Histogram: 2 0 1 2 0 0 3 0 0 2
Nghost:    35.4 ave 53 max 18 min
Histogram: 1 1 1 1 0 2 2 1 0 1
Neighs:    92.8 ave 190 max 0 min
Histogram: 2 0 1 0 2 2 0 1 0 2
FullNghs:  185.6 ave 406 max 0 min
Histogram: 2 0 2 1 1 0 1 1 1 1

Total # of neighbors = 1856
Ave neighs/atom = 29.9355
Neighbor list builds = 1
Dangerous builds = 0

# Equilibration

reset_timestep 0
velocity 		all create 0.0 1231 dist gaussian
fix 			1 all nve
fix 			2 all langevin 0 0 0.001 904297
thermo_style 	custom step etotal press vol temp density pxx pyy pzz
thermo	  	1000
timestep 		0.0003
run 			20000
Memory usage per processor = 5.04914 Mbytes
Step TotEng Press Volume Temp Density Pxx Pyy Pzz 
       0   -258.98482   -16.335355    860.61399            0   0.54518234    4.2240815   -52.460361  -0.76978563 
    1000   -258.98482   -18.467099    860.61399 2.6373086e-07   0.54518234     4.279785   -58.922537  -0.75854527 
    2000   -258.98482   -18.408935    860.61399  2.57724e-07   0.54518234    4.3413454   -58.764647  -0.80350341 
    3000   -258.98482   -18.351508    860.61399 2.5287401e-07   0.54518234    4.3926884   -58.599195  -0.84801847 
    4000   -258.98482   -18.296634    860.61399 2.4862843e-07   0.54518234    4.4353853   -58.435041  -0.89024532 
    5000   -258.98482   -18.244933    860.61399 2.4477753e-07   0.54518234    4.4695285   -58.274125  -0.93020235 
    6000   -258.98482   -18.196254    860.61399 2.4121733e-07   0.54518234    4.4957193    -58.11642   -0.9680618 
    7000   -258.98483   -18.150042    860.61399 2.3788626e-07   0.54518234     4.514807   -57.961022   -1.0039105 
    8000   -258.98483   -18.105594    860.61399 2.3474327e-07   0.54518234    4.5277382    -57.80672   -1.0377988 
    9000   -258.98483    -18.06221    860.61399 2.3175872e-07   0.54518234    4.5354502   -57.652305   -1.0697762 
   10000   -258.98483    -18.01927    860.61399 2.2891003e-07   0.54518234     4.538808   -57.496715   -1.0999023 
   11000   -258.98483   -17.976252    860.61399 2.2617938e-07   0.54518234    4.5385756   -57.339085   -1.1282465 
   12000   -258.98483    -17.93274    860.61399 2.2355234e-07   0.54518234    4.5354099   -57.178746   -1.1548854 
   13000   -258.98483   -17.888417    860.61399 2.2101696e-07   0.54518234    4.5298616   -57.015212   -1.1798998 
   14000   -258.98483   -17.843046    860.61399 2.1856326e-07   0.54518234    4.5223879   -56.848154   -1.2033727 
   15000   -258.98483   -17.796462    860.61399 2.1618277e-07   0.54518234     4.513367   -56.677366   -1.2253871 
   16000   -258.98484   -17.748553    860.61399 2.1386829e-07   0.54518234    4.5031093   -56.502744   -1.2460252 
   17000   -258.98484   -17.699254    860.61399 2.1161362e-07   0.54518234    4.4918689   -56.324264   -1.2653669 
   18000   -258.98484   -17.648533    860.61399 2.094134e-07   0.54518234    4.4798528   -56.141963   -1.2834899 
   19000   -258.98484   -17.596388    860.61399  2.07263e-07   0.54518234      4.46723   -55.955925   -1.3004684 
   20000   -258.98484   -17.542835    860.61399 2.0515839e-07   0.54518234    4.4541384   -55.766269   -1.3163738 
Loop time of 21.1757 on 10 procs for 20000 steps with 62 atoms

Pair  time (%) = 5.64706 (26.6677)
Neigh time (%) = 0 (0)
Comm  time (%) = 2.03111 (9.59171)
Outpt time (%) = 7.50698 (35.4509)
Other time (%) = 5.99055 (28.2897)

Nlocal:    6.2 ave 13 max 0 min
Histogram: 2 0 1 2 0 0 3 0 0 2
Nghost:    35.4 ave 53 max 19 min
Histogram: 1 1 1 1 0 2 2 1 0 1
Neighs:    93 ave 192 max 0 min
Histogram: 2 0 1 1 1 2 0 1 0 2
FullNghs:  186 ave 406 max 0 min
Histogram: 2 0 2 1 1 0 0 3 0 1

Total # of neighbors = 1860
Ave neighs/atom = 30
Neighbor list builds = 0
Dangerous builds = 0
unfix 1
unfix 2

######################################################################
# Deformation

run 0.
WARNING: No fixes defined, atoms won't move (../verlet.cpp:54)
Memory usage per processor = 5.04914 Mbytes
Step TotEng Press Volume Temp Density Pxx Pyy Pzz 
   20000   -258.98484    -17.55251     860.1396 2.0515839e-07   0.54548301    4.4565949   -55.797026   -1.3170998 
Loop time of 1.80244e-05 on 10 procs for 0 steps with 62 atoms

Pair  time (%) = 0 (0)
Neigh time (%) = 0 (0)
Comm  time (%) = 0 (0)
Outpt time (%) = 0 (0)
Other time (%) = 1.80244e-05 (100)

Nlocal:    6.2 ave 13 max 0 min
Histogram: 2 0 1 2 0 0 3 0 0 2
Nghost:    35.4 ave 53 max 19 min
Histogram: 1 1 1 1 0 2 2 1 0 1
Neighs:    93 ave 192 max 0 min
Histogram: 2 0 1 1 1 2 0 1 0 2
FullNghs:  186 ave 406 max 0 min
Histogram: 2 0 2 1 1 0 0 3 0 1

Total # of neighbors = 1860
Ave neighs/atom = 30
Neighbor list builds = 0
Dangerous builds = 0
reset_timestep 0

#Define Variables

variable 	yl equal "ylo+1"
variable 	yh equal "yhi-1"
variable 	yh2 equal "yhi*3"
variable 	y_step equal "yhi/2"

variable 	p3 equal "ly"
variable 	L0 equal ${p3}
variable 	L0 equal 20.774388695406
variable 	strain equal "(ly - v_L0)/v_L0"
variable 	p1 equal v_strain
variable 	p2 equal "-pyy/10"

# define groups

region   	1 block INF INF INF ${yl} INF INF
region   	1 block INF INF INF 1.61162482988567 INF INF
group		wall region 1
1 atoms in group wall
region	2 block INF INF ${yh} INF INF INF
region	2 block INF INF 20.3860135252917 INF INF INF
group		pull region 2
1 atoms in group pull
group		boundary union wall pull
2 atoms in group boundary
group		mobile subtract all wall pull
60 atoms in group mobile

# computes

compute peratom all stress/atom pair

# initial velocities

velocity	mobile create 0.01 887723
velocity	pull set 0.0 .5 0.0
velocity	mobile ramp vy 0 ${y_step} y 0 ${yh2} sum yes
velocity	mobile ramp vy 0 10.6930067626458 y 0 ${yh2} sum yes
velocity	mobile ramp vy 0 10.6930067626458 y 0 64.1580405758751 sum yes

# fixes

fix	f1 all nve
fix 	ftemp all langevin 25 25 0.05 904297
fix	f2 boundary setforce NULL 0.0 0.0

# Dump
dump 100 all custom 1000 stress_per_atom.txt id type x y z c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]
fix def1 all print 100 "${p1} ${p2}" file PE2_ss_nc10_cl1000.def1.txt screen no

thermo_style custom step temp pyy lx ly lz v_strain
thermo 10000
timestep 0.0005
run   50000
Memory usage per processor = 6.66095 Mbytes
Step Temp Pyy Lx Ly Lz strain 
       0    7.8956181    175.93613    7.6947242    20.774389    5.3808097 1.3681129e-15 
   10000    24.371339   -9630.1918    6.3750014    22.806131    5.8410614  0.097800348 
   20000    27.667952   -15903.107    5.0358511    25.508789    5.5483935   0.22789599 
   30000    26.446686   -151820.66    2.8828797    28.219945    3.9085843   0.35840073 
   40000    25.713185   -224272.78    2.5381237    30.770621    4.4009193    0.4811806 
   50000    38.778302   -90892.535    2.8497228    33.235559    3.9395652   0.59983331 
Loop time of 52.1821 on 10 procs for 50000 steps with 62 atoms

Pair  time (%) = 11.5983 (22.2266)
Neigh time (%) = 0.000626159 (0.00119995)
Comm  time (%) = 3.42308 (6.55987)
Outpt time (%) = 20.3336 (38.9666)
Other time (%) = 16.8265 (32.2458)

Nlocal:    6.2 ave 10 max 1 min
Histogram: 1 0 1 0 1 2 2 1 1 1
Nghost:    27.6 ave 40 max 16 min
Histogram: 2 0 2 1 0 1 0 2 1 1
Neighs:    71.1 ave 132 max 6 min
Histogram: 2 0 0 1 1 2 2 1 0 1
FullNghs:  142.2 ave 279 max 6 min
Histogram: 2 0 0 0 4 0 1 2 0 1

Total # of neighbors = 1422
Ave neighs/atom = 22.9355
Neighbor list builds = 26
Dangerous builds = 0

print "All Done"
All Done

LAMMPS dump files

The following dump files were produced during the simulation and can be used to make a movie. The stress-strain data was also dumped and can be used to make a stress-strain curve.


# Dump
dump 100 all custom 1000 stress_per_atom.txt id type x y z c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]
fix def1 all print 100 "${p1} ${p2}" file Stress_strain.def1.txt screen no

Creating a Movie in OVITO

  • Open Ovito.
  • Select "File/Open Local File," Choose the LAMMPS dump file you wish to animate, and "Open." For this tutorial select "stress_per_atom.txt".
MEAM polymer ani.gif
    • You must tell Ovito what each column of data in the dumpfile is. Often, Ovito can do this for you if you click "Auto-assign columns."
    • For the c_peratom columns choose stress tensor as the particle property. The 6 stress components are dumped in the order: XX, YY, ZZ, XY, XZ, YZ.
  • Use the Modifier List to add the desired modifications to your images and future animation. For this example, the "Color Coding" modifier was used to visualize the stress.
    • Under "Color Coding", choose YY as the property. This will animate the stress in the Y-direction (deformation direction).
  • Test the animation by using the animation tools below the graphics.
  • Click the "Render" tab.
  • Select "Complete animation", and under background select "transparent". Select save to file, choose output type (.gif for this example) and choose an output filename in "Output Filename."
  • Select the view that you would like to animate and click "Render Active Viewport" above "Render settings."
  • Now you should have an animated .gif file for every frame of your simulation.

Go Back

Personal tools
Namespaces

Variants
Actions
home
Materials
Material Models
Design
Resources
Projects
Education
Toolbox