1 /*
2
3 Copyright (C) 2012 Lucas K. Wagner
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18
19 */
20
21 #include "clark_updates.h"
build_excitation_list(Array3<Array1<int>> & occupation,int f)22 void Excitation_list::build_excitation_list(Array3 <Array1 <int> > & occupation,int f) { //(function,det,spin) (orb #)
23 int ndet=occupation.GetDim(1);
24 int ns=occupation.GetDim(2);
25 Array1 <int> nspin(ns);
26 for(int s=0; s< ns; s++) nspin(s)=occupation(0,0,s).GetDim(0);
27 excitations.Resize(ndet);
28 //Assume that the base determinant is the first one!
29 int base=0;
30 for(int d=0; d< ndet; d++) {
31 //Find the differences with the base determinant
32 vector <vector <int> > tot_missing(ns);
33 vector <vector <int> > tot_additional(ns);
34 //This search could likely be improved, but be careful with ordering,
35 for(int s=0; s< ns; s++) {
36 for(int e1=0; e1 < nspin(s); e1++) {
37 int o1=occupation(0,base,s)(e1);
38 bool found=false;
39 for(int e2=0; e2 < nspin(s); e2++) {
40 if(occupation(0,d,s)(e2)==o1) {
41 found=true;
42 break;
43 }
44 }
45 if(!found)tot_missing[s].push_back(o1);
46 int o2=occupation(0,d,s)(e1);
47 found=false;
48 for(int e2=0; e2 < nspin(s); e2++) {
49 if(occupation(0,base,s)(e2)==o2) {
50 found=true;
51 break;
52 }
53 }
54 if(!found) {
55 tot_additional[s].push_back(o2);
56 }
57
58 }
59 }
60 //The update formulation works in terms of replacing an occupied orbital with a
61 //virtual one. By default, QWalk converters order determinants by orbital number.
62 //These representations are equivalent up to a sign, which we determine here in
63 //a fairly general way
64
65 excitations[d].sign.Resize(ns);
66 excitations[d].sign=1;
67 for(int s=0; s< ns; s++) {
68 Array1<int> newocc=occupation(0,base,s);
69 int nex=tot_additional[s].size();
70 for(int e=0; e< nspin(s); e++) {
71 for(int ex=0; ex< nex; ex++) {
72 if(newocc[e]==tot_missing[s][ex])
73 newocc[e]=tot_additional[s][ex];
74 }
75 }
76 int count=0;
77 for(int e1=0; e1 < nspin(s); e1++) {
78 for(int e2=e1+1; e2 < nspin(s); e2++) {
79 if(newocc[e2]<newocc(e1)) count++;
80 if(occupation(0,d,s)(e2) < occupation(0,d,s)(e1)) count--;
81 }
82 }
83 // cout << "d " << d << " count " << count << endl;
84 if(count%2==1) excitations[d].sign(s)*=-1;
85 }
86
87
88 //tot_missing and tot_additional should be filled now.
89 excitations[d].g.Resize(ns);
90 excitations[d].e.Resize(ns);
91
92 for(int s=0;s < ns; s++) {
93 int nex_s=tot_missing[s].size();
94 if(nex_s != tot_additional[s].size()) {
95 error("In build_excitationscitation_list: different numbers of holes and particles");
96 }
97 excitations[d].g[s].Resize(nex_s);
98 excitations[d].e[s].Resize(nex_s);
99 for(int i=0; i< nex_s; i++) {
100 excitations[d].g[s][i]=tot_missing[s][i];
101 excitations[d].e[s][i]=tot_additional[s][i];
102 cout << "d " << d << " s " << s << " i " << i << " : " << excitations[d].g[s][i] << " -> " << excitations[d].e[s][i] << " sign " << excitations[d].sign(s) << endl;
103 }
104 }
105
106 }
107
108
109 //Now we build some of the lookup variables
110 int nex=excitations.GetDim(0);
111 alle.resize(ns);
112 allg.resize(ns);
113 for(int e=0; e< nex; e++) {
114 for(int s=0; s< ns; s++) {
115 for(int i=0; i< excitations(e).g(s).GetDim(0); i++) {
116 int g=excitations(e).g(s)(i);
117 bool found=false;
118 for(vector<int>::iterator gi=allg[s].begin(); gi!=allg[s].end(); gi++) {
119 if(g==*gi) {
120 found=true;
121 break;
122 }
123 }
124 if(!found) allg[s].push_back(g);
125 int ex=excitations(e).e(s)(i);
126 found=false;
127 for(vector<int>::iterator ei=alle[s].begin(); ei!=alle[s].end(); ei++) {
128 if(ex==*ei) {
129 found=true;
130 break;
131 }
132 }
133 if(!found) alle[s].push_back(ex);
134 }
135 }
136 }
137
138
139
140 //Array1 <Excitation> remap(nex);
141 remap.Resize(nex);
142 for(int e=0; e< nex; e++) {
143 remap(e).g.Resize(ns);
144 remap(e).e.Resize(ns);
145 for(int s=0; s< ns; s++) {
146 int n=excitations(e).g(s).GetDim(0);
147
148 remap(e).g(s).Resize(n);
149 remap(e).e(s).Resize(n);
150
151 int ng=allg[s].size();
152 int ne=alle[s].size();
153 for(int i=0; i< n; i++) {
154 for(int j=0; j< ng; j++) {
155 if(allg[s][j]==excitations(e).g(s)(i)) {
156 remap(e).g(s)(i)=j;
157 break;
158 }
159 }
160 for(int j=0; j< ne; j++) {
161 if(alle[s][j]==excitations(e).e(s)(i)) {
162 remap(e).e(s)(i)=j;
163 break;
164 }
165 }
166 }
167 }
168 }
169
170
171 }
172
173
174