Hello everyone,
I have simulated single layered graphene sheet with CVFF potential. The stress-strain curve for pristine graphene started from 0, but when I started adding -OH groups on the surface, the curve no longer started from 0. For, 100% functionalization, the curve started from 90 GPa and fractured at 140 GPa. Is it because, the O atoms are pulling the C atoms up, causing residual stresses, or am I calculating stresses incorrectly? But I am getting good results for pristine graphene. I am using per/atom stress command to calculate the stresses for C atoms. Here’s part of my input script:
I will greatly appreciate any suggestion. Thank you.
INPUT:
graphene (1-layer) (functionalized with -OH group)
strain in x-direction (zigzag)
Initialization
echo screen
log debug_grap1ohfull.log
units real
atom_style full
boundary p p p
Bonded Interaction
bond_style morse
angle_style harmonic
dihedral_style harmonic
improper_style cvff
Non-Bonded Interaction
pair_style lj/cut/coul/long 10.0 10.0 # cutoff1 = 2.5*sigma
read_data data.grap1ohfull
kspace_style ewald 1e-6
neighbor 5.0 bin
#neigh_modify delay 0 every 1 one 10000 check yes
neigh_modify delay 0 every 1 check yes
pair_modify mix arithmetic tail no
#------------ define groups ------------------xxx
group graphene type 1
group oxygen type 2
group ho type 3
compute tmpall all temp
compute mobile_temp all temp
variable temp equal “300” # Kelvin
velocity all create ${temp} 878234 temp mobile_temp dist gaussian mom yes rot yes
#velocity graphene set 1e-6 0.0 0.0 temp mobile_temp units box dist gaussian mom yes rot no
#------------------energy minimization-----------------------------xxx
#------------------NVE + Temp/rescale------------------------------xxx
minimize 1.0e-4 1.0e-6 100 1000
min_style cg
min_modify dmax 0.1
fix 2 all nve
fix 3 all temp/rescale 1 {temp} {temp} 1.0 0.5
fix 4 all langevin {temp} {temp} 100.0 2341456
timestep 1
thermo 100
#thermo_style custom step temp pe etotal press vol density evdwl ecoul elong epair ebond eangle edihed eimp emol # epair = evdwl+ecoul+elong;
thermo_style custom step temp pe etotal press vol density
thermo_modify lost error flush yes # emol = ebond+eangle+edihed+eimp;
run 100 # 500 ps
dump 1 all xyz 1 grap.xyz
#dump_modify 1 element C H O
write_data data.grap1ohfull_nve pair ii # Data after NVE+Langevin
unfix 2
unfix 3
unfix 4
#---------------------------NPT------------------------------------xxx
reset_timestep 0
fix 7 all npt temp {temp} {temp} 100.0 iso 0.0 0.0 1000.0 pchain 1 # Pressure unit in (atm)
#fix 7 all npt temp {temp} {temp} 100.0 x 0.0 0.0 1000.0 y 1.0 1.0 1000.0 z 1.0 1.0 1000.0 pchain 1 # Pressure unit in (atm)
timestep 1
thermo 100
#thermo_style custom step temp pe etotal press vol density evdwl ecoul elong epair ebond eangle edihed eimp emol # epair = evdwl+ecoul+elong;
thermo_style custom step temp pe etotal press vol density
thermo_modify lost error flush yes # emol = ebond+eangle+edihed+eimp;
run 300 # 1000 ps
dump 2 all xyz 1 grap.xyz
write_data data.grap1ohfull_npt pair ii # Data after NVE+temp/rescaling+Langevin+NPT
unfix 7
#---------------------------MD simulation------------------------------------xxx
#---------------------------Stress Calculation-------------------------------xxx
reset_timestep 0
compute peratom graphene stress/atom mobile_temp virial
compute p graphene reduce sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]
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]
#variable stress equal (c_p[1]+c_p[2]+c_p[3])/(3*2477.6588) # per-atom stress is negative of total pressure p
Calculation of stresses and strains # per-atom stress is sum of diagonal components in stress tensor
variable multiplication_factor equal 101.3251e-6 # stress converstion from atm to GPa
variable f equal {multiplication_factor}
variable sigmax equal (v_f)*(c_p[1]/{gvol})
variable sigmay equal (v_f)(c_p[2]/${gvol})
Applying axial and shear 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
fix print all print 100 "{s1} {p1} {p2} {p3} {p4}" file wgw.defo.txt screen no
Run the simulation
timestep 1
thermo 100
#thermo_style custom step temp v_strainx c_Utot
thermo_style custom step temp v_strainx v_sigmax v_sigmay v_mu v_lambda v_nu v_E lx v_C11 v_C12
run 300000 # 300 ps
write_data data.grap1ohfull_const pair ii
unfix 8
unfix 9
print “SIMULATION DONE!”