周期性体系¶

不建议使用pyscf来计算周期性体系

周期性体系一般指的就是晶体,由于其原子周期性重复的特征,输入文件和计算要求都和分子体系都有很大的不同,周期性体系以一个超晶胞为单位,我们需要指定原子坐标和晶格矢量.

输入文件¶

from pyscf.pbc import gto
cell = gto.M(
  atom = '''H 0 0 0; H 1 1 1''',
  basis = 'gth-dzvp',
  pseudo = 'gth-pbe',
  a = numpy.eye(3) * 2)

首先可以看到的是,这个gto模块不在是从pyscf中直接导入,而是从pyscf.pbc中导入,但是gto.M的用法和之前一致,依旧会自动完成build方法.

pyscf在处理周期性体系的时候往往喜欢使用赝势(pseudo)和赝势基组(gth-...),对轻元素,有gth-dzvp,gth-tzvp,gth-aug-dzvp...之类的基组,但是对重元素不适用,重元素的基组需要做一定的矫正,pyscf的pbc基组库里面支持重元素的basis只有gth-dzvp-molopt-sr和gth-szv-molopt-sr,量化计算的高斯型基组也能用.

同时,赝势基组和赝势要正确匹配,对赝势的书写要尽可能准确,例如,Cu匹配的gth-pbe赝势为gth-pbe-q11,不写q11这个后缀程序可能不能正确识别.至于后缀写什么,从pyscf的文件里面去查.但是很逆天的是,这个赝势和上面那个赝势基组是不配对的...,需要换一个赝势,同样,不同的赝势支持的元素也是有限的,还是需要从文件夹里面查,很重要的一点是,gth系列的基组只能和gth系列的赝势匹配 ,虽然你把赝势改成别的赝势也能跑.

与分子体系类似的是,我们可以指定晶胞的净自旋和净电荷,charge参数默认是0

还有一个我们要指定的参数就是晶格矢量a,默认是单位矩阵,有了晶格矢量,才能够重复描述整个周期性体系.

指定方式可以用矩阵形式代表三个向量(注意, 必须是要在右手坐标系下的坐标 ):

cell.a = [[2,0,0],[0,2,0],[0,0,2]]

k点¶

根据Bloch定理,周期性势场中粒子的波函数可以用Bloch波来表示:

\[ \phi_k(\vec{r})=\exp{(i\vec{k}\cdot \vec{r})}u_k(\vec{r}) \]

k是波矢量,k点可以认为是k矢量三个坐标方向的采样点,k点的个数越多,采样点越密,得到的结果越精确,相应的耗时也会增加.

k点网格可以使用晶体对象的make_kpts()方法指定,然后进行电子结构计算的时候需要相应输入:

kpts=cell.make_kpts([12,12,12])
kmf=scf.KRHF(cell,kpts=kpts)

如果不指定k点的话,程序默认在gamma点(布里渊区的原点)计算

k点采样等价于使用超胞进行计算(虽然我不是很理解这句话),例如k点网格为[2,2,2],就等价与使用含有8个原胞的超胞进行计算,此时的spin值为8个原胞内所有未成对电子的数目,对于一般的金属晶体, 通常都是闭壳层 (价电子可以自由移动从而进行配对),如果此时k点设置成为奇数个,很有可能导致总电子个数也为奇数个,spin就要设置成1才能跑,但是这违背了自旋为0的物理意义,也会加大计算时间(KRKS比KUKS要快的多),所以会有问题,通常我们都采用偶数个k点进行计算.

pyscf的pbc还有一个bug是,spin虽然指定的是等价超胞内的未配对原子个数,计算也是按照这样来的,但是程序会认为是你写的单胞内的自旋个数,然后不断抛出一大堆warning,这并不影响计算,接着算就完事了.

晶格参数优化¶

以Cu的晶格参数优化为例,我们尝试执行pyscf的周期性计算模块,通过不断计算大量的单点能拟合函数,然后求取E-a函数的极小点作为我们的优化结果.

在进行大规模优化的时候往往使用并行化计算方法:

from multiprocessing import Pool

# ... 在 if __name__ == '__main__': 块中
with Pool() as p:
    E = p.map(energy, x)
E = [e * 27.211386 for e in E]  # 转换为 eV

Pool是python的multiprocessing模块中的一个类,用于创建一个进程池执行并行计算任务,map方法用于将函数应用到可迭代对象的每个元素上,并返回一个结果列表.

但是应用到pyscf可能会有点问题,因为pyscf是默认全核并行的,两个并行之间可能会产生冲突,所以我们要换一种方式,每次计算都采用上一次计算的结果作为初始猜测,以此来降低计算时间.

Cu的FCC结构的晶格参数优化¶

使用密度矩阵继承的方法,由于每次改变晶格参数,其密度矩阵都大差不差,所以不妨在上一次迭代的基础上继续迭代以减少计算量,那么函数就需要传入两个参数以及传出两个参数.

代码为:

from pyscf.pbc import gto, dft
import numpy as np
import matplotlib.pyplot as plt

def energy(a0, prev_mf=None):
    cell = gto.M(
        atom=[['Cu', 0, 0, 0]],
        a=np.array([[0.5, 0.5, 0], [0, 0.5, 0.5], [0.5, 0, 0.5]])* a0,
        basis='sto-3g',
        ecp='lanl2dz',
        verbose=4,        # Reduce verbosity
        symmetry=True,
        spin=1,
        mesh=[10,10,10]
    )

    # Create Monkhorst-Pack k-point mesh
    kpts = cell.make_kpts([3, 3, 3])

    mf = dft.KUKS(cell, kpts=kpts).density_fit()
    mf.xc = 'pbe'
    mf.conv_tol = 1e-6    # Adjust convergence tolerance
    mf.max_cycle = 50     # Limit maximum number of SCF cycles

    if prev_mf is not None:
        # Use previous converged density matrix as initial guess
        dm = prev_mf.make_rdm1()
        mf.kernel(dm0=dm)
    else:
        mf.kernel()

    return mf.e_tot, mf

if __name__ == '__main__':
    x = np.linspace(3, 4, 10)
    E = []
    prev_mf = None
    for a0 in x:
        e_tot, prev_mf = energy(a0, prev_mf)
        E.append(e_tot * 27.211386)  # Convert to eV

    # Plotting code
    plt.figure(figsize=(10, 6))
    plt.scatter(x, E, color='black', label='data')
    plt.plot(x, E, '--', color='black')
    plt.xlabel('a constant (Å)')
    plt.ylabel('E_tot (eV)')
    plt.title('Optimization of Lattice Parameter')
    plt.legend(loc='best')
    plt.grid(True)
    plt.savefig('铜晶格FCC参数优化.png', dpi=300, bbox_inches='tight')

计算结果为:

alt text

最优晶格参数为3.12A左右,实验值为3.615A,计算结果偏小.但是有这个精度已经很感人了,误差大约为13%,还可以接受,如果想要进一步提升精度,则需要提高网格的密度,增加k点的数量,使用更高级的基组以及合适的赝势.