計算化學程式的實現：粒子數表象下波函式的表示
阿新 • • 發佈：2019-01-08
佔據數向量的表示
我們現在來考慮程式的實現。需要說明的是，這種計算並不是很輕鬆，所以一些簡單好用的程式語言是不適合拿來編寫量子化學的計算程式的。基本上，可以選擇的只有 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
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 自旋軌道，並且有
而總的
uint64_t 整數的個數為 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 的書寫方式左右顛倒，索引值是從左到右遞減的，最右側的軌道是第一個軌道。對於
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 位無符號整數的個數，即 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 purposes
std::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 uint64
std::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 string
int 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 string
int 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 string
int* 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 string
int* 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 level
int 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 operators
double 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 0
class 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 <--> end
std::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 <--> end
std::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: alpha
int 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: beta
int 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_a
int* SlaterDet::occ_list_a(int norb) const
{
return BasicOps::compute_occ_list(_strs, _nset, norb, nelec_a());
}
// return an array, size = nelec_b
int* SlaterDet::occ_list_b(int norb) const
{
return BasicOps::compute_occ_list((_strs + _nset), _nset, norb, nelec_b());
}
// return an array, size = norb - nelec_a
int* SlaterDet::vir_list_a(int norb) const
{
return BasicOps::compute_vir_list(_strs, _nset, norb, nelec_a());
}
// return an array, size = norb - nelec_b
int* 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: alpha
double 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: beta
double SlaterDet::compute_cre_des_sign_b(int p, int q)
{
return BasicOps::compute_cre_des_sign(p, q, (_strs + _nset), _nset);
}
// singly alpha-excited det
SlaterDet 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 det
SlaterDet 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)) ^