计算化学程序的实现：粒子数表象下波函数的表示

占据数矢量的表示

我们现在来考虑程序的实现。需要说明的是，这种计算并不是很轻松，所以一些简单好用的编程语言是不适合拿来编写量子化学的计算程序的。基本上，可以选择的只有 Fortran 以及 C/C++。在此文中使用 C++ 来实现，并且我会用到一些不常见的手段来提高运行速度。如果你觉得阅读这部分有困难，那么可以读完此部分之后再去读一些关于 C/C++ 语言的位运算之类的文章。

程序的第一步从构建波函数开始，而波函数则是由一个或多个占据数矢量构成。既然每个占据数矢量都可以拆分成 α 和 β 两部分，每部分都是一个由“占据”或者“非占据”两种状态组成的串，串上的每一个位置表示一个自旋轨道。那么很自然地，这种形式可以使用由 0 和 1 组成的序列来描述。但是应该使用怎样的数据类型呢？

最直接的想法是，既然每个位置的状态只有 0 和 1 两种，那么使用一个 bit（位） 来表示一个位置是最好的选择。不过计算机本身进行操作的最小内存单位是 byte（字节），等于 8 个 bit。计算机中所有的数据类型都至少是 1 byte 大小，即使是布尔类型也不例外。并且，任意长度的 bit 串是不可行的，这里存在内存对齐的问题。目前直接实现的 bit 串主要有以下两个：

• std::bitset<n>，需要在编译时就指定其大小。除非你想一遍遍地对不同的体系重新修改和编译程序，否则这并不是一个好的解决办法。
• std::vector<bool>，对于某些实现而言，这个类型是 bit 数组而非布尔类型的数组。但是此类型的访问速度很慢，无论是什么情况都不应该使用（除非内存占用确实很重要，重要到可以让你放弃性能）。

使用 bit 串存在的另一个问题是如何对其进行循环遍历，或者，至少需要一种方法能够高效地判断给定位置的状态。

另一个方法是使用数组，事实上很多现有的计算程序就是采用这种方式。数组形式的好处是可以方便直观地进行循环遍历。考虑到 int 类型的大小是 4 byte（32 位整数，对于 x86 以及 x86_64 处理器而言这是默认的情况），而实际需要表示的只有 0 和 1 两个数，所以即使真的要使用数组，也应该使用 int8_t（8 位整数）或者 char 这种最小的整数类型。但无论使用什么数据类型，数组形式都会浪费大量的内存空间，不适用于行列式数目非常多的情况。

这里我们使用一种折衷的实现方式：使用 64 位无符号整数 uint64_t 类型来表示一个长度为 64 的 0/1 序列。这种实现方式与文章 https://arxiv.org/abs/1311.6244 中的实现方式类似。对于有 norb 个轨道（此处指的是空间轨道，即 alpha 自旋轨道与 beta 自旋轨道都各有 norb 个）组成的体系，如果 norb 不大于 64，则需要 2 个 uint64_t 整数来表示这个占据数矢量，即 alpha 与 beta 各一个；否则，若 norb 大于 64，则需要多个 uint64_t 整数。简单来说，因为每 64 个轨道为一组，我们需要 nset 组 alpha 自旋轨道才能完全表示总共 norb 个 alpha 自旋轨道，并且有

nset=(norb+63)/64

而总的 uint64_t 整数的个数为 2×nset。这样我们就可以定义如下的类：

class SlaterDet{public:    // 省略的代码 ...private:    static int _nset; // 此处为静态成员变量    std::uint64_t* _strs; // 注意此处数组的长度为 2*_nset};

注意这里只使用了一个总长度为 2 * nset 的数组，其中前 nset 个是 alpha 序列，后 nset 个是 beta 序列。每个数都有 64 个 bit，表示 64 个自旋轨道，并且 1 表示有电子占据而 0 表示无电子占据。考虑到 x86 处理器是 little-endian，这里的每一个 uint64_t 数也采用从低位到高位逐位递进的方式。注意此处的“逐位”，我们的操作对象不再是 byte 而是 bit。例如，对序列 2ud0，其中 2 表示 alpha/beta 双占据，u 表示 alpha（或 spin up）占据，d 表示 beta（或 spin down）占据，0 表示此空间轨道无电子占据，这样其 alpha 序列为 1100 而 beta 序列为 1010，在内存中其形式为

alpha: ...0011 # 索引为 [63, 62, ..., 3, 2, 1, 0]，... 表示 60 个零beta: ...0101 # 索引为 [63, 62, ..., 3, 2, 1, 0]，... 表示 60 个零

简单来说就是与 2ud0 的书写方式左右颠倒，索引值是从左到右递减的，最右侧的轨道是第一个轨道。对于 nset>1 的情况也是类似的，只不过将整个序列切成了多个分块。例如，对于 norb=150，有 nset=3，此时索引值为

alpha[63, ..., 1, 0], alpha[127, ..., 64], alpha[191, ..., 128], beta[63, ..., 1, 0], beta[127, ..., 64], beta[191, ..., 128]

因为需要内存对齐，所以除非 norb 恰好是 64 的整数倍，否则总是会有多余的部分。这个例子中有 150 个轨道，所以索引值 >= 150 的位置全都忽略。显然，整个过程中需要清楚地知道轨道数 norb（由轨道数 norb 可以计算出使用的 64 位无符号整数的个数，即 2×nset）。考虑到实际计算中轨道数是确定的（所以 nset 也是确定的，开始计算之后就不会再改变），而且遍历过程中更常使用的是 nset，所以上面的实现中将 nset 作为静态成员变量，由所有的 SlaterDet 对象共享。必须说明的是，在创建第一个 SlaterDet 对象之前，必须先计算并设定好 nset 的数值。

对 64 位无符号整数的操作

在进一步介绍占据数矢量的操作之前，必须先说明对于 64 位无符号整数 std::uint64_t 类型的一些操作。注意，此类型需要使用 cstdint 头文件并且启用编译器的 C++11 支持。再次强调，此处的 std::uint64_t 元素是作为一个长度为 64 的 0/1 序列来使用的，并不是作为整数本身。

需要对占据数矢量进行的操作主要有计算电子数、确定电子占据轨道的位置、确定空轨道的位置这三种，在此基础上还可以组合出计算电子激发的个数、计算电子激发的位置等操作。可以定义如下的函数

/* --- file: basicops.h   --- author: LUO Zhen,               Institute of Theoretical and Computational Chemistry,               Nanjing University*/#ifndef BASIC_OPERATIONS_FUNCTIONS_HEADER#define BASIC_OPERATIONS_FUNCTIONS_HEADER#include <cstdint>#include <string>namespace BasicOps{    // Function to print a 64-bits string as a string for debug purposes    std::string int2bin(std::uint64_t x);    // Function to transform a 0/1 string to an uint64    std::uint64_t bin2int(const std::string& str_1010);    // Compute the number of set bits (= 1) in a 64-bits string    int popcount(std::uint64_t x);    // Compute the number of trailing (right) zeros in a 64-bits string    // eg. trailz(1) == 0, trailz(0xf000000000000000) == 60    int trailz(std::uint64_t x);    // Compute a list of occupied orbitals for a given string    int* compute_occ_list(const std::uint64_t* str, int nset, int norb, int nelec);    // Compute a list of occupied orbitals for a given string    int* compute_vir_list(const std::uint64_t* str, int nset, int norb, int nelec);    // Compare two strings and compute excitation level    int n_excitations(const std::uint64_t* str1, const std::uint64_t* str2, int nset);    // Compute orbital indices for a single excitation     int* get_single_excitation(const std::uint64_t* str1, const std::uint64_t* str2, int nset);    // Compute orbital indices for a double excitation     int* get_double_excitation(const std::uint64_t* str1, const std::uint64_t* str2, int nset);    // Compute sign for a pair of creation and desctruction operators    double compute_cre_des_sign(int p, int q, const std::uint64_t* str, int nset);}#endif // !BASIC_OPERATIONS_FUNCTIONS_HEADER

函数实现如下（假设上面的代码保存为头文件 basicops.h）：

/* --- file: basicops.cpp   --- author: LUO Zhen,               Institute of Theoretical and Computational Chemistry,               Nanjing University*/#include "basicops.h"// Function to print a 64-bits string as a string for debug purposesstd::string BasicOps::int2bin(std::uint64_t x){    std::string str(64, '0');    for(int pos = 0; pos < 64; ++pos)        if(x & (1ULL << pos)) str[63 - pos] = '1';    return str;}// Function to transform a 0/1 string to an uint64std::uint64_t BasicOps::bin2int(const std::string& str_1010){    std::uint64_t i = 0ULL;    for(int pos = 0; pos < str_1010.length(); ++pos)    {        i <<= 1;        if(str_1010[pos] == '1') i |= 1ULL;    }    return i;}// Compute the number of set bits (= 1) in a 64-bits stringint BasicOps::popcount(std::uint64_t x){#ifdef __INTEL_COMPILER    return _popcnt64(x);#else    const std::uint64_t m1 = 0x5555555555555555; //binary: 0101...    const std::uint64_t m2 = 0x3333333333333333; //binary: 00110011..    const std::uint64_t m4 = 0x0f0f0f0f0f0f0f0f; //binary:  4 zeros,  4 ones ...    const std::uint64_t m8 = 0x00ff00ff00ff00ff; //binary:  8 zeros,  8 ones ...    const std::uint64_t m16 = 0x0000ffff0000ffff; //binary: 16 zeros, 16 ones ...    const std::uint64_t m32 = 0x00000000ffffffff; //binary: 32 zeros, 32 ones    x = (x & m1) + ((x >> 1) & m1); //put count of each  2 bits into those  2 bits     x = (x & m2) + ((x >> 2) & m2); //put count of each  4 bits into those  4 bits     x = (x & m4) + ((x >> 4) & m4); //put count of each  8 bits into those  8 bits     x = (x & m8) + ((x >> 8) & m8); //put count of each 16 bits into those 16 bits     x = (x & m16) + ((x >> 16) & m16); //put count of each 32 bits into those 32 bits     x = (x & m32) + ((x >> 32) & m32); //put count of each 64 bits into those 64 bits     return x;#endif // INTEL_COMPILER}// Compute the number of trailing zeros in a 64-bits stringint BasicOps::trailz(std::uint64_t x){    int c = 64;    // Trick to unset all bits except the first one    x &= -(std::int64_t)x;    if(x) c--;    if(x & 0x00000000ffffffff) c -= 32;    if(x & 0x0000ffff0000ffff) c -= 16;    if(x & 0x00ff00ff00ff00ff) c -= 8;    if(x & 0x0f0f0f0f0f0f0f0f) c -= 4;    if(x & 0x3333333333333333) c -= 2;    if(x & 0x5555555555555555) c -= 1;    return c;}// Compute a list of occupied orbitals for a given stringint* BasicOps::compute_occ_list(const std::uint64_t* str, int nset, int norb, int nelec){    int* occ = new int[nelec];    int off = 0;    int occ_ind = 0;    for(int k = 0; k < nset; ++k)    {        int i_max = ((norb - off) < 64 ? (norb - off) : 64);        for(int i = 0; i < i_max; ++i)        {            int i_occ = (str[k] >> i) & 1;            if(i_occ)            {                occ[occ_ind] = i + off;                ++occ_ind;            }        }        off += 64;    }    return occ;}// Compute a list of occupied orbitals for a given stringint* BasicOps::compute_vir_list(const std::uint64_t* str, int nset, int norb, int nelec){    int* vir = new int[norb - nelec];    int off = 0;    int vir_ind = 0;    for(int k = 0; k < nset; ++k)    {        int i_max = ((norb - off) < 64 ? (norb - off) : 64);        for(int i = 0; i < i_max; ++i)        {            int i_occ = (str[k] >> i) & 1;            if(!i_occ)            {                vir[vir_ind] = i + off;                ++vir_ind;            }        }        off += 64;    }    return vir;}// Compare two strings and compute excitation levelint BasicOps::n_excitations(const std::uint64_t* str1, const std::uint64_t* str2, int nset){    int d = 0;    for(int p = 0; p < nset; ++p)        d += BasicOps::popcount(str1[p] ^ str2[p]);    return d / 2;}// Compute orbital indices for a single excitation int* BasicOps::get_single_excitation(const std::uint64_t* str1, const std::uint64_t* str2, int nset){    int* ia = new int[2];    for(int p = 0; p < nset; ++p)    {        std::uint64_t str_tmp = str1[p] ^ str2[p];        std::uint64_t str_particle = str_tmp & str2[p];        std::uint64_t str_hole = str_tmp & str1[p];        if(BasicOps::popcount(str_particle) == 1)            ia[1] = BasicOps::trailz(str_particle) + 64 * p;        if(BasicOps::popcount(str_hole) == 1)            ia[0] = BasicOps::trailz(str_hole) + 64 * p;    }    return ia;}// Compute orbital indices for a double excitation int* BasicOps::get_double_excitation(const std::uint64_t* str1, const std::uint64_t* str2, int nset){    int* ijab = new int[4];    int particle_ind = 2;    int hole_ind = 0;    for(int p = 0; p < nset; ++p)    {        std::uint64_t str_tmp = str1[p] ^ str2[p];        std::uint64_t str_particle = str_tmp & str2[p];        std::uint64_t str_hole = str_tmp & str1[p];        int n_particle = BasicOps::popcount(str_particle);        int n_hole = BasicOps::popcount(str_hole);        if(n_particle == 1)        {            ijab[particle_ind] = BasicOps::trailz(str_particle) + 64 * p;            ++particle_ind;        }        else if(n_particle == 2)        {            int a = BasicOps::trailz(str_particle);            ijab[2] = a + 64 * p;            str_particle &= ~(1ULL << a);            int b = BasicOps::trailz(str_particle);            ijab[3] = b + 64 * p;        }        if(n_hole == 1)        {            ijab[hole_ind] = BasicOps::trailz(str_hole) + 64 * p;            ++hole_ind;        }        else if(n_hole == 2)        {            int i = BasicOps::trailz(str_hole);            ijab[0] = i + 64 * p;            str_hole &= ~(1ULL << i);            int j = BasicOps::trailz(str_hole);            ijab[1] = j + 64 * p;        }    }    return ijab;}// Compute sign for a pair of creation and desctruction operatorsdouble BasicOps::compute_cre_des_sign(int p, int q, const std::uint64_t* str, int nset){    double sign;    int nperm;    int i;    int pg = p / 64;    int qg = q / 64;    int pb = p % 64;    int qb = q % 64;    if(pg > qg)    {        nperm = 0;        for(i = qg + 1; i < pg; ++i)            nperm += BasicOps::popcount(str[i]);        nperm += BasicOps::popcount(str[pg] << (64 - pb));        nperm += BasicOps::popcount(str[qg] >> (qb + 1));    }    else if(pg < qg)    {        nperm = 0;        for(i = pg + 1; i < qg; ++i)            nperm += BasicOps::popcount(str[i]);        nperm += BasicOps::popcount(str[qg] << (64 - qb));        nperm += BasicOps::popcount(str[pg] >> (pb + 1));    }    else    {        std::uint64_t mask;        if(p > q)            mask = (1ULL << pb) - (1ULL << (qb + 1));        else            mask = (1ULL << qb) - (1ULL << (pb + 1));        nperm = BasicOps::popcount(str[pg] & mask);    }    if(nperm % 2 == 1) // #elec between p and p is odd        sign = -1.0;    else // #elec between p and p is even        sign = 1.0;    return sign;}

这里使用了大量的位操作，好处是降低了循环的次数，并且充分利用了 CPU 的 Cache。如果对位操作不熟悉，可以参考文章 http://www.cppblog.com/biao/archive/2012/03/20/168357.html 或者查阅 C/C++ 语言的相关书籍。另外需要注意的是，对于大规模的计算程序，使用 new/delete 或者 malloc/free 来实现内存的分配与回收可能会影响并行的性能，这种情况下建议使用 intel Thread Building Blocks 提供的 scalable_malloc/scalable_free。

占据数矢量的最终实现

基于以上的分析和代码，我们可以实现最终的 SlaterDet 类：

/* --- file: slaterdet.h   --- author: LUO Zhen,               Institute of Theoretical and Computational Chemistry,               Nanjing University*/#ifndef SLATER_DETERMINANT_CLASS_HEADER#define SLATER_DETERMINANT_CLASS_HEADER#include <cstdint>#include <cstring>#include <iostream>#include <string>#include "basicops.h"// creation site: p, r// annihilation site: q, s// orbital index starts from 0class SlaterDet{public:    SlaterDet();    // construct from string "2ud0"    SlaterDet(const std::string& str_2ud0);    SlaterDet(const SlaterDet& det);    SlaterDet(SlaterDet&& det);    SlaterDet& operator=(const SlaterDet& det);    SlaterDet& operator=(SlaterDet&& det);    ~SlaterDet();    void swap(SlaterDet& det);    // number of alpha/beta string    static int nset();    // set the static parameter: _nset    static void set_nset(int n);    std::uint64_t* alpha() const;    std::uint64_t* beta() const;    friend bool operator==(const SlaterDet& det1, const SlaterDet& det2);    friend bool operator!=(const SlaterDet& det1, const SlaterDet& det2);    friend bool operator<(const SlaterDet& det1, const SlaterDet& det2);    // alpha string, start <--> end    std::string String_a(int norb) const;    // beta string, start <--> end    std::string String_b(int norb) const;    // count the number of electrons before site pos: alpha    int count_elec_a(int pos) const;    // count the number of electrons before site pos: beta    int count_elec_b(int pos) const;    // number of alpha electrons    int nelec_a() const;    // number of beta electrons    int nelec_b() const;    // number of electrons, neleca + nelecb    int nelec() const;    // whether stra[pos] is occupied    bool occupied_a(int pos) const;    // whether strb[pos] is occupied    bool occupied_b(int pos) const;    // occupied orbitals: alpha. return int[nelec_a], remember to delete[]!    int* occ_list_a(int norb) const;    // occupied orbitals: beta. return int[nelec_b], remember to delete[]!    int* occ_list_b(int norb) const;    // empty orbitals: alpha. return int[norb - nelec_a], remember to delete[]!    int* vir_list_a(int norb) const;    // empty orbitals: beta. return int[norb - nelec_b], remember to delete[]!    int* vir_list_b(int norb) const;    // compute sign for a pair of creation and desctruction operators: alpha    double compute_cre_des_sign_a(int p, int q);    // compute sign for a pair of creation and desctruction operators: beta    double compute_cre_des_sign_b(int p, int q);    // singly alpha-excited det    SlaterDet Excited_a_pq(int p, int q) const;    // doubly alpha-excited det    SlaterDet Excited_a_pqrs(int p, int q, int r, int s) const;    // singly beta-excited det    SlaterDet Excited_b_pq(int p, int q) const;    // doubly beta-excited det    SlaterDet Excited_b_pqrs(int p, int q, int r, int s) const;private:    static int _nset;    std::uint64_t* _strs;};#endif // !SLATER_DETERMINANT_CLASS_HEADER

如之前所述，为了尽量减少内存占用并且减少计算次数，这里将 nset 作为 SlaterDet 类的静态成员变量，轨道数 norb 作为一些成员函数的参数。因为实际的程序中，轨道数是开始计算时就已经确定的，没有必要将其作为成员变量（或者作为静态成员变量也可以。但是实际上 norb 决定了 nset，没必要将两者都作为静态成员；nset 在成员函数中经常被使用，将其设定为成员变量是更好的选择，这样可以避免很多重复的计算）。这个类的实现是

/* --- file: slaterdet.cpp   --- author: LUO Zhen,               Institute of Theoretical and Computational Chemistry,               Nanjing University*/#include "slaterdet.h"int SlaterDet::_nset = 0;SlaterDet::SlaterDet(){    _strs = nullptr;}SlaterDet::SlaterDet(const std::string& str_2ud0){    _strs = new std::uint64_t[_nset * 2];    std::memset(_strs, 0, (sizeof(std::uint64_t) * _nset * 2));    std::uint64_t* stra;    std::uint64_t* strb;    int iset = (str_2ud0.length() + 63) / 64;    for(int i = 0; i < iset; ++i)    {        stra = _strs + i;        strb = _strs + _nset + i;        std::string subStr = str_2ud0.substr(i * 64, 64);        for(auto i = subStr.crbegin(); i != subStr.crend(); ++i)        {            *stra <<= 1;            *strb <<= 1;            switch(*i)            {            case '2':                *stra |= 1ULL;                *strb |= 1ULL;                break;            case 'u':                *stra |= 1ULL;                break;            case 'd':                *strb |= 1ULL;                break;            default:                break;            }        }    }}SlaterDet::SlaterDet(const SlaterDet& det){    _strs = new std::uint64_t[_nset * 2];    std::memcpy(_strs, det._strs, (sizeof(std::uint64_t) * _nset * 2));}SlaterDet::SlaterDet(SlaterDet&& det){    _strs = det._strs;    det._strs = nullptr;}SlaterDet& SlaterDet::operator=(const SlaterDet& det){    if(this != &det)    {        _strs = new std::uint64_t[_nset * 2];        std::memcpy(_strs, det._strs, (sizeof(std::uint64_t) * _nset * 2));    }    return *this;}SlaterDet& SlaterDet::operator=(SlaterDet&& det){    if(this != &det)    {        _strs = det._strs;        det._strs = nullptr;    }    return *this;}SlaterDet::~SlaterDet(){    delete[] _strs;}void SlaterDet::swap(SlaterDet& det){    std::swap(_strs, det._strs);}int SlaterDet::nset(){    return _nset;}void SlaterDet::set_nset(int n){    _nset = n;}std::uint64_t* SlaterDet::alpha() const{    return _strs;}std::uint64_t* SlaterDet::beta() const{    return (_strs + _nset);}bool operator==(const SlaterDet& det1, const SlaterDet& det2){    for(int i = 0; i < (SlaterDet::nset() * 2); ++i)    {        if(det1._strs[i] != det2._strs[i])            return false;    }    return true;}bool operator!=(const SlaterDet& det1, const SlaterDet& det2){    return !(det1 == det2);}bool operator<(const SlaterDet& det1, const SlaterDet& det2){    for(int i = 0; i < det1._nset; ++i)    {        if(det1._strs[i] < det2._strs[i])            return true;        else if(det1._strs[i] > det2._strs[i])            return false;    }    return false; // det1 == det2}// alpha string, start <--> endstd::string SlaterDet::String_a(int norb) const{    std::string stra(norb, '0');    int index;    for(int iset = 0; iset < _nset; ++iset)    {        std::uint64_t a = *(_strs + iset);        for(int i = 0; i < 64; ++i)        {            index = iset * 64 + i;            if(index < norb)            {                if(a & (1ULL << i))                    stra[index] = '1';            }            else            {                break;            }        }    }    return stra;}// beta string, start <--> endstd::string SlaterDet::String_b(int norb) const{    std::string strb(norb, '0');    int index;    for(int iset = 0; iset < _nset; ++iset)    {        std::uint64_t b = *(_strs + _nset + iset);        for(int i = 0; i < 64; ++i)        {            index = iset * 64 + i;            if(index < norb)            {                if(b & (1ULL << i))                    strb[index] = '1';            }            else            {                break;            }        }    }    return strb;}// count the number of electrons before site pos: alphaint SlaterDet::count_elec_a(int pos) const{    int iset = (pos + 63) / 64;    int off = pos % 64;    std::uint64_t tail = *(alpha() + iset - 1);    int n = 0;    for(int i = 0; i < iset - 1; ++i)        n += BasicOps::popcount(*(alpha() + i));    n += BasicOps::popcount(tail << (64 - off));    return n;}// count the number of electrons before site pos: betaint SlaterDet::count_elec_b(int pos) const{    int iset = (pos + 63) / 64;    int off = pos % 64;    std::uint64_t tail = *(beta() + iset - 1);    int n = 0;    for(int i = 0; i < iset - 1; ++i)        n += BasicOps::popcount(*(beta() + i));    n += BasicOps::popcount(tail << (64 - off));    return n;}int SlaterDet::nelec_a() const{    int _nelec_a = 0;    for(int i = 0; i < _nset; ++i)        _nelec_a += BasicOps::popcount(*(_strs + i));    return _nelec_a;}int SlaterDet::nelec_b() const{    int _nelec_b = 0;    for(int i = 0; i < _nset; ++i)        _nelec_b += BasicOps::popcount(*(_strs + _nset + i));    return _nelec_b;}int SlaterDet::nelec() const{    return nelec_a() + nelec_b();}bool SlaterDet::occupied_a(int pos) const{    return ((*(alpha() + (pos / 64))) & (1ULL << (pos % 64)));}bool SlaterDet::occupied_b(int pos) const{    return ((*(beta() + (pos / 64))) & (1ULL << (pos % 64)));}// return an array, size = nelec_aint* SlaterDet::occ_list_a(int norb) const{    return BasicOps::compute_occ_list(_strs, _nset, norb, nelec_a());}// return an array, size = nelec_bint* SlaterDet::occ_list_b(int norb) const{    return BasicOps::compute_occ_list((_strs + _nset), _nset, norb, nelec_b());}// return an array, size = norb - nelec_aint* SlaterDet::vir_list_a(int norb) const{    return BasicOps::compute_vir_list(_strs, _nset, norb, nelec_a());}// return an array, size = norb - nelec_bint* SlaterDet::vir_list_b(int norb) const{    return BasicOps::compute_vir_list((_strs + _nset), _nset, norb, nelec_b());}// compute sign for a pair of creation and desctruction operators: alphadouble SlaterDet::compute_cre_des_sign_a(int p, int q){    return BasicOps::compute_cre_des_sign(p, q, _strs, _nset);}// compute sign for a pair of creation and desctruction operators: betadouble SlaterDet::compute_cre_des_sign_b(int p, int q){    return BasicOps::compute_cre_des_sign(p, q, (_strs + _nset), _nset);}// singly alpha-excited detSlaterDet SlaterDet::Excited_a_pq(int p, int q) const{    SlaterDet det(*this);    *(det._strs + (p / 64)) ^= (1ULL << (p % 64));    *(det._strs + (q / 64)) ^= (1ULL << (q % 64));    return det;}// doubly alpha-excited detSlaterDet SlaterDet::Excited_a_pqrs(int p, int q, int r, int s) const{    SlaterDet det(*this);    *(det._strs + (p / 64)) ^= (1ULL << (p % 64));    *(det._strs + (q / 64)) ^= (1ULL << (q % 64));    *(det._strs + (r / 64)) ^= (1ULL << (r % 64));    *(det._strs + (s / 64)) ^= (1ULL << (s % 64));    return det;}// singly beta-excited detSlaterDet SlaterDet::Excited_b_pq(int p, int q) const{    SlaterDet det(*this);    *(det._strs + _nset + (p / 64)) ^= (1ULL << (p % 64));    *(det._strs + _nset + (q / 64)) ^= (1ULL << (q % 64));    return det;}// doubly beta-excited detSlaterDet SlaterDet::Excited_b_pqrs(int p, int q, int r, int s) const{    SlaterDet det(*this);    *(det._strs + _nset + (p / 64)) ^= (1ULL << (p % 64));    *(det._strs + _nset + (q / 64)) ^= (1ULL << (q % 64));    *(det._strs + _nset + (r / 64)) ^= (1ULL << (r % 64));    *(det._strs + _nset + (s / 64)) ^= (1ULL << (s % 64));    return det;}

产生算符与湮灭算符都可以用异或操作来实现：在位置 p 上的产生和湮灭算符都是只在 p 位置是 1 而其他位置全为 0 的 std::uint64_t 整数，将这个数与原序列进行异或操作。若原序列 p 位置为 0，则异或操作将此位置变为 1，相当于在此位置产生一个电子；反之，若原序列 p 位置为 1，异或操作将此位置变为 0，相当于在此位置湮灭一个电子。关于产生和湮灭算符作用之后的符号（相位）的计算，可以参考 Trygve Helgaker 等人所著的 MOLECULAR ELECTRONIC-STRUCTURE THEORY 一书。