hi every one

i want to simulate cu-water.and compute thermal conductivity.i have compared results with the others.the results are good but there is a problem.in articles result is about 7/5… but my answer is 7.5*e-014.thus i changed variable dt to 20*e+14(it was 02) and the answer became 7.5*e-005…is changing the variable dt resonable?? does anybody know what should i do to give a resonable answer?

thanks

best wishes

variable T equal 298

variable V equal vol

variable dt equal 200000000000 # .00000000002

variable x equal 23.41

variable y 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

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 3 units box

create_atoms 3 region Cu

#set group oxygen charge -1.040 #???

#set group hydrogen charge .520 #???

#set group cu charge 0.000

pair_style hybrid lj/cut/coul/long 0.1521 3.157 eam lj/cut 5 # .583 #2.8 # 3.157 # 7.5

pair_coeff 1 1 lj/cut/coul/long 0.0460 0.4000 #H-H epsilon sigm # 108.0e-21 32.0e-11

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

bond_coeff 1 450 0.9572 #O-H

angle_coeff 1 55 104.52 #H-O-H

# ++++++++++++++++setting+++++++++++++++++++++

neighbor 2.0 bin

neigh_modify delay 0 every 1 check yes

min_modify dmax 0.01

minimize 1.0e-8 10.e-4 100000 300000

timestep .0000000000000002 # ${dt}

thermo $d

velocity all create 298 4928459 rot yes dist gaussian #23482341

fix 1 hydrogen shake 1e-6 500 0 m 1.0 a 1 #for hydrogen

fix 12 water npt temp 298 298 100.0 iso 0.0 0.0 1000.0

# ---------- Relaxation -----------------------------------------

# minimization : avoid atoms overlapping

#min_style fire

#thermo_style custom step etotal enthalpy pe press ke

#thermo_modify flush yes

run 1200

reset_timestep 0

#------------------------dump--------------------------------

#dump 1 all custom 10000 dump.equilibrium. id type x y z vx vy vz

# settings

# Green-Kubo viscosity calculation

# Define distinct components of symmetric traceless stress tensor

variable pxy equal pxy

variable pxz equal pxz #-press

variable pyz equal pyz

fix SS all ave/correlate $s $p $d &

v_pxy v_pxz v_pyz type auto file S0St.dat ave running

variable scale equal 1.0/$t*vol*s*dt
variable scale equal {convert}/(${kB}*$T)*$V*s*{dt}

variable v11 equal trap(f_SS[3])*{scale}
variable v22 equal trap(f_SS[4])*{scale}
variable v33 equal trap(f_SS[5])*${scale}

#*******************************

*${scale}*

variable k33 equal trap(f_JJ[5])

*{scale} variable k22 equal trap(f_JJ[4])*{scale}*

variable k11 equal trap(f_JJ[3])*#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}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/$V

variable k11 equal trap(f_JJ[3])

variable k33 equal trap(f_JJ[5])

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 400

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”