change_box command with reaxff potential

Dear All,

I’m trying to do a static tensile simulation by Reax force field with change_box command. However, the stress-strain curves that I got were abnormal. Because the tensile stress could be roughly 100 times larger than tensile stress I got by using fix deform command. Meanwhile, my minimizer can’t convergence to the value that I desire. Could you give me some suggestions, and my input script is as follows.

Yuan Xin


units real
dimension 3
boundary p p p
atom_style charge
neighbor 2.5 bin
neigh_modify every 2 delay 10 check yes

read_restart poly.restart.720000

restart 10000 poly.restart

ReaxFF potentials

pair_style reax 10.0 0 1 1.0e-6
pair_coeff * * ffield.reax 1 2 3

group pmma type 1 2 3

group up_1 id >= 909
group down_1 id <= 1135
group chain_1 intersect up_1 down_1

group up_2 id >= 2044
group down_2 id <= 2270
group chain_2 intersect up_2 down_2

group up_3 id >= 3179
group down_3 id <= 3405
group chain_3 intersect up_3 down_3

group up_4 id >= 4314
group down_4 id <= 4540
group chain_4 intersect up_4 down_4

group up_5 id >= 5449
group down_5 id <= 5675
group chain_5 intersect up_5 down_5

group 3chains union chain_1 chain_2 chain_3
group 5chains union chain_1 chain_2 chain_3 chain_4 chain_5

compute reax all pair reax
variable eb equal c_reax[1]
variable ea equal c_reax[2]
variable elp equal c_reax[3]
variable emol equal c_reax[4]
variable ev equal c_reax[5]
variable epen equal c_reax[6]
variable ecoa equal c_reax[7]
variable ehb equal c_reax[8]
variable et equal c_reax[9]
variable eco equal c_reax[10]
variable ew equal c_reax[11]
variable ep equal c_reax[12]
variable efi equal c_reax[13]
variable eqeq equal c_reax[14]

store final cell length for strain calculations

variable px equal “-pxx1.013/10"
variable py equal "-pyy
variable pz equal “-pzz1.013/10"
variable time equal "step

compute s all stress/atom
compute strea all reduce sum c_s[1]
variable stressa equal c_strea/vol
compute streb all reduce sum c_s[2]
variable stressb equal c_streb/vol
compute strec all reduce sum c_s[3]
variable stressc equal c_strec/vol


timestep 0
thermo 100
thermo_style custom step v_strain temp etotal eangle ke pe vol lx ly lz fmax fnorm press v_px v_py v_pz
dump 1 all atom 500 dump_tensile.lammpstrj
dump_modify 1 scale no
dump 2 all xyz 500
dump_modify 2 element C H O
dump 3 chain_1 xyz 500
dump_modify 3 element C H O
dump 4 3chains xyz 500
dump_modify 4 element C H O
dump 5 5chains xyz 500
dump_modify 5 element C H O

min_style sd
minimize 0.0 1.0e-3 100000 100000

for loop

variable i loop 1000
label forloop

fix 1 all box/relax x 0 nreset 10
fix 2 all box/relax y 0 nreset 10

change_box all z delta 0.0 0.04 remap units box
min_style cg
minimize 0.0 1.0e-3 100000 100000
next i
jump in.txt forloop

From looking at your script, you seem to be trying to do the following:

  1. Read in a PMMA structure
  2. do for 1000 iterations
  3. Deform: Increase the length in z direction by 0.04 angstroms
  4. Relaxation: find local energy minimum and find Pxx=Pyy=0

This would be fine and dandy if energy relaxation and stress relaxation were completely robust operations in LAMMPS. They are not. They require a lot of experimentation with the settings and carefully monitoring. Moreover, your PMMA sample is probably an amorphous glass. That creates a lot of issues associated with the fact that the energy landscape is very rough with many local minima with large differences in energy.

I suggest first testing if you can successfully relax the undeformed sample. Once you can do that, you can experiment with very small strains, building up to large strains.

Finally comments:

  1. You should combine the two fix box/relax commands in to one.
  2. PMMA at 300K may behave very differently than PMMA at zero temperature. Consult the literature on MD simulation of polymer glasses in tension at finite temperature.