units real
dimension 3
boundary p p p
neighbor 3.0 bin
neigh_modify delay 0 every 1 check yes page 1000000 one 20000
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
pair_style lj/cut/coul/long 12.25 18.0
pair_modify mix geometric
special_bonds lj/coul 0.0 0.0 0.5
kspace_style pppm 0.00001

read_data "system.data"

pair_coeff 1 1 0.066 3.5
pair_coeff 2 2 0.066 3.5
pair_coeff 3 3 0.03 2.5
bond_coeff 1 268.0 1.529
bond_coeff 2 340.0 1.09
angle_coeff 1 58.35 112.7
angle_coeff 2 33.0 107.8
angle_coeff 3 37.5 110.7
dihedral_coeff 1 1.3 -0.05 0.2 0.0
dihedral_coeff 2 0.0 0.0 0.3 0.0
dihedral_coeff 3 0.0 0.0 0.3 0.0
set type 1 charge -0.18
set type 2 charge -0.12
set type 3 charge 0.06

#-----------------------------variable---------------------

variable    V  equal vol
variable    kB equal 1.3806504e-23
variable    A2m equal 1.0e-10
variable    fs2s equal 1.0e-15
variable    kcal2j equal 4186.0/6.02214e23
variable    convert equal ${kcal2j}*${kcal2j}/${fs2s}/${A2m}

#------------------------------------minimize-----------------

variable dt equal 1.0
variable T equal 600
timestep ${dt}

dump mydump1 all custom 5000 traj_eq1_min.lammpstrj id mol type x y z ix iy iz
thermo_style custom step temp ke pe etotal press vol density 
thermo 1000


min_style sd
minimize 1e-4 1e-6 10000 10000

undump mydump1
write_data system_after_eq1_min.data

#----------------------nvt------------------------------------

dump mydump2 all custom 5000 traj_eq2_nvt.lammpstrj id mol type x y z ix iy iz

velocity all create $T 456783 mom yes rot yes dist gaussian

fix fxnvt all nvt temp $T $T 100

run 500000

undump mydump2

write_data system_after_eq2_nvt.data

unfix fxnvt

#----------------------------npt--------------------------------

dump mydump3 all custom 5000 traj_eq2_npt.lammpstrj id mol type x y z ix iy iz

fix fxnpt all npt temp $T $T 100 iso 29.607 29.607 1000

run 500000

undump mydump3

write_data system_after_eq3_npt.data

unfix fxnpt

#----------------------------------nve-----------------------

fix fxnve all nve

dump mydump4 all custom 5000 traj_nve.lammpstrj id mol type x y z ix iy iz

run 500000

undump mydump4

write_data system_after_eq4_nve.data

#----------------------------------------product------------------


reset_timestep 0

variable s equal 10
variable p equal 1000
variable d equal $p*$s
variable w equal 250
variable r equal $w*$d

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]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol

fix JJ all ave/correlate $s $p $d c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running

variable scale equal ${convert}/${kB}/$T/$T/$V*$s*${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}

dump tc_dump all custom 10000 tc_dump.lammpstrj id type x y z

thermo ${d}
thermo_style custom step v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 press vol temp

run ${r}

variable k equal (v_k11+v_k22+v_k33)/3.0

print "average conductivity: $k [W/mK] @ $T K."