############################################################################################ ## ## ## sample LAMMPS script for EMD (Green-Kubo method) [not fully developed but works] ## written by Afnan Mostafa ## ## ############################################################################################ dimension 3 units metal processors * * 2 boundary p p p neighbor 0.3 bin neigh_modify every 1 delay 0 check yes #define variables variable T equal 300.0 variable V equal lx*ly*10.05 # 3.35 (1st layer) + gap (0) + 3.35 (2nd layer) variable dt equal 0.0001 variable Nevery equal 3500 variable Nrepeat equal 100 variable Nfreq equal 350000 #convert from LAMMPS metal units to SI units variable kB equal 1.3806504e-23 variable eV2J equal 1.602177e-19 variable A2m equal 1.0e-10 variable ps2s equal 1.0e-12 variable convert equal ${eV2J}*${eV2J}/${ps2s}/${A2m} #define lattice atom_style atomic read_data graphene_3_layers_ab_10x10.data #define interatomic potentials pair_style airebo 3 1 0 pair_coeff * * CH.airebo C C C #group hydro type 4 #group carbon type 1 2 3 compute 1 all stress/atom NULL compute 2 all displace/atom #compute 3 carbon coord/atom cutoff 1.73 group all #compute 4 hydro coord/atom cutoff 1.3 group carbon #variable coords equal "c_3+c_4" dump mini all custom 5000 minimize.xyz id type x y z c_1[1] c_1[2] c_1[3] c_1[4] c_1[5] c_1[6] dump_modify mini sort id #define thermo variable settings timestep ${dt} #### BOND FORMATION #### fix noCOM all recenter INIT INIT INIT units box thermo_style custom step time temp press pxx pyy pzz pe ke etotal vol lx ly lz f_noCOM thermo 1000 thermo_modify lost error # ======= minimization ======= min_style cg minimize 1e-8 1e-8 400000 400000 #initialize velocities velocity all create $T 23121 mom yes rot yes dist gaussian units box # ======= thermal equilibration ======= fix NVE all nve fix ts all temp/rescale 1 $T $T $(100.0*dt) 1.0 run 120000 unfix ts unfix NVE fix 1 all nvt temp $T $T $(100.0*dt) drag 2.0 run 120000 unfix 1 fix 1 all nve run 50000 unfix 1 unfix noCOM undump mini dump emd all xyz 100000 emd.xyz dump_modify emd sort id fix 1 all nve #thermal conductivity calculation reset_timestep 0 thermo ${Nfreq} compute myKE all ke/atom compute myPE all pe/atom compute myStress all stress/atom NULL virial compute flux all heat/flux myKE myPE myStress variable Jx equal c_flux[1]/$V variable Jy equal c_flux[2]/$V variable Jz equal c_flux[3]/$V variable t equal step*${dt}/1000 fix JJ all ave/correlate ${Nevery} ${Nrepeat} ${Nfreq} & c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.data ave running variable scale equal ${convert}/${kB}/$T/$T/$V*${Nevery}*${dt} variable k11 equal trap(f_JJ[3])*${scale} variable k22 equal trap(f_JJ[4])*${scale} variable k33 equal trap(f_JJ[5])*${scale} variable k equal (v_k11+v_k22)/2.0 thermo_style custom step temp press v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 v_k vol pe ke fnorm lx ly fix print_tc all print ${Nfreq} "$t ${k11} ${k22} ${k33} $k" file "./tc_on_the_fly.txt" screen no title "#time_in_ns kxx kyy kzz k" restart 8750000 tmp.restart run 87500000 print "average thermal conductivity: $k [W/mK] @ $T K" write_data post_emd.data