Nemd 方法计算水的热导率

1.模型简介：两侧硅壁面中间为水

2.in文件： 

 1.变量设置：

variable n index   2 3 4 5 
log graphene-h2o.log
variable       T equal 353
variable       t equal 293
variable	   rc equal 10
variable       V equal vol
variable       dt equal  1e-3
variable    kb equal 8.63e-5
variable     p equal 50000         # correlation length
variable     s equal 10          # sample interval
variable     d equal $p*$s       # dump interval
variable    nd      equal 10                
variable    run     equal ${nd}*$d              
variable	  tdamp	 equal   ${dt}*100
variable	  pdamp	 equal   ${dt}*1000
variable  lattice_high  equal 1.8075
variable  lattice_low   equal 2.7156 
variable  z_hezi        equal     146.044696756
variable  Dscale_high     equal   ${lattice_high}/${z_hezi}
variable  Dscale_low      equal   ${lattice_low}/${z_hezi}

echo          screen
units          metal
dimension      3
boundary     p p p 
atom_style  full
bond_style     harmonic
angle_style    harmonic
read_data      si-h2o.data
################################################################
variable Z1 	 equal zlo
variable Z2 	 equal zhi
variable   Len   equal   ${Z2}-${Z1}
variable   B1    equal   ${Z1}+3*${lattice_high}              #fix1   
variable   B2    equal   ${Z1}+15*${lattice_high}             #bath1
variable   B26   equal   ${Z2}-15*${lattice_high}             #bath2
variable   B27   equal   ${Z2}-3*${lattice_high}             #fix2

2.region和势函数设置

region         wall1   block       INF INF INF INF  INF     ${B1}        units box side in
region         hot    block         INF INF INF INF ${B1}   ${B2}      units box  side in
 
region         cold   block          INF INF INF INF    ${B26}   ${B27}   units box  side in
region         wall4   block         INF INF INF  INF   ${B27}    INF   units box  side in
region         wall union 2 wall1 wall4 units box side in
##############################################################
 pair_style  hybrid  lj/cut/tip4p/long   2 3 1 1 0.150 10 8.5   lj/cut 10    sw
pair_coeff * *      sw     Si.sw    Si  NULL NULL 
pair_coeff  1 2      lj/cut                       0.121  2.631
pair_coeff  2  2     lj/cut/tip4p/long      0.006739    3.166  
pair_coeff  *  3      lj/cut                       0      0
  
bond_coeff      1 450  1
angle_coeff      1 55   109.47

kspace_style   pppm/tip4p 1.0e-5
set type 2 charge -0.8476  
set type 3 charge 0.4238
 
mass  1  28.086000  # si
mass   2  15.999400 # o*
mass   3   1.007970  # hw
#############################################################
group    water  type  3 2
group    si  type 1
group    wall  region wall
group    wall1 region wall1
group   wall4 region wall4
group   hot  region hot
group  cold  region cold
group  heat  subtract si hot cold wall
group   other union water heat hot cold 
 

3.计算部分设置：

 
timestep            ${dt}
neighbor	          2 bin
neigh_modify    every 2 delay 10 check yes 

#fix                     003 si spring/self 30
#fix_modify        003 energy yes

fix              202 si  setforce 0 0 0
min_style         cg
minimize         1e-7 1e-7 10000  10000
unfix  202
######驰预阶段
compute         Twater water temp
compute_modify Twater dynamic/dof yes extra/dof 0
compute         Thot hot temp
compute        Tcold cold temp
compute        Twall wall temp
compute    Theat heat temp
thermo           1000


dump         1 all custom 1000 water1.lammpstrj id type x y z
velocity       all create $T 1565155 mom yes rot yes dist gaussian units box
fix              water  water shake 1e-4 100 0 b 1 a 1
 
fix              02 wall  setforce 0 0 0
velocity      wall set 0 0 0

fix             10     water nvt temp $t $t ${tdamp}
fix             110   hot nvt temp $t $t ${tdamp}
fix             1110 cold nvt temp $t $t ${tdamp}
fix             101   heat nve
thermo_style    custom step etotal pe ke  temp press density c_Twater c_Theat  c_Tcold c_Thot c_Twall 
run            400000
unfix           10
 unfix           110
unfix           1110
unfix          101

fix             1     water nve
fix             121 heat nve
fix             11   hot nvt temp $t $t ${tdamp}
fix             111 cold nvt temp $t $t ${tdamp}
thermo_style    custom step etotal pe ke  temp press density c_Twater c_Theat c_Tcold c_Thot c_Twall 
run           400000
unfix           1
 unfix           11
unfix           111
unfix     121
#################### NVE ###########################
reset_timestep 0
fix                 6 other nve 
 timestep            1e-3
#fix 66  water efield        0     0     1    ##########################z方向施加1v/A的电场############
fix             hot  hot  langevin  $T $T  0.1 598034 tally yes
fix             cold  cold  langevin  $t $t  0.1 59804 tally yes
thermo          $d
compute         myKE1 water ke/atom
variable        temp1 atom c_myKE1/1.0/${kb}
compute         myKE2 si ke/atom
variable        temp2 atom c_myKE2/1.5/${kb}

compute         layer all chunk/atom bin/1d z lower ${Dscale_high}  units reduced
fix	      7 all ave/chunk $s $p $d layer v_temp1  v_temp2   temp density/mass density/number   file si-profile.txt

compute         layers all chunk/atom bin/1d z lower ${Dscale_low} units reduced
fix	      7 all ave/chunk $s $p $d layers   v_temp1 v_temp2    temp density/mass density/number   file water-profile.txt
 
thermo_style    custom step etotal pe ke  temp press density c_Twater c_Theat c_Tcold c_Thot c_Twall  f_hot f_cold
run              ${run}    

3.计算结果；0.60 W/mK

欢迎关注公众号------硕博科研小助手