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