[交流] 【转帖】热导率计算的in文件 已有12人参与

看到有不少人在找热导率计算方面的in文件，我就贡献三个in文件吧，仅供参考。     同时，附件里贴出了我的计算结果。EMD的输出结果（compute heat/flux command+compute tc command的计算结果）中， “ac.dat”（见附件中的"ac.wmf"）是热流自相关函数（我已经修改了compute_tc.cpp，目前输出的是normalized HCACF，但结果中给出的还是没有归一化的热流自相关函数，但形状和归一化的是一样的，请注意！）随m的变化，"tc.dat"(见附件中的"tc.wmf"是热导率随m的变化（m的涵义请参看热导率计算的Green Kubo离散化公式,见附件"Comparison of atomic-level simulation methods for computing thermal conductivity”中的(9)式），"tc_time.dat"(见附件中的"tc_time.wmf"是热导率随时间的变化；NEMD的结果中，"temp.profile"(见附件中的"temp_distribution.wmf"是z方向上的温度分布，"thermal_conductivity.dat“(见附件中的"NEMD.wmf"表示热导率随时间的变化，两者都包含了fix heat command 和fix thermal/conductivity command的计算结果。 1.用compute heat/flux command+compute tc command得到热流自相关函数和热导率（EMD方法） # MD simulation of Ar thermal conductivity # Initialization units             lj dimension         3 newton            on boundary          p   p   p atom_style        atomic neighbor          0.3    bin neigh_modify      check  yes lattice           fcc   0.844 region            box  block -4  4  -4  4  -4  4  units lattice create_box     1  box create_atoms     1  box mass      1  1.0 velocity          all  create  0.71 458127641 mom yes  rot yes dist gaussian units box # LJ potential ********************************************************* pair_style        lj/cut 2.8 pair_coeff        1  1   1.0                  1.0                     #  LJ parameters for Ar-Ar fix               temp all  temp/berendsen 0.71 0.71 0.000466 fix               nve  all  nve thermo_style      custom step temp etotal vol thermo_modify     lost warn thermo            100 # Run timestep          0.000466 run               200000 reset_timestep    0 # -------------- Flux calculation in nve --------------- compute     myKE all ke/atom compute     myPE all pe/atom compute     myStress all stress/atom virial variable          factor_ac equal 1.0 variable          factor_tc equal 1.3806504e-23*sqrt(1.67e-21/6.633e-26)/3.405e-10^2 compute      jflux all heat/flux myKE myPE myStress compute           tc all tc c_thermo_temp c_jflux v_factor_ac v_factor_tc iso first 10000 900000 100000 fix               tc_out  all  ave/time  1  1  1  c_tc   file  tc_time.dat thermo_style      custom  step  temp restart           100000   restart.*             run             1000000 2. 用fix thermal/conductivity command得到温度梯度，进而得到热导率(NEMD方法) # MD simulation of Ar thermal conductivity # Initialization units             lj dimension         3 newton            on boundary          p   p   p atom_style        atomic neighbor          0.3    bin neigh_modify      check  yes lattice           fcc   0.844 region            box  block -4  4  -4  4  -4  4  units lattice create_box     1  box create_atoms     1  box region            up1    block  INF INF  INF INF  -0.5  -0.25  units lattice region            up2    block  INF INF  INF INF   0.5   0.75  units lattice region            up  union 2 up1 up2 region            down1  block  INF INF  INF INF  -3.5  -3.25  units lattice region            down2  block  INF INF  INF INF   3.5   3.75  units lattice region            down union 2 down1 down2 mass      1  1.0 velocity          all  create  0.71 458127641 mom yes  rot yes dist gaussian units box # Tersoff potential ********************************************************* pair_style        lj/cut 2.8 pair_coeff        1  1   1.0                  1.0                     #  LJ parameters for Ar-Ar fix               temp all  temp/berendsen 0.71 0.71 0.0466 fix               nve  all  nve compute           ke  all  ke/atom variable          temp atom  c_ke/(1.5*1.0) fix               temp_profile    all    ave/spatial  1  100000  100000  z  lower  0.25      v_temp  file  temp.profile    units  lattice compute           up_temp    all  temp/region up compute           down_temp  all  temp/region down variable          delta_temp   equal  c_up_temp-c_down_temp fix               delta_out  all  ave/time  1  100000  100000  v_delta_temp   file  delta_temp.dat thermo_style      custom step temp etotal vol thermo_modify     lost warn thermo            100 # Run timestep          0.000466 run               100001 unfix             temp fix               heat_swap   all  thermal/conductivity  10  z  32 fix               e_exchange  all  ave/time  10  10000  100000  f_heat_swap  file  e_exchange.dat    variable          thermal_conductivity equal f_e_exchange/(0.000466*10.0*4.0*f_delta_out)*1.3806504e-23/3.405e-10/3.405e-10*sqrt(1.67e-21/6.633e-26)*6.0/8.0 # 以上variable命令需要特别注意，因为我所模拟的系统，盒子边长Lx=Ly=Lz，热导率计算公式经过推导变成为e_exchange/(4.0*t*L*delta_T)， # 为了不在in文件里给L赋值，我修改了fix_thermal_conductivity.cpp文件(见附件)，将e_exchange修改成了 e_exchange += force->mvv2e * (all[0].value - all[1].value) / (domain->zprd); 同时在end_of_step() 里添加了一句 “   e_exchange = 0.0;“,详见附件中的fix_thermal_conductivity.cpp，这样所得的e_exchange曲线基本上是一条水平曲线，而不是用原来的fix thermal/conductivity command所得到的斜向上的曲线，请注意！！！ # 所以才出现以上variable的表达式。 # 请看明白后再做计算，免得算出错误的结果！！！ fix               thermal_conductivity_out  all  ave/time  100000  1   100000  v_thermal_conductivity   file  thermal_conductivity.dat                    # Run run               10000000 3. 用fix heat command建立温度梯度，进而得到热导率（NEMD方法） # MD simulation of Ar thermal conductivity # Initialization units             lj dimension         3 newton            on boundary          p   p   p atom_style        atomic neighbor          0.3    bin neigh_modify      check  yes lattice           fcc   0.844 region            box  block -4  4  -4  4  -4  4  units lattice create_box     1  box create_atoms     1  box region            up1    block  INF INF  INF INF  -0.5  -0.25  units lattice region            up2    block  INF INF  INF INF   0.5   0.75  units lattice region            up  union 2 up1 up2 region            down1  block  INF INF  INF INF  -3.5  -3.25  units lattice region            down2  block  INF INF  INF INF   3.5   3.75  units lattice region            down union 2 down1 down2 region            hot  block  INF INF  INF INF     0.0   0.25  units lattice group             hot  region  hot region            cold  block  INF INF  INF INF   -4.0  -3.75  units lattice group             cold  region  cold mass      1  1.0 mass0             6.633e-26 epsilon0          1.67e-21 sigma0            3.405e-10 velocity          all  create  0.71 458127641 mom yes  rot yes dist gaussian units box # Tersoff potential ********************************************************* pair_style        lj/cut 2.8 pair_coeff        1  1   1.0                  1.0                     #  LJ parameters for Ar-Ar fix               temp all  temp/berendsen 0.71 0.71 0.0466 fix               nve  all  nve compute           ke  all  ke/atom variable          temp atom  c_ke/(1.5*1.0) fix               temp_profile    all    ave/spatial  1  100000  100000  z  lower  0.25      v_temp  file  temp.profile    units  lattice compute           up_temp    all  temp/region up compute           down_temp  all  temp/region down variable          delta_temp   equal  c_up_temp-c_down_temp fix               delta_out  all  ave/time  1  100000  100000  v_delta_temp   file  delta_temp.dat thermo_style      custom step temp etotal vol thermo_modify     lost warn thermo            100 # Run timestep          0.000466 run               100001 unfix             temp fix               hot   all  heat  1   50   region hot fix               cold  all  heat  1  -50   region cold variable          thermal_conductivity equal 50.0*0.5*1.67e-21/3.405e-10/sqrt(6.633e-26/1.67e-21)/((4.0*8.0*8.0*8.0/0.844)^(1.0/3.0)*3.405e-10*2.0*f_delta_out*1.67e-21/1.3806504e-23)*6.0/8.0 fix               thermal_conductivity_out  all  ave/time  100000  1   100000  v_thermal_conductivity   file  thermal_conductivity.dat                    # Run run               10000000 【转自】http://www.mdbbs.org/thread-26451-1-1.html

回复此楼