Dear users,
I’m stuck at trying to equilibrate 3 graphene sheets (for a total of 11808 atoms).
------------------------ INITIALIZATION ----------------------------
units metal
dimension 3
boundary p p p
newton on
--------------------------- BOND TYPE ------------------------------
atom_style molecular
bond_style none
-------------------------- FORCE FIELD -----------------------------
pair_style airebo 3.0 1 1
pair_coeff * * CH.airebo C
---------------------------- SETTINGS ------------------------------
group 1atoms id 1:3936
group 2atoms id 3937:7872
group 3atoms id 7873:11808
-------------------------- MINIMIZATION ----------------------------
variable n loop 400
label minloop
fix 1 all box/relax iso 0.0 vmax 1.0e-10 nreset 100
fix 2 setforce 1atoms 0.0 0.0 0.0
fix 3 setforce 2atoms 0.0 0.0 0.0
fix 4 setforce 3atoms 0.0 0.0 0.0
min_style cg
min_modify dmax 2.0 line quadratic
minimize 1.0e-10 1.0e-10 10000 100000
unfix 1
min_style cg
min_modify dmax 2.0 line quadratic
minimize 1.0e-10 1.0e-10 10000 100000
next n
jump SELF minloop
This led the system to shrink and bonds to form between atoms of adjacent layers; it was suggested in previous posts that this may be due to wrapped atom coordinates but I didn’t quite understand what was meant by it. Moreover, the sheets were quite deformed at the edges while the middle section was unaffected.
A similar outcome was obtained when aveforce commands were exchanged for setforce ones, although a lower potential energy was obtained with aveforce (-22000 compared to -10500). When neither aveforce nor setforce commands were issued, the system still appeared the same as described above.
Increasing vmax of the fix box/relax to more than 1.0e-5 produced high negative pressures and it was at 1.0e-10 that the system closely approached the 0 bar target set by box/relax.
Then, I read that rotational motion should also be fixed to prevent atoms from flying off the layers so I added these commands in between the minimizations:
fix 5 allatoms rigid group 3 1atoms 2atoms 3atoms force 13 off off off torque 13 off off off
fix 6 allatoms momentum 1 linear 1 1 1 angular rescale
timestep 0.0005
run 50
unfix 5
unfix 6
The result from loop n=1 read the following:
Step Temp Press Volume TotEng Lx Ly Lz Pxx Pyy Pzz Fmax Fnorm
0 0 -18.166496 48000000 -296.24299 400 400 300 -19.508143 -19.544627 -15.446719 0.01175742 1.0431687
10000 0 -35.502223 34992000 -507.43287 360 360 270 -35.055259 -35.037493 -36.413916 0.022450616 1.991915
10000 0 1.735124e-15 34992000 -507.43287 360 360 270 1.6136653e-13 -2.602686e-14 -1.301343e-13 0.023323621 2.1226287
10200 0 -3.9907853e-14 34992000 -507.43287 360 360 270 3.6437605e-14 -4.6848349e-14 -1.0931281e-13 0.023323621 2.1226287
10200 0 -35.502223 34992000 -507.43287 360 360 270 -35.055259 -35.037493 -36.413916 0.022450616 1.991915
10433 0 -27.615845 34992000 -756.32193 360 360 270 -41.406116 -41.44138 -3.7791408e-05 3.4501836e-08 3.061151e-06
Stopping criterion = linesearch alpha is zero
Results from other loops were similar.
Any help would be greatly appreciated.
Michael