hi every one
i’m computing the thermal conductivity of water-cu by Green-kubo formula.in order to articles the obtained answer is’nt reasonable and it shows weired number :4.5e-011 or 6.1e-012…
the answer should be around 6.5 and i tried changing neighboring .npt,nve but it didn’t work.i have used the codes in lammps manual.
thank you for helping me
regard
units real
variable T equal 298
variable V equal vol
variable dt equal 0.001
variable x equal 23.41
variable y equal 23.41
variable z equal 23.41
variable rho equal 0.6
variable t equal 20
variable rc equal 2.5
variable p equal 200 # correlation length
variable s equal 2 # 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 kCal2J equal 4186.0/6.02214e23
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal {kCal2J}*{kCal2J}/{fs2s}/{A2m} #thermal conductivity
variable convert equal {atm2Pa}*{atm2Pa}{fs2s}*{A2m}{A2m}*{A2m} #*${kCal2J}
#set up problem
dimension 3
echo screen
boundary p p p
newton on
atom_style full
bond_style harmonic #hybrid harmonic
angle_style harmonic #hybrid harmonic
kspace_style pppm 1.0e-5
read_data water.data
group hydrogen type 1
group water type 1 2
group cu type 3
group oxygen type 2
lattice fcc 3.615 #Cu lattice constant
region Cu sphere 0 0 0 7 units box
create_atoms 3 region Cu
set group oxygen charge -0.834
set group hydrogen charge .520
set group cu charge 0.000
pair_style hybrid lj/cut/coul/long 0.1521 10 eam lj/cut 5
pair_coeff 1 1 lj/cut/coul/long 0.0460 0.4000 #H-H epsilon sigm #
pair_coeff 1 2 lj/cut/coul/long 0.0836 1.7753 #O-H epsilon sigma
pair_coeff 1 3 lj/cut 0.6589 0.2117 #H-Cu epsilon sigma
pair_coeff 2 2 lj/cut/coul/long 0.1521 3.157 #O-O epsilon sigma # 0 0
pair_coeff 2 3 lj/cut 1.198 1.587 #O-Cu epsilon sigma
pair_coeff 3 3 eam cu.eam #Cu-Cu
for cu-cu bond sigma=.227 epsilon(Lj)=.583 ev # sigma=2.34 epsilon=9.4512 kcal/mol … cu eam cut off= 4.95 Ang
pair_modify mix arithmetic compute no
bond_coeff 1 450 0.9572 #O-H
angle_coeff 1 55 104.52 #H-O-H
++++++++++++++++setting+++++++++++++++++++++
neighbor 0.3 bin #bin #
neigh_modify delay 0 every 1 check yes
timestep ${dt}
thermo $d
velocity all create 298 458127641 rot yes dist gaussian
fix 12 water shake 1e-4 200 0 b 1 a 1
#fix 12 water npt temp 298 298 100.0 iso 0.0 0.0 1000.0
---------- Relaxation -----------------------------------------
unfix 12
min_style cg#fire#
min_modify dmax 0.01
minimize 1.0e-17 1.0e-17 0 90000
run 800
reset_timestep 0
#------------------------dump--------------------------------
dump 1 all xyz 10000 final.xyz
settings
Green-Kubo viscosity calculation
Define distinct components of symmetric traceless stress tensor
#*******************************#thermal conductivity
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom 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/$Vs*{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}
thermo_style custom step temp press v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 #v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
run 6000
variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print “average conductivity: $k[W/mK] @ T K, {ndens} /A^3”