ASE原子模拟环境入门之晶体和能带结构-测试狗·科研服务

预存

Document

当前位置：文库百科 › 文章详情

ASE原子模拟环境入门之晶体和能带结构

来源：时间：2022-12-16 16:41:29 浏览：5135次

这章，我们计算晶体的性质。

建立体结构

ASE提供了三种框架来建立体结构：

ase.build.bulk()可以建立单质的体结构和少量化合物的点阵类型和点阵常数，但自定义有限。
ase.spacegroup.crystal()根据典型的晶体学信息(如空间群、晶格参数和基)创建原子。
ase.lattice从晶格和基底显式地创建原子。

习题：使用ase.build.bulk()来得到一个银单质的原胞，然后将它可视化。

已知银单质为FCC结构，所以这个功能可能会得到一个FCC原胞。但是我们总会希望能够验证我们所得到的结果。你能认出它的FCC结构吗？你可以，例如，使用ASE GUI来重复结构和识别A-B-C堆叠。ASE还应该能够验证它确实是一个FCC原胞，并告诉我们选择了什么晶格常数:


print(atoms.cell.get_bravais_lattice())

ASE中的周期结构用atoms.cell和atoms.pbc来表示。晶胞是一个cell对象，它使用三个向量来表示晶体的晶格。pbc是一个由三个布尔值组成的数组，它表示系统的每个方向上是否具有周期性。

体结构计算

对于周期性DFT计算，我们通常使用一些k点来自对布里渊区取样。包括GPAW和Aims在内的许多计算软件都接受kpts关键字，它可以是一个元组，如(4,4,4)。在GPAW中，平面波模式非常适合于较小的周期系统。使用平面波模式，我们还应该设置一个平面波截断(单位为eV):


from gpaw import GPAW, PW
calc = GPAW(mode=PW(600), kpts=(8, 8, 8),setups={'Ag': '11'}, ...)

这里我们使用了setups关键字来制定我们想要的是11个电子的PAW数据集，而不是默认的17个电子的数据集，这使得计算更快。

（原则上，我们应该确保k点和平面波截断都收敛——也就是说，编写一个循环来测试不同的采样，以确定两者都足够好，可以准确地描述我们想要的量。）

习题：使用GPAW运行一个银的体结构的单点计算。使用calc.write('Ag.gpw')保存它的基态，按照GPAW自己的格式。

态密度

保存了基态之后，我们可以重新加载它给ASE以提取态密度：


import matplotlib.pyplot as plt
from ase.dft.dos import DOS
from gpaw import GPAW
calc = GPAW('groundstate.gpw')
dos = DOS(calc, npts=500, width=0)
energies = dos.get_energies()
weights = dos.get_dos()
plt.plot(energies, weights)
plt.show()

调用width=0的DOS类意味着ASE很好地使用线性四面体插值方法计算DOS，这种方法花费时间，但给出了更好的表示。们也可以给它一个非零宽度，比如默认值0.1 (eV)。在这种情况下，它将使用一个简单的高斯模糊宽度，但我们需要更多的k点来得到相同质量的图。注意，能量轴的零点是费米能量。

习题：画出（银的）DOS。你可能记得银原子包含10个d电子和1个s电子。（DOS谱中）哪一部分来自s电子？哪一部分来自于d电子？

你大概知道，过渡金属的d轨道更加地局域在核的附近，而s电子更加地离域。

在体相结构中，s轨道有很多重叠，因此在很宽的能量范围内分裂成很宽的能带。d态的重叠要少得多，所以也分裂得少：它们形成一个窄带，并且具有很高的DOS，因为d电子的数量是s电子的10倍。

所以来回答这个问题：d能带占据了-6.2eV到-2.6eV之间的一个很高的、窄的能带的绝大部分。而所有超出这个区间的能带都是由更宽的s带能构成。

贵金属Cu、Ag和Au的特征，是d带被完全占据。即，整个d带在费米能级（能量=0）以下。如果我们计算任何其他的过渡金属，费米能级应该在d带之内的某个地方。

注意：我们可以计算s，p和d投影的DOS来更加确定地看到哪个能带具有什么特性。这时候我们需要查看GPAW的文档，或者其他计算软件的文档。我们不在这里做这些事情。

能带结构

我们计算一下银的能带结构。

首先我们需要建立一个能带路径。图像搜索引擎可以给我们展示一些参考图表。我们可以从Exciting和GPAW计算出来的能带结构中看到布里渊区路径“W L Γ X W K”。ASE知道这些字母，可以帮助我们图形化倒易晶胞：


lat = atoms.cell.get_bravais_lattice()
print(lat.description())
lat.plot_bz(show=True)

一般来说，ase.lattice模块提供了BravaisLattice类，用于分别表示3D和2D中的14+5个布拉菲格点。这些类知道高对称k点和标准布里渊区路径（使用AFlow约定）。

习题：为（W L G X W K）构建一个能带路径。你可以使用path = atom.cell.bandpath(...)—可以查阅cell文档来看这里需要哪些参数。这给我们一个BandPath对象。

你可以print()这个能带路径对象来查看它的基本信息，或者使用它的write()方法来将能带路径保存到一个json文件中，比如path.json。然后使用下面的命令将它可视化：


$ ase reciprocal path.json

一旦我们确定我们有一个合理数量的k点的良好路径，我们就可以运行能带结构的计算。如何触发一个带结构的计算取决于我们使用的是哪个计算软件，所以我们通常会参考该计算软件的文档(ASE总有一天会提供快捷方式，使这与普通计算软件协调更容易):


calc = GPAW('groundstate.gpw')
atoms = calc.get_atoms()
path = atoms.cell.bandpath(<...>)
calc.set(kpts=path, symmetry='off', fixdensity=True)

我们在这里告诉GPAW使用我们的k点路径，而不是使用对称性约化的k点，并固定电子密度。

然后我们就可以触发一个新的计算，这是一个非自洽计算，然后提取并保存能带结构：


atoms.get_potential_energy()
bs = calc.band_structure()
bs.write('bs.json')

同样，ASE命令行工具提供了一个有用的命令来从文件中绘制带结构:


$ ase band-structure bs.json

习题：计算，保存，然后画出银的能带图，使用K点路径[WLΓXWK]

你可能需要放大一些视野，一次看到全部。图将费米能级显示为虚线（但是并没有像DOS图那样定义为零）。观察能带结构，我们可以看到大部分由d轨道构成的复杂的缠结，以及少量由s态构成的低于d能带的能级（在Γ点）和高于d能带的能级（穿过费米能级）。

状态方程

通过状态方程的计算，可以找到最优格参数，并计算体积模量。这意味着在一个范围内采样能量和晶格常数，以获得最小值以及曲率，这给了我们体积模量。

ASE的在线文档已经提供了如何使用经验EMT势的教程:https://wiki.fysik.dtu.dk/ase/tutorials/eos/eos.html

习题：运行EOS教程。

复杂晶体和晶胞优化

对于简单的FCC结构，我们只有一个参数a, EOS的拟合告诉我们需要知道的一切。

对于更复杂的结构，我们首先需要一个更高级的框架来构建原子，比如ase.spacegroup.crystal()函数。

以下内容会帮助我们如何构建金红石结构，使我们省去了查找原子基（basis）和其他晶体结构信息的麻烦。金红石是二氧化钛的一种常见矿物形式。

习题：构建和可视化金红石结构。

我们来优化结构。除了原子位置，我们必须优化晶格常数，它是以两个长度a和c为特征的四角晶格。

优化晶格常数需要能量对晶格常数的导数，也就是应力张量。atoms.get_stress()计算并返回应力作为一个6元向量(Voigt形式)。使用它需要附加的计算软件支持应力张量。GPAW的平面波模式做到了这一点。

ase.constraints.ExpCellFilter允许同时优化晶胞和位置。它通过将自由度暴露给优化器来实现这一点，就像它们是额外的位置一样——因此起到了一种过滤器的作用。我们用它来包裹atoms变量:


from ase.optimize import BFGS
from ase.constraints import ExpCellFilter
opt = BFGS(ExpCellFilter(atoms), ...)
opt.run(fmax=0.05)

习题：使用GPAW的平面波模式优化金红石的原胞。你可能要使用至少500eV的平面波截断。优化后的晶格常数a和c是多少？

习题：计算金红石的能带结构。它与搜索引擎搜到的金红石的能带结构符合吗？

习题答案

Ag基态：


from ase.build import bulk
from gpaw import GPAW, PW
atoms = bulk('Ag')
calc = GPAW(mode=PW(350), kpts=[8, 8, 8], txt='gpaw.bulk.Ag.txt',
            setups={'Ag': '11'})
atoms.calc = calc
atoms.get_potential_energy()
calc.write('bulk.Ag.gpw')

Ag DOS：


import matplotlib.pyplot as plt
from gpaw import GPAW
from ase.dft.dos import DOS
calc = GPAW('bulk.Ag.gpw')#energies, weights = calc.get_dos(npts=800, width=0)
dos = DOS(calc, npts=800, width=0)
energies = dos.get_energies()
weights = dos.get_dos()
ax = plt.gca()
ax.plot(energies, weights)
ax.set_xlabel('Energy [eV]')
ax.set_ylabel('DOS [1/eV]')
plt.savefig('dos.png')
plt.show()

Ag能带结构：


from gpaw import GPAW
calc = GPAW('bulk.Ag.gpw')
atoms = calc.get_atoms()
path = atoms.cell.bandpath('WLGXWK', density=10)
path.write('path.json')
calc.set(kpts=path, fixdensity=True, symmetry='off')
atoms.get_potential_energy()
bs = calc.band_structure()
bs.write('bs.json')

金红石晶胞优化：


from ase.constraints import ExpCellFilter
from ase.io import write
from ase.optimize import BFGS
from ase.spacegroup import crystal
from gpaw import GPAW, PW
a = 4.6
c = 2.95
# Rutile TiO2:
atoms = crystal(['Ti', 'O'], basis=[(0, 0, 0), (0.3, 0.3, 0.0)],
                spacegroup=136, cellpar=[a, a, c, 90, 90, 90])
write('rutile.traj', atoms)
calc = GPAW(mode=PW(800), kpts=[2, 2, 3],
            txt='gpaw.rutile.txt')
atoms.calc = calc
opt = BFGS(ExpCellFilter(atoms), trajectory='opt.rutile.traj')
opt.run(fmax=0.05)
calc.write('groundstate.rutile.gpw')
print('Final lattice:')
print(atoms.cell.get_bravais_lattice())

金红石能带结构：

calc = GPAW('groundstate.rutile.gpw')

atoms = calc.get_atoms()

path = atoms.cell.bandpath(density=7)

path.write('path.rutile.json')

calc.set(kpts=path, fixdensity=True,

         symmetry='off')

atoms.get_potential_energy()

bs = calc.band_structure()

bs.write('bs.rutile.json')

计算狗·模拟计算

上一篇：ASE原子模拟环境入门之结构优化: H2O

下一篇：ASE原子模拟环境入门之纳米颗粒

相关文章

基础理论丨一文了解XPS（概念、定性定量分析、分析方法、谱线结构）

2020-05-03

手把手教你用ChemDraw 画化学结构式：基础篇

2021-06-19

晶体结构可视化软件 VESTA使用教程（下篇）

2021-01-22

【科研干货】电化学表征：循环伏安法详解（上）

2019-10-25

【科研干货】电化学表征：循环伏安法详解（下）

2019-10-25

Zeta电位的基本理论、测试方法和应用

2020-08-24

项目推荐/Project

第一性原理-COHP

第一性原理-反应能垒

热门文章/popular

基础理论丨一文了解XPS（概念、定性定量分析、分析方法、谱线结构）

手把手教你用ChemDraw 画化学结构式：基础篇

晶体结构可视化软件 VESTA使用教程（下篇）

【科研干货】电化学表征：循环伏安法详解（上）

电化学实验基础之电化学工作站篇（二）三电极和两电极体系的搭建和测试

【科研干货】电化学表征：循环伏安法详解（下）

微信扫码分享文章

举报有奖

TEL: 191-3608-6524

如：在网络上恶意使用“测试狗”等相关关键词误导用户点击、恶意盗用测试狗商标、冒称官方工作人员等情形，请您向我们举报，经查实后，我们将给予您奖励。

举报内容：

200

上传附件：

文件格式不正确，请重新上传文件格式不正确，请重新上传文件格式不正确，请重新上传

文件格式：jpg、jpeg、png、gif、tif、doc、docx、ppt、pptx、xls、xlsx、pdf、zip、rar

联系方式

姓名

电话

提交意见

意见反馈

Suggestions

您可以在此留下您宝贵的意见，您的意见或问题反馈将会成为我们不断改进的动力。

意见类型

测试服务

网站功能

财务报账

其他类型

意见内容

200

联系方式

姓名

电话

提交意见