Dear all LAMMPS users,

I am new user of lammps and I am trying to find the effect of tensile strain on graphene thermal conductivity. The program runs well without applying the tensile strain.

When I applied tensile strain, I am getting error as “Divide by 0 in variable formula” I find that in the denominator the value is getting zero. In the The input script is given below. In finding the Variable Z the denominator f_5 is getting zero after 50000 run i.e after 2nd equilibration run. The v_tdiff(without the application of strain) is attached here.

I want to know how the f_5 value is getting zero in case of tensile strain which is related to v_tdiff and it calculated as f_2[96][3]-f_2[1][3]. As I am new user of lammps

f_2[96][3]-f_2[1][3] how this getting calculated which is getting stored in profile.mp which is attached here(without application of strain).

Please refer the input script below.

variable T equal 300

variable V equal vol

variable dt equal 0.001

#variable p equal 1000 # correlation length

#variable s equal 5 # sample interval

#variable d equal $p*$s # dump interval

# convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann

variable eV2J equal 1.602e-19

variable A2m equal 1.0e-10

variable ps2s equal 1.0e-12

variable convert equal {eV2J}*{eV2J}/{ps2s}/{A2m}

variable alat equal 2.46

# setup problem

units metal

boundary p p p

atom_style atomic

read_data Read_data_1st

mass 1 12.0

pair_style tersoff

pair_coeff * * /home/kpark39/programs/lammps-1Feb14/potentials/SiC_1994.tersoff C

neighbor 2.0 bin

neigh_modify delay 5

timestep ${dt}

variable tmp equal lx

variable lo equal ${tmp}

variable strain equal (lx-v_lo)/v_lo

variable p1 equal “-pxx/10000”

variable p2 equal “-pyy/10000”

# 2nd equilibration run

compute ke all ke/atom

# units of Kb = J/K, units of ke = eV or J, so temp becomes K

variable temp atom c_ke*{eV2J}/(1.5*{kB})

fix 1 all nve

compute cc1 all chunk/atom bin/1d x lower 5

fix 2 all ave/chunk 10 100 1000 cc1 v_temp file profile.mp #units box

fix 3 all thermal/conductivity 10 x 190 swap 2

fix 4 all temp/berendsen 300 300 100

variable tdiff equal f_2[96][3]-f_2[1][3]

thermo_style custom step temp epair etotal f_3 v_tdiff

thermo 1000

dump 2 all cfg 1000 dump.config.*.cfg mass type xs ys zs vx vy vz x y z v_temp

run 50000

unfix 4

# thermal conductivity calculation

fix 5 all ave/time 1 1 1000 v_tdiff ave running file ave.dat

variable Z equal f_3*{eV2J}*1000/2/(110*3.4)/{A2m}/{dt}/100000/{ps2s}/f_5

fix 7 all deform 1 x erate 0.0001 units box

fix def all print 1000 “{strain} {p2} ${p1}” file Graphene_tensile_data_300K.txt

thermo_style custom step temp epair etotal f_3 v_tdiff f_5 v_Z

thermo 1000

dump 7 all custom 5000 Graphene_thermal_deform*.lammpstrj id type x y z

run 100000

ave.dat (1.57 KB)

profile.mp (833 KB)