计算化学程序的实现：哈密顿矩阵元的计算

数学推导

本文是上一篇文章 http://blog.csdn.net/luozhen07/article/details/78701311 的继续。

量子化学计算的核心过程是构建体系的哈密顿矩阵并对角化，从而得到能量（特征值）和对应的波函数（特征矢）。波函数由一个或多个占据数矢量线性组合而成：前者称为“单行列式波函数”，譬如最基础的限制性 Hartree-Fock 自洽场（Hartree-Fock SCF）方法以及密度泛函（Density Functional Theory）方法得到的波函数；后者称为“多行列式波函数”，例如多组态自洽场（Multi-Configuration SCF, MCSCF）方法以及这里要讨论的组态相互作用（Configuration Interaction, CI）方法得到的波函数。用符号表示，为

|Ψ⟩=∑kckϕk

哈密顿矩阵元可以表示为

HMN=⟨ΨM|H^|ΨN⟩=∑i,jc∗M,icN,j⟨ϕi|H^|ϕj⟩

即，最终需要计算两个行列式之间的哈密顿量。考虑哈密顿算符的形式

H^=∑pqhpqEpq+12∑pqrsgpqrsepqrs

其中，单激发算符为

Epq=a†pαaqα+a†pβaqβ

双激发算符为

epqrs=a†pαa†rαasαaqα+a†pαa†rβasβaqα+a†pβa†rαasαaqβ+a†pβa†rβasβaqβ

并且 hpq 与 gpqrs 分别为单、双电子积分。考虑

⟨ϕi|ϕj⟩=δij

可见，只有当两个行列式的差异不超过两次激发（即，通过不超过两次电子移动，可以让左右两个行列式完全相同）时，哈密顿量才有非零值。这部分的具体讨论，可以参考下面的代码或者搜索“Slater-Condon 规则”。

程序实现

这里的程序主要是用于计算行列式（或者说“占据数矢量”）之间的哈密顿矩阵元 ⟨ϕi|H^|ϕj⟩。在上面的数学推导中，我们知道此计算还需要单、双电子积分。在这部分代码中，我们使用了名为 Integral 的类来保存积分，并且 Integral::h1e(p,q) 函数返回 p 轨道与 q 轨道之间的单电子积分，Integral::h2e(p,q,r,s) 函数返回 p、q、r 与 s 轨道之间的双电子积分。除了计算哈密顿矩阵元，其他计算主要有计算电子激发个数和激发的位置。

/* --- file: ci.h   --- author: LUO Zhen,               Institute of Theoretical and Computational Chemistry,               Nanjing University*/#ifndef CONFIGURATION_INTERACTION_FUNCTIONS_HEADER#define CONFIGURATION_INTERACTION_FUNCTIONS_HEADER#include "slaterdet.h"#include "integral.h"#include "basicops.h"// excitations are always from the second parameter (det_j) to the first one (det_i).namespace CI{    // number of excitations from det_j to det_i, alpha    int n_excitations_a(const SlaterDet& det_i, const SlaterDet& det_j);    // number of excitations from det_j to det_i, beta    int n_excitations_b(const SlaterDet& det_i, const SlaterDet& det_j);    // number of excitations from det_j to det_i, alpha + beta    int n_excitations(const SlaterDet& det_i, const SlaterDet& det_j);    // compute orbital indices for a single excitation: alpha, return int[2]    int* get_single_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j);    // compute orbital indices for a single excitation: beta, return int[2]    int* get_single_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j);    // compute orbital indices for a double excitation: alpha, return int[4]    int* get_double_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j);    // compute orbital indices for a double excitation: beta, return int[4]    int* get_double_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j);    // compute Hamiltonian matrix element H_ij = <det_i|H|det_j>    double H_ij(const SlaterDet& det_i, const SlaterDet& det_j, const Integral& integrals, int norb);}#endif // !CONFIGURATION_INTERACTION_FUNCTIONS_HEADER

这些函数的实现如下：

/* --- file: ci.cpp   --- author: LUO Zhen,               Institute of Theoretical and Computational Chemistry,               Nanjing University*/#include "ci.h"// number of excitations from det_j to det_i, alphaint CI::n_excitations_a(const SlaterDet& det_i, const SlaterDet& det_j){    return BasicOps::n_excitations(det_i.alpha(), det_j.alpha(), SlaterDet::nset());}// number of excitations from det_j to det_i, betaint CI::n_excitations_b(const SlaterDet& det_i, const SlaterDet& det_j){    return BasicOps::n_excitations(det_i.beta(), det_j.beta(), SlaterDet::nset());}// number of excitations from det_j to det_i, alpha + betaint CI::n_excitations(const SlaterDet& det_i, const SlaterDet& det_j){    return n_excitations_a(det_i, det_j) + n_excitations_b(det_i, det_j);}// compute orbital indices for a single excitation: alpha, return int[2]int* CI::get_single_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j){    return BasicOps::get_single_excitation(det_i.alpha(), det_j.alpha(), SlaterDet::nset());}// compute orbital indices for a single excitation: beta, return int[2]int* CI::get_single_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j){    return BasicOps::get_single_excitation(det_i.beta(), det_j.beta(), SlaterDet::nset());}// compute orbital indices for a double excitation: alpha, return int[4]int* CI::get_double_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j){    return BasicOps::get_double_excitation(det_i.alpha(), det_j.alpha(), SlaterDet::nset());}// compute orbital indices for a double excitation: beta, return int[4]int* CI::get_double_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j){    return BasicOps::get_double_excitation(det_i.beta(), det_j.beta(), SlaterDet::nset());}// compute Hamiltonian matrix element H_ij = <det_i|H|det_j>double CI::H_ij(const SlaterDet& det_i, const SlaterDet& det_j, const Integral& integrals, int norb){    double h1 = 0.0; // single particle    double h2 = 0.0; // double particle    int n_excit_a = CI::n_excitations_a(det_i, det_j);    int n_excit_b = CI::n_excitations_b(det_i, det_j);    if(n_excit_a + n_excit_b < 3) // ignore >= 3 excitations    {        int nelec_a = det_i.nelec_a();        int nelec_b = det_i.nelec_b();        if(n_excit_a == 0 && n_excit_b == 0) // i == j, diagonal term        {            int* occ_list_a = det_i.occ_list_a(norb);            int* occ_list_b = det_i.occ_list_b(norb);            // single particle            // h1 = SUM_(p in occ_a){h1e(p,p)} + SUM_(p in occ_b){h1e(p,p)}            for(int idx1 = 0; idx1 < nelec_a; ++idx1) // alpha            {                int p = occ_list_a[idx1];                h1 += integrals.h1e(p, p);            }            for(int idx1 = 0; idx1 < nelec_b; ++idx1) // beta            {                int p = occ_list_b[idx1];                h1 += integrals.h1e(p, p);            }            // double particle            // h2 = 0.5 * (SUM(p,r in occ_a){h2e(p,p,r,r) - h2e(p,r,r,p)}            //            +SUM(p,r in occ_b){h2e(p,p,r,r) - h2e(p,r,r,p)}            //            +SUM(p in occ_a, r in occ_b){h2e(p,p,r,r)}            //            +SUM(p in occ_b, r in occ_a){h2e(p,p,r,r)})            for(int idx1 = 0; idx1 < nelec_a; ++idx1)            {                int p = occ_list_a[idx1];                // SUM(p,r in occ_a){h2e(p,p,r,r) - h2e(p,r,r,p)}                for(int idx2 = 0; idx2 < nelec_a; ++idx2)                {                    int r = occ_list_a[idx2];                    h2 += (integrals.h2e(p, p, r, r) - integrals.h2e(p, r, r, p));                }                // SUM(p in occ_a, r in occ_b){h2e(p,p,r,r)}                for(int idx2 = 0; idx2 < nelec_b; ++idx2)                {                    int r = occ_list_b[idx2];                    h2 += integrals.h2e(p, p, r, r);                }            }            for(int idx1 = 0; idx1 < nelec_b; ++idx1)            {                int p = occ_list_b[idx1];                // SUM(p,r in occ_b){h2e(p,p,r,r) - h2e(p,r,r,p)}                for(int idx2 = 0; idx2 < nelec_b; ++idx2)                {                    int r = occ_list_b[idx2];                    h2 += (integrals.h2e(p, p, r, r) - integrals.h2e(p, r, r, p));                }                // SUM(p in occ_b, r in occ_a){h2e(p,p,r,r)}                for(int idx2 = 0; idx2 < nelec_a; ++idx2)                {                    int r = occ_list_a[idx2];                    h2 += integrals.h2e(p, p, r, r);                }            }            h2 *= 0.5;            delete[] occ_list_b;            delete[] occ_list_a;        }        else if(n_excit_a == 1 && n_excit_b == 0) // a -> a        {            int* ia = CI::get_single_excitation_a(det_i, det_j);            int* occ_list_a = det_i.occ_list_a(norb);            int* occ_list_b = det_i.occ_list_b(norb);            // single particle            // h1 = sign(p,q) * h1(p,q)            int p = ia[1];            int q = ia[0];            double sign = BasicOps::compute_cre_des_sign(p, q, det_j.alpha(), SlaterDet::nset());            h1 = sign * integrals.h1e(p, q);            // double particle            // h2 = sign(p,q) * (SUM_(r in occ_a){h2e(p,q,r,r) - h2e(p,r,r,q)}            //                  +SUM_(r in occ_b){h2e(p,q,r,r)})            for(int idx1 = 0; idx1 < nelec_a; ++idx1)            {                int r = occ_list_a[idx1];                h2 += (integrals.h2e(p, q, r, r) - integrals.h2e(p, r, r, q));            }            for(int idx2 = 0; idx2 < nelec_b; ++idx2)            {                int r = occ_list_b[idx2];                h2 += integrals.h2e(p, q, r, r);            }            h2 *= sign;            delete[] occ_list_b;            delete[] occ_list_a;            delete[] ia;        }        else if(n_excit_a == 0 && n_excit_b == 1) // b -> b        {            int* ia = CI::get_single_excitation_b(det_i, det_j);            int* occ_list_a = det_i.occ_list_a(norb);            int* occ_list_b = det_i.occ_list_b(norb);            // single particle            // h1 = sign(p,q) * h1(p,q)            int p = ia[0];            int q = ia[1];            double sign = BasicOps::compute_cre_des_sign(p, q, det_j.beta(), SlaterDet::nset());            h1 = sign * integrals.h1e(p, q);            // double particle            // h2 = sign(p,q) * (SUM_(r in occ_b){h2e(p,q,r,r) - h2e(p,r,r,q)}            //                  +SUM_(r in occ_a){h2e(p,q,r,r)})            for(int idx1 = 0; idx1 < nelec_b; ++idx1)            {                int r = occ_list_b[idx1];                h2 += (integrals.h2e(p, q, r, r) - integrals.h2e(p, r, r, q));            }            for(int idx2 = 0; idx2 < nelec_a; ++idx2)            {                int r = occ_list_a[idx2];                h2 += integrals.h2e(p, q, r, r);            }            h2 *= sign;            delete[] occ_list_b;            delete[] occ_list_a;            delete[] ia;        }        else if(n_excit_a == 2 && n_excit_b == 0) // a,a -> a,a        {            int* ijab = CI::get_double_excitation_a(det_i, det_j);            int p = ijab[0];            int r = ijab[1];            int s = ijab[2];            int q = ijab[3];            // make sure s < q and p < r            if(q < s) std::swap(q, s);            if(r < p) std::swap(p, r);            double sign = BasicOps::compute_cre_des_sign(p, r, det_i.alpha(), SlaterDet::nset());            sign *= BasicOps::compute_cre_des_sign(s, q, det_j.alpha(), SlaterDet::nset());            h2 += sign * (integrals.h2e(p, s, r, q) - integrals.h2e(p, q, r, s));            delete[] ijab;        }        else if(n_excit_a == 0 && n_excit_b == 2) // b,b -> b,b        {            int* ijab = CI::get_double_excitation_b(det_i, det_j);            int p = ijab[0];            int r = ijab[1];            int s = ijab[2];            int q = ijab[3];            // make sure s < q and p < r            if(q < s) std::swap(q, s);            if(r < p) std::swap(p, r);            double sign = BasicOps::compute_cre_des_sign(p, r, det_i.beta(), SlaterDet::nset());            sign *= BasicOps::compute_cre_des_sign(s, q, det_j.beta(), SlaterDet::nset());            h2 += sign * (integrals.h2e(p, s, r, q) - integrals.h2e(p, q, r, s));            delete[] ijab;        }        else // a,b -> a,b        {            int* ia = CI::get_single_excitation_a(det_i, det_j);            int* jb = CI::get_single_excitation_b(det_i, det_j);            int p = ia[0];            int q = ia[1];            int r = jb[0];            int s = jb[1];            double sign = BasicOps::compute_cre_des_sign(q, p, det_j.alpha(), SlaterDet::nset());            sign *= BasicOps::compute_cre_des_sign(s, r, det_j.beta(), SlaterDet::nset());            h2 += (sign * integrals.h2e(p, q, r, s));            delete[] jb;            delete[] ia;        }    }    return h1 + h2;}

多组态波函数哈密顿矩阵元的计算

由之前的数学推导可知，对于多组态波函数哈密顿矩阵元，只需要对左矢和右矢分别遍历，计算每一组行列式的矩阵元与对应系数的乘积并求加和。定义如下的行列式集合类

class DetSet{public:    // 省略的代码 ...private:    std::vector<double> _coeffs;    std::vector<SlaterDet> _dets;};

则哈密顿矩阵元可以通过如下方式计算

double H_ij(const DetSet& set_i, const DetSet& set_j, const Integral& integrals, int norb){    double h = 0.0;#pragma omp parallel for reduction(+:h)    for(int idx1 = 0; idx1 < set_i.size(); ++idx1)    {        for(int idx2 = 0; idx2 < set_j.size(); ++idx2)        {            double coeff_product = set_i.coeff(idx1) * set_j.coeff(idx2);            h += (coeff_product * CI::H_ij(set_i.det(idx1), set_j.det(idx2), integrals, norb));        }    }    return h;}