Minimization of a Bi-layer Graphene structure

Hi All,

I am trying to simulate a bi-layer graphene structure (lammps version: 25 Feb 2016). I am using a combination of tersoff (intralayer) and LJ (interlayer) potential. I am simulating it in two scenario:
1) Initial configuration is perfect AA stacking (both layer are exactly on top of each other, please visit the link for image)
2) Initial configuration is SP stacking (see the image attached, please visit the link for image)

The minimization leads to the different structure depending the initial configuration. In case (1), it gives a structure in AB stacking with both the graphene layers are flat. In case (2), however, the the minimized structure is not a perfect AB stacking and have corrugations in both layer (i.e. out of plane transverse displacement). The case (1) makes sense, but the case (2) seems completely non-sense (atleast to me); because corrugation should increase the bending energy, hence a both layer should remain flat in AB stacking in both the cases. To see where the simulation are going, I checked the thermo output for the case (2), some steps of the thermo output are shown below:

Step TotEng LJ_E Est1 Est2
610 -61497.056 -481.48763 -30507.784 -30507.784
620 -61497.056 -481.48763 -30507.784 -30507.784
630 -61497.056 -481.48763 -30507.784 -30507.784
640 -61497.056 -481.48763 -30507.784 -30507.784
650 -61497.056 -481.48763 -30507.784 -30507.784
660 -61497.121 -480.7168 -30508.224 -30508.181
670 -61497.682 -478.71814 -30509.511 -30509.453
680 -61497.958 -477.9428 -30509.996 -30510.02
690 -61498.168 -477.18335 -30510.466 -30510.518
700 -61498.255 -477.01739 -30510.616 -30510.621
710 -61498.301 -476.94745 -30510.681 -30510.673
720 -61498.329 -476.9662 -30510.682 -30510.681

Here, TotEng is the total energy of the system, LJ_E is interlayer Lj interaction energy, Est1 and Est2 are the individual intralayer energy for bottom layer and top layer respectively. We can see that , the simulation are going correct till step 650, upto which lammps is trying to minimize the LJ interaction energy (the system achieved the AB stacking while the both layer are flat). The step 650, is the onset of the corrugation in the Graphene layers. After step 650, the LJ energy start increasing, and Est1 and Est2 start decreasing, hence decreasing the total energy of the system. But this does not make sense, as bending the layer should cost more energy, moreover why does the different initial configuration leading to different final configuration during the minimization. Any suggestion or comment where I am making a mistake is highly appreciated.

Thank you for time.

Best,
Vinit

Here is link for the images (Because of size limitation)

https://drive.google.com/open?id=1wBjR15ucGI_HlAaw9Hp_EGK6ivNJHFx5

Here is my input file:
# --------------------MATERIAL PROPERTIES --------------------
#10.1103/PhysRevB.83.235428
variable epC1C2 equal 0.0046
variable sigC1C2 equal 3.276

variable massC equal 12.0107

# ------------------- INITIALIZATION ------------------------
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array

newton on
timer full

# GR-2 setup
read_data ${inread}
replicate 2 24 1 #(3.1 x 3.1 x 3.1)
change_box all boundary p p p

# Group
group GR type 1 2
group GR1 type 1 # Bottom layer
group GR2 type 2 # top layer

# Set mass
mass 1 \{massC\} mass 2 {massC}

# Force field
pair_style hybrid tersoff tersoff lj/sf 10.0
pair_coeff * * tersoff 1 C.tersoff C NULL # C1
pair_coeff * * tersoff 2 C.tersoff NULL C # C2 tersoff-LJ
pair_coeff 1 2 lj/sf \{epC1C2\} {sigC1C2} # C1-C2

# Neighbor updating criteria
neighbor 5.0 bin
neigh_modify every 10

#-------------------- Minimization ------------------------------

compute 2 GR1 pe/atom pair
compute 3 GR1 reduce sum c_2
compute 4 GR2 pe/atom pair
compute 5 GR2 reduce sum c_4

compute 6 GR2 group/group GR1
variable fx equal c_6[1]
variable fy equal c_6[2] ############## computing LJ energy and forces on top layer
variable fz equal c_6[3]

variable E equal "2*c_6" # total LJ energy

variable Est1 equal "c_3-c_6"
variable Est2 equal "c_5-c_6" # individual Energy in each layer

displace_atoms GR2 move 0.0 0.709 3 units box # y= 2.836 --> AA stacking, y= 0.709 --> SP stacking

thermo 10
thermo_style custom step etotal v_E v_fx v_fy v_fz v_Est1 v_Est2
dump d1 all atom 10 confinit0.lammpstrj

min_style cg
minimize 0.0 1e-4 50000 100000
undump d1

print "-------------------------Mini done--------"

Hi All,

Any thoughts or comments on this:

I am trying to simulate a bi-layer graphene structure (lammps version: 25 Feb 2016). I am using a combination of tersoff (intralayer) and LJ (interlayer) potential. I am simulating it in two scenario:
1) Initial configuration is perfect AA stacking (both layer are exactly on top of each other, please visit the link for image)
2) Initial configuration is SP stacking (see the image attached, please visit the link for image)

The minimization leads to the different structure depending the initial configuration. In case (1), it gives a structure in AB stacking with both the graphene layers are flat. In case (2), however, the the minimized structure is not a perfect AB stacking and have corrugations in both layer (i.e. out of plane transverse displacement). The case (1) makes sense, but the case (2) seems completely non-sense (atleast to me); because corrugation should increase the bending energy, hence a both layer should remain flat in AB stacking in both the cases. To see where the simulation are going, I checked the thermo output for the case (2), some steps of the thermo output are shown below:

Step TotEng LJ_E Est1 Est2
610 -61497.056 -481.48763 -30507.784 -30507.784
620 -61497.056 -481.48763 -30507.784 -30507.784
630 -61497.056 -481.48763 -30507.784 -30507.784
640 -61497.056 -481.48763 -30507.784 -30507.784
650 -61497.056 -481.48763 -30507.784 -30507.784
660 -61497.121 -480.7168 -30508.224 -30508.181
670 -61497.682 -478.71814 -30509.511 -30509.453
680 -61497.958 -477.9428 -30509.996 -30510.02
690 -61498.168 -477.18335 -30510.466 -30510.518
700 -61498.255 -477.01739 -30510.616 -30510.621
710 -61498.301 -476.94745 -30510.681 -30510.673
720 -61498.329 -476.9662 -30510.682 -30510.681

Here, TotEng is the total energy of the system, LJ_E is interlayer Lj interaction energy, Est1 and Est2 are the individual intralayer energy for bottom layer and top layer respectively. We can see that , the simulation are going correct till step 650, upto which lammps is trying to minimize the LJ interaction energy (the system achieved the AB stacking while the both layer are flat). The step 650, is the onset of the corrugation in the Graphene layers. After step 650, the LJ energy start increasing, and Est1 and Est2 start decreasing, hence decreasing the total energy of the system. But this does not make sense, as bending the layer should cost more energy, moreover why does the different initial configuration leading to different final configuration during the minimization. Any suggestion or comment where I am making a mistake is highly appreciated.

Thank you for time.

Best,
Vinit

Here is link for the images (Because of size limitation)

https://drive.google.com/open?id=1wBjR15ucGI_HlAaw9Hp_EGK6ivNJHFx5

Here is my input file:
# --------------------MATERIAL PROPERTIES --------------------
#10.1103/PhysRevB.83.235428
variable epC1C2 equal 0.0046
variable sigC1C2 equal 3.276

variable massC equal 12.0107

# ------------------- INITIALIZATION ------------------------
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array

newton on
timer full

# GR-2 setup
read_data ${inread}
replicate 2 24 1 #(3.1 x 3.1 x 3.1)
change_box all boundary p p p

# Group
group GR type 1 2
group GR1 type 1 # Bottom layer
group GR2 type 2 # top layer

# Set mass
mass 1 \{massC\} mass 2 {massC}

# Force field
pair_style hybrid tersoff tersoff lj/sf 10.0
pair_coeff * * tersoff 1 C.tersoff C NULL # C1
pair_coeff * * tersoff 2 C.tersoff NULL C # C2 tersoff-LJ
pair_coeff 1 2 lj/sf \{epC1C2\} {sigC1C2} # C1-C2

# Neighbor updating criteria
neighbor 5.0 bin
neigh_modify every 10

#-------------------- Minimization ------------------------------

compute 2 GR1 pe/atom pair
compute 3 GR1 reduce sum c_2
compute 4 GR2 pe/atom pair
compute 5 GR2 reduce sum c_4

compute 6 GR2 group/group GR1
variable fx equal c_6[1]
variable fy equal c_6[2] ############## computing LJ energy and forces on top layer
variable fz equal c_6[3]

variable E equal "2*c_6" # total LJ energy

variable Est1 equal "c_3-c_6"
variable Est2 equal "c_5-c_6" # individual Energy in each layer

displace_atoms GR2 move 0.0 0.709 3 units box # y= 2.836 --> AA stacking, y= 0.709 --> SP stacking

thermo 10
thermo_style custom step etotal v_E v_fx v_fy v_fz v_Est1 v_Est2
dump d1 all atom 10 confinit0.lammpstrj

min_style cg
minimize 0.0 1e-4 50000 100000
undump d1

print "-------------------------Mini done--------"