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