1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: ras/determinants.cc
4 // Copyright (C) 2013 Toru Shiozaki
5 //
6 // Author: Shane Parker <shane.parker@u.northwestern.edu>
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // This program is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program. If not, see <http://www.gnu.org/licenses/>.
23 //
24
25 #include <iostream>
26 #include <iomanip>
27 #include <cassert>
28
29 #include <src/util/combination.hpp>
30 #include <src/ci/ras/determinants.h>
31
32 using namespace std;
33 using namespace bagel;
34
RASDeterminants(const int norb1,const int norb2,const int norb3,const int nelea,const int neleb,const int max_holes,const int max_particles,const bool mute)35 RASDeterminants::RASDeterminants(const int norb1, const int norb2, const int norb3, const int nelea, const int neleb,
36 const int max_holes, const int max_particles, const bool mute)
37 : ras_{{norb1, norb2, norb3}}, max_holes_(max_holes), max_particles_(max_particles) {
38
39 if ( nelea < 0 || neleb < 0) throw runtime_error("nele < 0");
40
41 // check that nbit__ is big enough
42 if ( max_particles_ >= nbit__ || std::max(nelea,neleb) >= nbit__)
43 throw logic_error("Inappropriate value for \"nbit__\". Must be greater than nele AND max_particles");
44
45 // Construct spaces and with them, a list of strings
46 if (!mute) cout << " o Restricted Active Spaces:" << endl;
47 if (!mute) cout << " - RAS1 -> " << ras_[0] << endl;
48 if (!mute) cout << " - RAS2 -> " << ras_[1] << endl;
49 if (!mute) cout << " - RAS3 -> " << ras_[2] << endl << endl;
50
51 if (!mute) cout << " o Constructing all possible strings with up to " << max_holes_ << " holes and " << max_particles_ << " particles" << endl;
52 {
53 list<shared_ptr<const RASString>> alpha;
54 list<shared_ptr<const RASString>> beta;
55
56 for (int nholes = 0; nholes <= max_holes_; ++nholes) {
57 const int nele1 = norb1 - nholes;
58 for (int nparticles = 0; nparticles <= max_particles_; ++nparticles) {
59 const int nele3 = nparticles;
60 const int nele2a = nelea - (nele1 + nele3);
61 const int nele2b = neleb - (nele1 + nele3);
62
63 if (nele1 >= 0 && nele3 <= norb3) {
64 if (nele2a >= 0 && nele2a <= norb2) {
65 auto astring = make_shared<RASString>(nele1, norb1, nele2a, norb2, nele3, norb3);
66 if (astring->size() > 0)
67 alpha.push_back(astring);
68 }
69
70 if (nele2b >= 0 && nele2b <= norb2) {
71 auto bstring = make_shared<RASString>(nele1, norb1, nele2b, norb2, nele3, norb3);
72 if (bstring->size() > 0)
73 beta.push_back(bstring);
74 }
75 }
76 }
77 }
78 if (!alpha.empty())
79 alphaspaces_ = make_shared<CIStringSet<RASString>>(alpha);
80 if (!beta.empty())
81 betaspaces_ = make_shared<CIStringSet<RASString>>(beta);
82 }
83
84 size_ = 0;
85
86 if (alphaspaces_ && betaspaces_) {
87 if (alphaspaces_->size()*betaspaces_->size() > 0) {
88 if (!mute) cout << " - alpha strings: " << alphaspaces_->size() << endl;
89 if (!mute) cout << " - beta strings: " << betaspaces_->size() << endl << endl;
90
91 if (!mute) cout << " o Constructing alpha and beta displacement lists" << endl;
92 construct_phis_<0>(alphaspaces_, phia_, phia_ij_);
93 if (!mute) cout << " - alpha lists: " << phia_->size() << endl;
94 construct_phis_<1>(betaspaces_, phib_, phib_ij_);
95 if (!mute) cout << " - beta lists: " << phib_->size() << endl;
96
97 if (!mute) cout << " o Constructing pairs of allowed string spaces" << endl;
98
99 for (int nholes = 0; nholes <= max_holes_; ++nholes) {
100 for (int nha = nholes; nha >= 0; --nha) {
101 const int nhb = nholes - nha;
102 for (int npart = 0; npart <= max_particles_; ++npart) {
103 for (int npa = npart; npa >= 0; --npa) {
104 const int npb = npart - npa;
105
106 auto block = make_shared<const CIBlockInfo<RASString>>(space<0>(nha, npa), space<1>(nhb, npb), size_);
107 blockinfo_.push_back(block);
108 if (!block->empty()) size_ += block->size();
109 }
110 }
111 }
112 }
113 }
114 if (!mute) cout << " - size of restricted space: " << size_ << endl;
115 }
116 }
117
spin_adapt(const int spin,const bitset<nbit__> alpha,const bitset<nbit__> beta) const118 pair<vector<tuple<bitset<nbit__>, bitset<nbit__>, int>>, double> RASDeterminants::spin_adapt(const int spin, const bitset<nbit__> alpha, const bitset<nbit__> beta) const {
119 vector<tuple<bitset<nbit__>, bitset<nbit__>, int>> out;
120
121 // bit pattern for doubly occupied orbitals
122 bitset<nbit__> common = (alpha & beta);
123
124 bitset<nbit__> alpha_without_common = alpha ^ common;
125 bitset<nbit__> beta_without_common = beta ^ common;
126
127 // alpha pattern without highest spin orbitals
128 vector<int> salpha_array = bit_to_numbers(alpha_without_common);
129 vector<int> ualpha_array;
130 if (salpha_array.size() < spin) throw logic_error("Something is wrong? Determinants::spin_adapt");
131 for (int i = 0; i != spin; ++i) {
132 ualpha_array.push_back(salpha_array.back());
133 salpha_array.pop_back();
134 }
135 bitset<nbit__> salpha = numbers_to_bit(salpha_array);
136 bitset<nbit__> ualpha = numbers_to_bit(ualpha_array);
137 bitset<nbit__> common_plus_alpha(common.to_ulong() + ualpha.to_ulong());
138
139 // number of unpaired alpha orbitals (minus Ms)
140 const int nalpha = salpha.count();
141
142 // a vector of number that specify open orbitals
143 vector<int> open = bit_to_numbers(salpha^beta_without_common);
144 assert((salpha^beta_without_common) == (salpha|beta_without_common));
145
146 // take a linear combination to make a vector singlet coupled.
147 // TODO for the time being, we just leave Ms highest orbitals and singlet-couple other orbitals
148 assert(nalpha*2 == open.size());
149 int icnt = 0;
150 do {
151 bitset<nbit__> ialpha = common_plus_alpha;
152 bitset<nbit__> ibeta = common;
153 for (int i =0; i!=nalpha; ++i) ialpha.flip(open[i]);
154 for (int i=nalpha; i!=open.size(); ++i) ibeta.flip(open[i]);
155
156 // sign is always compensated by moving alpha to the left and beta to the right
157 // our convention is (aaaa)(bbbb)|0> due to the alpha-beta string algorithm
158 const double sign = 1.0;
159
160 out.push_back(make_tuple(ibeta, ialpha, sign));
161 ++icnt;
162 } while (boost::next_combination(open.begin(), open.begin()+nalpha, open.end()));
163
164 // scale to make the vector normalized
165 const double factor = 1.0/sqrt(static_cast<double>(icnt));
166 return {out, factor};
167 }
168