科学网

 找回密码
  注册
科学网 标签 DFT

tag 标签: DFT

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

[转载]减少DFT算法模式中实变换时的频域乘法的次数
SciteJushi 2016-5-13 08:17
原载 http://blog.sina.com.cn/s/blog_729a92140102wc7i.html 本短文用“实变换时”表明:限于Tpwp中,特别地假定被处理信号和滤波器冲激响应,都是实数序列。不管Tpwp之外的变换,因为使用DFT的条件和方法不一样。这里,不涉及人们从文献中可能见到的奇偶时间点序列分组、FIR以及时间延迟之类的问题,给定的一个被处理序列x和滤波器序列F(高通或低通),必须靠“移位方向参量为空值”另外指明它们都是频域区间 的结果仍等于v。 以博文《减少小波包变换的DFT算法模式中的DFT的次数和点数》(2015-10-26)为基础,减少了部分复数乘法,但是,新增了复共轭、向量倒序的处理内容,改变了那里的“频域数据重复”与频域的乘积运算的顺序。 用这两个函数子程序,覆盖掉以前用的旧子程序,重新运行博文《冗余小波包合成以及降噪试验和能量守恒定标》(2015-12-27)的最后三个降噪实验。可见在Scilab-5.3.3的命令窗中显示了相同的数字(format(‘v’,10) 格式)结果。但是,有些遗憾,没有见到实验的总耗时量明显减少,甚至反而可见高达百分之十几的增加,这可能与引入的其它处理以及系统工作平台中整体的解释、管理等有关。据此类经验,比较算法的性能时,居士希望尽量交代清楚方案、条件、环境。 新浪赛特居士SciteJushi-2016-05-13。 图片 1. 减少了紧凑DFT算法模式中实变换时的频域乘法的次数的程序段 附 . 最近,搜索小波的最新研究进展,在《百度百科》中,注意到了中国已有研究成果“比小波分析创始人Mallat提出的著名Mallat算法快一个数量级以上,而且效果和质量远优于Mallat算法”。看到这很振奋又好奇,然而还未搜到演示资料或方法细述,故特将其复制为留念于下。原来,还有了国际小波分析应用研究中心。虽然,姓名,是居士在博客中多次提到(因个人认为,翻译得已很不错)的小波名著“十讲”的译者,也是前面(2016-04-20)的实变换的滤波器构造中曾用过那篇中文参考文献的作者。但是,居士真不记得当年离开小波的研究时,对这个姓名已有印象,实际上,除去博文已明确提到的部分,至今仍几乎未知其具体研究(注意,在此无褒贬之意)。 : 李建平(电子科技大学教授) 简介?编辑 1986年获重庆大学应用数学理学学士,1989年获西安交通大学软件工程、计算数学双硕士学位,1998年获重庆大学计算机科学理论工学博士学位,1999年至2001年到香港浸会大学做博士后研究,2002年至2003年到法国普罗旺斯大学做访问学者,2005年至2006年到加拿大圭尔夫大学、多伦多大学做访问学者. 教学情况 现为电子科技大学计算机科学与工程学院、示范性软件学院副院长,国际小波分析应用研究中心主任,国际学术进展 International Progress on Wavelet Active Media Technology and Information Processing主编,国际学术期刊International Journal of Wavelet Multiresolution and Information Processing副主编,兼任国家科学技术奖励评审委员、国家自然科学基金项目评审委员、中华人民共和国公安部技术顾问等十几个学术、社会职务,是2004年国际计算机学术大会、第三届国际小波分析及其应用学术大会(ICWAA)、第二届智能体媒介技术国际学术大会(ICAMT)主席,是国际上小波分析与信息处理研究领域十分活跃的知名专家。 科研情况 他在国际上独立提出并系统建立了“小波变换的加速方法”、“矢量积小波变换理论”、“基于小波分析的电子签名系统”等系列理论与方法,该理论比小波分析创始人 Mallat提出的著名Mallat算法快一个数量级以上,而且效果和质量远优于Mallat算法,这为小波理论在信号处理、信息分析等许多方面的应用提供了先进的算法,为小波分析的实时处理和产品化提供了理论基础,为小波分析开辟了广阔的应用前景。他在国际上首次提出“基于‘三大特征’(机器特征、文档特征、人体特征)的信息最安全传输的模型与方法”,为当今研究热点的网络与信息安全做出了重要贡献。他先后主持国家863高技术项目、国家自然科学基金项目等30项。他在国内外著名学术期刊上发表重要论文150篇,被国际三大检索机构SCI、EI、ISTP等检索收录论文48篇,出版学术专著15部,其中2部被多次修订重印。他主持研制的“小波指纹加密系统”、“分布式网络监控系统”等高技术产品产生了广泛的经济效益和社会影响
789 次阅读|1 个评论
7th International School and Workshop on Time-Dependent DFT
gcshan 2016-2-1 16:29
7th International School and Workshop on Time-Dependent Density-Functional Theory Dear Colleagues, Apologies for possible cross-posting. This is to announce the CECAM / Psi-k funded 7th International School and Workshop on Time-Dependent Density-Functional Theory: Prospects and Applications which will be held September 12th - September 23rd, 2016, at the beautiful Centro de Ciencias de Benasque Pedro Pascual (Benasque, Spain; http://benasque.org/ ) Please see http://www.benasque.org/2016tddft/ for details, and read more below! Summary: This School Workshop is the seventh of a very successful series that started in 2004. The positive response to this first event, also held in the “Centro de Ciencias de Benasque Pedro Pascual”, from August 28th to September 12th 2004, encouraged the organization of the sequels in the same place with an approximate periodicity of two years. The purpose has in all occasions been to (1) make a very intense introduction to both the theory, the practice, and the numerical implementation of time-dependent density-functional theory (TDDFT), mainly (but not exclusively) oriented to young scientists willing to initiate or strengthen their knowledge and skills on TDDFT, followed by (2) a workshop on the subject in which all the main aspects are to be covered by the leading experts. All the students of the school are expected to participate in the workshop, in order to learn about the state-of-the-art of the subject, after being exposed to the fundamentals. Since TDDFT is a rapidly evolving field of Science, the precise content of both school and workshop have changed over the years –although the format of the events has been largely unaltered. In all occasions there has been a very large number of applicants for the school, that has increased every edition to become more than 150. This is not only a testimony of the strong pulse of the scientific field itself, but also of the good quality of the school. Since we want to maximize the learning experience of the students via a close interaction with the teachers (and also due to the logistic limitations imposed by the hands-on tutorial), the participants of the school will not exceed 40. The total number of participants in the full School Workshop has been close to 100 in all occasions. It is worth mentioning that participants came from all over the world, making this series of schools and workshops a truly global event. Location/Timing: This event will take place at the “Centro de Ciencias de Benasque Pedro Pascual”, Benasque, Spain ( http://benasque.org/ ), from September 12th -- 23rd, 2016 (travel to Benasque on the 11th). Benasque is a beautiful town in the heart of the Pyrenees. The school will take place from September 12th to September 19th, while the workshop will be from September 20th to September 23rd. Participants The call for participation will be mainly directed to students and scientists specialized on computational physics, quantum chemistry and biophysics. We call for students willing to participate at the School (and attend the workshop), and for scientists willing to present a contributed talk or poster at the Workshop. We will limit the number of students to the school to 40 and participants to the workshop to less than 100, in order to ensure a maximum interaction between all the scientists participating. Attendance of graduate students and post-docs will be strongly encouraged through the inclusion of short contributed talks and a poster session. Furthermore, we will award to PhD students who present an outstanding poster short oral presentations. Applications/Support: All persons who wish to participate should fill out the application form at: http://benasque.org/2016tddft/nbsp ; In the comments section, please indicate if you wish to participate in the Summer Summer School or in the Workshop (or in both). For participants coming from the USA, please check the following address for support: http://www.mcc.uiuc.edu/ School only: As we have a very limited number of places for the school, students will be selected from among an open pool of applicants who have demonstrated a strong interest in computational sciences, applied to chemistry, physics, materials science and biology.
个人分类: Information|6 次阅读|0 个评论
[转载]减少小波包变换的DFT算法模式中DFT/IDFT的次数和点数
SciteJushi 2015-10-26 09:33
原载 http://blog.sina.com.cn/s/blog_729a92140102vwr9.html 在小波(包)和滤波器的普通理论中,常见的是Z变换和傅立叶积分变换。为应用方便,基于子带分解塔、周期化共轭滤波器组和Wickerhauser等的最佳基思想,居士建立有限长复向量空间的小波包分割与覆盖以及多分辨分析。其重要意义,不受限制于,连续时间小波函数和无穷时间序列滤波。在博客中,已包括了许多创新,不明现连续时间的基函数(多分辨分析和小波的纯粹数学理论的主要构造目标)、Z变换域分析以及滤波器阵的反混叠条件(anti-aliasing condition)等内容,而反复强调了关键词“基向量”。因为其尺度、时间和频率都是有限的离散量,所以这种小波包变换,与离散傅立叶变换的联系,更直接,工整自然,当然,其理论基础需要Z变换和连续频率。 以《随机设置小波包变换的DFT算法模式》(2015-05-13)为基础,还可使离散傅立叶变换(从通常的数字信号处理教科书即可知其所需性质),与小波包变换的算法实现之间,更紧密地结合。例如: (1).底层(相应于常规的“卷积、滤波、下抽样、插0”处理)的子程序可输入和输出频域数据,供调用程序缓存和重复利用,以免相邻尺度的运算过程中的反复的DFT和IDFT处理。尤其在信号重建中,一般并不需要中间经过的节点上的时域数据。高通支路输出与低通支路输出的相加,在时域或频域均可实现。 (2).用频域半周期折叠平均,代替时域隔2下抽样 ( 2- 抽样,二个点中取一舍一) ,使获得时域系数的IDFT的点数减半,回避下抽样时要抛弃的一半。 (3).用频域数据“重复”,代替时域插0变抽样率,绕过了插零,可使求频域系数的DFT的点数减半。 新旧两套方案对照,如图片1.的右半部的同名文件所示。 这种调整,不可避免地要改变数字误差,就不妨顺便用矩阵乘法代替用求和函数sum做非DFT模式的分解。调用程序,需较大改动,变得更复杂些,尤其在计算所谓的双向树时,才能兼容地管理时域数据和频域数据以及非DFT模式,但是,这不影响此前讨论的应用和测试程序运行及判断。 重新运行《用小波包变换试验多孔滤波器组的零频点移位》(2015-09-09)和《用双重噪声化相谱的多孔滤波器组测验小波包双向树》(2015-09-25)的测验程序,其状况如图片1.左上半部所示。 第二测验程序的耗时仍很长,其中约40%的测试不使用DFT算法模式。仅改动其中的语句“if abs(Wp.RightShift)5e5, Wp.RightShift= ;”即强迫全部用DFT算法模式处理,可重命名函数和运行测试,则耗时从2286秒降至了568秒。 新浪赛特居士SciteJushi-2015-10-26。 图片 1. 更紧凑DFT算法模式及其三个测验程序的情况
736 次阅读|0 个评论
原子电子结构计算程序atom-0.1.0
yyhxwhd 2014-9-21 22:59
atom-0.1.0 atom-0.1.0.tar.gz 首先声明:这个程序的原始版本不俺写的,但是在此基础上我加入了一些个人的想法 包括:1.修改了之前的CMake ,重新按照在SIESTA Makefile的格式进行了修改 2. 加入了SIESTA fdf 库函数,方便标签读取 3. 加入了文件IO操作子过程,io.f90 4. 加入了可以扩展的数学库 现在程序可以计算中性原子的电子波函数,对应的本征值,势函数等 这个程序可以作为初学者学习电子结构基本理论的实例,当然程序还有许多地方需要修改,最终目标是将其修改成可以产生赝势的程序包。如果有谁对这一块比较感兴趣,可以一起讨论讨论.
个人分类: Program;Script|3082 次阅读|0 个评论
[转载]用基向量从随机设置小波包变换树上捕捉DFT负频率的幽灵
SciteJushi 2014-3-27 15:53
原载 http://blog.sina.com.cn/s/blog_729a92140101oxbx.html 在Matlab-R2011a中,用 =wavefun('bior2.2',10); subplot(2,1,1); plot(x,w); axis tight; subplot(2,1,2); plot(x,w1); axis tight; 看到国际标准 JPEG2000的小波的“连续时间函数”,像感觉到分形曲线的味道。至于得到的bior3.1、bior3.3、bior3.5的函数波形,如果是居士自己搞的,那么很可能自认为它们是错误的。不过,尖叫“Matlab牛”后,居士还是偏爱“基向量”,并以为,本短文对简化、理清文献中的连续时间小波包基函数、其自然排序、其频率排序等问题的长期缠绕,有重要意义。 无先验知识以明确被分析数据与尺度函数的关系,用Tpwp的“小波包基向量”,无需工具箱生成函数(Building Wavelet Packets, Wavelet Packet Atoms),即可方便地确定小波包变换树上任意位置的频率中心(功率或能量谱中心,区别正频率、负频率和直流)。 常用三个索引指标,尺度、频率、位置,标示一个小波包。但是,频率指标值越大,不同于1兆赫兹大于1千赫兹这样的事情,也不必意味着相应的小包中的高频能量(功率)就更大、或在低频时更小。实际上,尺度指标,常与水平(level)、深度(depth)混用,也并不直接就是小波函数概念中的“尺度值”。频率指标是,把“基本 小波和尺度函数 (或序列) ”以适当方式“运算、组合”而 形成某个“波的包裹” 的过程的一种编码。在离散变换中,可以简明地看到“高通、低通”滤波器的组合过程,但是注意,“下抽样”使频谱搬移了、使高低频率的位置转换了! 普通离散小波分解中,第一级分解后,得到的“细节系数序列”中的“直流部分”,对应原信号中的“高频成分”。 在某个尺度上,如果它的节点数(频率指标的个数)不超过被处理序列的长度N的一半,那么简单用功率谱的中心、波形过零率等物理直观真实的正频率概念,区别小波包,是自然的。 然而,在最大的尺度上,节点数(频率指标的个数)等于N,那么这些普通的“频率分辨观念”就可能显得不方便。例如,实数序列的离散傅里叶变换DFT中正频率部分(0,fs/2)中,只有N/2-1点,难道,小波包分解能分辨出约N个“正”频率点吗? 这时,小波包基向量,存在配对关系:幅度谱相同,则其定义的频率中心就相同,然而,相谱不同、有不相等的频率索引指标。这类似于DFT中的正负频率、共轭对称、幅相谱等的状况和问题。单个节点内,不再存在,可以用于分辨相位的“位置指标”。 如果遵从那些小波工具箱的处理,用滤波器和信号序列的长度,来限制分解深度,那么就遇不到这种状况。 在《随机设置小波包变换及其优选基的随机降秩矩阵》(2014-03-19)的基础上,居士做试验程序PwpRandNegf.sce,附于本文末。运行 clear ; LessLevel = 0 ; rand ( 'seed' , 1e9 ) ; exec ( 'PwpRandNegf.sce' ) ; 可得图片 1。 假设信号序列长度,随机地取为16、32、64、128,相应的最大分解深度分别为4、5、6、7。在最大分解尺度上,频率分辨率达到极限。抽出最大尺度上的所有小波包基向量,做两两比较,检验它们的FFT幅度谱的关系。 两基向量,其规范化幅度谱之间的误差,RMSE,er0,小于千之一,则视为是配对的。对于正交变换,附带地确认:两基向量为等腰三角形的直角边,即不可能是因为两向量本身相同所以幅度谱才相同,如若不然,程序报告错误。 从数据序列长度的二分之一中,减去(倒记数)相配的对数,得图片1.中的第一幅曲线图。结果始终为1,与滤波器、小波包变换的随机设置、序列长度都无关。 余下一对(故结果恒为1):涉及直流成分,和高频向量。980次实验的直流基向量偏离理想值的误差(RMSE)的累加总和,TestDC = 0.000038,很小。高频向量,被归入某对的计数结果,TestHF = 0。 配对向量的幅度谱之间的实际误差er0的总和,示于图片1.中的第二幅曲线图。“总和”用的数值的数目,与序列长度有关,但结果仍都很小。 在一个稳定的变换中,尺度参数、频率参数都相同而只是位置参数不同的小波包基向量,它们也具有相同的幅度谱,但与图片1.反应的情形不同。例如,没有直流基向量;配对向量,经历的数字运算过程更相近,所以它们的幅度谱之间的数字误差更小。 运行 clear ; LessLevel = 1 ; rand ( 'seed' , 1e9 ) ; exec ( 'PwpRandNegf.sce' ) ; 可得图片 2,反应“次最大”尺度上的基向量的情况。无单纯的直流基向量,所以TestDC = 125.14948,很大;TestHF = 980,即,等于随机实验次数,所以高频向量也总是成对的。 由于,必须覆盖相同的物理正频率带,所以, LessLevel =0 时的向量对,与 LessLevel =1 时的向量对,不可能有实质上不同的正频率带宽度,然而频率中心须错位。这种错位,改变了正、负频率中心的距离,提供了一种频率信息“分辨能力”。在 LessLevel =0 时,配对向量,可来自前面的尺度上的不同频率指标块。为此转换问题,居士杜撰“负频率幽灵”一词,供有兴趣者参考。 程序的Matlab版的处理,约更快1倍。 附程序: PwpRandNegf-2014-03-27.zip // PwpRandNegf.sce // for the concept: wavelet packet basis vector // test the random settings of wavelet packet vectors // and the ghost/compare of negative frequency of DFT basis // in order to run with fixed seed, use one of following lines: // clear; LessLevel=0; rand('seed',1e9); exec('PwpRandNegf.sce'); // clear; LessLevel=1; rand('seed',1e9); exec('PwpRandNegf.sce'); // in order to run with variable seed, use one of following lines: // clear; LessLevel=0; exec('PwpRandNegf.sce'); // clear; LessLevel=1; exec('PwpRandNegf.sce'); // note: TpwpSubs.bin, TpwpSubs_E01.bin, MatOrtWlts.bin are // in Scilab current directory // reference: PwpRandPsbm.sce, uploaded, 2014-03-19 // in Scilab-5.3.3,Baiyu Tang( tang.baiyu@gmail.com ) // last revised,2014-Mar xdel(winsid()); // kill all figures mode(0); ieee(1); // clear; // LessLevel=0; // rand('seed',1e9); // ---start timing date1=getdate(); date1=date1( ); disp( ); tic(); load('TpwpSubs.bin'); // load function subroutines load('TpwpSubs_E01.bin'); // extended set of subroutines for PreSet Basis load('MatOrtWlts.bin'); // load the filter cases TestDC=0; // DC vector error, all, total TestHF=0; // count high frequency vectors in pairs TestCount=zeros(98,10); // count matrix TestError=zeros(98,10); // error matrix // ---compute for ii=1:98; // filter case index. 98 filter cases for jj=1:10; // random test index. 10 tests per filter-case // ---random length and decomposition depth of test signal Wp.DecL=4+round(3*rand(1,1,'uniform'));// maximal level, 4,5,6 or 7 Wp.DatL=2^Wp.DecL; // data length, 16,32,64 or 128 Wp.DecL=Wp.DecL-LessLevel; // less-level decomposition CIndx0= *(Wp.DecL+1); // coefficient index, Last-Level // ---random settings Set0=rand(1,9,'normal')*1e6; Wp.RightShift=Set0(1); Wp.ExchangeR_D=Set0(2); Wp.FlipFirst=Set0(3); Wp.MatchFilter=Set0(4); Wp.LowT_zero=fix(Set0(5)); // integer Wp.HighT_shift=fix(Set0(6)); // integer Wp.FlipTime=Set0(7); Wp.PeriodFirst=Set0(8); Wp.AlterL4fSign=Set0(9); // ---get filters if Wp.ExchangeR_D0 then =WltFilters(ii); else =WltFilters(ii); end if length(hcoef(:))2 then if length(rhcoef(:))1 then hcoef=rhcoef; rhcoef= =TpwpAllGH(hcoef,rhcoef,Wp.DatL,Wp.DecL,... Wp.MatchFilter,Wp.PeriodFirst,Wp.LowT_zero,Wp.HighT_shift,Wp.FlipTime); if Wp.AlterL4fSign0 then // forced sign. usually not used hlen=max(round(length(hcoef)/2),round(length(rhcoef)/2))/2; if fix(hlen)==hlen then pgcoef=-pgcoef; // rpgcoef=-rpgcoef; end end =PsbMatrix(CIndx0,Wp.DecL,Wp.DatL,phcoef,pgcoef,Wp.RightShift); // ---match, compare x1=tM0(1,:)-1/sqrt(Wp.DatL); // DC error, realization and idea TestDC=TestDC+sqrt(mean(x1.*x1)); // DC RMSE accumulation TestCount(ii,jj)=Wp.DatL/2; // half the length of vectors and signals tM0f=zeros(Wp.DatL,Wp.DatL); for k1=1:Wp.DatL; x0=fft(tM0(k1,:)); // fast Fourier Transform, DFT x1=abs(x0); // magnitude spectrum tM0f(k1,:)=x1/sqrt(sum(x1.*x1)); // normalized with power end for k1=1:(Wp.DatL-1); // check the pairs of base vectors k12=0; // count vectors in one group for k2=(k1+1):Wp.DatL; x1=tM0f(k1,:)-tM0f(k2,:); er0=sqrt(mean(x1.*x1)); // magnitude RMSE if er00.001 then // equal if ii81 then // orthogonal cases x1=tM0(k1,:)-tM0(k2,:); er1=sqrt(sum(x1.*x1)); // vector error, l_2 norm if er11.4 then // confirm difference, two disp( ); error('Strange Vector Pair.'); end end if (k1==(Wp.DatL/2+1))|(k2==(Wp.DatL/2+1)) then TestHF=TestHF+1; // count, high frequency end k12=k12+1; TestCount(ii,jj)=TestCount(ii,jj)-1; // count, pairs TestError(ii,jj)=TestError(ii,jj)+er0; // RMSE accumulation,total end end if k121 then disp( ); error('The group has too many vectors for PwpRandNegf.sce .'); end end end end // ---display TestDC, TestHF, subplot(2,1,1); plot2d('nn',TestCount, rect= ); set(gca(),'tight_limits','on'); set(gca(),'box','on'); set(gca(),'grid', ); xlabel('Filters: sym1-sym35,coif1-coif5,dmey,db4-db42,18 biorthogonal cases'); ylabel('Length by 2 minus Magnitude-Match Count'); AuLabel1= ; xstring(1,-4+%eps,AuLabel1); Note1= ; if LessLevel0 then title('PwpRandNegf.sce: One Level Less than Maximal Decomposition'); xstring(1,2,Note1( )); else title('PwpRandNegf.sce: Tpwp Transform Trees Hold Negative Frequency of DFT ?'); xstring(1,1.5,Note1( )); end subplot(2,1,2); plot2d('nl', TestError+%eps, rect= ); set(gca(),'tight_limits','on'); set(gca(),'box','on'); set(gca(),'grid', ); xlabel('Filters: sym1-sym35,coif1-coif5,dmey,db4-db42,18 biorthogonal cases'); ylabel('Error of Magnitude of Matched Pairs'); xstring(1,2e-16,AuLabel1); if LessLevel0 then title('PwpRandNegf.sce: One Level Less than Maximal Decomposition'); xstring(1,1e-5,Note1( )); else title('PwpRandNegf.sce: Tpwp Transform Trees Hold Negative Frequency of DFT ?'); xstring(1,1e-5,Note1( )); end // ---timing Time_In_second=toc(), // --clear; LessLevel=0; rand('seed',1e9); exec('PwpRandNegf.sce'); // !Date and Time: 2014 3 25 8 41 23 ...... ! // TestDC = // 0.0000380 // TestHF = // 0. // Time_In_second = // 209.516 // --clear; LessLevel=1; rand('seed',1e9); exec('PwpRandNegf.sce'); // !Date and Time: 2014 3 25 8 46 55 ...... ! // TestDC = // 125.14948 // TestHF = // 980. // Time_In_second = // 194.25 // ***** The Version in Matlab ***** // clear; LessLevel=0; rng(1e9); PwpRandNegf; // Date and Time: 2014 3 25 8 57 ...... // TestDC = // 4.6382e-005 // TestHF = // 0 // Time_In_second = // 84.647 // clear; LessLevel=1; rng(1e9); PwpRandNegf; // Date and Time: 2014 3 25 8 59 ...... // TestDC = // 124.64 // TestHF = // 980 // Time_In_second = // 81.178 新浪赛特居士SciteJushi-2014-03-27。 图片 1。用小波包基向量从变换树上捕捉DFT中负频率的影子 图片 2。用同样的处理检验在次最大尺度上小波包基向量
1125 次阅读|0 个评论
[转载][Repost]Explained: The Discrete Fourier Transform
qhhuangscut 2014-1-28 07:25
Source: http://web.mit.edu/newsoffice/2009/explained-fourier.html Larry Hardesty, MIT News Office Science and technology journalists pride themselves on the ability to explain complicated ideas in accessible ways, but there are some technical principles that we encounter so often in our reporting that paraphrasing them or writing around them begins to feel like missing a big part of the story. So in a new series of articles called Explained, MIT News Office staff will explain some of the core ideas in the areas they cover, as reference points for future reporting on MIT research. In 1811, Joseph Fourier, the 43-year-old prefect of the French district of Isère, entered a competition in heat research sponsored by the French Academy of Sciences. The paper he submitted described a novel analytical technique that we today call the Fourier transform, and it won the competition; but the prize jury declined to publish it, criticizing the sloppiness of Fourier’s reasoning. According to Jean-Pierre Kahane, a French mathematician and current member of the academy, as late as the early 1970s, Fourier’s name still didn’t turn up in the major French encyclopedia the Encyclopædia Universalis. Now, however, his name is everywhere. The Fourier transform is a way to decompose a signal into its constituent frequencies, and versions of it are used to generate and filter cell-phone and Wi-Fi transmissions, to compress audio, image, and video files so that they take up less bandwidth, and to solve differential equations, among other things. It’s so ubiquitous that “you don’t really study the Fourier transform for what it is,” says Laurent Demanet, an assistant professor of applied mathematics at MIT. “You take a class in signal processing, and there it is. You don’t have any choice.” The Fourier transform comes in three varieties: the plain old Fourier transform, the Fourier series, and the discrete Fourier transform. But it’s the discrete Fourier transform, or DFT, that accounts for the Fourier revival. In 1965, the computer scientists James Cooley and John Tukey described an algorithm called the fast Fourier transform, which made it much easier to calculate DFTs on a computer. All of a sudden, the DFT became a practical way to process digital signals. To get a sense of what the DFT does, consider an MP3 player plugged into a loudspeaker. The MP3 player sends the speaker audio information as fluctuations in the voltage of an electrical signal. Those fluctuations cause the speaker drum to vibrate, which in turn causes air particles to move, producing sound. An audio signal’s fluctuations over time can be depicted as a graph: the x-axis is time, and the y-axis is the voltage of the electrical signal, or perhaps the movement of the speaker drum or air particles. Either way, the signal ends up looking like an erratic wavelike squiggle. But when you listen to the sound produced from that squiggle, you can clearly distinguish all the instruments in a symphony orchestra, playing discrete notes at the same time. That’s because the erratic squiggle is, effectively, the sum of a number of much more regular squiggles, which represent different frequencies of sound. “Frequency” just means the rate at which air molecules go back and forth, or a voltage fluctuates, and it can be represented as the rate at which a regular squiggle goes up and down. When you add two frequencies together, the resulting squiggle goes up where both the component frequencies go up, goes down where they both go down, and does something in between where they’re going in different directions. The DFT does mathematically what the human ear does physically: decompose a signal into its component frequencies. Unlike the analog signal from, say, a record player, the digital signal from an MP3 player is just a series of numbers, each representing a point on a squiggle. Collect enough such points, and you produce a reasonable facsimile of a continuous signal: CD-quality digital audio recording, for instance, collects 44,100 samples a second. If you extract some number of consecutive values from a digital signal — 8, or 128, or 1,000 — the DFT represents them as the weighted sum of an equivalent number of frequencies. (“Weighted” just means that some of the frequencies count more than others toward the total.) The application of the DFT to wireless technologies is fairly straightforward: the ability to break a signal into its constituent frequencies lets cell-phone towers, for instance, disentangle transmissions from different users, allowing more of them to share the air. The application to data compression is less intuitive. But if you extract an eight-by-eight block of pixels from an image, each row or column is simply a sequence of eight numbers — like a digital signal with eight samples. The whole block can thus be represented as the weighted sum of 64 frequencies. If there’s little variation in color across the block, the weights of most of those frequencies will be zero or near zero. Throwing out the frequencies with low weights allows the block to be represented with fewer bits but little loss of fidelity. Demanet points out that the DFT has plenty of other applications, in areas like spectroscopy, magnetic resonance imaging, and quantum computing. But ultimately, he says, “It’s hard to explain what sort of impact Fourier’s had,” because the Fourier transform is such a fundamental concept that by now, “it’s part of the language.”
个人分类: 科研天地|2983 次阅读|0 个评论
[转载]Vasp计算的过程遇到的问题
plgongcat 2013-11-20 10:05
vasp计算的过程遇到的问题 (2012-12-14 10:59:37) 转载 ▼ 标签: 杂谈 分类: Vasp 最近在学vasp,这篇文章是百度文库找到的,看了不错,转载一把。 另外附上 vasp 程序 ,linux中下载后无须安装即可使用。单机中可能会出现内存溢出问题,可以放机群上使用。 01 、第一原理计算的一些心得 ( 1 )第一性原理其实是包括基于密度泛函的从头算和基于 Hartree-Fock 自洽计算的从头算,前者以电子密度作为基本变量(霍亨伯格 - 科洪定理),通过求解 Kohn-Sham 方程,迭代自洽得到体系的基态电子密度,然后求体系的基态性质;后者则通过自洽求解 Hartree-Fock 方程,获得体系的波函数,求基态性质; 评述: K-S 方程的计算水平达到了 H-F 水平,同时还考虑了电子间的交换关联作用。 ( 2 )关于 DFT 中密度泛函的 Functional ,其实是交换关联泛函 包括 LDA , GGA ,杂化泛函等等 一般 LDA 为局域密度近似,在空间某点用均匀电子气密度作为交换关联泛函的唯一变量,多数为参数化的 CA-PZ 方案; GGA 为广义梯度近似,不仅将电子密度作为交换关联泛函的变量,也考虑了密度的梯度为变量,包括 PBE,PW,RPBE 等方案, BLYP 泛函也属于 GGA ; 此外还有一些杂化泛函, B3LYP 等。 ( 3 )关于赝势 在处理计算体系中原子的电子态时,有两种方法,一种是考虑所有电子,叫做全电子法,比如 WIEN2K 中的 FLAPW 方法 ( 线性缀加平面波 ) ;此外还有一种方法是只考虑价电子,而把芯电子和原子核构成离子实放在一起考虑,即赝势法,一般赝势法是选取一个截断半径,截断半径以内,波函数变化较平滑,和真实的不同,截断半径以外则和真实情况相同,而且赝势法得到的能量本征值和全电子法应该相同。 赝势包括模守恒和超软,模守恒较硬,一般需要较大的截断能,超软势则可以用较小的截断能即可。另外,模守恒势的散射特性和全电子相同,因此一般红外,拉曼等光谱的计算需要用模守恒势。 赝势的测试标准应是赝势与全电子法计算结果的匹配度,而不是赝势与实验结果的匹配度,因为和实验结果的匹配可能是偶然的。 ( 4 )关于收敛测试 ( a ) Ecut ,也就是截断能,一般情况下,总能相对于不同 Ecut 做计算,当 Ecut 增大时总能变化不明显了即可;然而,在需要考虑体系应力时,还需对应力进行收敛测试,而且应力相对于 Ecut 的收敛要比总能更为苛刻,也就是某个截断能下总能已经收敛了,但应力未必收敛。 ( b ) K-point ,即 K 网格,一般金属需要较大的 K 网格,采用超晶胞时可以选用相对较小的 K 网格,但实际上还是要经过测试。 ( 5 )关于磁性 一般何时考虑自旋呢?举例子,例如 BaTiO3 中, Ba 、 Ti 和 O 分别为 +2 , +4 和 -2 价,离子全部为各个轨道满壳层的结构,就不必考虑自旋了;对于 BaMnO3 中,由于 Mn+3 价时 d 轨道还有电子,但未满,因此需考虑 Mn 的自旋,至于 Ba 和 O 则不必考虑。其实设定自旋就是给定一个原子磁矩的初始值,只在刚开始计算时作为初始值使用,具体的可参照磁性物理。 ( 6 )关于几何优化 包括很多种了,比如晶格常数和原子位置同时优化,只优化原子位置,只优化晶格常数,还有晶格常数和原子位置分开优化等等。 在 PRL 一篇文章中见到过只优化原子位置,晶格常数用实验值的例子( PRL100, 186402(2008) );也见到过晶格常数先优化,之后固定晶格常数优化原子位置的情况;更多的情况则是 Full geometry optimization 。 一般情况下,也有不优化几何结构直接计算电子结构的,但是对于缺陷形成能的计算则往往要优化。 ( 7 )关于软件 软件大致分为基于平面波的软件,如 CASTEP 、 PWSCF 和 ABINIT 等等,计算量大概和体系原子数目的三次方相关;还有基于原子轨道线性组合的软件 (LCAO) ,比如 openmx , siesta , dmol 等,计算量和体系原子数目相关,一般可模拟较多原子数目的体系。 VASP 是使用赝势和平面波基组,进行从头量子力学分子动力学计算的软件包,它基于 CASTEP1989 版开发。 VAMP/VASP 中的方法基于有限温度下的局域密度近似(用自由能作为变量)以及对每一 MD 步骤用有效矩阵对角方案和有效 Pulay 混合求解瞬时电子基态。这些技术可以避免原始的 Car-Parrinello 方法存在的一切问题,而后者是基于电子、离子运动方程同时积分的方法。离子和电子的相互作用超缓 Vanderbilt 赝势 (US-PP) 或投影扩充波 (PAW) 方法描述。两种技术都可以相当程度地减少过渡金属或第一行元素的每个原子所必需的平面波数量。力与张量可以用 VAMP/VASP 很容易地计算,用于把原子衰减到其瞬时基态中。 02 、 VASP 程序的亮点 : 1. VASP 使用 PAW 方法或超软赝势,因此基组尺寸非常小,描述体材料一般需要每原子不超过 100 个平面波,大多数情况下甚至每原子 50 个平面波就能得到可靠结果。 2. 在平面波程序中,某些部分代码的执行是三次标度。在 VASP 中,三次标度部分的前因子足可忽略,导致关于体系尺寸的高效标度。因此可以在实空间求解势的非局域贡献,并使正交化的次数最少。当体系具有大约 2000 个电子能带时,三次标度部分与其它部分可比,因此 VASP 可用于直到 4000 个价电子的体系。 3.VASP 使用传统的自洽场循环计算电子基态。这一方案与数值方法组合会实现有效、稳定、快速的 Kohn-Sham 方程自洽求解方案。程序使用的迭代矩阵对角化方案( RMM-DISS 和分块 Davidson )可能是目前最快的方案。 4. VASP 包含全功能的对称性代码,可以自动确定任意构型的对称性。 5. 对称性代码还用于设定 Monkhorst-Pack 特殊点,可以有效计算体材料和对称的团簇。 Brillouin 区的积分使用模糊方法或四面体方法。四面体方法可以用 Blöchl 校正去掉线性四面体方法的二次误差,实现更快的 k 点收敛速度。 03 、 VASP5.2 的新功能: 1. 大规模并行计算需要较少的内存。 2. 加入新的梯度校正泛函 AM05 和 PBEsol ;用标准 PBEPOTCAR 文件提供新泛函;改善了单中心处理。 3. 离子位置和格矢中加入有限差分,从而得到二阶导,用于计算原子间力常数和声子(需要超晶胞近似),和弹性常数。计算中自动考虑对称性。 4. 离子位置和静电场中加入线性响应,从而得到二阶导,用于计算原子间力常数和声子(需要超晶胞近似), Born 有效电荷张量,静态介电张量(电子和离子贡献),内应变张量,压电张量(电子和离子贡献)。线性响应只能用于局域和半局域泛函。 5. 精确的非局域交换和杂化泛函: Hartree-Fock 方法;杂化泛函,特别是 PBE0 和 HSE06 ;屏蔽交换;(实验性的)简单模型势 GW-COHSEX ,用于经验的屏蔽交换内核;(实验性的)杂化泛函 B3LYP 。 6. 通过本征态求和计算含频介电张量:使用粒子无关近似,或通过 GW 的随机相近似。可用于局域,半局域,杂化泛函,屏蔽交换,和 Hartree-Fock 。 7. 完全含频 GW ,速度达到等离子极点模型:单发 G0W0 ;在 G 和 W 中迭代本征矢直至自洽;(实验性的)迭代 G (也可以选 W )本征矢的自洽 GW ;(实验性的)对相关能使用 RPA 近似的 GW 总能量;用 LDA 计算 G 和 W 的顶点校正(局域场效应),仅能用于非自旋极化的情况;(实验性的) W 的多体顶点校正,仅能用于非自旋极化的情况。 8. 实验性的功能:用 TD-HF 和 TD- 杂化泛函求解 Cassida 方程(仅能用于非自旋极化的 Tamm-Dancoff 近似); GW 顶点的 Bethe-Salpeter (仅能用于非自旋极化的 Tamm-Dancoff 近似)。 1 、 V ASP 能够进行哪些过程的计算?怎样设置 ? 我们平时最常用的研究方法是做单点能计算,结构优化、从头计算的分子动力学和电子结构相关性质的计算。 一般我们的研究可以按照这样的过程来进行 如果要研究一个体系的最优化构型问题可以首先进行结构弛豫优化,然后对优化后的结构进行性质计算或者单点能计算。 如果要研究一个体系的热力学变化过程可以首先进行分子动力学过程模拟,然后在某个温度或压强下进行性质计算或者单点能计算。 如果要研究一个体系的热力学结构变化可以首先在初始温度下进行 NVT 计算,然后进行分子动力学退火,然后在结束温度下进行性质计算研究。 2 、 什么是单点能计算 (single pointenergy) ?如何计算? 跟其它软件类似, VASP 具有单点能计算的功能。也就是说,对一个给定的固定不变的结构(包括原子、分子、表面或体材料)能够计算其总能,即静态计算功能。 单点能计算需要的参数最少,最多只要在 KPOINTS 文件中设置一下合适的 K 点或者在 INCAR 文件中给定一个截断能 ENCUT 就可以了。还有一个参数就是电子步的收敛标准的设置 EDIFF ,默认值为 EDIFF=1E-4 ,一般不需要修改这个值。 具体来说要计算单点能,只要在 INCAR 中设置 IBRION=-1 也就是让离子不移动就可以了。 3 、 什么是结构优化 (structureoptimization) ?如何计算? 结构优化又叫结构弛豫( structurerelax ),是指通过对体系的坐标进行调整,使得其能量或内力达到最小的过程,与动力学退火不同,它是一种在0K下用原子间静力进行优化的方法。可以认为结构优化后的结构是相对稳定的基态结构,能够在实验之中获得的几率要大些(当然这只是理论计算的结果,必须由实验来验证)。 一般要做弛豫计算,需要设置弛豫收敛标准,也就是告诉系统收敛达成的判据( convergencebreak condition) ,当系统检测到能量变化减小到一个确定值时例如 EDIFFG=1E-3 时视为收敛中断计算,移动离子位置尝试进行下一步计算。 EDIFFG 这个值可以为负,例如 EDIFFG=-0.02 ,这时的收敛标准是当系统发现所有离子间作用力都小于给定的数值,如0 . 02 eV/A 时视为收敛而中断。 弛豫计算主要有两种方式:准牛顿方法 (quasi-Newton RMM-DIIS) 和共轭梯度法 (CG) 两种。准牛顿方法计算速度较快,适合于初始结构与平衡结构(势能面上全局最小值)比较接近的情况,而 CG 方法慢一些,找到全局最小的可能性也要大一些。选择方法为 IBRION=1 时为准牛顿方法而 IBRION=2 时为 CG 方法。 具体来说要做弛豫计算,设置 IBRION=1 或者 2 就可以了,其它参数根据需要来设置。 NSW 是进行弛豫的最大步数,例如设置 NSW=100 ,当计算在 100 步之内达到收敛时计算自动中断,而 100 步内没有达到收敛的话系统将在第 100 步后强制中止(平常计算步数不会超过 100 步,超过 100 步可能是计算的体系出了问题)。参数通常可以从文献中发现,例如收敛标准 EDIFFG 等。 有的时候我们需要一些带限制条件的弛豫计算,例如冻结部分原子、限制自旋的计算等等。冻结部分原子可以在 POSCAR 文件中设置 selectivedynamic 来实现。自旋多重度限制可以在 INCAR 中以 NUPDOWN 选项来设置。另外 ISIF 选项可以控制弛豫时的晶胞变化情况,例如晶胞的形状和体积等。 费米面附近能级电子分布的 smearing 是一种促进收敛的有效方法,可能产生物理意义不明确的分数占据态情况,不过问题不大。在 INCAR 文件中以 ISMEAR 来设置。一般来说 K 点只有一两个的时候采用 ISMEAR=0 ,金属体材料用 ISMEAR=1 或 2, 半导体材料用 ISMEAR=-5 等等。不过有时电子步收敛速度依然很慢,还需要设置一些算法控制选项,例如设置 ALGO=Very_Fast ,减小真空层厚度,减少 K 点数目等。 弛豫是一种非常有效的分析计算手段,虽然是静力学计算但是往往获得一些动力学得不到的结果。 4 、 vasp 的分子动力学模拟 vasp 做分子动力学的好处,由于 vasp 是近些年开发的比较成熟的软件,在做电子 scf 速度方面有较好的优势。缺点:可选系综太少。 尽管如此,对于大多数有关分子动力学的任务还是可以胜任的。主要使用的系综是 NVT 和 NVE 。 一般做分子动力学的时候都需要较多原子,一般都超过 100 个。 当原子数多的时候, k 点实际就需要较少了。有的时候用一个 k 点就行,不过这都需要严格的测试。通常超过 200 个原子的时候,用一个 k 点,即 Gamma 点就可以了。 INCAR : EDIFF 一般来说,用 1E-4 或者 1E-5 都可以,这个参数只是对第一个离子步的自洽影响大一些,对于长时间的分子动力学的模拟,精度小一点也无所谓,但不能太小。 IBRION=0 分子动力学模拟 IALGO=48 一般用 48 ,对于原子数较多,这个优化方式较好。 NSW=1000 多少个时间步长。 POTIM=3 时间步长 , 单位 fs, 通常 1 到 3. ISIF=2 计算外界的压力 . NBLOCK=1 多少个时间步长,写一次 CONTCAR , CHG 和 CHGCAR , PCDAT. KBLOCK=50 NBLOCK*KBLOCK 个步长写一次 XDATCAR. (个离子步写一次 PCDAT. ) ISMEAR=-1 费米迪拉克分布 . SIGMA =0.05 单位 : 电子伏 NELMIN=8 一般用 6 到 8, 最小的电子 scf 数 . 太少的话 , 收敛的不好 . LREAL=A APACO=10 径向分布函数距离 , 单位是埃 . NPACO=200 径向分布函数插的点数 . LCHARG=F 尽量不写电荷密度 , 否则 CHG 文件太大 . TEBEG=300 初始温度 . TEEND=300 终态温度。不设的话,等于 TEBEG. SMASS=-3 NVEensemble;-1 用来做模拟退火。大于 0NVT 系综。正确: SMASS=1,2,3 是没有区别的。都是 NVTensemble 。 SMASS 只要是大于 0 就是 NVT 系综。 CONTCAR 是每个离子步之后都会写出来的,但是会用新的把老的覆盖 CHG 是在每 10 个离子步写一次,不会覆盖 CHGCAR 是在任务正常结束之后才写的。 5 、收敛判据的选择 结构弛豫的判据一般有两中选择:能量和力。这两者是相关的,理想情况下,能量收敛到基态,力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大,力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的,能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。 到底是以能量为收敛判据,还是以力为收敛判据呢?关心能量的人,觉得以能量为判据就够了;关心力相关量的人,没有选择,只能用力作为收敛标准。对于超胞体系的结构优化,文献大部分采用 Gamma 点做单点优化。这个时候即使采用力为判据( EDIFFG=-0.02 ),在做静态自洽计算能量的时候,会发现,原本已经收敛得好好的力在不少敏感位置还是超过了结构优化时设置的标准。这个时候,是不是该怀疑对超胞仅做 Gamma 点结构优化的合理性呢?是不是要提高 K 点密度再做结构优化呢。 在我看来,这取决于所研究的问题的复杂程度。我们的计算从原胞开始,到超胞,到掺杂结构,到吸附结构,到反应和解离。每一步都在增加复杂程度。结构优化终点与初始结构是有关的,如果遇到对初始结构敏感的优化,那就头疼了。而且,还要注意到,催化反应不仅与原子本身及其化学环境有关,还会与几何构型有关。气固催化反应过程是电子的传递过程,也是分子拆分与重新组合的过程。如果优化终点的构型不同,可能会导致化学反应的途径上的差异。仅从这一点来看,第一性原理计算的复杂性,结果上的合理性判断都不是手册上写的那么简单。 对于涉及构型敏感性的结构优化过程,我觉得,以力作为收敛判据更合适。而且需要在 Gamma 点优化的基础上再提高 K 点密度继续优化,直到静态自洽计算时力达到收敛标准的。 6 、结构优化参数设置 结构优化,或者叫弛豫,是后续计算的基础。其收敛性受两个主要因素影响:初始结构的合理性和弛豫参数的设置 初始结构 初始结构包括原子堆积方式,和自旋、磁性、电荷、偶极等具有明确物理意义的模型相关参数。比如掺杂,表面吸附,空位等结构,初始原子的距离,角度等的设置需要有一定的经验积累。 DFT 计算短程强相互作用(相对于范德华力),如果初始距离设置过远(如超过 4 埃),则明显导致收敛很慢甚至得到不合理的结果。 比较好的设置方法可以参照键长。比如 CO 在 O 顶位的吸附,可以参照 CO2 中 C-O 键长来设置(如增长 20% )。也可以参照文献。记住一些常见键长,典型晶体中原子间距离等参数,有助于提高初始结构设置的合理性。实在不行,可以先在小体系上测试,然后再放到大体系中算。 弛豫参数 弛豫参数对收敛速度影响很大,这一点在计算工作没有全部铺开时可能不会觉察到有什么不妥,反正就给 NSW 设置个“无穷大”的数,最后总会有结果的。但是,时间是宝贵的,恰当的设置 3 小时就收敛的结果,不恰当的设置可能要一个白天加一个黑夜。如果你赶文章或者赶着毕业,你就知道这意味这什么。 结构优化分电子迭代和离子弛豫两个嵌套的过程。电子迭代自洽的速度,有四个响很大的因素:初始结构的合理性, k 点密度,是否考虑自旋和高斯展宽( SIGMA );离子弛豫的收敛速度,有三个很大的影响因素:弛豫方法( IBRION ) , 步长( POTIM )和收敛判据( EDIFFG ) . 一般来说,针对理论催化的计算,初始结构都是不太合理的。因此一开始采用很粗糙的优化( EDIFF=0.001 , EDIFFG=-0.2 ),很低的 K 点密度 (Gamma) ,不考虑自旋就可以了,这样 NSW60 的设置就比较好。其它参数可以默认。 经过第一轮优化,就可以进入下一步细致的优化了。就我的经验, EDIFF=1E-4,EDIFFG=-0.05 ,不考虑自旋, IBRION=2 ,其它默认, NSW=100; 跑完后可以设置 IBRION= 1 ,减小 OPTIM (默认为 0.5 ,可以设置 0.2 )继续优化。 优化的时候让它自己闷头跑是不对的,经常看看中间过程,根据情况调节优化参数是可以很好的提高优化速度。这个时候,提交两个以上的任务排队是好的方式,一个在调整的时候,下一个可以接着运行,不会因为停下当前任务导致机器空闲。 无论结构优化还是静态自洽,电子步的收敛也常常让新手头痛。如果电子步不能在 40 步内收敛,要么是参数设置的问题,要么是初始模型太糟糕(糟糕的不是一点点)。 静态自洽过程电子步不收敛一般是参数设置有问题。这个时候,改变迭代算法( ALGO ),提高高斯展宽( SIGMA 增加) , 设置自洽延迟( NELMDL )都是不错的方法。对于大体系比较难收敛的话,可以先调节 AMIN,BMIX 跑十多步,得到电荷密度和波函数,再重新计算。实在没办法了,可以先放任它跑 40 步,没有收敛的迹象的话,停下来,得到电荷密度和波函数后重新计算。一般都能在 40 步内收敛。 对于离子弛豫过程,不调节关系也不大。开始两个离子步可能要跑满 60 步(默认的),后面就会越来越快了。 总的说来,一般入门者,多看手册,多想多理解,多上机实践总结,比较容易提高到一个熟练操作工的水平。 如果要想做到“精确打击”,做到能在问题始发的时候就立刻采取有效措施来解决,就需要回归基础理论和计算方法上来了。 7 、优化结果对初始结构和“优化路径”的依赖 原子吸附问题不大,但是小分子吸附,存在初始构型上的差异。 slab 上水平放置,还是垂直放置,可能导致收敛结果上的差异。根据 H-K 理论,理想情况下,优化得到的应该是全局最小,但在数值计算的时候可能经常碰到不是全局最小的情况。实际操作中发现,多个不同初始结构优化收敛后在能量和结构上存在一定差异。 为了加快收敛速度,特别是对于表面 - 分子吸附结构,初始放松约束,比如 EDIFF=1E-3 , EDIFFG=-0.3,NSW=30 可能是很好的设置。但是下面的情况应当慎重: EDIFF=1E-3; EDIFFG=-0.1; !或者更小 NSW=500 ;!或者更大 电子步收敛约束较小,而离子步约束偏大,离子步数又很多,这种情况下,可能导致的结果是结构弛豫到严重未知的区间。 再在这个基础上提高约束来优化,可能就是徒劳的了——结果不可逆转的偏向不正常的区间。 好的做法,是对初始结构做比较松弛的约束,弛豫离子步 NSW 应该限制在一个较小的数值内。 EDIFF=1E-3 的话, EDIFFG 也最好是偏大一些,如 -0.3 而不是 -0.1. 这样可以在较少的步数内达到初步收敛。 对于远离基态的初始结构,一开始在非常松弛的约束下跑若干离子步,时间上带来的好处是很大的。对于 100 个原子的体系用 vasp 做 Gamma 点优化,如果一开始就是正常优化( EDIFF=1E-4,EDIFFG=-0.02 )设置,开始十个离子步可能都要花上几个小时。如果这个时候才发现输入文件有错误,那下午的时间就白费了,顺便带上晚上机器空转。 所以,我习惯的做法,是在初始几步优化后,会用 xcrysden 检查一下 XDATCAR 中的数据,用 xdat2xyz.pl 生成 movie.xyz ,然后看看弛豫过程是不是按照设想的那样。后续过程跑完一个收敛过程,就再检查一下 movie.xyz 。如此这般,才放心的展开后续计算。 8 、目的导向的结构优化 结构优化到这个阶段,是高级的了。为了得到特定结构,或者为了验证某些猜想,需要设计合理的初始结构,然后在这个基础上小心优化,比如 POTIM=0.1 跑几步看看,然后修改优化参数。 我遇到过的一件跟结构优化关系很大的算例是 CeO2 氧空位结构电子局域的问题 (http://emuch.net/bbs/viewthread.php?tid=2954558) 。按照一般方式(从优化好的 bulk 建 slab 模型,然后优化)得到一个 O 空位留下的两个电子均匀局域到 O 次外层三个 Ce 原子上,得到空位形成能 2.34eV. 经高人指点后,调节空位附近 O 原子位置,打破对称性后重新优化,两个电子完美的局域到两个 Ce 原子上了。并且空位形成能降低到 2.0XeV 。从这个例子可以看到,结构优化存在不少技巧的,这些技巧建立在研究者对模拟对象的物理意义的理解上。对物理图像的直观深入理解,才能做好模型预设,在此引导下才可能有目的的优化出不比寻常的结果。 目前第一性原理理论中的交换关联泛函部分包含经验参数。考虑这一点对优化结果的影响也很有意思。比如有专家提到, DFT+U 参数对某些结构的收敛终态构型有影响。构型的变化可能影响表面反应过程。基于这一点,一个好的计算研究可能就出来了。 真实过程总是复杂多变的。无论何种模拟,估计都可以找到一些试验现象来验证。但是到底应该如何评判模拟结果,如何从第一性原理研究中得出有意义的结论需要很好的洞察力。这样的模拟不见得就必须建立的试验的基础上,完全凭空设计的模型有可能更能优美的解释本质。 9 、在单机上计算态密度好像不会出问题。我先谈一下我的看法: 第一个 WARNING ,可以在 INCAR 文件中设置 NGX,NGY 和 NGZ 的值,设置的值要足够大,就可以消除这个 warning 。设置多大合适呢?这就要用到编译 vasp 时,同时也编译得到的 makeparam 小程序, makeparam 可以帮助你预先检查你设置的文件是否正确,以及某些参数的值是否合适。要得到合适的 NGX,NGY,NGZ 以及 NBANDS, 先在 INCAR 中不设置这些参数的值,然后运行 makeparamparam.inc ,其中 param.inc 是包含了输出结果的文件,在 param.inc 文件中你可以看到这些参数的值,以及计算大概需要多少的内存。然后把 param.inc 文件中的 NGX,NGY,NGZ 和 NBANDS 的值拷贝到 INCAR 文件中。 第二个是计算态密度时,我个人的做法是,一般把 KPOINTS 文件中的 k 点增多,然后把 INCAR 文件中的 ISTART=1,ICHARG=11 ,当然还设置 RWIGS 。最后把静止自洽计算得到的 CHG 和 CHGCAR 文件拷贝到当前目录下。从我在单机上的计算来看,没有 WAVECAR 文件也是可以计算态密度的。我想你出现的这个问题,可能是你 cluster 上计算时,每个节点上的 CHGCAR 和 WAVECAR 文件不一致造成的。 第三个是当 k 点数增加了,会出现一个 WARING ,要把此 WARNING 消失掉,在 INCAR 文件中设置 NELMDL ,它的值小于等于默认值(默认值好像是 -5, 你可以设为 -6) 。没有 cluster 的系统用来计算,也没有这样的经历,我仅从在单机上的计算经验来谈,有错还请包涵 10 、如何用 VASP 计算铁磁、反铁磁和顺磁 顺磁,意味进行 non-spin polarized 的计算,也就是 ISPIN=1 。 铁磁,意味进行 spin-polarized 的计算, ISPIN=2 ,而且每个磁性原子的初始磁矩设置为一样的值,也就是磁性原子的 MAGMOM 设置为一样的值。对非磁性原子也可以设置成一样的非零值(与磁性原子的一样)或零,最后收敛的结果,非磁性原子的 local 磁矩很小,快接近 0 ,很小的情况,很可能意味着真的是非磁性原子也会被极化而出现很小的 local 磁矩。 反铁磁,也意味着要进行 spin-polarized 的计算, ISPIN=2 ,这是需采用反铁磁的磁胞来进行计算,意味着此时计算所采用的晶胞不再是铁磁计算时的最小原胞。比如对铁晶体的铁磁状态,你可以采用 bcc 的原胞来计算,但是在进行反铁磁的 Fe 计算,这是你需要采用 sc 的结构来计算,计算的晶胞中包括两个原子,你要设置一个原子的 MAGMOM 为正的,另一个原子的 MAGMOM 设置为负,但是它们的绝对值一样。因此在进行反铁磁的计算时,应该确定好反铁磁的磁胞,以及磁序,要判断哪种磁序和磁胞是最可能的反铁磁状态,那只能是先做好各种可能的排列组合,然后分别计算这些可能组合的情况,最后比较它们的总能,总能最低的就是可能的磁序。同样也可以与它们同铁磁或顺磁的进行比较。了解到该材料究竟是铁磁的、还是顺磁或反铁磁的。 亚铁磁,也意味要进行 spin-polarized 的计算, ISPIN=2 ,与反铁磁的计算类似,不同的是原子正负磁矩的绝对值不是样大。非共线的磁性,那需采用专门的 non-collinear 的来进行计算,除了要设置 ISPIN , MAGMOM 的设置还需要指定每个原子在 x,y,z 方向上的大小。这种情况会复杂一些。 举个例子来说,对于 Mn-Cu(001)c(2x2) 这种体系,原胞里面有 2 个 Mn 原子,那么你直接让两个 Mn 原子的 MAGMOM 的绝对值一样,符号相反就可以了,再加上 ISPIN=2 。这样就可以实现进行反铁磁的计算了 11 、 vasp 在计算磁性的时候, oszicar 中得到的磁矩和 outcar 中得到各原子磁矩之和不一致在投稿的是否曾碰到有审稿人质疑,对于这个不一致你们一般是怎么解释的了? 答: OSZICAR 中得到的磁矩是 OUTCAR 中最后一步得到的总磁矩是相等的。总磁矩和各原子的磁矩 (RMT 球内的磁矩 ) 之和之差就是间隙区的磁矩。因为有间隙区存在 , 不一致是正常的。 12 、磁性计算应该比较负责。你应该还使用别的程序计算过磁性,与 vasp 结果比较是否一致,对磁性计算采用的程序有什么推荐。 ps :由于曾使用 vasp 和 dmol 算过非周期体系磁性,结构对磁性影响非常大,因此使用这两个程序计算的磁性要一致很麻烦。还不敢确定到底是哪个程序可能不可靠。 答:如果算磁性,全电子的结果更精确,我的一些计算结果显示磁性原子对在最近邻的位置时, PAW 与 FPLAW 给出的能量差不一致,在长程时符合的很好。虽然并没有改变定性结论。感觉 PAW 似乎不能很好地描述较强耦合。我试图在找出原因,主要使用 exciting 和 vasp 做比较。计算磁性推荐使用 FP-LAPW, FP-LMTO, FPLO 很吸引人 ( 不过是商业的 ) ,后者是 O(N) 算法。 13 、 vasp 学习笔记 POTCAR 的建立 POTCAR 将要告诉 vasp 计算的系统中所包含的各种元素的赝势 pesudopotential , vasp 本身就带有比较完善的赝势包,我们需要做的就是选择我们需要具体哪种赝势,然后把相应的文件拷贝形成我们具体的 POTCAR 文件。我们以 GaAs 为例。 1 )赝势的选择: vasp 的赝势文件放在目录 ~/vasp/potentials 下,可以看到该目录又包含五个子目录 potpot_GGA potpaw potpaw_GGA potpaw_PBE ,其中每一个子目录对应一种赝势形式。 赝势按产生方法可以分为 PP(standard pesudopotential ,其中大部分是 USPP,ultrasoft pesudopotential) 和 PAW(projector augmented wavemethod) 。按交换关联函数的不同又可以有 LDA(local density approximation) 和 GGA(generalized gradient approximation) ,其中 GGA 之下又可以再分为 PW91 和 PBE 。 以上各个目录对应起来分别是 pot== PP, LDA ; pot_GGA == PP, GGA ;potpaw == PAW, LDA ; potpaw_GGA ==PAW, GGA, PW91 ; potpaw_PBE == PAW , GGA,PBE 。选择某个目录进去,我们还会发现对应每种元素往往还会有多种赝势存在。这是因为根据对截断能量的选取不同还可以分为 Ga,Ga_s,Ga_h ,或者根据半芯态的不同还可以分为 Ga,Ga_sv,Ga_pv 的不同。 一般推荐选取 PAW_PBE 。其中各个元素具体推荐哪种形式的赝势可以参考 vaspworkshop 中有关赝势部分的 ppt 。当然自己能测试之后在选择是最好不过的了,以后再聊。 2 ) .POTCAR 的建立: 选好哪一种赝势之后,进入对应的目录,你会看到里边有这么几个文件, POTCAR.Z PSCTR.Z V_RHFIN.Z WS_FTP.LOG 。我们需要的是第一个。把它解压,如 zcatPOTCAR.Z Ga 。对 As 元素我们也可以类似得到一个 As 文件。用 cp 命令或者 mv 命令把这两个文件都移到我们的工作目录里。然后再用 cat 命令把这两个文件合并在一起,如 cat GaAs POTCAR ,这样就得到了我们需要的 POTCAR 。同理,有多个元素的 POTCAR 也可以这样产生。这里需要注意的是,记住元素的排列顺序,以后在 POSCAR 里各个元素的排列就是按着这里来的。 3 ) .POTCAR 里的信息: 如果你想看 POTCAR 长什么样,可以用 vimPOTCAR 命令,进去后可以用上下键移动光标。想出来的时候,可以敲入 :q! 就可以。具体的 vim 的命令可以在网上查到。一般我会看 POTCAR 里的截断能量为多大,用 grep-in enmax POTCAR 。 据说 B3LYP 的赝势计算比较准,我在 MS 上面测试过,好像 DOS 和能带图的计算确实比较准。不过不知道 vasp 有没有类似的赝势包。 hybridfunctional 的计算,并不需要特定的 hybrid functional 的赝势。大部分就是基于 GGA-PBE 的赝势来做,也就是芯电子与价电子的交换关联作用,以及芯电子与芯电子的交换关联作用还是基于 GGA-PBE 的,只是将价电子与价电子的交换关联作用通过 hybridfunctional 交换关联来描述。 谢谢老师的解答。那具体操作是不是像网上写的那样,使用 GGA 的赝势,设置 GGA =B3 ,然后更改 POTCAR 里面的 LEXCH=B3 就行了。我试过了,可以跑,不过结果没做详细的分析。 14 、 VASP 中所有能量的物理意义及它们之间的区别,让你彻底搞清楚 VASP 的所有能量 (一)首先我们应明白,固体的结合能就是固体的内能 E (结合) =U (内能), 原因如下: 一般情况都把孤立原子的能量作为能量参考点。前段时间有个同学问 VASP 中得出的绝对能量是相对于什么的,其实就是相对孤立原子得。 (二)其次我们根据自由能与内能之间的关系 F=U-TS 而且我们都知道 VASP 的所有计算都是在绝对 0 度下的情况, T=0 代入上式,有 F=U 。所以结合就等于内能等于自由能。肯定有 Freeenergy TOTEN=energy without entropy 恒成立 ... 这时候肯定有人会说不对啊,可以看 VASP 手册,候博的参考书作证,肯定不对得。 现在我告诉你确实它们二者确实有区别,区别在下面的情况 ( 1 )当我们用 ISMEAR=-5 时,费米能这儿没有展宽,它算出来的就是完全在绝对 0 度的能量。 Freeenergy TOTEN=energy without entropy 恒成立。 ( 2 )有时为了在数学上处理的方便,为了更容易积分,我们也用 ISMEAR ! =-5 (! = 是不等于的意思)的方法,这个时候费米能这儿有一定的展宽。此时,我们容易想到,有展宽不就是相当有一定的熵值吗?所以这个时候虽然算的是绝对 0 度的情况,但是有一定的熵值(我们应明白,这个熵值不是由一定的温度带来的,而是数学处理的结果)。所以在 SMEAR ! =-5 的方法我们会发现 Freeenergy TOTEN 和 energywithout entropy 有一定的差别。此时 energywithout entropy 是 Freeenergy TOTEN 在 SIGMA 趋于 0 的极限。 注意 : ( 1 )有人在算单个原子的能量时会发现单个原子的能量虽然很小但并不是 0, 但是按我上面的推导 , 固体中的结合能是相对孤立体系的能量而来的 , 所以单个原子得到的 TOTEN 肯定是 0 啊 , 原因在于我们的 POTCAR 不可能绝对合理 , 而且我们也知道计算单个原子的能量就是为了检测赝势 , 单原子得到的 TOTEN 越小说明赝势越好。但一般不会正好是 0. 对这个说法我还存在点疑问,写在了最后面。 ( 2 )如果你注意的话, energywithout entropy 与 Freeenergy TOTEN 在 SIGMA 趋于 0 也不是完全相等,但是也会发现它们之间的差别在 10E-3 左右,原因在于计算机求积分、求极限不能像我们人一样达到任意的精度。 15 、 VASP 中过渡态计算设置的一点体会 计算过渡态先要摆正心态,不急于下手。步骤如下: ( 1 )做模型,初态 IS 和终态 FS ,分别结构优化到基态; ( 2 )线形插入 images: nebmake.pl POSCAR.IS POSCAR.FSN N 为 image 个数。 ( 3 ) nebmovie.pl, 生成 movie.xyz 。用 Xcrysden --xyz movie.xyz 反复观看动画,仔细检查过程的合理性。这里要提醒, POSCAR.IS 和 POSCAR.FS 中原子坐标列表的顺序必须对应。 ( 4 )写 INCAR, 选 IOPT 。注意,最好忘记 vasp 自带的 NEB ,而全部改用包含 vtstool 的 vasp.IBRION=3,POTIM=0 关闭 vasp 自带的 NEB 功能。 ( 5 )过渡态计算第一个离子步最耗时,也最容易出问题,也是模型设计合理性检验的首要环节。所以可以选小一些的 ENCUT ,可以不用考虑自旋 (ISPIN=1) ,也不用考虑 DFT+U 。而且用最快最粗糙的算法( IOPT=3, 其他默认)。 ( 6 )带 vtstool 的 vasp-ClNEB(NEB) 过渡态计算 ICHAIN=0 作为入口,这个也是默认的。 LCLIMB=TRUE 也是默认的。如果不要 climbimage, 可以设置 LCLIMB= False. ( 7 )收敛判据 EDIFFG0 。过渡态计算要以力为收敛判据,而不是能量。一般 EDIFFG=-0.05 就可以接受, -0.02 或者 -0.01 更好。但是作为开始的过渡态计算,可以设置很宽的收敛条件,如 EDIFFG=-1. ( 8 )初步过渡态收敛后 , 修改 INCAR 中的优化器( IOPT ),并修改相应参数(参考 vtstool 官方论坛), EDIFFG 改小(如 -0.05 ),然后运行 vfin.pl ,这个脚本自动帮你准备在原来的基础上继续运行新的过渡态计算(完成 cp CONTCAR POSCAR, 保留电荷密度和波函数的操作)。 ( 9 )过渡态如何验算虚频呢? 比如一个 6 层原子层的 slab 上表面吸附小分子。 slab 下部 3 层原子是固定的。验算虚频的时候,是不是还是固定下面三层原子,然后按照一般频率计算方法来算虚频?这样的话,可以移动的原子数在 20 数量级上,考虑三个自由度,及其组合,就有很多很多可能了。请问该怎么设置这样的过渡态虚频计算呢? 16 、关于概念的问题做个讨 论 (一)关于结合能。你说 “ 结合能是定义为相距无穷远的原子结合形成一定结构的物质所放出的能量 ” 你和我说的没区别,我说的是结合能是相对于 “ 孤立原子做参考点的 ” ,也就是它与周围任何原子没有相互作用,和你所说的相距无穷远一回事,我这个好像没有任何错误。 =============== 这里你说的是没有错误,但是我觉得有必要先澄清一下。 (二)关于单点能。你说 “ 它是第一性原理计算直接得到的能量,或者说是赝能,是一个空间点阵平均每阵点上采用赝势计算所得到的能量,其中包含了结合能的贡献,但是更多的,也包含了靠近芯区附近的电子在采用赝势近似下的能量,这一部分能量既不是原子芯区附近电子能量的真实反应,也不会影响化学键性质,不会对结合能有所贡献 ” 。 我赞同你的大部分观点,也提出你说的几点错误,单点能准确的来说它包含了所有哈密顿的量,而且这儿的单点能不是你所说的 “ 平均每个原子的能量 ” ,而是你计算的整个原胞的能量。但是这个能量有一个参考点。你可以看候博得,也可以看我回的下一个贴子,至于 “ 影响不影响成键之类的内容 ” 固体力学上已经说的很清楚了。 ============== 从你的回复中,我可以知道你肯定没有学过晶体学或者空间群理论,你应该看看晶体学国际表中对于阵点的定义,阵点并不是每个原子,这里你的理解有问题,阵点是一个抽象点,一个晶体中包含所有对称性的可以仅通过平移来构造整个晶体的结构所占据的位置就是一个阵点,换句话说,一个阵点,就是一个满足平移对称性的原子集团,且该集团内部的位置满足该晶体结构的全部对称性,而且它不仅仅是 “ 原胞 ” (三)你说 “ FreeenergyTOTEN 是体系总能,要减去阵点上分布的原子的能量再除以平均原子数才是结合能(当然,这个和你的计算脚本的设计有关),而且这还没有考虑不加展宽时没有被计算到的能带的因素 ” 我不赞同你后面说的几点。 FreeenergyTOTEN 从字面意思上我们也知道它的结果是自由能,你可以说它是总能,因为根据我上面的推导,它们至少在数值是相等得。不在于你把它说成什么,你就是把它说成总能,其实它还是等于结合能,等于自由能,等于内能。至于除不除原子总数在于你想得到的是平均每个原子的还是总体系的,这在于个人。考虑不考虑展宽,那要看 ISMEAR 等于几,做几个实例就会感觉到它考虑没考虑了。 =============== 这个能量确切的说应该是叫做考虑电子振动熵的体系总自由能,当不考虑展宽的时候,它是等于总能的,如果你读过 Vasp 的代码,就知道 TOTEN 在 vasp 的计算中就是总能,这个和结合能不是一个概念,还包含有非成键部分的贡献,至于内能的定义,如果你阅读过塞兹的现代固体理论,或者 Pauling 的书,或者读过 Morse 当时提出 morse 势的那篇文献,就应该知道,固体物理中所使用的内能,指的是离子实的动能和原子的 “ 结合能 ” 之和 —— 这里结合能之所以要打引号,是因为按定义,是要形成稳定结构或者亚稳态结构时才能称之为结合能,内能定义并没有考虑粒子芯区附近电子能量的影响,正如你所说,是 “ 相对于 “ 孤立原子做参考点的 ”” ,在 0K 下已经不考虑动能,因此就应是总能减去孤立原子的能量和才行,至于结合能,则是稳定状态下结构的这个能量。 (四)你说 “ 是否考虑展宽和结合能的定义没有关系 ” 我也没说和它的定义有什么关系啊,但是由于数学处理带来的误差,它对结合能的结果有一定的影响啊。 (五)对单个原子的结合能的计算应该只计算 Gamma 点的能量,且用削除简并。 你说得第五点我不懂,我算单个原子的能量时一直是按三个表面的构造方法来算的,也没想过什么简并。还望有高手给我帮助,第五点怎么理解。 ============= 简单的操作是计算单个原子能量只考虑 Gamma 点,然后三边都设置在 10A 以上,且不相等 至于原因,应该去查量子力学的书,记得本科的时候,老师都会讲到的。 根据 F=U-TS , E (结合) =E ( bulk ) -nE( 单个离散 ) 。而 E ( bulk )就是 U ,通常我们取 E (离散)为参考点。也就是把 E (离散)看为 0 ,这样推出来,在 T=0 时, F=E (结合)。 它是包含你所说得那些,而且就像我在前面的贴子中说得,它就是总得 H 得到的能量,但是它有一个参考,而这个参考就是离散原子的能量 17 、用 vasp 软件研究表面小分子的吸附解离遇到几个问题。 1 关于反应物和产物的结构,请问你们是如何构建的?(参考文献吗?)能不能介绍一下相关的经验,个人认为好的开端是成功的一半,所以需要更多的经验帮助。 答: . 据说现在流行的过渡态算法多用 CINEB , vasp 自带的 NEB 也可以,但是镜像受力优化和势垒精确度以及收敛速度略逊于 CINEB , lz 可以 google 到其官方网站了解。初始和终了的结构是局部能量最低的构型( vasp 做弛豫收敛的构型),至于 lz 所关心的是从具体结构之间的过渡的化,一是猜测,二是文献吧。基本出发点是能量趋向降低的构型。 2 关于结构的优化,能量收敛和力收敛的标准是否应该比默认值高,有利于过渡态的搜寻? 答: DFT 理论计算能量可能更准确,对力的计算由于是对能量再求导,涉及数值算法原因可能不准,所以收敛标准基本采用默认即可,有时适当提高(一个量级以内)也是可以的。在合理的精度范围内讨论出合理的结果就达到目的了,这是师兄告我的。基本上收敛精度高于默认值去跑的话,那就 god bless you 了。 3 关于结构虚频的处理,通过学习了解了消除虚频的方法 -- 将对应原子的坐标加到原坐标上。由于 vasp 的振动模式不能可视化,所以请问如何找到对应的原子?另外关于虚频,大家是如何处理的呢。 答: . 虚频没有关注过。只知道在鞍点位置才应该有一个虚频,属于鞍点具有的特性吧。为何要消除呢?根据这个好像才能按速率理论算出反应率的吧,这个我没有经验。只关注过势垒。 4 关于 images 的问题,最近利用脚本产生了 2 、 4 的 POSCAR 的文件,发现只有 2 的能跑起来,我的服务器是 4 个核的。看别人发的贴:指定 4 个 cpu 同时计算,应该是 1 个 cpu 一个 image 。为什么我设置 4 个 image 跑不起来呢? 答: image 数量越多,相互之间受力需要同时满足收敛准则,速度肯定会降低,但是这样找到的能量最低路径也会与实际更加符合。没错,每个 image 事实上是需要相当于单独计算的构型的 cpu 的数目的,所以 4 个核的话,平均一个 image 一个 cpu ,实在是不一定跑得动。如果体系原子数较多, k 点网格较密,能量截断值再一高,估计经过很久才可能有结果。 5 、在计算 dos 时,为什么会出现 WARNING:stress and forces are notcorrect 。这是那些原因造成的,计算的是多铁物质。还有就是在优化时, TOTAL-FORCE 这项怎么全是 0 啊,原子这个没有关系的 答:用优化后的 CONTCAR 算 band 或 dos 都会出现这个提示位置没有变化啊,这是怎么回事啊? 6 、我在做能带计算,发现自洽计算和非自洽计算中, OUTCAR 文件中都可以找到一个费米能级的大小,但是现在迷茫的是,不知道在画能带时减去的费米能级的大小应该取哪个结果文件中得值?是取自洽计算的结果文件中的 E-fermi ,还是取非自洽的结果文件中的 E-fermi 的大小?? 答:算能带的那次计算,指的是自洽(求解 WAVECAR 和 CHGCAR 的那一步)的那步计算,还是非自洽(直接读取自洽的电荷和波函数,进行非自洽计算)的那次计算啊?我对这个说法的标准不是很理解,希望您能再指导一下,我觉得画能带图时,费米能级的取值是否正确是很关键的,取自洽计算得到的 E-fermi 转自:http://blog.sina.com.cn/s/blog_b364ab230101e9dp.html
个人分类: 老生经验之谈|4976 次阅读|0 个评论
DFT+U理论和参数U的确定
热度 4 ChemiAndy 2012-11-29 23:19
最近在做DFT+U计算,有些收获分享一下。本人数学功底先天不足,难免有理解错误的地方,欢迎指正。 首先是对+U必要性的感悟。DFT(LDA和GGA)对一般体系的结构、能带计算还是很令人满意的。这些一般体系是指前3周期的体系和纯金属体系。但是,对于包含d电子和f电子的体系,特别是过渡金属氧化物或氮化物,在“金属/绝缘体”的定性判断上常常出错。LDA和GGA往往把一些绝缘体/半导体的gap算的太小,甚至最高占据轨道/能带(HOMO)在Fermi面之上,即成金属了。这些最高占据轨道往往是来自这些金属原子的d电子和f电子。那么为什么DFT处理d/f电子会出现较大误差呢?我认为这些d/f电子的特殊性在于它们不仅离原子核较远,轨道空间伸展形状奇特(比如穿透效应,即穿透更加外层的s,p轨道)而且多数会发生自旋极化,即原应自旋相反的电子对发生自旋跃迁,导致金属原子神奇地出现了净的磁矩。如果单胞内所有金属原子的净磁矩方向相同,则单胞总体产生净磁矩,材料呈现铁磁性;如果单胞内相邻金属原子的净磁矩相反,则单胞没有总体净磁矩,但存在一上一下的磁矩排列,材料呈现为反铁磁性。正是d/f电子这些神奇的特性,使得DFT计算的过渡金属氧化物或氮化物能带结构频频失误。为什么呢?原因在于DFT计算电子-电子之间的排斥作用所用的交换-相关泛函,是基于单粒子近似发展起来的。单粒子近似是把自旋配对的电子对看成单个粒子,这样Kohn-Sham轨道数量为总电子数的一半。但是当处理d/f电子时,不得不把电子对分开为自旋向上和自旋向下的两个粒子,即开壳层计算open shell。对开壳层体系使用单粒子近似下的交换-相关泛函时,其中的相关泛函就不能完全考虑d/f电子的相关效应。d/f电子的相关效应,就是自旋电子与自旋电子之间的排斥,和自旋电子与内层电子之间的排斥。这些排斥作用使得轨道/能带之间分割较远,轨道/能带较窄。但是DFT的相关泛函对这些排斥考虑的不够,结果轨道/能带过宽,轨道与轨道相互接近甚至重叠。这就是为什么DFT把绝缘体计算成了金属。一般把这些过度金属半导体/绝缘体体系称为强相关体系strongly correlated system。 怎么办呢?可能你已经直觉到了,直接加上一项,把d/f轨道上的自旋电子间的排斥作用作为(u/2*int(n_up*n_down))一个能量项加进去。这个就是Hubbard U model: E_DFT+U = E_DFT + U/2*int(n_up*n_down) 其中第一项E_DFT仍为正常的DFT计算。n_up, n_down分别代表自旋向上和自旋向下的电荷密度。二者相乘表示了它们之间的排斥作用。完整的Hubbard model还要考虑其中的重复计算,减去一个项,称为J项。+J项的计算刚刚在QE中实现。如果你追求结果的定量准确性,可以尝试玩玩+J。 这个+U项排斥作用前面有个系数U/2,也就是说这个排斥能并非直接的库仑排斥作用,而是有个系数。这个系数并非一个经验参数,而是代表了能量对d/f轨道上电子密度变化的二级导系数【1】。你可能会困惑,这有点像GGA?是的,不过GGA是通用的交换-相关能对密度的二阶展开的泛函。而这里的U是纯DFT基态总能量对特定体系d/f电子密度变化的二阶导。这有点像是个二次方项的收缩函数,即你们这些自旋作用的电子在轨道上的占据必须往特定方向收缩,否则会有能量惩罚。 上图中,不加U时左边的最高占据轨道和右边的空轨道之间的gap很窄,是不正确的;+U后变宽,更加接近实验结果。需要注意的是,+U是对电子结构的改善,对原子几何、晶体结构的影响可能很有限。此时只要不涉及电子性质,一般的GGA, LDA的结果是可以接受的。 那么这个+U能量惩罚项是人为的经验项,还是真实的物理存在呢?Cococcioni曾在【1】 PRB71,035105 中试图用图1来说明Hubbard U能量项的本质。他说真实的体系与环境交换电子,其能量变化应该是交换前后的两个状态的线性组合,这样的能量变化应该是线性的。但LDA和GGA对部分(非整数的)占据的Kohn-Sham轨道的自作用的不正确处理导致能量变化是非线性的。这二者之间的能量差就是Hubbard U能量项的来源。从这个解释可以看出, (1) +U能量惩罚项是物理存在的,并非经验校正。参数U的确定也不应该是经验的,而是必须考察DFT与真实势能面的差。 (2) 这个二级导是体系环境变化而变化的。这就导致必须对每个体系,甚至体系的不同状态(比如键长变化)计算相应的参数U。 注意如果体系内有对称性不相同的金属位点,需要对每个位点上的金属原子计算其相应的参数U,不管它们是否为同种金属! 获得合适的参数U以后,才能使用DFT+U模型精确确定自旋极化和能带、带隙等问题。(当然,更精确的模型还要考虑旋轨耦合spin-orbital coupling。DFT+U并不包括旋轨耦合效应) 怎么确定某个体系参数U呢?因为U是纯DFT基态总能量对这个体系d/f电子密度变化的二阶导,那么很自然的一个想法,我们对密度加以扰动,求能量的变化(响应),可以得到数值二阶导(微扰法)。不过,密度并非一个理想的微扰项,因为一个金属原子(Hubbard位点)涉及到多个电子。怎么办呢?目前我看到的文献中共有2种策略。 第一种是Cococcioni的线性响应法。它向单个Hubbard位点施加一个势alpha,引起该位点的密度(占据数)重新分配,然后看密度对此alpha势变化的响应。此时响应有两种情况,1是纯DFT能量对势alpha的响应(不+U, 即考虑相关效应前)称为bare response(X0);2是+U能量对加势和占据数变化后总的响应,称为self consistent response(X)。二者的线性响应系数的倒数之差即为U: U = 1/X0 - 1/X 此法原始文献即PRB71. 上图中,红色为bare response,其线性系数为X0=-0.32,数据从自洽前第一个循环计算得到的Hubbard site上(即带d/f电子的金属原子)的占据数得到;蓝色为self-consistent response, 线性系数为X=-0.15,数据从自洽完成后Hubbard site的占据数得到. U = 1/X0 - 1/X = 3.47eV。 (图片来源:教程2) 在这个教程【3】里有很详细的例子: 点击打开并下翻到 July, 23 Tursday: LSDA, collinear and noncollinear magnetism; spin-orbit coupling; DFT+U 在其中的Labs 1. Matteo Cococcioni,Examples of spin-polarized and DFT+U calculations 注意上面这种方法基于一个假定,即势alpha足够小,此时得到的U近似等于alpha=0时的U,即U_0。这是线性响应理论。即体系能量实际对势alpha的响应并不是线性的,但我们假定扰动足够小的情况下假定这种响应是线性的。上述方法对单个Hubard位点施加势alpha,但是在周期性体系中,任何扰动也都是周期性的。怎么办呢? Cococcioni提出的方法是延伸单胞为超胞,做响应计算;然后再次延伸为更大超胞,做响应计算。。。这样随着盒子越来越大,势在其中的作用越来越小,直到盒子无限大,则响应系数收敛。这种方法并不需要计算实际超胞,而是对密度响应进行外推。 Marzari提供了另一种不使用超胞外推即可获得无扰动的U_0的自洽方法(self-consistent U)(【2】 原始文献PRL97,103001(2006) )注意到尽管我们需要的是纯DFT能量对alpha的响应,但是上面的所有计算其实都是按照DFT+U的模型计算,只不过,进行alpha势扰动的时候,设置U为一个很小的值(称为初始U值,U_in),以获得纯DFT基态能对alpha的响应,从而得到U(称为U_out)。 Marzari并不让初始U_in近似为0, 而是设置它等于一系列值(0.5, 1.0, 1.5, 2.0, 2.5, 3.0 ...),然后运用上述方法去计算相应的U_out,然后将上述(U_in, U_out)数据系列在线性相关区拟合直线,可外推到U_in=0时的U_scf。如下图所示。 上图中,对给定的U_in = 0.5, 1.0, 1.5, 2.0, 2.5, 3.0做DFT+U计算,结束后可根据Hubbard位点上占据数得到一个U_out = Uscf - U_in / m。然后外推U_in=0时,U_scf=3.56即为正确的参数U。其中的U_0是使用第一种方法在势alpha扰动近似为0得到的近似U_0值,并非alpha=0的U值。(图片来源:教程2) 第一种超胞法(称为linear response法)外推和U_in/U_out外推法得到的参数U应该是一致的,但第二种方法更加优雅。 这个方法在这里有很好的教程【4】: Calculating the Hubbard U 其中对两种方法都有相应的例子和Python脚本。注意两种方法的切换要修改其中的Variable.py文件,设定hubbylist或者alphalist。其中所用的例子是6重态MnO的cluster,并非固体。 虽然参数U的确定很麻烦,但是还是值得的。从PRL 97, 103001 (2006)中可以看到,GGA+U后算FeO的平衡键长、简谐振动频率,非谐频率的结果,直接逼近CCSD(T),当时哥就震惊了。 在使用+U model的过程中,我一直在思考这个U参数的本质,即U大一点意味着什么,小一点意味着什么?我们想,既然它是DFT能量能量对电荷密度变化的二阶导,那么有点类似与振动频率。在分子中,如果一根键的强度很高,那么振动频率就很高;如果键强度低,振动就弱,比如C=C双键就比C-C单键的振动频率高。类比的话,应该说U代表了d/f电荷密度随外部势(比如其它原子核对它的作用)的变化而变化的难易程度。如果U大,那么不容易变化;如果U小,则容易变化。在不使用+U模型时,那就是说U=0,这时d/f电子最容易受外部势的影响而左右摇摆,从宏观上来看,就是分布比较均匀,能带很宽,而能带之间的带系带隙变宽。一旦+U以后,立即受到束缚,从而能带变窄而带隙变宽。这和+U的实际作用正好一致。 进一步,Hubbard原子上的U参数与材料的金属性-非金属性的有一种联系:材料的金属性越强,U越小;材料非金属性越强,U越大。这从Cococcioni博士论文中的数据可以看出来:在金属态中,铁Fe的U为2.2, 铁磁FeO中为4.29, 而绝缘体Fe2SiO4中为 4.9/4.6. 所以,如果你得到的U很小,很可能意味着你的体系的金属性很强。当然,这并不意味着你可以不使用U,U虽小仍然可能大大改善能带的分布。 【1】参数U确定的线性响应法:PRB 71, 035105 (2005) 【2】参数U确定的In-Out自洽法:PRL 97, 103001 (2006) 【3】教程1:UCSB09, Matteo Cococcioni: DFT+U calculations : Lecture slides , examples , instructuction to examples 【4】教程2: Heather Kulik: Calculating the Hubbard U 原发于小木虫: http://emuch.net/bbs/viewthread.php?tid=5237955 后记: 强相关体系材料的能带和band gap问题,北大蒋洪教授有篇2012年的《化学进展》讲的非常系统,深刻:《 带隙问题:第一原理电子能带理论研究现状》,化学进展,Vol.24 No.6, pp910. 看了后茅塞顿开,受益良多。 蒋老师指明LDA/GGA带隙问题的原因是 标准方法对长程相关的考虑不足,或者说是电子密度的自作用误差导致的。 解决带隙问题的两大主流方法是屏蔽杂化泛函和格林函数方法。蒋老师 在一个讲座视频中称+U方法只对d电子的相关效应有较好修正。推荐阅读。
个人分类: 科普文章|20816 次阅读|5 个评论
[转载]VASP包含范德瓦尔斯力的计算
dwd0826 2012-9-7 15:31
转自: http://cms.mpi.univie.ac.at/vasp/vasp/vdW_DF_functional_Langreth_Lundqvist_et_al.html vdW-DF functional of Langreth and Lundqvist et al.The vdW-DF proposed by Dion et al. is a non-local correlation functional that approximately accounts for dispersion interactions . In VASP the method is implemented using the algorithm of Roman-Perez and Soler which transforms the double real space integral to reciprocal space and reduces the computational effort. Several propsed versions of the method can be used : the original vdW-DF , vdW-DF with exchange functionals optimised for the correlation part , and the vdW-DF2 of Langreth and Lundqvist groups . N.B. : This feature has been implemented by J. Klimeš. If you make use of the vdW-DF functionals presented in this section, we ask that you cite the following paper: J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83 , 195131 (2011). Correlation functionals The method is invoked by setting LUSE_VDW = .TRUE. Moreover, the PBE correlation correction needs to be removed since only LDA correlation is used in the functionals. This is done by setting AGGAC = 0.0000 The two tags above need to be used for all the following functionals. Exchange functionals To use the different exchange functionals, the GGA tag needs to be set appropriately. The original version of Dion et al uses revPBE which can be set by GGA = RE More accurate exchange functionals for the vdW correlation functional have been proposed in and . They can be used by setting GGA = OR for optPBE, GGA = BO PARAM1 = 0.1833333333 PARAM2 = 0.2200000000 for the optB88 functional , or GGA = MK PARAM1 = 0.1234 PARAM2 = 1.0000 for the optB86b functional . For the vdW-DF2 functional the rPW86 exchange functional is used: GGA = ML moreover, the vdW functional needs to be changed to the vdW2 correlation which requires only a change of a parameter: Zab_vdW = -1.8867 An overview of the performance of the different approaches can be found for example in for gas phase clusters and in for solids. Important remarks : The method needs a precalculated kernel which is distributed via the VASP download portal ( VASP - src - vdw_kernel.bindat ) and on the ftp server ( vasp5/src/vdw_kernel.bindat ). If VASP does not find this file, the kernel will be calculated. This, however, is rather demanding calculation. The kernel needs to be either copied to the VASP run directory for each calculation or can be stored in a central location and read from there. The location needs to be set in routine PHI_GENERATE. This does not work on some clusters and the kernel needs to be copied into the run directory in such cases. Currently the evaluation of the vdW energy term is not done fully within the PAW method but the sum of the pseudo-valence density and partial core density is used. This approximation works rather well, as is discussed in , and the accuracy generally increases when the number of valence electrons is increased or when harder PAW datasets are used . For example, for adsorption it is recommended to compare the adsorption energy obtained with standard PAW datasets and more-electron POTCARs for both PBE calculation and vdW-DF calculation to asses the quality of the results. The optimisation of the cell (ISIF=3 and higher) is currently not possible. The spin polarised calculations are possible, but strictly speaking the non-local vdW correlation is not defined for spin-polarised systems. For spin-polarised calculation the non-local vdW correlation energy is evaluated on the sum of the spin-up and spin-down densities. The evaluation of the vdW energy requires some additional time. Most of it is spent on performing FFTs to evaluate the energy and potential. Thus the additional time is determined by the number of FFT grid points in the calculation, basically size of the simulation cell. It is almost independent on the number of the atoms in the cell. Thus the relative cost of the vdW-DF method depends on the ``filling" of the cell and increases with the amount of vacuum in the cell. The relative increase is high for isolated molecules in large cells, but small for solids in smaller cells with many k-points. This feature has been implemented by J. Klimeš. If you make use of the vdW-DF functionals presented in this section, we ask that you cite the following paper: J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83 , 195131 (2011). http://cms.mpi.univie.ac.at/vasp/vasp/DFT_D2_method_Grimme.html#tab:grimme DFT-D2 method of Grimme LVDW= .TRUE. | .FALSE. (Available as of VASP.5.2.11) Default: LVDW=.FALSE. Popular density functionals are unable to describe correctly van der Waals interactions resulting from dynamical correlations between fluctuating charge distributions. A pragmatic method to work around this problem has been given by the DFT-D approach , which consists in adding a semi-empirical dispersion potential to the conventional Kohn-Sham DFT energy: (90) In the DFT-D2 method of Grimme , the van der Waals interactions are described via a simple pair-wise force field, which is optimized for several popular DFT functionals. The dispersion energy for periodic systems is defined as (91) where the summations are over all atoms and all translations of the unit cell , the prime indicates that for , is a global scaling factor, denotes the dispersion coefficient for the atom pair , is a position vector of atom after performing translations of the unit cell along lattice vectors. In practice, terms corresponding to interactions over distances longer than a certain suitably chosen cutoff radius contribute only negligibly to and can be ignored. The term is a damping function (92) whose role is to scale the force field such as to minimize contributions from interactions within typical bonding distances. Combination rules for dispersion coefficients and vdW radii are (93) and (94) respectively. The global scaling parameter has been optimized for several different DFT functionals such as PBE ( ), BLYP ( ), and B3LYP ( ). The DFT-D2 method can be activated by setting LVDW=.TRUE. Optionally, the forcefield parameters can be controlled using the following flags (the default values are listed): VDW_RADIUS = 30.0 cutoff radius () for pair interactions VDW_SCALING = 0.75 global scaling factor VDW_D = 20.0 damping parameter VDW_C6 = ,... parameters ( ) for each species defined in POSCAR VDW_R0 = ,... parameters () for each species defined in POSCAR The default values for VDW_C6 and VDW_R0 are compiled in Tab. 2 . As the potential energy, interatomic forces as well as stress tensor are corrected by adding contribution from the forcefield, simulations such as the atomic and lattice relaxations, molecular dynamics, and vibrational analysis can be performed. The number of atomic pairs contributing to and the estimated vdW energy are written in OUTCAR (check lines following the expression 'Grimme's potential'). The forces and stresses written in OUTCAR contain the vdW correction but the corrected energy should be read from OSZICAR (energies in OUTCAR do not contain the vdW term). IMPORTANT NOTE: The defaults for VDW_C6 and VDW_R0 are defined only for the first five rows of periodic table of elements (see Tab. 2 ) - if the system contains other elements the user must provide the corresponding parameters. Table 2: Parameters used in the empirical force-field of Grimme . Element C R Element C R Jnm mol Jnm mol H 0.14 1.001 K 10.80 1.485 He 0.08 1.012 Ca 10.80 1.474 Li 1.61 0.825 Sc-Zn 10.80 1.562 Be 1.61 1.408 Ga 16.99 1.650 B 3.13 1.485 Ge 17.10 1.727 C 1.75 1.452 As 16.37 1.760 N 1.23 1.397 Se 12.64 1.771 O 0.70 1.342 Br 12.47 1.749 F 0.75 1.287 Kr 12.01 1.727 Ne 0.63 1.243 Rb 24.67 1.628 Na 5.71 1.144 Sr 24.67 1.606 Mg 5.71 1.364 Y-Cd 24.67 1.639 Al 10.79 1.716 In 37.32 1.672 Si 9.23 1.716 Sn 38.71 1.804 P 7.84 1.705 Sb 38.44 1.881 S 5.57 1.683 Te 31.74 1.892 Cl 5.07 1.639 I 31.50 1.892 Ar 4.61 1.595 Xe 29.99 1.881
13871 次阅读|0 个评论
DFT计算新算法:编织最不寻常的电子舞
热度 1 GoogleMIT 2012-7-18 20:18
APS的观点评论最近刊出了一篇关于 DFT 计算的viewpoint,作者是 Neepa T. Maitra , Department of Physics and Astronomy, Hunter College and the City University of New York, 695 Park Avenue, New York, NY 10065, USA。 DFT (density-functional theory)方法是计算材料学和量子化学中计算电子结构的一种重要方法。它最初于1964年被提出,已经被应用了许多年,相关的文献和研究成果恐怕已经无法计数 。而 TD-DFT (time dependentdensity-functional theory)是DFT的扩展,主要用于处理非平衡态的电子与受外场束缚的电子 ,但是什么才是正确的方法通过“弯曲”有吸引力的原子核位势去描述不相互作用的电子(或者用文中生动的方法,叫做编排电子的舞蹈),这些电子恰恰又正好包含了真实的时间相关系统的准确信息。这个问题困扰了很久。然而限制看来,除了少数双电子系统由于计算水平限制之外,其他大部分系统的问题恐怕都可以得到解决。最近的 Phys. Rev. Lett 上,NYU的 James Ramsden 等人提出了一种新的算法来解决这个问题,其实也不是很陌生,它就是含时的Kohn-Sham potential,他们已经将此种方法应用到含有运动电子的半导体材料中 Exact Density-Functional Potentials for Time-Dependent Quasiparticles.pdf ,之前TDDFT忽略了一些 Kohn-Sham potential中重要的细节。 原文链接: http://physics.aps.org/articles/v5/79 References: P. Hohenberg and W. Kohn, “Inhomogeneous Electron Gas,” Phys. Rev.136, B864 (1964) . W. Kohn and L. J. Sham, “Self-Consistent Equations Including Exchange and Correlation Effects,” Phys. Rev.140, A1133 (1965) . E. Runge and E. K. U. Gross, “Density-Functional Theory for Time-Dependent Systems,” Phys. Rev. Lett.52, 997 (1984) . Fundamentals of Time-Dependent Density Functional Theory, edited by M. A. L. Marques, N. T. Maitra, F. Nogueira, E. K. U. Gross, and A. Rubio, Lecture Notes in Physics Vol. 837 (Springer, Berlin, 2012) . J. D. Ramsden and R. W. Godby, ”Exact Density-Functional Potentials for Time-Dependent Quasiparticles,” Phys. Rev. Lett.109, 036402 (2012) .
个人分类: 论文评论等|794 次阅读|2 个评论
Zn真的有+3价的氧化态吗?--评一篇JACS文章
热度 22 ChemiAndy 2012-6-29 13:05
请先看: 胡新根:你相信金属锌(Zn)有+3价氧化态吗? 卤素,F, Cl, Br, I具有强烈的吸电子能力。其吸电子能力可以通过计算它们的电子亲和能得到,即X吸收一个电子形成X-所需的能量。Cl具有卤素中最高的电子亲和能,3.6eV。但是,你能不能设计一个具有更高电子亲和能的物质?这种物质对工业中所需要的强氧化过程非常有用。 可以。BO2和AuF6就是这样的具有更高电子亲和能的物质,称为superhalogen,“高级卤素”,其形成BO2(-)和AuF6(-)阴离子的电子亲和能分别为4.5和8.4。 还能不能找到更高电子亲和能的物质?能,通过把F, BO2, 和AuF6这样的superhalogen和过渡金属阳离子结合形成配合物,它们的电子亲和能会更高,称为hyperhalogen,“超级卤素”。这就是弗吉尼亚理工的Puru Jena用电子亲和能玩出来的概念。他和他的研究生Devleena Samanta不断设计新的具有高电子亲和能的各种配体。 设计能够稳定存在的“超级卤素”的技巧,关键是要让“高级卤素”和极度缺电子的过渡金属阳离子结合,而且金属阳离子氧化态越高越好。过渡金属由于d轨道不满,容易形成多种高配位数的配合物。 他们最近发表在JACS上的一篇文章,“Zn的+3价氧化态”,不仅设计出了电子亲和能高达9.4的超级卤素Zn(AuF6)3,而且理论上预测了一种稳定存在的非常不寻常的Zn的氧化态,+3(作者强调为氧化态,或者配位态,而非阳离子价态)。我们知道Zn的电子组态时3d10,4s2,因此,一般只有+2价;但是,既然同族的Hg被证实可以有+4价,为什么Zn不可以有+3价?(Hg的+4价也是理论预测存在,并在长达20年后被实验证实存在的。)他们在这篇文章中,通过纯的DFT理论计算,预测了超级卤素Zn(AuF6)3的稳定存在。具体是不是真的存在。。。who knows? 纯理论计算发JACS还是比较难的,尤其是纯的B3LYP/6-311+g*,除非这个计算背后的理论或者概念足够重要。 不过,我认为这篇文章的对Zn的+3价的确定并非无懈可击,甚至可以说有重大缺陷,令人很难信服。这篇预测一种新型Zn氧化态配合物的文章背后涉及两个基本的计算理论问题:第一,怎么证明你设计的,或者通俗说你用GaussView搭出来的一堆原子组成的化合物真的能够稳定存在,尤其是一个从来没有人观测到过的Zn(III)价态?第二,你怎么确定你的那堆配合物中Zn真的是Zn(III)氧化态?因为对于非离子型化合物,确定原子化合价是比较困难的。 对于第一个问题,即怎么验证一个假想中的化合物能在实际中稳定存在,作者的回答还是比较好的。确定假想化合物热力学能够稳定存在,需要分析这个假想化合物各种不同的可能组成形式(构象),和各种不同的可能的分解形式,然后说明哪种状态下能量最低:是处于化合态时的能量低,还是处于解离态能量低。比如, H3明显是热力学不稳定的,因为分别优化计算H3, H2和H得到其能量,可知H3--H2 + H反应是放热的。因此H3只是假象中的化合物,实际中即使短时间碰到一起也会立即分解。 (注:有些化合物是动力学稳定的,即它虽然是热力学不稳定的高能态,但解离途径存在能垒,导致化合物可以短暂呆在这个状态;此文说的稳定是指热力学稳定,即只考虑产物态势能,而不考虑过渡态的势能) 因此,在这篇文章中,作者对中性ZnX3的各种解离产物进行结构优化,并比较各种解离路线是放热还是吸热的。如果所有解离路线都是吸热的,那么OK,说明ZnX3是能够稳定存在的。这对于ZnF3不难,但是对于Zn(AuF6)3就不容易了,要考虑Zn(AuF6)2 + AuF4 + F2这样的路线,和ZnF2 + Au2F10 + AuF4 + F2 这样的路线,等等。文章中称这种解离能为片段化能Fragmentation energies. 文章附件对各种体系罗列了总共数十种可能的解离方式,十分细致。他们发现,ZnF3不稳定,分解成ZnF2和(1/2)F2单质要放热;而中性Zn(BO2)3虽然能够稳定存在,但其中Zn是+2价的。只有Zn(AuF6)3这个东东可以稳定存在且其中Zn为3配位的。 对于第二个问题,即他们怎么确定配合物Zn(AuF6)中Zn的氧化价态,作者给出的理由并不充分,关键处含混其词。基本上来说,他们对优化好的结构做了NBO布局分析何轨道占据分析,然而却又指出NBO分析并不能确定无误地给出氧化价态的结果,最后给出了一个自己的判别条件,认为只要3个配体在优化中没有塌缩到一起的倾向,就能称之为3价氧化态存在。这实际上还是从几何结构稳定性上分析,认为3配位等于3价氧化态。不过,这3个配体却并非C3对称,而是C2对称,即有两个配体更接近一些。对于这一点,作者甚至回避了解释。并总结说明,预测了一个“配体不会聚集”的化合物,以强调其3配位态的稳定性。但是,3配位态等于3价氧化态吗?这个化合物中是真的存在3根配位键,还是2根配位键加一个很弱的分子间作用?我认为这是本文暗含的一个逻辑缺陷。 结论是,这是一篇值得看看的文章。如果有兴趣,可以接着讨论如何确定其中Zn的氧化态到底是否Zn(III);或者,沿着作者的思路,看看能不能和其它过渡态金属结合,预测新的,有更高电子亲和能的配合物。 参考文献: (1)JACS原文: Zn in the +3 Oxidation State DOI: 10.1021/ja3029119 (2)Phys.org 报道1: Enhancing the chemistry of zinc (3)Phys.org 报道2: Researchers discover a new class of highly electronegative chemical species (3)Angewandte通信: Hyperhalogens: Discovery of a New Class of Highly Electronegative Species 后续:Zn真的有Zn(III)价的氧化态吗?--JACS刊登反驳文章 http://bbs.sciencenet.cn/home.php?mod=spaceuid=38740do=blogquickforward=1id=590937
个人分类: 科学进展|15245 次阅读|32 个评论
[FP-LMTO方法] LmtART 7.04 的安装
热度 5 wonvein 2010-9-25 21:24
LmtART使用Fortran 90编写,用完全势线性muffin-tin(FP-LMTO)轨道方法计算电子结构。它在密度泛函理论的框架内进行键结构,总能量,和力的计算。 详细介绍见: 官方主页: http://www.fkf.mpg.de/andersen/docs/lmtoart_programs.html 官方主页: http://www.physics.ucdavis.edu/~mindlab/MaterialResearch/Scientific/index_lmtart.htm 最近用到lmtART这个软件,安装时遇到了问题,发现这个软件并不提供Makefile文件,而且是f和f90混编。为了编译LmtART,我只好自己写了Makefile文件,可以顺利编译。 其安装步骤为: 1 Edit the file man_lmtsetup.f (lines 119-122) and specify the path to the scratch and atomdat directories. Also check that other items match your computer settings. NAMEDENS='/public/home/wangwei/Apps/LmtART/den.' ! atomic density data NAMEATOM='/public/home/wangwei/Apps/LmtART/rat.' ! input atomic files NAMESTRC='/public/home/wangwei/Apps/LmtART/str.' ! structure files NAMELMTO='/public/home/wangwei/Apps/LmtART/lmt.' ! LMTO atom files 2. make Makefile for LmtART: # Makefile for LmtART v7.04 # Written by Wei Wang # 2010/07/26 #------------------------------------------------------------------------------- # Fortran compiler #------------------------------------------------------------------------------- FC=ifort FFLAGS = -O3 AR = ar #------------------------------------------------------------------------------- # Suffix rules #------------------------------------------------------------------------------- .SUFFIXES: .o .f .f.o: $(FC) $(FFLAGS) -c $ .SUFFIXES: .o .f90 .f90.o: $(FC) $(FFLAGS) -c $ #------------------------------------------------------------------------------- # Source files #------------------------------------------------------------------------------- SRC_mod_F77= mod_atoms.f mod_common.fmod_hubbard.fmod_models.f mod_qmc.f mod_supra.f\ mod_dimart.fmod_init.f mod_phonons.fmod_response.fmod_work.f SRC_mod_F90=mod_cls_l_lib.f90mod_cls_l_time.f90 mod_cls_l.f90 SRC_main=man_main.f SRC_F77=bnd_allocate.fchi_hubpar.f del_dtosint.f dmf_lmtqup.f frc_energy.f imp_hybint.f ini_makettr.f lib_ranw.f mdl_init.f opt_optpar.f pot_multasa.f rat_shletot.f \ bnd_bndpar.f chi_intchi.f del_dtospar.f dmf_lmtrat.f frc_forces.f imp_hyblev.f ini_readini.f lib_ratmom.f mdl_mixdmf.f opt_optpsi.f pot_multftr.f rat_shlfiles.f \ bnd_conv.f chi_inthiq.f del_dynbnd.f dmf_lmtrea.f frc_intvar.f imp_impfock.fini_readpnt.f lib_relharm.f mdl_modelhyb.fopt_optptb.f pot_multrho.f rat_shlmat.f \ bnd_coreny.f chi_lmtchi.f del_dynsym.f dmf_mixdmf.f frc_lmtvar.f imp_impfsig.fini_readscf.f lib_simq.f mdl_readmdl.f opt_opttet.f pot_potfiles.frat_siggrf.f \ bnd_fndfrm.f chi_lmtvec.f del_epibnd.f dmf_ordbnd.f frc_pulend.f imp_imphub1.fini_readstr.f lib_sphfun.f mdl_scfhub1.f opt_setopt.f pot_print2.f rat_slater.f \ bnd_ftrmat.f chi_magmat.f del_episym.f dmf_ordvec.f frc_pulint.f imp_impirat.fini_rpr48.f lib_sphharm.f mdl_scfirat.f opt_symopt.f pot_vcftr.f rat_storage.f \ bnd_getenr.f chi_magpar.f del_lmtvec.f dmf_setdos.f ftb_allocate.fimp_impqmcf.fini_storage.f lib_splin3.f mdl_scfqmcf.f phn_ctrlinf.f pot_vcoul.f rat_vcoul.f \ bnd_intmat.f chi_mttrd.f del_nmtfiles.fdmf_setfat.f ftb_denmat.f imp_impratf.fini_test.f lib_splkff.f mdl_scfratf.f phn_deltet.f pot_vexch.f rat_vexch.f \ bnd_intpar.f chi_mttrs.f del_psivec.f dmf_setfrs.f ftb_energy.f imp_impsbmf.flib_broy6.f lib_stnmat.f mdl_scfsbmf.f phn_dynmat.f pot_vxftr.f rho_emader.f \ bnd_lmtbnd.f chi_orbpar.f del_rintpar.f dmf_setgrf.f ftb_ftbbnd.f imp_setdos.f lib_calc.f lib_timel.f mdl_setgrf.f phn_epitet.f pot_vzero.f rho_mixro1.f \ bnd_lmtdft.f chi_printq.f del_rtosint.f dmf_sethyb.f ftb_ftbfat.f imp_sigrat.f lib_cinv.f lib_transcub.f mdl_storage.f phn_init.f qmc_allocate.frho_mixrob.f \ bnd_lmtfat.f chi_setchi.f del_sgmpar.f dmf_sgmdat.f ftb_ftbhub.f imp_sunhub1.flib_cmpdiag.f lib_wigmat.f phn_readlrt.f qmc_ctrlinf.f rho_orbmag.f \ bnd_lmtmat.f chi_sethiq.f del_strvec.f dmf_sgmmas.f ftb_ftbmft.f imp_sunirat.flib_conjgrad.f man_artinp.f phn_runtask.f qmc_getgrf.f rho_rencor.f \ bnd_lmtpar.f chi_setjey.f del_symdps.f dmf_sigmaw.f ftb_ftbpar.f imp_sunsbmf.flib_corsch.f man_artout.f phn_setepi.f qmc_init.f rho_renrov.f \ bnd_nmtpar.f chi_symchi.f del_symdrho.f dmf_tosint.f ftb_ftbstr.f imp_vertex.f lib_csplines.f man_atoms.f phn_setphn.f qmc_qmcpra.f rho_rhofiles.f \ bnd_orbpar.f cls_chem1.f dmf_allocate.fdmf_tospar.f ftb_hubint.f ini_blowmts.flib_cubharm.f man_electrons.fphn_spectra.f qmc_qmcyif.f rho_rhoful.f \ bnd_ovrpar.f cls_chem2.f dmf_denmat.f dmf_ttrint.f ftb_mixden.f ini_ctrlinf.flib_delrad.f man_ftbscf.f phn_storage.f qmc_readqmc.f str_dstr.f \ bnd_placeny.f cls_shellxnot.fdmf_energy.f dpt_delpot.f ftb_symden.f ini_getchi.f lib_deriv1.f man_impscf.f phn_symnuc.f qmc_setgrf.f str_hstr.f \ bnd_potpar.f del_allocate.f dmf_fndfrm.f dpt_delvhub.f hop_alpcon.f ini_getdos.f lib_det.f man_lmtchn.f plz_allocate.fqmc_storage.f str_sitegen.f \ bnd_seny.f del_d2intpar.f dmf_ftbfat.f dpt_dmultftr.fhop_findirr.f ini_getfat.f lib_dgemm.f man_lmtscf.f plz_plipar.f rat_agfmat.f str_strmsh.f \ bnd_setdos.f del_d2lmat.f dmf_ftbima.f dpt_dmultrho.fhop_gethtb.f ini_getfrs.f lib_dilog.f man_lmtsetup.f plz_plsbnd.f rat_allocate.fstr_vecgen.f \ bnd_setenr.f del_d2lmto.f dmf_ftbmas.f dpt_dvcftr.f hop_readhop.f ini_getgrf.f lib_drsub.f plz_plssym.f rat_cgfmat.f sup_ctrlinf.f \ bnd_seteny.f del_d2tosend.f dmf_ftbqup.f dpt_dvcoul.f hop_scrcon.f ini_gethiq.f lib_eigen1c.f man_mdlchn.f odf_allocate.fplz_plwcrd.f rat_clsgor.f sup_elifun.f \ bnd_setfat.f del_d2tosint.f dmf_ftbrat.f dpt_dvexch.f hop_scrstr.f ini_gethyb.f lib_erf.f man_mdlsetup.f odf_ftbodf.f plz_plzbnd.f rat_clsmat.f sup_init.f \ bnd_setfrs.f del_d2vint.f dmf_ftbrea.f dpt_dvxftr.f hop_stralp.f ini_getodf.f lib_forcesym.f man_models.f odf_lmtodf.f plz_plzend.f rat_clssav.f sup_omegaq.f \ bnd_strgnt.f del_delbnd.f dmf_grfbnd.f dpt_gradpot.f hub_hubpar.f ini_getopt.f lib_formt.f man_phead.f odf_odfbnd.f plz_plzfiles.frat_clsyin.f sup_phndos.f \ bnd_sum.f del_deleny.f dmf_grfexp.f dpt_potfiles.fhub_hubpot.f ini_groups.f lib_gradfun.f man_phnchn.f odf_odfint.f plz_plzinl.f rat_ctrlinf.f sup_readdyn.f \ bnd_symtos.f del_delexp.f dmf_grfgrp.f dpt_print3.f hub_mixhub.f ini_init.f lib_hubmat.f man_phnsetup.f odf_odfmat.f plz_plzint.f rat_energy.f sup_readepi.f \ bnd_tosend.f del_delgnt.f dmf_grflev.f dpt_scrpot.f hub_msbmesh.f ini_makeatm.flib_inverse.f man_phonons.f odf_odfpar.f plz_plzmat.f rat_getads.f sup_readsup.f \ bnd_tosint.f del_delmat.fdmf_grfloc.f dro_dbroy4n.f hub_readhub.f ini_makeenv.flib_lsqmom.f man_qmcchn.f odf_odfpsi.f plz_plzmto.f rat_getgrf.f sup_storage.f \ bnd_tospar.f del_delmto.f dmf_grfmat.f dro_delrho.f hub_readrep.f ini_makefft.flib_mklegw.f man_qmc.f odf_odfttr.f plz_plzmts.f rat_init.f sup_ttrint.f \ bnd_ttrint.f del_dhubpar.f dmf_grfpar.f dro_drofiles.fhub_rhohub.f ini_makegnt.flib_morefun.f man_qmcsetup.f odf_setodf.f plz_plzpar.f rat_matatm.f sup_widths.f \ chi_allocate.fdel_dinlmat.f dmf_grfwgt.f dro_gradrho.f imp_agfirat.f ini_makegrp.flib_order.f man_ratchn.f odf_symodf.f plz_plzpin.f rat_mixrat.f ttr_fermicof.f \ chi_chifiles.fdel_dintmat.f dmf_hubmat.f dro_gradrps.f imp_combrep.f ini_makehan.flib_pade.f man_ratsetup.f opt_allocate.fplz_plzsym.f rat_ratden.f ttr_fermiint.f \ chi_chimat.f del_dintpar.f dmf_imphyb.f dro_magdro.f imp_crfhub1.f ini_makeplw.flib_paulim.f man_supchn.f opt_ftbopt.f plz_polarz.f rat_ratmesh.f ttr_mcttrint.f \ chi_chipar.f del_dlmtpar.f dmf_implev.f dro_mixbrd.f imp_crfirat.f ini_makerad.flib_pcoefs.f man_supra.f opt_lmtopt.f plz_screen.f rat_ratscf.f ttr_ttrvel.f \ chi_chitet.f del_dnmtpar.f dmf_impmod.f dro_mixdps.f imp_crfqmcm.f ini_makescf.flib_prattpols.fman_supsetup.f opt_optdhk.f pot_exchcorr.frat_readrat.f \ chi_delmsh.f del_dovrpar.f dmf_lmtfat.f dro_mixdro.f imp_crfsbmf.f ini_makesmt.flib_pzeros.f mdl_allocate.f opt_opthan.f pot_gga91.f rat_setads.f \ chi_ftbchi.f del_dpotpar.f dmf_lmtima.f dro_spldps.f imp_fftsub.f ini_makesym.flib_qd.f mdl_ctrlinf.f opt_optint.f pot_gga96.f rat_setgrf.f \ chi_ftbvec.f del_dtosend.f dmf_lmtmas.f dro_spldro.f imp_hybfun.f ini_maketei.flib_radsch.f mdl_getgrf.f opt_optmat.f pot_lsda.f rat_shells.f SRC_F90=cls_angular1.f90 cls_exactdiag.f90 cls_fockvec1.f90cls_greenfun1.f90cls_hamilton2.f90cls_main1.f90 cls_subrtn.f90qmc_run.f90 qmc_sampling_PC_diag.f90 \ cls_angular2.f90cls_fill1.f90 cls_fockvec2.f90cls_greenfun2.f90cls_impurity.f90 cls_main2.f90lib_rapx.f90 qmc_run_switch.f90qmc_sampling_PC_random.f90 \ cls_bath.f90 cls_fill2.f90 cls_fun.f90cls_hamilton1.f90cls_l_diag.f90 cls_main.f90 qmc_fourier.f90 qmc_sampling.f90 OBJ_mod_F77 = $(SRC_mod_F77:.f=.o) OBJ_mod_F90 = $(SRC_mod_F90:.f90=.o) OBJ_F77 = $(SRC_F77:.f=.o) OBJ_F90 = $(SRC_F90:.f90=.o) OBJ_main = $(SRC_main:.f=.o) OBJ=$(OBJ_mod_F77) $(OBJ_mod_F90) $(OBJ_main) $(OBJ_F77) $(OBJ_F90) EXE = lmtart lmtart: $(OBJ) $(FC) $(FFLAGS) -o $(EXE) $(OBJ) $(LIB) clean: rm -f *.o *.mod *~ fort.* ifc* *.log $(EXE)
个人分类: 量化软件|9656 次阅读|4 个评论
MIT给大家解释离散傅立叶变换(DFT)
毛宁波 2009-12-1 01:43
我们大家知道,傅立叶变换是我们很多科学领域的重要的数学工具。有人说没有傅立叶变换就没有我们现代科学一些新分支或者很多学科就可能不太完善,我觉得这句话一点没有错。我们在地球物理的信号,无线电信号,现代的通讯技术、光学、声学甚至音乐中都大量应用傅立叶变换,尤其是离散的傅立叶变换。美国麻省理工学院(MIT)最近在她的网站主页上用非数学的语言给普通民众大谈离散傅立叶变换(DFT),这篇文章很通俗,推荐给大家共享,希望大家多多讨论。 网站链接: http://web.mit.edu/newsoffice/2009/explained-fourier.html Explained: The Discrete Fourier Transform The theories of an early-19th-century French mathematician have emerged from obscurity to become part of the basic language of engineering. Larry Hardesty, MIT News Office November 25, 2009 Science and technology journalists pride themselves on the ability to explain complicated ideas in accessible ways, but there are some technical principles that we encounter so often in our reporting that paraphrasing them or writing around them begins to feel like missing a big part of the story. So in a new series of articles called Explained, MIT News Office staff will explain some of the core ideas in the areas they cover, as reference points for future reporting on MIT research. In 1811, Joseph Fourier, the 43-year-old prefect of the French district of Isre, entered a competition in heat research sponsored by the French Academy of Sciences. The paper he submitted described a novel analytical technique that we today call the Fourier transform, and it won the competition; but the prize jury declined to publish it, criticizing the sloppiness of Fouriers reasoning. According to Jean-Pierre Kahane, a French mathematician and current member of the academy, as late as the early 1970s, Fouriers name still didnt turn up in the major French encyclopedia the Encyclopdia Universalis. Now, however, his name is everywhere. The Fourier transform is a way to decompose a signal into its constituent frequencies, and versions of it are used to generate and filter cell-phone and Wi-Fi transmissions, to compress audio, image, and video files so that they take up less bandwidth, and to solve differential equations, among other things. Its so ubiquitous that you dont really study the Fourier transform for what it is, says Laurent Demanet, an assistant professor of applied mathematics at MIT. You take a class in signal processing, and there it is. You dont have any choice. The Fourier transform comes in three varieties: the plain old Fourier transform, the Fourier series, and the discrete Fourier transform. But its the discrete Fourier transform, or DFT, that accounts for the Fourier revival. In 1965, the computer scientists James Cooley and John Tukey described an algorithm called the fast Fourier transform, which made it much easier to calculate DFTs on a computer. All of a sudden, the DFT became a practical way to process digital signals. Summing together only three discrete frequencies can produce a much more erratic composite. The Fourier transform provides a way to decompose signals into their constituent frequencies. To get a sense of what the DFT does, consider an MP3 player plugged into a loudspeaker. The MP3 player sends the speaker audio information as fluctuations in the voltage of an electrical signal. Those fluctuations cause the speaker drum to vibrate, which in turn causes air particles to move, producing sound. An audio signals fluctuations over time can be depicted as a graph: the x-axis is time, and the y-axis is the voltage of the electrical signal, or perhaps the movement of the speaker drum or air particles. Either way, the signal ends up looking like an erratic wavelike squiggle. But when you listen to the sound produced from that squiggle, you can clearly distinguish all the instruments in a symphony orchestra, playing discrete notes at the same time. Thats because the erratic squiggle is, effectively, the sum of a number of much more regular squiggles, which represent different frequencies of sound. Frequency just means the rate at which air molecules go back and forth, or a voltage fluctuates, and it can be represented as the rate at which a regular squiggle goes up and down. When you add two frequencies together, the resulting squiggle goes up where both the component frequencies go up, goes down where they both go down, and does something in between where theyre going in different directions. The DFT does mathematically what the human ear does physically: decompose a signal into its component frequencies. Unlike the analog signal from, say, a record player, the digital signal from an MP3 player is just a series of numbers, representing very short samples of a real-world sound: CD-quality digital audio recording, for instance, collects 44,100 samples a second. If you extract some number of consecutive values from a digital signal 8, or 128, or 1,000 the DFT represents them as the weighted sum of an equivalent number of frequencies. (Weighted just means that some of the frequencies count more than others toward the total.) The application of the DFT to wireless technologies is fairly straightforward: the ability to break a signal into its constituent frequencies lets cell-phone towers, for instance, disentangle transmissions from different users, allowing more of them to share the air. The application to data compression is less intuitive. But if you extract an eight-by-eight block of pixels from an image, each row or column is simply a sequence of eight numbers like a digital signal with eight samples. The whole block can thus be represented as the weighted sum of 64 frequencies. If theres little variation in color across the block, the weights of most of those frequencies will be zero or near zero. Throwing out the frequencies with low weights allows the block to be represented with fewer bits but little loss of fidelity. Demanet points out that the DFT has plenty of other applications, in areas like spectroscopy, magnetic resonance imaging, and quantum computing. But ultimately, he says, Its hard to explain what sort of impact Fouriers had, because the Fourier transform is such a fundamental concept that by now, its part of the language.
个人分类: 美国麻省理工学院见闻|8010 次阅读|1 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-6-3 09:57

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部