请选择 进入手机版 | 继续访问电脑版

分子模拟论坛 Molecular Simulation Forums

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 879|回复: 10

[原创]热导率计算的in文件

[复制链接]

11

主题

117

帖子

251

积分

中级会员

Rank: 3Rank: 3

积分
251
发表于 2010-6-7 14:39:00 | 显示全部楼层 |阅读模式
本帖最后由 zhxlhdd2008 于 2010-10-28 16:23 编辑
看到有不少人在找热导率计算方面的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
unitslj
dimension3
newtonon
boundaryppp
atom_styleatomic
neighbor0.3bin
neigh_modifycheckyes
latticefcc0.844
regionboxblock -44-44-44units lattice
create_box1box
create_atoms1box
mass11.0
velocityallcreate0.71 458127641 mom yesrot yes dist gaussian units box
# LJ potential *********************************************************
pair_stylelj/cut 2.8
pair_coeff111.01.0#LJ parameters for Ar-Ar
fixtemp alltemp/berendsen 0.71 0.71 0.000466
fixnveallnve
thermo_stylecustom step temp etotal vol
thermo_modifylost warn
thermo100
# Run
timestep0.000466
run200000
reset_timestep0
# -------------- Flux calculation in nve ---------------
computemyKE all ke/atom
computemyPE all pe/atom
computemyStress all stress/atom virial
variablefactor_ac equal 1.0
variablefactor_tc equal 1.3806504e-23*sqrt(1.67e-21/6.633e-26)/3.405e-10^2
computejflux all heat/flux myKE myPE myStress
computetc all tc c_thermo_temp c_jflux v_factor_ac v_factor_tc iso first 10000 900000 100000
fixtc_outallave/time111c_tcfiletc_time.dat
thermo_stylecustomsteptemp
restart100000restart.*
run1000000
2. 用fix thermal/conductivity command得到温度梯度,进而得到热导率(NEMD方法)
# MD simulation of Ar thermal conductivity
# Initialization
unitslj
dimension3
newtonon
boundaryppp
atom_styleatomic
neighbor0.3bin
neigh_modifycheckyes
latticefcc0.844
regionboxblock -44-44-44units lattice
create_box1box
create_atoms1box
regionup1blockINF INFINF INF-0.5-0.25units lattice
regionup2blockINF INFINF INF0.50.75units lattice
regionupunion 2 up1 up2
regiondown1blockINF INFINF INF-3.5-3.25units lattice
regiondown2blockINF INFINF INF3.53.75units lattice
regiondown union 2 down1 down2
mass11.0
velocityallcreate0.71 458127641 mom yesrot yes dist gaussian units box
# Tersoff potential *********************************************************
pair_stylelj/cut 2.8
pair_coeff111.01.0#LJ parameters for Ar-Ar
fixtemp alltemp/berendsen 0.71 0.71 0.0466
fixnveallnve
computekeallke/atom
variabletemp atomc_ke/(1.5*1.0)
fixtemp_profileallave/spatial1100000100000zlower0.25v_tempfiletemp.profileunitslattice
computeup_tempalltemp/region up
computedown_tempalltemp/region down
variabledelta_tempequalc_up_temp-c_down_temp
fixdelta_outallave/time1100000100000v_delta_tempfiledelta_temp.dat
thermo_stylecustom step temp etotal vol
thermo_modifylost warn
thermo100
# Run
timestep0.000466
run100001
unfixtemp
fixheat_swapallthermal/conductivity10z32
fixe_exchangeallave/time1010000100000f_heat_swapfilee_exchange.dat
variablethermal_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
[color=]_conductivity.cpp,这样所得的e_exchange曲线基本上是一条水平曲线,而不是用原来的fix thermal/conductivity command所得到的斜向上的曲线,请注意!!!
# 所以才出现以上variable的表达式。
# 请看明白后再做计算,免得算出错误的结果!!!
fixthermal_conductivity_outallave/time1000001100000v_thermal_conductivityfilethermal_conductivity.dat
# Run
run10000000
3. 用fix heat command建立温度梯度,进而得到热导率(NEMD方法)

# MD simulation of Ar thermal conductivity
# Initialization
unitslj
dimension3
newtonon
boundaryppp
atom_styleatomic
neighbor0.3bin
neigh_modifycheckyes
latticefcc0.844
regionboxblock -44-44-44units lattice
create_box1box
create_atoms1box
regionup1blockINF INFINF INF-0.5-0.25units lattice
regionup2blockINF INFINF INF0.50.75units lattice
regionupunion 2 up1 up2
regiondown1blockINF INFINF INF-3.5-3.25units lattice
regiondown2blockINF INFINF INF3.53.75units lattice
regiondown union 2 down1 down2
regionhotblockINF INFINF INF0.00.25units lattice
grouphotregionhot
regioncoldblockINF INFINF INF-4.0-3.75units lattice
groupcoldregioncold
mass11.0
#mass06.633e-26
#epsilon01.67e-21
#sigma03.405e-10
velocityallcreate0.71 458127641 mom yesrot yes dist gaussian units box
# Tersoff potential *********************************************************
pair_stylelj/cut 2.8
pair_coeff111.01.0#LJ parameters for Ar-Ar
fixtemp alltemp/berendsen 0.71 0.71 0.0466
fixnveallnve
computekeallke/atom
variabletemp atomc_ke/(1.5*1.0)
fixtemp_profileallave/spatial1100000100000zlower0.25v_tempfiletemp.profileunitslattice
computeup_tempalltemp/region up
computedown_tempalltemp/region down
variabledelta_tempequalc_up_temp-c_down_temp
fixdelta_outallave/time1100000100000v_delta_tempfiledelta_temp.dat
thermo_stylecustom step temp etotal vol
thermo_modifylost warn
thermo100
# Run
timestep0.000466
run100001
unfixtemp
fixhotallheat150region hot
fixcoldallheat1-50region cold
variablethermal_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
fixthermal_conductivity_outallave/time1000001100000v_thermal_conductivityfilethermal_conductivity.dat
# Run
run10000000

回复

使用道具 举报

80

主题

214

帖子

514

积分

高级会员

Rank: 4

积分
514
发表于 2010-6-7 14:46:00 | 显示全部楼层
收下了 ,受益匪浅

回复

使用道具 举报

33

主题

131

帖子

301

积分

中级会员

Rank: 3Rank: 3

积分
301
发表于 2010-6-7 18:20:00 | 显示全部楼层
本帖最后由 senorcengou 于 2010-6-7 19:01 编辑
lz好人啊!大赞
不过有些地方有点不明白,比如 variabletemp atomc_ke/(1.5*1.0)的意思应该是求得单个原子的温度吧?表达式里面的1.5应该是Ke=3/2*N*Kb*T中的3/2,那么1.0是什么呢?
回复

使用道具 举报

11

主题

117

帖子

251

积分

中级会员

Rank: 3Rank: 3

积分
251
 楼主| 发表于 2010-6-7 22:59:00 | 显示全部楼层
本帖最后由 zhxlhdd2008 于 2010-6-7 23:07 编辑
1.0也就是update.cpp里赋值的force->boltz这个常数(force.h里有定义)
回复

使用道具 举报

4

主题

9

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2010-6-8 00:02:00 | 显示全部楼层
这个算出来的热导率波动的厉害么,我自己算出来的都波动的比较厉害,对于最后热导率应该在哪段范围取值很是头疼
回复

使用道具 举报

33

主题

161

帖子

361

积分

中级会员

Rank: 3Rank: 3

积分
361
发表于 2010-6-8 01:52:00 | 显示全部楼层
在那输出(dump, thermo)的结果中哪一个是热导率???
回复

使用道具 举报

33

主题

161

帖子

361

积分

中级会员

Rank: 3Rank: 3

积分
361
发表于 2010-6-8 09:26:00 | 显示全部楼层
本帖最后由 szd419 于 2010-6-8 10:30 编辑
besides, i can't find the compute tc command from the lastest manual.
-------
Problem is solved.
I find the script you wrote and posted in this website.
Thanks from me and the other guys
回复

使用道具 举报

11

主题

117

帖子

251

积分

中级会员

Rank: 3Rank: 3

积分
251
 楼主| 发表于 2010-6-8 09:29:00 | 显示全部楼层
本帖最后由 zhxlhdd2008 于 2010-6-8 09:44 编辑
compute tc command是我自己写的,参看http://simulation.5d6d.com/thread-26298-1-3.html
回复

使用道具 举报

11

主题

117

帖子

251

积分

中级会员

Rank: 3Rank: 3

积分
251
 楼主| 发表于 2010-6-8 09:43:00 | 显示全部楼层
本帖最后由 zhxlhdd2008 于 2010-6-8 10:21 编辑
计算结果见主题帖中的附件。
回复

使用道具 举报

31

主题

503

帖子

1043

积分

金牌会员

Rank: 6Rank: 6

积分
1043
发表于 2010-6-8 09:48:00 | 显示全部楼层
学习了 呵呵 谢谢分享
回复

使用道具 举报

Archiver|手机版|小黑屋|分子模拟论坛  

GMT+8, 2019-7-16 14:48 , Processed in 0.129691 second(s), 28 queries .

Powered by Discuz! X3.3

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表