1 /*
2
3 Copyright (C) 2007 Zachary Helms
4 with further modifications by Lucas K. Wagner
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License along
17 with this program; if not, write to the Free Software Foundation, Inc.,
18 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
19
20 */
21
22 #include "Program_options.h"
23 #include "Determinant_select_method.h"
24 #include "qmc_io.h"
25 #include "System.h"
26 #include "MatrixAlgebra.h"
27 #include "ulec.h"
28 /*!
29 Read the "words" from the method section in the input file
30 via doinput() parsing, and store section information in private
31 variables orbs, resolution, and minmax.
32 Set up MO_matrix and Sample_point objects for wavefunction from input
33 */
read(vector<string> words,unsigned int & pos,Program_options & options)34 void Determinant_select_method::read(vector <string> words,
35 unsigned int & pos,
36 Program_options & options)
37 {
38 cout << "rrreeead " << endl;
39
40 sys=NULL;
41 allocate(options.systemtext[0], sys);
42
43 vector <vector < string> > orbgroup_txt;
44 pos=0;
45 vector < string> dummy;
46 while( readsection(words,pos,dummy,"ORB_GROUP")) {
47 orbgroup_txt.push_back(dummy);
48 }
49 cout << "aaaahree " << endl;
50 orbital_groups.Resize(orbgroup_txt.size());
51 for(int i=0; i< orbital_groups.GetDim(0); i++) {
52 int norb_this=orbgroup_txt[i].size();
53 orbital_groups[i].Resize(norb_this);
54 for(int j=0; j< norb_this; j++) {
55 orbital_groups[i][j]=atoi(orbgroup_txt[i][j].c_str())-1;
56 }
57 }
58
59 cout << "here " << endl;
60
61
62 vector <string> orbtext;
63 if(!readsection(words, pos=0, orbtext, "ORBITALS")){
64 error("Need ORBITALS");
65 }
66 allocate(orbtext, sys, mymomat);
67
68 vector <string> occupied_up, virtu_up, occupied_down, virtu_down;
69 if(!readsection(words, pos=0, occupied_up, "OCCUPIED_UP"))
70 error("Need OCCUPIED section");
71 if(!readsection(words, pos=0, virtu_up, "VIRTUAL_UP"))
72 error("Need VIRTUAL section");
73 if(!readsection(words, pos=0, occupied_down, "OCCUPIED_DOWN"))
74 error("Need OCCUPIED section");
75 if(!readsection(words, pos=0, virtu_down, "VIRTUAL_DOWN"))
76 error("Need VIRTUAL section");
77 nocc.Resize(2); nvirt.Resize(2);
78 nocc(0)=occupied_up.size();
79 nocc(1)=occupied_down.size();
80 nvirt(0)=virtu_up.size();
81 nvirt(1)=virtu_down.size();
82 occ.Resize(2,max(nocc(0),nocc(1)));
83 virt.Resize(2,max(nvirt(0),nvirt(1)));
84 for(int i=0; i< nocc(0); i++) occ(0,i)=atoi(occupied_up[i].c_str())-1;
85 for(int i=0; i< nocc(1); i++) occ(1,i)=atoi(occupied_down[i].c_str())-1;
86
87 for(int i=0; i< nvirt(0); i++) virt(0,i)=atoi(virtu_up[i].c_str())-1;
88 for(int i=0; i< nvirt(1); i++) virt(1,i)=atoi(virtu_down[i].c_str())-1;
89 cout << "done read " << endl;
90 mywalker=NULL;
91 sys->generateSample(mywalker);
92
93 mymomat->buildLists(orbital_groups);
94 cout << " done build " << endl;
95 }
96
97 //----------------------------------------------------------------------
showinfo(ostream & os)98 int Determinant_select_method::showinfo(ostream & os) {
99 os << "Determinant_select method " << endl;
100 return 1;
101 }
102
103 #include "Generate_sample.h"
104
105 //----------------------------------------------------------------------
run(Program_options & options,ostream & output)106 void Determinant_select_method::run(Program_options & options, ostream & output) {
107 int nsamples=1000;
108 for(int g=0; g< orbital_groups.GetDim(0); g++) {
109 int norb=orbital_groups[g].GetDim(0);
110 Array1 <Array1 <doublevar> > r_sample(nsamples);
111 generate_mo_sample(mywalker,sys,mymomat,g,nsamples,r_sample);
112 Array2 <doublevar> vals(nsamples,norb);
113 Array2 <doublevar> tmp_val(norb,2);
114 Array1 <doublevar> orb_norm(norb,0.0);
115 for(int s=0;s < nsamples;s++) {
116 mywalker->setElectronPos(0,r_sample(s));
117 mymomat->updateVal(mywalker,0,g,tmp_val);
118 for(int i=0; i< norb; i++) {
119 vals(s,i)=tmp_val(i,0);
120 cout << s << i << " vals " << vals(s,i) << endl;
121 orb_norm(i)+=vals(s,i)*vals(s,i);
122 }
123 }
124 for(int i=0; i< norb; i++)
125 orb_norm(i)/=nsamples;
126
127 Array4 <doublevar> voverlap(norb,norb,norb,norb);
128 voverlap=0.0;
129 for(int e1=0; e1 < nsamples; e1++) {
130 cout << "sample " << e1 << endl;
131 for(int e2=e1+1; e2< nsamples; e2++) {
132 mywalker->setElectronPos(0,r_sample(e1));
133 mywalker->setElectronPos(1,r_sample(e2));
134 mywalker->updateEEDist();
135 Array1 <doublevar> ree(5);
136 mywalker->getEEDist(0,1,ree);
137 doublevar rij=1.0/ree(0);
138 //cout << "rij " << rij << endl;
139
140 for(int i=0; i< norb; i++) {
141 for(int j=0; j< norb; j++) {
142 for(int k=0; k< norb; k++) {
143 for(int l=0; l< norb; l++) {
144 voverlap(i,j,k,l)+=vals(e1,i)*vals(e2,j)*rij*vals(e1,k)*vals(e2,l);
145 }
146 }
147 }
148 }
149 }
150 }
151 Array2 <int> lookup_occ=occ, lookup_virt=virt;
152 for(int s=0; s< 2; s++) {
153 for(int i=0; i< occ.GetDim(1); i++) {
154 for(int j=0; j< norb; j++) {
155 if(occ(s,i)==orbital_groups(g)(j))
156 lookup_occ(s,i)=j;
157 }
158 }
159 for(int i=0; i< virt.GetDim(1); i++) {
160 for(int j=0; j< norb; j++) {
161 if(virt(s,i)==orbital_groups(g)(j))
162 lookup_virt(s,i)=j;
163 }
164 }
165 }
166 doublevar renorm=1.0/(nsamples*(nsamples-1)/2.0);
167 for(int s1=0; s1 < 2; s1++) {
168 for(int s2=s1; s2 < 2; s2++) {
169 output << "----------Between channels " << s1 << " and " << s2 << endl;
170 for(int i=0; i< nocc(s1); i++) {
171 for(int j=0; j< nocc(s2); j++) {
172 for(int k=0; k< nvirt(s1); k++) {
173 for(int l=0; l< nvirt(s2); l++) {
174 int oi=lookup_occ(s1,i);
175 int oj=lookup_occ(s2,j);
176 int ok=lookup_virt(s1,k);
177 int ol=lookup_virt(s2,l);
178 //cout << "i " << i << " j " << j << " k " << k << " l " << l <<endl;
179 //cout << "oi " << oi << " oj " << oj << " ok " << ok << " ol " << ol << endl;
180 output << occ(s1,i)+1 << " " << occ(s2,j)+1 << " " << virt(s1,k)+1 << " "
181 << virt(s2,l)+1 << " "
182 << voverlap(oi,oj,ok,ol)*renorm/sqrt(orb_norm(oi)*orb_norm(oj)*orb_norm(ok)*orb_norm(ol)) << endl;
183 }
184 }
185 }
186 }
187 }
188 }
189
190
191 }
192 }
193