本文详细介绍 VASP 官方 wiki《DFT+DMFT calculation》内容。VASP可通过接口与外部DMFT代码(如TRIQS/solid_dmft)协同工作通过求解杂质模型获得自能，进行电荷自洽的CSC DFT+DMFT循环，输出可经最大熵方法解析延拓至实频得到局域谱函数。

 

 

 

 

 

 

前置准备

 

VASP    6.5.1       

TRIQS    大于3.0.0    

通过conda安装

conda create -n triqsconda activate triqsconda install -c conda-forge triqs

或借助mamba

conda install -c conda-forge mambamamba install -c conda-forge triqs

 

DFTTools  

CTHYB      

MaxEnt 

  Soild-dmf

 conda install -c conda-forge triqs_dft_tool triqs_cthybpip3 install  solid_dmf  triqs-maxen

Wannier90    3.1.0  

 

 

 

VASP  基态 SCF

 输入文件准备

IPOSCAR、

 NiO4.17  0.0000000000  0.5000000000  0.5000000000  0.5000000000  0.0000000000  0.5000000000  0.5000000000  0.5000000000  0.0000000000 Ni O  11 Direct  0.0000000000  0.0000000000  0.0000000000  0.5000000000  0.5000000000  0.5000000000

INCAR 

SYSTEM  = NiOISMEAR  = -5          # 四面体展宽法SIGMA   = 0.01EDIFF   = 1.E-10      # 严格收敛标准NELMIN  = 35NBANDS  = 32          # 略多于默认值LORBIT  = 14          # 优化投影器LMAXMIX = 4# 绘制 DOS 用的能量范围EMIN    = -6EMAX    = 18NEDOS   = 5001# 投影到 Ni d 轨道和 O p 轨道LOCPROJ = 1 : d : Pr  # 关键标签：为 Ni 创建 d 轨道投影

KPOINTS：金属体系至少 12×12×12，绝缘体 8×8×8，必须做收敛测试

执行VASP计算后，投影会存储在文件vaspgamma.h5和LOCPROJ中。

 

将VASP输出转换为TRIQS输入

准备一个配置文件plo.cfg：

[General]DOSMESH = -10103001
[Group 1]SHELLS = 1EWINDOW = -1010BANDS = 514
[Shell 1] # Ni d shellLSHELL = 2IONS = 1CORR = TrueTRANSFORM = 1.0  0.0  0.0  0.0  0.0            0.0  1.0  0.0  0.0  0.0            0.0  0.0  1.0  0.0  0.0            0.0  0.0  0.0  1.0  0.0            0.0  0.0  0.0  0.0  1.0

转换脚本converter.py

from triqs_dft_tools.converters.vasp import VaspConverterimport triqs_dft_tools.converters.plovasp.converter as plo_converter
# Generate and store PLOsplo_converter.generate_and_output_as_text('plo.cfg', vasp_dir='./')
# run the converterConverter = VaspConverter(filename = 'vasp')Converter.convert_dft_input()

直接运行python converter.py 运行python脚本，会生成可用于DMFT计算的输入文件vasp.h5。

 

 使用py4vasp 绘制 VASP 的 DFT DOS 与 VASP 投影轨道的 DOS 图：

import py4vaspimport matplotlib.pyplot as pltimport numpy as np
# get VASP DOScalc = py4vasp.Calculation.from_file('scf/vaspout.h5')dft_dos = calc.dos.to_dict()
# load projector (PLO) DOS from TRIQSplo_dos_Ni = np.loadtxt('scf/pdos_0_0.dat')
fig = plt.figure(figsize=(7,4), dpi=150)
plt.plot(dft_dos['energies'], dft_dos['total'], label=r'VASP total DOS')
# sum PLO DOS for all orbitalsplt.plot(plo_dos_Ni[:,0],np.sum(plo_dos_Ni[:,1:],axis=1), label=r'PLO Ni$_{d}$')
plt.xlim(-10,4)plt.ylim(0,12)plt.xlabel(r'$\epsilon - {\text{E}_\text{F}}$ (eV)')plt.ylabel(r'states/eV')plt.legend()plt.show()
图片

DMFT计算

准备dmft_config.toml文件：

[general]seedname = "vasp"          # 对应 vasp.h5jobname = "dmft_os"beta = 20                  # 电子温度：1/eV，对应室温一半U = 8.0                    # Hubbard U (eV)J = 1.0                    # Hund 交换 (eV)dc = true                  # 开启双计数修正n_iter_dmft = 20[solver]solver = "ctseg"          # CT-QMC 分段求解器length_cycle = 500n_warmup_cycles = 1e+4n_cycles_tot = 1e+7       # 总蒙特卡洛步数perform_tail_fit = true    # 高频尾巴拟合fit_max_moment = 4fit_min_w = 20fit_max_w = 30[advanced]dc_fixed_value = 59.0

运行 solid_dmft

mpirun -n 32 solid_dmft

计算完成后，会生成以jobname命名的结果目录，核心看两个文件：

  1. observables_imp0.dat：每轮迭代的关键可观测量，重点看G(beta/2)，该值趋近于 0，说明体系打开能隙，符合 NiO 莫特绝缘体的特性

  2. conv_obs0.dat：收敛判据，重点看δGimp，该值持续下降并稳定在 1e-4 以下，说明 DMFT 迭代收敛

收敛后的自能会保存在结果目录的vasp.h5文件中，为后续 CSC 计算提供初始猜解。

绘制结果

import numpy as npimport matplotlib.pyplot as pltfrom h5 import HDFArchivefrom triqs.gf import *
# ---------- 工具函数 ----------def mesh_to_np_arr(mesh):    from triqs.gf import MeshImTime, MeshReFreq, MeshImFreq    ifisinstance(mesh, MeshReFreq):        mesh_arr = np.linspace(mesh.w_min, mesh.w_max, len(mesh))    elif isinstance(mesh, MeshImFreq):        mesh_arr = np.linspace(mesh(mesh.first_index()).imag, mesh(mesh.last_index()).imag, len(mesh))    elif isinstance(mesh, MeshImTime):        mesh_arr = np.linspace(0, mesh.beta, len(mesh))    else:        raise AttributeError('input mesh must be either MeshReFreq, MeshImFreq, or MeshImTime')    return mesh_arr
def plot_sigma_iw(S_iw, ax, color, label='', marker='-o', subtract=True):    mesh = mesh_to_np_arr(S_iw.mesh)    mid = len(mesh)//2    if subtract:        ax[0].plot(mesh[mid:], S_iw.data[mid:,0,0].real-S_iw.data[-1,0,0].real, marker, color=color, label=label)    else:        ax[0].plot(mesh[mid:], S_iw.data[mid:,0,0].real, marker, color=color, label=label)    ax[1].plot(mesh[mid:], -1*S_iw.data[mid:,0,0].imag, marker, color=color, label=label)    ax[0].set_xlabel(r"$i \omega_n$ (eV)")    ax[1].set_xlabel(r"$i \omega_n$ (eV)")    ax[0].set_ylabel(r"$Re \Sigma (i \omega_n)$  (eV)")    ax[1].set_ylabel(r"$- Im \Sigma (i \omega_n)$  (eV)")
# ---------- 读取数据并画图 ----------with HDFArchive('./vasp.h5', 'r') as ar:    S_iw = ar['DMFT_results/last_iter/Sigma_freq_0']
fig, ax = plt.subplots(1, 2, figsize=(14, 5), dpi=100)
# 画出 up_0 和 up_2 轨道plot_sigma_iw(S_iw['up_0'], ax, color='C0', label='t2g')plot_sigma_iw(S_iw['up_2'], ax, color='C1', label='eg')
ax[0].legend()plt.tight_layout()plt.savefig('sigma_iw.png', dpi=150)plt.show()
图片

官网结果

图片

 

CSC DFT+DMFT循环

准备新的文件夹和NCAR、POSCAR、KPOINTS、POTCAR、plo.cfg、CHGCAR、WAVECAR和dmft_config.toml。 

在dmft_config.toml中添加

[dft]plo_cfg = "plo.cfg"dft_code = "vasp"projector_type = "plo"n_iter = 4n_cores = 32mpi_env = "default"mpi_exe = "/path/to/mpirun"dft_exec = "/path/to/bin/vasp_std"

在[general部分修改

csc = truen_iter_dmft_first = 2n_iter_dmft_per = 2n_iter_dmft = 20
load_sigma = truepath_to_sigma = "/path/to/dmft-os/U8.0-J1.0-beta20-qmc1e+7/vasp.h5"

在INCAR中添加

 ICHARG  = 5       AMIX    = 0.08    LSYNCH5 = True

用于读入hdf5文件，并设置环境允许同时读写 h5 文件

export HDF5_USE_FILE_LOCKING=FALSE

然后继续运行solid_dmft即可，可实现VASP的调用。

Ni-d 从 MaxEnt 投影的局部谱函数

图片