博文
用量子化学软件Gaussian计算周期性体系
|
Gaussian软件提供了PBC关键词来计算周期性体系。该部分内容主要来自Rice大学G. E. Scuseria等的贡献。G. E. Scuseria在2003年发表的JCP文章，用HSE方法在Gaussian软件包里进行了PBC计算。之后还有一系列相关工作，虽然该方法在Gaussian里用的并不多，但在Crystal里似乎已经成为一个重要的方法。
Gaussian进行周期性体系计算的例子网络上很少，特别是关于能带的计算。我这里做了一个测试，在此记录一下。
我们用一个简单的二维石墨烯体系为例，也就是Gaussian官网例子里的第二个。
Gaussian的输入文件如下：
%chk=graphene.chk |
其中iop(5/103)是控制K点能量本征值的选项，在iops列表里可以找到。
这里采用了较低的lsda/sto-3g计算级别，可参考G09软件包下tests文件里的test0474和test0475。
输出的graphene.log文件里会有这些行：
Using k point mesh (from -Pi to Pi): |
说明对1024个K点，进行了46×46的划分。
最后会输出各个K点的前10个本征值。
SCF Done: E(RS-VWN) = -74.8752722793 A.U. after 3 cycles |
将占据轨道和非占据轨道的能量本征值提取出来，即可得到能带。但Gaussian有一个缺点就是不能指定k point path，只能均匀mesh。为了获得指定高对称点之间的能量本征值，需要从这些mesh points里提取出那些位于高对称点之间连线上的点。
由于VASP可以灵活的设定k point path，为了和Gaussian的计算结果对比，我们也选取均匀mesh的方法。
VASP计算中，C用的是PAW_PBE基组，INCAR所有参数均采用默认。
考虑到二者所用的方法和基组并不相同，最后得到的能带结构的一致性还是比较好的，均在k=14左右找到了狄拉克点。
但没发现Gaussian里如何输出Fermi energy，所以狄拉克点能量在0 eV 以下，这可以跟VASP的结果对比来看。
其中k path在倒格子中的路径为
默认的Gaussian的k path不是这样的，为了对比进行了顺序调整。
从Gaussian log文件中提取能带数据的脚本为
extract_eig.sh.gz。(可用gzip -d解压)
https://wap.sciencenet.cn/blog-485752-1345468.html
上一篇：磁化率和磁矩的单位
