LAMMPS Polymer

From EVOCD
Revision as of 15:08, 1 December 2011 by Dzhuk (Talk | contribs)

Jump to: navigation, search

Contents

Abstract

In this tutorial Molecular dynamics simulation in LAMMPS is used to show what happens to a polymer chain at a certain temperature after some time. Chain's movement is caused by a molecular forces between atoms in the chain and by temperature of the chain. This factors are modeled in LAMMPS in order to show the behavior of this polymer chain. The chain was previously prepared in MATLAB. It contains 100 atoms with bound between neighbor atoms. During the simulation, firstly, the chain was left for do deform freely under the molecular forces and a random atomic fluctuation due to the temperature. After that the chain was minimized to find it's minimal energy condition. The tutorial explains basic LAMMPS commands and shows how to make a movie using AtomEye and ImageJ.

Author(s): Mark A. Tschopp, Dmitry I. Zhuk.

Corresponding Author: Mark Tschopp

Input File

Figure 1. Not-deformed polymer chain.

The polymer chain of 100 atoms was specially prepared in MATLAB. The atom's Z coordinate does not varies much, all of them are within 2Å. The distance between atoms is about 1.5Å. Basically, the chain goes from left upper to right lower corner of the box. At the right you can see the picture of the initial condition of the chain.

The data file is shown below and available for download here: https://icme.hpc.msstate.edu/mediawiki/images/e/e1/PE_cl100.txt

# Model for PE

     100     atoms
      99     bonds
      98     angles
      97     dihedrals

         1     atom types
         1     bond types
         1     angle types
         1     dihedral types

    0.0000   158.5000 xlo xhi
    0.0000   158.5000 ylo yhi
    0.0000   100.0000 zlo zhi

Masses

         1          14.02

Atoms

         1         1         1    5.6240    5.3279   51.6059
         2         1         1    7.4995    7.4810   50.2541
         3         1         1    8.2322    8.0236   51.2149
         4         1         1    9.6108    9.9075   51.7682
         5         1         1   11.5481   11.3690   50.4167
         6         1         1   12.9409   13.4562   50.2481
         7         1         1   14.4708   14.8569   50.0868
         8         1         1   16.1916   16.4790   50.5665
         9         1         1   17.1338   17.6853   51.8189
        10         1         1   19.1109   19.4000   50.3869

...

Bonds

         1         1         1         2
         2         1         2         3
         3         1         3         4
         4         1         4         5
         5         1         5         6
         6         1         6         7
         7         1         7         8
         8         1         8         9
         9         1         9        10
        10         1        10        11
 
...

Angles

         1         1         1         2         3
         2         1         2         3         4
         3         1         3         4         5
         4         1         4         5         6
         5         1         5         6         7
         6         1         6         7         8
         7         1         7         8         9
         8         1         8         9        10

...

Dihedrals

         1         1         1         2         3         4
         2         1         2         3         4         5
         3         1         3         4         5         6
         4         1         4         5         6         7
         5         1         5         6         7         8
         6         1         6         7         8         9
         7         1         7         8         9        10

Here, the first section defines the numbers of atoms, bonds, angles and dihedrals. Lower you can see the types and the box sizes. Below that are the simulation box sizes.

Then comes the Atoms section. The first column is the atom number, then comes the atom type and some other information not used in this tutorial. Last three columns is the atom's x, y and z coordinates correspondingly.

Figure 2. Angles between atoms defined in LAMMPS.
Figure 2. Dihedral angles between atoms defined in LAMMPS.


In the Angle section listed the angles between atoms. First column is the angle number, second - angle ID, then comes the atom's number between which the angle is defined. Example of the distribution of the atoms can be seen on the picture to the left. Dihedrals are a little bit more complicated. To define a diherdal four atoms are needed. Syntax are pretty much similar to the angles section. The only difference is that one more atom is needed to define a dihedral. Basically, the dihedral angle is the angle between the planes formed by 2 groups of 3 neighbor atoms. You see the example on the picture to the right.











LAMMPS Script

Below is the script used for the actual simulation. This input script was run using the November 2010 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 data file.

#####################################################
#                                                   #
#                                                   #
# Filename: in.deform.polychain.txt                 #
# Author: Mark Tschopp, 2010                        #
#                                                   #
# The methodology outlined here follows that from   #
# Hossain, Tschopp, et al. 2010, Polymer.  Please   #
# cite accordingly. The following script requires   #
# a LAMMPS data file containing the coordinates and #
# appropriate bond/angle/dihedral lists for each    #
# united atom.                                      #
#                                                   #
# Execute the script through:                       #
# lmp_exe < in.deform.polychain.txt                 #
#                                                   #
#####################################################

# VARIABLES
variable fname index PE_cl100.txt
variable simname index PE_cl100

# Initialization
units		real
boundary	f f f
atom_style	molecular
log 		log.${simname}.txt
read_data	${fname}

# Dreiding potential information
neighbor	0.4 bin
neigh_modify	every 10 one 10000
bond_style      harmonic
bond_coeff	1 350 1.53
angle_style     harmonic
angle_coeff	1 60 109.5
dihedral_style	multi/harmonic
dihedral_coeff	1 1.73 -4.49 0.776 6.99 0.0
pair_style	lj/cut 10.5
pair_coeff	1 1 0.112 4.01 10.5

compute csym all centro/atom fcc
compute peratom all pe/atom 



#####################################################
# Equilibration (Langevin dynamics at 5000 K)

velocity 	all create 5000.0 1231
fix		1 all nve/limit 0.05
fix		2 all langevin 5000.0 5000.0 10.0 904297
thermo_style	custom step temp 
thermo          10000
timestep	1
run		1000000
unfix 1
unfix 2
write_restart 	restart.${simname}.dreiding1


#####################################################
# Define Settings
compute eng all pe/atom 
compute eatoms all reduce sum c_eng 

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

dump 		1 all cfg 6 dump.comp_*.cfg id type xs ys zs c_csym c_peratom fx fy fz

reset_timestep 0 
fix 1 all nvt temp 500.0 500.0 100.0
thermo 20 
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
min_style cg
minimize 1e-25 1e-25 500000 1000000 


print "All done"

In general, this script does equilibration and minimization to the polymer chain. Polymer chain data file named 'PE_cl100.txt' should be the same directory. Line "dump 1 all cfg 6 dum..." used to output information during simulation can be moved before the equilibration part of the script to output the process of equilibration. Please note that you need to change a time step for a dump command otherwise it will give too many or very few dump files. It may be reasonable to choose a timestep of 10,000 for equilibration and 6 for minimization. That will give about 100 files for each step.

# VARIABLES
variable fname index PE_cl100.txt
variable simname index PE_cl100

# Initialization
units		real
boundary	f f f
atom_style	molecular
log 		log.${simname}.txt
read_data	${fname}

This part is just opens the data file, defines the boundary conditions, units, logfile's name, etc.

# Dreiding potential information
neighbor	0.4 bin
neigh_modify	every 10 one 10000
bond_style      harmonic
bond_coeff	1 350 1.53
angle_style     harmonic
angle_coeff	1 60 109.5
dihedral_style	multi/harmonic
dihedral_coeff	1 1.73 -4.49 0.776 6.99 0.0
pair_style	lj/cut 10.5
pair_coeff	1 1 0.112 4.01 10.5

In this section goes the information about the bonds, angles, dihedrals in a chain. Bond_style and bond_coeff defines the type on the force field between atoms and a magnitude of this fields. "1" here corresponds to the second column of the "Bonds" section of the data file. Thus, every atom pair with "1" in the second column will be having such properties during the simulation. Similar goes to the angles and dihedrals. Angle_* and dihedral_* lines defines the angles and dihedral angles between atoms in the polymer chain. Pair_* commands used to define pair potentials between atoms that are within a cutoff distance. More about this commands and a parameters can be found at SANDIA website:

http://lammps.sandia.gov/doc/Section_commands.html#cmd_5

compute csym all centro/atom fcc
compute peratom all pe/atom 

This commands are used to define which properties are to be calculated during the simulation. For example, "pe/atom" simply means to calculate the potential energy for each atom. Information about other possible properties to calculate can be found here

# Equilibration (Langevin dynamics at 5000 K)

velocity 	all create 5000.0 1231
fix		1 all nve/limit 0.05
fix		2 all langevin 5000.0 5000.0 10.0 904297
thermo_style	custom step temp 
thermo          10000
timestep	1
run		1000000
unfix 1
unfix 2
write_restart 	restart.${simname}.dreiding1


<pre>
# Equilibration (Langevin dynamics at 5000 K)

velocity 	all create 5000.0 1231
fix		1 all nve/limit 0.05
fix		2 all langevin 5000.0 5000.0 10.0 904297
thermo_style	custom step temp 
thermo          10000
timestep	1
run		1000000
unfix 1
unfix 2
write_restart 	restart.${simname}.dreiding1

The equilibration step allows the chain to deform freely under the temperature driven movements of the atoms. Velocity command add the temperature to the chain. Fix command performs the check of the system every time step and updates the velocity and positions of the atoms. Thermo commands defines thermo output during the simulation. Run command stars the dynamic transformation.



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

dump 		1 all cfg 6 dump.comp_*.cfg id type xs ys zs c_csym c_peratom fx fy fz

reset_timestep 0 
fix 1 all nvt temp 500.0 500.0 100.0
thermo 20 
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
min_style cg
minimize 1e-25 1e-25 500000 1000000 

The other step in the program is the minimization. It finds the minimum energy condition for such configuration. The parameters for minimize command includes: stopping tolerance for energy, stopping tolerance for force, max iterations of minimizer and max number of force/energy evaluations. Where stopping tolerance for energy being the first parameter and the max number of force/energy evaluations the last one.


LAMMPS Logfile

Here is the logfile produced by LAMMPS during the simulation. Note that the temperature during the equilibration does not concave and just randomly changes over time. On the other hand, the potential energy during the minimization lowers over time until it reaches the minimum for this configuration within a tolerance.

read_data	${fname}
read_data	PE_cl100.txt
  1 = max bonds/atom
  1 = max angles/atom
  1 = max dihedrals/atom
  orthogonal box = (0 0 0) to (158.5 158.5 100)
  2 by 2 by 1 processor grid
  100 atoms
  99 bonds
  98 angles
  97 dihedrals
  2 = max # of 1-2 neighbors
  2 = max # of 1-3 neighbors
  4 = max # of 1-4 neighbors
  6 = max # of special neighbors

# Dreiding potential information
neighbor	0.4 bin
neigh_modify	every 10 one 10000
bond_style      harmonic
bond_coeff	1 350 1.53
angle_style     harmonic
angle_coeff	1 60 109.5
dihedral_style	multi/harmonic
dihedral_coeff	1 1.73 -4.49 0.776 6.99 0.0
pair_style	lj/cut 10.5
pair_coeff	1 1 0.112 4.01 10.5

compute csym all centro/atom fcc
compute peratom all pe/atom 



#####################################################
# Equilibration (Langevin dynamics at 5000 K)

velocity 	all create 5000.0 1231
fix		1 all nve/limit 0.05
fix		2 all langevin 5000.0 5000.0 10.0 904297
thermo_style	custom step temp 
thermo          10000
timestep	1
run		1000000
Memory usage per processor = 5.96847 Mbytes
Step Temp 
       0         5000 
   10000    4503.8382 
   20000    4389.5108 
   30000    4137.4417 
   40000    5341.7194 
   50000    4458.9062 
   60000    5377.3337 
   70000    5004.3444 
   80000    4338.4843 
   90000    4691.7776 
  100000    4733.5682 
  110000    4701.5107 
  120000    4855.5683 
  130000    4477.8402 
  140000    4538.7404 
  150000    4588.4541 
  160000    4862.2936 
  170000    4684.8557 
  180000    4485.3563 
  190000    4759.6888 
  200000     5381.965 
  210000    5137.9384 
  220000    4812.9015 
  230000    4628.1747 
  240000    4645.6596 
  250000    4913.3042 
  260000    5404.8637 
  270000    4420.5497 
  280000     5357.551 
  290000    5009.1386 
  300000    4610.7359 
  310000    4752.3123 
  320000    4358.3474 
  330000    4961.3549 
  340000    4317.1413 
  350000    4154.0022 
  360000    5339.4602 
  370000     4460.511 
  380000    4997.8213 
  390000    4354.7411 
  400000    4797.0648 
  410000     4830.892 
  420000    5087.6876 
  430000    4723.8778 
  440000    4629.8021 
  450000    4976.1032 
  460000    4609.4548 
  470000    5019.8023 
  480000    4446.1004 
  490000    4862.2646 
  500000     5026.155 
  510000     4431.281 
  520000    4391.5373 
  530000    4214.6362 
  540000    4846.1485 
  550000    4807.4181 
  560000    4488.3661 
  570000    5360.1736 
  580000    5071.4535 
  590000     4463.498 
  600000    4555.1718 
  610000    5114.4069 
  620000    4699.2292 
  630000    5142.6987 
  640000    4479.5147 
  650000    4671.0304 
  660000    4564.5161 
  670000    4265.3884 
  680000    4814.3015 
  690000    4964.3923 
  700000    5121.2814 
  710000    5230.1069 
  720000    4297.1403 
  730000    4433.7983 
  740000    5176.7691 
  750000    4783.2431 
  760000    4792.0254 
  770000    4318.1008 
  780000    5279.5916 
  790000    4970.3985 
  800000    4163.3936 
  810000    5014.5128 
  820000    4423.9304 
  830000    4860.6778 
  840000    4502.2809 
  850000    4738.7186 
  860000    4585.4628 
  870000    4851.0205 
  880000    4747.5203 
  890000    4511.0158 
  900000    4841.8463 
  910000    5075.7912 
  920000    5070.3709 
  930000    4343.2652 
  940000    4743.7859 
  950000    4020.8883 
  960000    4671.5891 
  970000    4446.4685 
  980000     4448.264 
  990000    4507.5559 
 1000000    4873.6706 
Loop time of 36.1922 on 4 procs for 1000000 steps with 100 atoms

Pair  time (%) = 4.03491 (11.1486)
Bond  time (%) = 8.30473 (22.9462)
Neigh time (%) = 1.95553 (5.40319)
Comm  time (%) = 19.297 (53.3182)
Outpt time (%) = 0.00102639 (0.00283595)
Other time (%) = 2.59897 (7.18103)

Nlocal:    25 ave 39 max 5 min
Histogram: 1 0 0 0 0 1 0 0 1 1
Nghost:    45.75 ave 62 max 25 min
Histogram: 1 0 0 0 0 1 0 1 0 1
Neighs:    267.75 ave 425 max 44 min
Histogram: 1 0 0 0 1 0 0 0 1 1
FullNghs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0

Total # of neighbors = 1071
Ave neighs/atom = 10.71
Ave special neighs/atom = 5.88
Neighbor list builds = 100000
Dangerous builds = 100000
unfix 1
unfix 2
write_restart 	restart.${simname}.dreiding1
write_restart 	restart.PE_cl100.dreiding1


#####################################################
# Define Settings
compute eng all pe/atom 
compute eatoms all reduce sum c_eng 

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

dump 		1 all cfg 6 dump.comp_*.cfg id type xs ys zs c_csym c_peratom fx fy fz

reset_timestep 0 
fix 1 all nvt temp 500.0 500.0 100.0
thermo 20 
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
min_style cg
minimize 1e-25 1e-25 500000 1000000 
WARNING: Resetting reneighboring criteria during minimization
Memory usage per processor = 6.81045 Mbytes
Step PotEng Lx Ly Lz Press Pxx Pyy Pzz eatoms 
       0    1237.2398        158.5        158.5          100    16.181003    10.716166    11.032966    26.793875    1237.2398 
      20    93.175603        158.5        158.5          100    23.387346    19.859621     23.64755    26.654868    93.175603 
      40    63.587362        158.5        158.5          100    25.830375    21.672133    27.047328    28.771665    63.587362 
      60    54.741088        158.5        158.5          100    26.066912    21.364531      27.7035    29.132705    54.741088 
      80    48.917554        158.5        158.5          100     25.87727    20.908929    27.620402    29.102479    48.917554 
     100    45.293817        158.5        158.5          100    26.386696    21.757562    27.905252    29.497275    45.293817 
     120    41.459225        158.5        158.5          100    25.974718    20.775913    27.911315    29.236928    41.459225 
     140    37.876111        158.5        158.5          100    25.799883    20.794051    27.553625    29.051972    37.876111 
     160    36.701952        158.5        158.5          100    25.852641    20.951326    27.544375    29.062221    36.701952 
     180     35.07739        158.5        158.5          100    26.004914    21.057098    27.774282    29.183362     35.07739 
     200    33.877666        158.5        158.5          100    25.961315     21.11421    27.590146     29.17959    33.877666 
     220     31.12208        158.5        158.5          100    23.995326    19.308908     24.96515    27.711921     31.12208 
     240    28.809216        158.5        158.5          100    25.831303    21.277521    27.270668    28.945719    28.809216 
     260     26.02543        158.5        158.5          100    25.951859     21.58113    27.157947    29.116498     26.02543 
     280    23.070199        158.5        158.5          100    26.005868    21.781819    27.222861    29.012922    23.070199 
     300    20.013665        158.5        158.5          100    25.453848    20.749404    26.676606    28.935533    20.013665 
     320    18.738385        158.5        158.5          100    25.954677    21.603359    27.136183    29.124488    18.738385 
     340    16.749844        158.5        158.5          100    25.789221    21.322696     27.03715    29.007817    16.749844 
     360      15.8834        158.5        158.5          100    25.985101    21.659604    27.170167    29.125532      15.8834 
     380    14.987185        158.5        158.5          100    25.977784    21.397522    27.357437    29.178393    14.987185 
     400    14.517165        158.5        158.5          100     25.95897     21.55348    27.296561    29.026868    14.517165 
     420      13.7945        158.5        158.5          100    26.006995    21.535186    27.281645    29.204153      13.7945 
     440    12.512098        158.5        158.5          100    25.840513     21.28497    27.217437    29.019132    12.512098 
     460    12.290892        158.5        158.5          100    25.924928    21.357401    27.298119    29.119263    12.290892 
     480    10.775174        158.5        158.5          100    25.830121    21.201044    27.158886    29.130434    10.775174 
     500    10.218356        158.5        158.5          100    25.614842      21.0218    26.893606     28.92912    10.218356 
     520    9.8614258        158.5        158.5          100     25.93082    21.345743    27.320707    29.126011    9.8614258 
     540    9.1390877        158.5        158.5          100    25.893538     21.32844    27.283388    29.068787    9.1390877 
     560    8.5527184        158.5        158.5          100    25.954932    21.387239    27.355923    29.121632    8.5527184 
     580    8.1607915        158.5        158.5          100    25.865475     21.24017    27.270322    29.085934    8.1607915 
     600    7.8012324        158.5        158.5          100    25.950545    21.404499    27.375776    29.071362    7.8012324 
     620    7.6152706        158.5        158.5          100    25.937275    21.370367    27.342438    29.099019    7.6152706 
     640    7.3593621        158.5        158.5          100    25.940706    21.344493    27.415647    29.061979    7.3593621 
     656    7.2144567        158.5        158.5          100    25.915345    21.326291    27.308608    29.111136    7.2144567 
Loop time of 0.712466 on 4 procs for 656 steps with 100 atoms

Minimization stats:
  Stopping criterion = linesearch alpha is zero
  Energy initial, next-to-last, final = 
         1237.23980569      7.21445672663      7.21445672663
  Force two-norm initial, final = 1373.33 1.10527
  Force max component initial, final = 285.91 0.235955
  Final line search alpha, max atom move = 6.31527e-09 1.49012e-09
  Iterations, force evaluations = 656 3308

Pair  time (%) = 0.0280486 (3.93683)
Bond  time (%) = 0.0298776 (4.19355)
Neigh time (%) = 0.00120121 (0.168599)
Comm  time (%) = 0.05711 (8.01581)
Outpt time (%) = 0.359286 (50.4284)
Other time (%) = 0.236944 (33.2568)

Nlocal:    25 ave 41 max 4 min
Histogram: 1 0 0 0 1 0 0 0 1 1
Nghost:    47.25 ave 65 max 26 min
Histogram: 1 0 0 0 0 1 1 0 0 1
Neighs:    287.25 ave 507 max 30 min
Histogram: 1 0 0 1 0 0 0 1 0 1
FullNghs: 573 ave 992 max 47 min
Histogram: 1 0 0 0 1 0 0 1 0 1

Total # of neighbors = 2292
Ave neighs/atom = 22.92
Ave special neighs/atom = 5.88
Neighbor list builds = 51
Dangerous builds = 0


print "All done"
All done


LAMMPS dump files

Also, the dump files have been produced during the simulation. The can be used to produce a movie.


Movie

First download AtomEye.

  • Go to AtomEye website.
  • Click on Download, this will take you to the raw binary files.
  • Download the A.i686 version by right-clicking on the link and "Save Target As..." to one of your directories.
  • Rename the downloaded file as A by typing "cp A.i686 A". A will be the executable.
  • Using UNIX, run "chmod A 755" on the file to change this to an executable.
  • To test, you need to save a CFG file as well, such as cnt8x3.cfg. Try running using "./A cnt8x3.cfg".

Next download ImageJ.

  • Go to ImageJ website.
  • Download the appropriate version for Windows, LINUX, or MAC.

This assumes that you already have AtomEye and ImageJ downloaded and installed.

  • Visualize the dumpfile in AtomEye by typing the following command, "/A dump.tensile_0.cfg" (UNIX).
  • Use the AtomEye options to select how you want to visualize deformation. In this example, the centrosymmetry parameter was used to show only atoms in a non-centrosymmetric environment (see Fig. 2).
    • Use Alt+0 to activate centrosymmetric (csym) view.
    • Adjust threshold, or set of atoms to view, by using Shift+T. This will allow creation of a set for the current parameter (in this case, csym). Please note that you need to adjust both lower and higher thresholds unless the atoms from following images that exceeds maximum value for the first one will be not shown. You can make it 5 or 10.
    • Make atoms with values outside of the threshold invisible by using Ctrl+A.
  • Press 'y' within AtomEye to produce an animation script.
  • The folder "Jpg" now contains snapshots of all dumpfiles.
  • Open ImageJ
  • Drag the folder Jpg into ImageJ
    • Select "Convert to RGB" to keep the color from the AtomEye images.
    • Choose "yes" to load a stack.
  • Adjust the size as needed (Image/Adjust/Size)
  • Adjust frame rate as desired (Image/Stacks/Tools/Animation Options)
  • Save as Animated Gif file


Figure 4. The equilibration process followed by minimization.

Go Back

Acknowledgments

The authors would like to acknowledge funding for this work through the Department of Energy.

Personal tools
Namespaces

Variants
Actions
home
Materials
Material Models
Design
Resources
Projects
Education
Toolbox