Hello everyone,

I did MD simulation of a bilayer graphene in LAMMPS. The stress-strain behavior of BLG didn’t seem to deviate from the single sheet. The fracture stress of single sheet graphene was 167 GPa. And the bilayer fracture stress are between 162-168 GPa. I am using CVFF potential for the simulation. When I created my data file, I created two layers as two types of C atoms and set up an interaction between type 1 and type 2 C atoms using lj/cut. The parameters for lj/cut are same as I used between 1-1 interaction or 2-2 interaction. If I change the lj/cut parameters, then the bilayer strength changes a lot. But that seems nonphysical to me. Can someone help me with what I am doing wrong? Thank you very much. Here’s a part of my code:

# Non-Bonded Interaction

pair_style hybrid lj/cut/coul/long 10.0 lj/cut 10.0 # cutoff1 = 2.5*sigma

read_data data.grapbi2

pair_coeff 1 1 lj/cut/coul/long 0.1479999981 3.6170487995

pair_coeff 2 2 lj/cut/coul/long 0.1479999981 3.6170487995

pair_coeff 1 2 lj/cut 0.1479999981 3.6170487995

compute peratom all stress/atom tmpall virial

# Storing pxx, pyy, pzz and pxy as variables

variable p1 equal c_p[1]

variable p2 equal c_p[2]

variable p3 equal c_p[3]

variable p4 equal c_p[4]

# Calculation of stresses and strains # per-atom stress is sum of diagonal components in stress tensor

variable prefactor equal 101.325*1e-6 # stress converstion from atm to GPa
variable f equal {prefactor}
variable sigmax equal (v_f)*(c_p[1]/{gvol})
variable sigmay equal (v_f)*(c_p[2]/${gvol})

compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]

# Applying axial strain by fix deform

fix 8 all nvt temp {temp} {temp} 5.0 drag 1

fix 9 all deform 1 x erate ${erate1} units box remap x flip no

Regards,

Baig