Polymer simulation

In this blog, I will be showing how to calculate the Radius of gyration and Radial distribution function(rdf) for a polymer chain.

 Input script

# 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)

####################################################

# dump for Equilibration

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

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

# undump 1 as the dump process is going to continue in next steps

undump 1

#####################################################

# Define Settings

compute eng all pe/atom 

compute eatoms all reduce sum c_eng 

#####################################################

# Minimization

dump 1 all cfg 6 dump.comp_*.cfg mass 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"

#as we are further going to perform simulation undump 1 and unfix 1 is being done

undump 1

unfix 1

###################################################################################

# for radius of gyration 

###################################################################################

#to store rg values in rg.in file

log rg.txt

#time step is being reset so that rg is calculated from the begining

reset_timestep 0

#to calculate rg for equlibration process the following codes are used

compute gyr all gyration

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 c_gyr[1] c_gyr[2] c_gyr[3] c_gyr[4] c_gyr[5] c_gyr[6]

thermo          10000

timestep 1

run 1000000

unfix 1

unfix 2

write_restart restart.${simname}.dreiding1

#to calculate rg for energy minimisation process the following codes are used

reset_timestep 0 

#compute gyrm all gyration

fix 1 all nvt temp 500.0 500.0 100.0

thermo 20 

thermo_style custom step c_gyr[1] c_gyr[2] c_gyr[3] c_gyr[4] c_gyr[5] c_gyr[6]

min_style cg

minimize 1e-25 1e-25 500000 1000000 

unfix 1

###################################################################################

#for rdf

# Provide an initial maxwellian distribution of velocity corresponding to temperature 5000K

velocity all create 5000 198728 dist gaussian

reset_timestep 0

#Perform an NVE integration with this initial position and velocity distribution

timestep 1 #<Time step in ps. So this is 0.001ps or 1 femto second>

compute rdfs all rdf 100 1 1 cutoff 10.0

fix 1 all ave/time 1 10000 10000 c_rdfs[1] c_rdfs[2] c_rdfs[3] file rdf.txt mode vector

run 1000000

unfix 1

#End

Visit LAMMPS Manual

Original script and structure file

Simulation

Colouring of atoms done based on centrosymmetry

Explanation

I have edited the original script and added codes for radius of gyration and rdf. Radius of gyration values can be found in rg.txt file and rdf values can be found in rdf.txt file after running the script.

If you have any doubts or sugessions, mention them in comment box. Follow the blog for further updates.