攻略：第一性原理计算结构优化

前言
结构优化的过程是获得合理结构的过程。所谓合理的结构，大多数场景下为在某些限制下能量更低的结构。这些限制，可以是结构维度的限制（比如限制为二维材料）、晶格的限制（比如施加了应变、放置于衬底上）、元素分布的限制（比如进行吸附时要求吸附物的空间分布）、对称性的限制（比如双层二维材料不同的堆垛方式导致体系对称性不同）。若使用不合理的结构进行后续的计算，比如体系能量、电子结构、磁态的计算结果就会出现较大的误差甚至是致命错误。
一、对结构优化影响较大的计算参数
一些解释：
1）弛豫：在本文里即为结构优化的意思。
2）电子步：自洽迭代时一次迭代称作一个电子步。
3）离子步：一次自洽计算称作一个离子步。自洽计算通常需要进行多步迭代，即一个离子步包含多个电子步。

4）分步优化：在结构优化时，可以先使用较低精度的参数组合，快速获得一个较为合理的结构；然后再使用更高精度的参数组合进一步优化，直至获得合理的结构。这里的分步即是多次结构优化计算。所谓的较为合理的结构，其反面对应的是离合理情况较远的的结构，常见的有：未弛豫的间隙原子掺杂、空位原子、表面小分子吸附、手动搭建的新结构等。特别是对于构建了缺陷的大体系，使用分步优化甚至可以节省一半以上的计算耗时。关于分步优化的进一步讨论见下文。
5）k-spacing：在倒空间中在 ka，kb，kc 三个倒格矢上相邻 k 点的距离，由 KPOINTS 文件控制。
6）如何获得初始结构 POSCAR：
从现有的材料数据库中导出cif 文件，然后使用 VESTA 或者VASPKIT 把 cif 文件转换为 POSCAR 格式的文件。常用数据库有：
https://materialsproject.org/#search/materials/
http://www.aflowlib.org/
https://materials.springer.com/periodictable#
基于现有结构搭建。有许多软件都可以用于建模，常用的有Materials Studio (简称 MS) 。使用 MS 建好模型后，导出cif 文件，然后使用VESTA 或者 VASPKIT 把cif 文件转换为POSCAR 格式的文件。
7）结构优化的过程见下列简单的流程图 (简单绘画，一些不当的细节请忽略)：VASP 会对当前结构进行自洽迭代获得该结构的电荷密度，计算该结构的能量和原子受力，然后判断是否达到结构优化收敛判据。若没达到，则生成新结构，使用新结构进行自洽计算 … 如此往复，直至达到收敛判据或者达到设定的最大离子步数(NSW 控制)，计算退出。如何判断弛豫是否完成见下文 “结构优化需要收集的信息” 。
8）备注：个人只使用过负的 EDIFFG，在此不考虑正数 EDIFFG 的情况。

PREC、ENCUT、EDIFF、KPOINTSPREC、ENCUT、EDIFF、KPOINTS 这四个参数会显著影响电子迭代的速度、以及得到的结果的精度。自洽迭代的精度会进一步影响各个原子上受力的精度。
1、对于大体系或者离合理结构较远的体系的结构优化，强烈建议进行分步优化。个人常用的分步优化方案见下文 “个人常用分步优化方案” 。
2、

ISPIN ISPIN 参数确定自洽计算是否考虑自旋极化。若考虑自旋极化，自旋向上和自旋向下的电子将被当作不同的对象进行处理，此时电子步耗时会增加 (有时还会导致自洽迭代步数增加，甚至导致自洽迭代不收敛) 。对于暂不确定是否有磁性的体系，建议在结构优化时打开ISPIN =2，并设置合理的 MAGMOM 值。若进行分步优化，我自己在第一步低精度优化时会考虑关掉自旋极化 (即 ISPIN = 1) (对于小体系，我可能会在此时直接打开 ISPIN=2 ) ，而在后面高精度结构优化时打开ISPN=2。

ALGO、IMIX(及 IMIX 相关参数) ALGO、IMIX(及 IMIX 相关参数) 控制电子自洽迭代的算法，既影响每一步电子步耗时，也影响自洽收敛所需电子步数。关于这一参数的具体讨论可参考本人的帖子《加快磁性材料电子迭代收敛经验小结》。

LREALLREAL 同样会影响电子迭代速度和计算精度。贴一段手册中对这一参数的描述：We recommend to use the real-space projection scheme for systems containing more than 20 atoms.
1）根据手册描述，若体系原子大于20 个，可以考虑使用 LREAL=Auto。
2）若关注的性质需要误差很小，比如极小的带隙 (meV 范围)，磁各向异性能 (多在μeV 范围内)，建议使用LREAL=.FALSE. 。

3）同样，若做分步优化，在第一步粗精度优化时，若体系原子大于 20 个，个人会考虑使用LREAL=Auto。

EDIFFGEDIFFG，负值为每个原子上的最大受力。若体系中每个原子受力小于该值，则判断为达到结构优化要求，完成计算。若做分步优化，在第一步粗精度优化时，个人常用EDIFFG=-0.05 ~ -0.08，最后的精优化常用 EDIFFG=-0.01。

NSWNSW，结构优化中最大离子步数。若计算 NSW 步离子步后，结构优化依然不收敛，则退出计算。若做分步优化，在第一步粗精度优化时，根据个人经验，NSW 最大设成200 即可。第一步粗优化跑完200 步依然不收敛也没有关系，此时我们的目的 (获得较合理的结构) 已经达到了，可以考虑进行下一步的精优化。

ISYMISYM，控制自洽迭代时是否需要考虑对称性。结构优化时多用ISYM=0 或 ISYM=2。设置 ISYM=2 打开对称性的话，VASP 中会涉及如下操作：
1、确定结构的对称性：结构所属的布拉维晶格类型、点群和空间群。
2、布拉维格子类型在OUTCAR 中“LATTPY”关键字后面打印出来 (阅读 VASP 源代码可知)，见下面第一张图。
3、在“LATTPY” 关键字后面几十行，打印了体系所属点群信息 (其中包括 static configuration 和 dynamic configuration，我目前还不确定这两个具体都有什么区别) ，见下面第二张图。
4、如果体系打开了自旋极化 (ISPIN=2)，还会给出考虑了磁态后所属的点群。MAGMOM 的设置会影响这一对称性。要提一句的是，MAGMOM 的正负值不影响磁态所属点群，其绝对值才影响磁态所属点群。以CrI3 结构为例，CrI3 属于D3d 点群。若MAGMOM = 6*0 2*3，此时磁态属于 D3d 点群。若把其中一个Cr 原子初始磁矩设为设为负数，即MAGMOM = 6*0 3 -3，此时磁态依然属于D3d 点群。若改变一下 Cr 原子初始磁矩大小，比如MAGMOM=6*0 3 3.5，此时磁态属于 D3 点群，见下面第三张图。

5、确定考虑了对称性后用于计算的不可约 K 点，查看 IBZKPT 文件可见。此时用于计算的 K 点个数减少，电子步耗时将减少，计算速度加快。
若还未确定合理结构所属对称性：在进行分步优化，个人建议在第一步粗优化时不要加上对称性 (即设置 ISYM=0) 。关于对称性的更具体讨论见下文。

IBRIONIBRION，控制产生新结构的算法。若进行分步优化，个人在第一步粗优化时使用IBRION=2。若认为当前结构离合理结构比较接近，个人会考虑换用 IBRION=1。IBRION 参数的说明从手册中摘抄如下：
IBRION=1：对于初始结构接近局域极小时是最好的选择。
IBRION=2：比较可靠，弛豫困难时，推荐使用。
IBRION=3：对非常糟糕的初始结构进行弛豫时可考虑使用。

ISIFISIF，在产生新结构时，决定晶胞和原子如何变化。手册中有一个表格说明了不同 ISIF 值对应什么意思，此处不再摘抄下来。我个人常用的 ISIF 解释如下：
1) ISIF=2：保持晶格常数不变，仅弛豫原子。施加双轴应变进行结构优化，或者弛豫零维体系时，就需要使用 ISIF=2。一些复杂体系弛豫时，为了加快结构优化，第一步粗优化也可以考虑使用 ISIF=2 弛豫后，再在后面的精优化中使用 ISIF=3。
2) ISIF=3：在弛豫原子时，也弛豫晶格。在弛豫晶格时，可以通过修改 VASP 源代码 constr_cell_relax.F 文件，实现一些 VASP 手册中没有提到的结构优化方式，比如固定某个晶格的优化，具体讨论见下文。在计算二维材料时，在设置了足够大的真空层后，个人喜欢使用 ISIF=3，并在优化时固定住 c 方向晶格 (真空层位于 c 方向) 。

POTIMPOTIM，(简单理解为) 在产生新结构时，控制原子移动的距离。VASP 会根据原子上的受力以及 POTIM 计算原子移动的距离。POTIM 越大，原子移动的距离越大。
1）在作分步优化时，对于极度不合理的结构，第一步粗优化个人有时会考虑设置增大 POTIM 值 (如 POTIM=0.2)，加快原子移动到较合理的位置。
2）若结构经多次弛豫依然无法收敛，此时结构可能离合理结构非常接近，但一直在势能极小值点附近振荡，此时将 POTIM 值调小，使用 IBRION=1 会有助于结构的收敛。
二、几种常见情况及讨论

对称性的影响对称性的影响：在结构优化过程中，对称性极为重要。寻找合理结构的过程，一方面也是确定结构对称性的过程。对称性主要影响三个方面：
1、计算速度：
考虑了对称性后，只需要使用简约区内的 k 点进行积分计算，电子步耗时缩短。同时，原子受力也会受到对称性的影响，原子只能在对称性允许的方向上移动，保证结构的对称性不会降低 (在弛豫过程中，结构对称性可以升高，但不会降低) 。一般情况下，给体系加上合理的对称性会使得体系更快收敛到合理的结构。但若加上不合理的对称性，会导致以下两个问题。在个人常用分步优化方案中说明了如何寻找结构合理的对称性。
2、体系的合理性：
在搭建新结构时，以某纯平面二维材料中掺入过渡金属 (TM) 为例。如果在初始搭建模型时，TM 与二维材料处于同一平面，那么在弛豫时，(即使不加对称性 ISYM=0) ，TM 在垂直平面方向受力始终为0，弛豫后整个材料依然会保持纯平面的结构。但如果 TM 稍微突出二维材料平面，弛豫后 TM 及其附近的二维材料出现较大起伏。两种结构的能量对比如下表，此时后者体系能量明显低于前者。

往体系中掺杂间隙原子 (比如上例) ，或者从已有结构做元素替换成新结构 (比如从 MoS2 出发把 Mo 替换成其他 TM，或者把一部分的 S 替换成其他元素) ，个人建议在构建初始结构时，手动移动某些关键原子，尽量使得结构只具有 P1 对称性，避免初始结构处于势能面上某些特殊点。这样往往可以获得能量更低的结构。
3、对电子结构(主要是对与能带相关的量)的影响：
以石墨烯为例。如果关注石墨烯布里渊区高对称点 K 点处的性质。如果在结构优化时，没有找到正确的对称性 (D6h)，而 VASP 给确定成了 C2h 对称性。C2h 对称性的石墨烯晶格常数 a ≠ b，格矢 a 和 b的夹角可能与 120°存在一定偏差，这会导致该结构倒空间中的高对称点 K 点与真正的 K 点 (D6h 对称性的石墨的 K 点) 有偏离，进一步导致与 K 点有关的量存在计算误差。而这个误差，是可以通过给结构加上正确的对称性完全避免的。

个人常用分步优化方案个人常用分步优化方案：对于基于现有结构产生的新结构 (比如构造缺陷，元素替换，原子位置调整) 和纯粹手搭的全新结构，常常需要多次优化才能获得合理结构。以两步优化为例，个人常用优化方案如下：
Step1：手动微调一些原子，使其只具有 P1 对称性 (比如在 MS 中以 0.01 的精度找到的对称性为 P1 )。
PREC = M
ISPIN = 1
EDIFF = 1E-04 (若体系较小或者认为离合理结构较近，考虑使用 1E-05)
EDIFFG = -0.05 (若体系较小或者认为离合理结构较近，考虑使用 -0.03；若体系很大，考虑使用 -0.08)
LREAL = Auto (原子数目大于 20 使用 Auto，小于 20 使用 .FALSE.)
NSW = 200
IBRION = 2
ISIF = 3 (有的时候使用 ISIF = 2 也是不错的选择)
ISMEAR = 0, SIGMA = 0.1
ISYM = 0
K-spacing = 0.03-1 ~ 0.04-1
Step2：将 step1 获得的结构导入 MS 中，按照 0.05 ~ 0.1 的精度寻找对称性并加上该对称性。该对称性可能会是体系合理的对称性。
PREC = N
ISPIN = 2 (若认为体系没有磁性，可使用 ISPIN = 1)
EDIFF = 1E-06
EDIFFG = -0.01
LREAL = Auto (原子数目大于 20 使用 Auto，小于 20 或者想进一步提高计算精度，则使用 .FALSE.)
NSW = 200
IBRION = 2
ISIF = 3
ISMEAR = 0, SIGMA = 0.1 (见下面讨论)
ISYM = 2
K-spacing = 0.02-1 ~ 0.025-1
POTIM (可以不调整，若体系多次优化依然不收敛，使用更小的 POTIM 值)
Step3：在 step1 的基础上，可以大致判断出体系为金属还是半导体，判断方法见本人帖子《加快磁性材料电子迭代收敛经验小结》。根据体系是半导体还是金属，调整 ISMEAR 和 SIGMA。
半导体：ISMEAR=-5
金属：ISMEAR=1，SIGAM=0.1。
可以继续使用 ISMEAR=0，但是需要根据 EENTRO 设置合理的 SIGMA 值。对于金属，根据 step1 中最后的 EENTRO 值调整 SIGMA 值，尽可能使得 EENTRO/atom 在 1meV ~ 2meV。(SIGMA 越大，EENTRO 最大)。
Step4：每一步需要收集相应的信息，见下文。如果 POSCAR 和 CONTCAR 的晶格常数差别较大，优化过程中使用的平面波对 CONTCAR 的描述可能已经有较大的偏差，(即使在 step2 中弛豫收敛了) 需要考虑使用 CONTCAR 进一步优化。个人常用的标准为 1% (即 POSCAR 和 CONTCAR 晶格常数只差的最大值)。、

已有结构的优化已有结构的优化：对于已知对称性的结构，可以直接加上其对称性，使用合理的晶格常数进行优化。

固定某个晶格的优化固定某个晶格的优化：通过修改 VASP 源代码 constr_cell_relax.F 文件，实现一些 VASP 手册中没有提到的结构优化方式，比如固定特定晶格的优化 (目前个人认为也能够实现固定特定晶格夹角的优化，未来若实现了这一功能，会发布详细的讨论及说明) 。如何实现固定特定晶格的优化，可以参考刘锦程博士的博文：
http://blog.wangruixing.cn/2019/05/05/constr/ 。
需要提醒的是：
1、OPTCELL 文件中 0 和 1 之间是没有空格的。
2、刘锦程博士的博文中还讨论了一种情况：将 OPTCELL 中上三角矩阵元设为 0 可以使得在弛豫过程中晶胞晶格不发生转动。但个人认为这么做可能存在一定的风险。具体讨论如下：
简单起见，假设晶格 c 轴垂直于 ab 轴平面。同时假设 a 轴沿笛卡尔坐标系 x 轴正方向（即 L_ax > 0），b 轴 x 方向分量 L_bx > 0 且 y 方向分量 L_by > 0。示意图如下：

格矢变换矩阵为：其中 F 矩阵为 FCELL 矩阵，L 矩阵为 A_OLD 矩阵，其中正数标记为绿色。

在 L[sub]ax[/sub]、L[sub]bx[/sub]、L[sub]by[/sub] 均大于 0 的情况下，在 F 矩阵元取不同值时，等式右边晶格增量矩阵中 a 晶格在 x 和 y 方向的增量：F[sub]ax[/sub]*L[sub]ax[/sub]+F[sub]ay[/sub]*L[sub]bx[/sub] 和 F[sub]ay[/sub]*L[sub]by[/sub] ，可能是正值也可能是负值，一共有四种情况，见下图。

若将F[sub]ay[/sub]矩阵元置零，则 a 晶格在 x 和 y 方向的增量分别从 F[sub]ax[/sub]*L[sub]ax[/sub]+F[sub]ay[/sub]*L[sub]bx[/sub] 和 F[sub]ay[/sub]*L[sub]by[/sub] 退化为 F[sub]ax[/sub]*L[sub]ax[/sub] 和 0。这时除了 y 方向增量被置零，x 方向的增量也受到了影响。若 F[sub]ax[/sub] > 0，对应 F[sub]ax[/sub]*L[sub]ax[/sub]> 0，此时若 F[sub]ay[/sub] < 0，是有可能导致 F[sub]ax[/sub]*L[sub]ax[/sub]+F[sub]ay[/sub]*L[sub]bx[/sub]0) 对应的能量。为最后给出的体系的能量
2、EENTRO：OUTCAR 中最后一个 EENTRO 值。为最后给出的体系的 “熵” 值。
3、require：搜索 OUTCAR 中是否打印了【reached required accuracy - stopping structural energy minimization】这一句话。若有，则说明弛豫收敛了，显示 “reached” 。
4、pressure：OUTCAR 中最后一个 external pressure 对应的值。对于不加压计算，其绝对值通常小于 1。
5、time：OUTCAR 中 Elapsed time 对应的时间。为计算耗时。
6、MAGNE：OSZICAR中最后一个 mag 对应的能量。为最后给出的体系的总磁矩。
7、Edisp：OUTCAR 中最后一个 Edisp 对应的能量。为考虑了 vdw 校正，最后给出的 vdw 校正能量。若不设置 vdw 校正，则输出 “--” 。
8、Pullay：OUTCAR 中最后一个 Pullay stress 对应的值。不加压计算，该值为 0。
9、sta_sym、dyn_sym、mag_sym：分别为 static configuration、dynamic configuration、magnetic configuration 对应的点群。
10、diff(%)：三个格矢三个分量以及各个晶格常数 (POSCAR - CONTCAR)/POSCAR 值。在 Total 一栏中的值若大于 1%，我个人会考虑再次优化。
11、subt(A)：三个格矢三个分量以及各个晶格常数 (POSCAR - CONTCAR) 值。

若体系弛豫是打开了自旋极化 (ISPIN=2)，还需要搜集整理磁矩相关的信息。
1、mag_tot：OSZICAR中最后一个 mag 对应的能量。为最后给出的体系的总磁矩。即为上一个截图中的 MAGNE。
2、ISPIN：若设置了 ISPIN=1，或者每个原子磁矩绝对值都小于 0.01，则该值为 1。
3、atom_nam：元素符号。
4、atom_mag：OUTCAR 中最后一个离子步结束后给出的原子磁矩。

下图为个人将上述信息简要记录到 word 文件中所用的模板：

上述两个截图均由脚本生成，脚本稍后会上传至 qq 群 (432953524，DFT计算之家-砥砺前行)。
四、其他讨论
结构优化给出的态密度是否可用：阅读手册 DOSCAR 的介绍部分，可以看到在弛豫时，DOSCAR 包含的是平均 DOS。手册不建议使用弛豫时的 DOSCAR。“For dynamic simulations and relaxations, an averaged DOS and an averaged integrated DOS is written to the file.… Mind: For relaxations, the DOSCAR is usually useless.”。个人认为，如果弛豫时离子步只跑了很少的步数 (比如 20 步以内)，结构没有发生较大变化，即弛豫过程体系电子结构没有发生较大变化，可以通过弛豫输出的 DOSCAR 获得体系的 DOS，从而大概判断体系为金属还是半导体，带隙大小。
结构优化给出的能量是否可用：如果在弛豫过程中，晶格变化较小 (比如晶格常数相差 5% 以内，晶格夹角相差 5% 以内)，结构没有发生较大的变化，可以使用弛豫给出的体系能量做一些快速判断 (比如吸附能，铁磁反铁磁能量差)。更进一步，弛豫在 20 步以内快速收敛，且晶格变化非常小 (1% 以内)，结构优化给出的数据与后续静态计算给出的数据几乎是一样的。下表为某体系吸附小分子的吸附自由能比较，可以看到分别使用弛豫和自洽计算的能量获得的吸附自由能相差在 1meV 以下。需要进一步说明的是，若要获得更加准确的能量，建议增加 K 点进行自洽计算。

五、总结
1、结构优化的重要性无论如何强调都不为过。
2、一定要花时间做各种检查：计算前检查初始文件，计算中检查弛豫过程，计算后检查结构及各种指标。
3、一定要花时间记录好弛豫后的信息，特别是对于磁性体系，强烈建议记录好结构的磁矩信息。之后的计算，包括单点能计算，能带计算，声子谱计算，都应该检查自洽结束后输出的磁矩信息，确保体系收敛到了弛豫给出的磁态。
本文作者为中科大无缺。
本文系学术之友微信公众号与材料人联合推出的计算攻略专栏。

【更多攻略】材料人APP持续推出了第一性原理、分子动力学系列课程，包括vasp、lammps、cp2k、qe、gromacs等软件的教学视频，永久回看，可开发票。
相关课程介绍
掌握vasp使用及固体物理与表面计算就在6月北京
想要系统学习LAMMPS分子动力学模拟技术与应用？看这套视频
视频教程：晶格振动模拟和计算 (phonopy 软件入门)
基于vasp的动力学模拟分析讲座即将开讲
七小时讲座带你入门利用QE开展计算化学工作
电化学计算系列之电池材料计算线上小班即将开课（最后一次课）
#计算化学#利用QE开展第一性原理讲座今天开讲

如有需要，欢迎下载材料人APP观看

下雨心会烦吧 · 2020-9-24 10:49:07

请问Step1-3中提到的利用MS手动调P1对称性是如何进行的呢谢谢

攻略：第一性原理计算结构优化

1 个回复