1 //------------------------------------------------------------------------
2 //src/Backflow_wf_data.cpp
3
4 #include "Qmc_std.h"
5 #include "qmc_io.h"
6 #include "Backflow_wf_data.h"
7 #include "Backflow_pf_wf_data.h"
8 #include "Wavefunction_data.h"
9 #include "Backflow_pf_wf.h"
10 #include <algorithm>
11
12
read(vector<string> & words,System * sys)13 void Pfaffian_keeper::read(vector <string> & words, System * sys) {
14 coef_eps=1e-10;
15 unsigned int pos=0;
16 vector <string> strpfwt;
17 vector <vector <string> > statevec;
18 unsigned int startpos=pos;
19 //unsigned int tmp;
20
21 vector <string> npairsstr;
22 if( mpi_info.node==0 )
23 cout<<"Pfaffian"<<endl;
24
25 readsection(words, pos, npairsstr, "NPAIRS");
26 if(npairsstr.size() != 3)
27 error("NPAIRS must have 3 elements");
28 npairs.Resize(3);
29 npairs(0)=atoi(npairsstr[0].c_str());
30 npairs(1)=atoi(npairsstr[1].c_str());
31 npairs(2)=atoi(npairsstr[2].c_str());
32 assert(npairs(0)>=npairs(1));
33 if(npairs(0)!=sys->nelectrons(0))
34 error("Number of upup pairs must be equal to number of spin up electrons in pfaffian matrix");
35 if(npairs(1)!=sys->nelectrons(1))
36 error("Number of down-up pairs must be equal to number of spin down electrons in pfaffian matrix");
37 if(npairs(2)>sys->nelectrons(0)+sys->nelectrons(1))
38 error("Too many upaired orbitals in pfaffian matrix");
39 if( (npairs(0)+ npairs(1)+ npairs(2))%2!=0)
40 error("need even number of elements in pffafian matrix");
41
42 pos=startpos;
43 if (readsection(words, pos, strpfwt, "PFWT")){
44 npf=strpfwt.size();
45 pfwt.Resize(npf);
46 for(int pf=0; pf < npf; pf++)
47 pfwt(pf)=atof(strpfwt[pf].c_str());
48 }
49 else {
50 npf=1;
51 pfwt.Resize(1);
52 pfwt(0)=1.0;
53 }
54
55 if(haskeyword(words, pos=startpos, "OPTIMIZE_PFWT")) {
56 optimize_pfwt=1;
57 }
58 else optimize_pfwt=0;
59
60 if(haskeyword(words, pos=startpos, "CHECK_PFWT_SIGN")) {
61 check_pfwt_sign=1;
62 }
63 else check_pfwt_sign=0;
64
65 vector <string> str_order;
66 pos=startpos;
67 int max=0;
68 int npairstotal=npairs(0)+ npairs(1)+ npairs(2);
69 if(!readsection(words, pos, str_order,"ORBITAL_ORDER")){
70 //error("need SINGLET_ORDER in pfaffian section");
71 if( mpi_info.node==0 )
72 cout << "Assuming the defaul configuration of i-th orbital per i-th pfaffian "<<endl;
73 order_in_pfaffian.Resize(npf);
74 max=npf-1;
75 for(int pf=0;pf<npf;pf++){
76 order_in_pfaffian(pf).Resize(npairstotal,npairstotal);
77 order_in_pfaffian(pf)=0;
78 for(int i=0;i<npairstotal;i++){
79 for(int j=i+1;j<npairstotal;j++){
80 order_in_pfaffian(pf)(i,j)=order_in_pfaffian(pf)(j,i)=pf;
81 }
82 }
83 }
84 }
85 else {
86 if (int(str_order.size())!=npairstotal*(npairstotal-1)*npf/2)
87 error("supply the right amount of values in the ORBITAL_ORDER section");
88
89 order_in_pfaffian.Resize(npf);
90 int counter=0;
91 for(int pf=0;pf<npf;pf++){
92 order_in_pfaffian(pf).Resize(npairstotal,npairstotal);
93 order_in_pfaffian(pf)=0;
94 for(int i=0;i<npairstotal;i++){
95 for(int j=i+1;j<npairstotal;j++){
96 order_in_pfaffian(pf)(i,j)=atoi(str_order[counter].c_str())-1;
97 order_in_pfaffian(pf)(j,i)=order_in_pfaffian(pf)(i,j);
98 if( order_in_pfaffian(pf)(i,j)>max )
99 max=order_in_pfaffian(pf)(i,j);
100 counter++;
101 }
102 }
103 }
104 }
105 if( mpi_info.node==0 ){
106 cout << "ORBITAL_ORDER"<<endl;
107 for(int pf=0;pf<npf;pf++){
108 for(int i=0;i<npairstotal;i++){
109 for(int j=0;j<npairstotal;j++){
110 if((i<npairs(0)+npairs(1) || j<npairs(0)+npairs(1)) && (i!=j))
111 cout << order_in_pfaffian(pf)(i,j)+1<<" ";
112 else
113 cout << order_in_pfaffian(pf)(i,j)<<" ";
114 }
115 cout<<endl;
116 }
117 cout <<endl;
118 }
119
120 }
121
122
123 nsfunc=max+1;
124 if( mpi_info.node==0 )
125 cout <<"Number of functions needed: "<<nsfunc<<endl;
126
127 occupation.Resize(nsfunc);
128 occupation_pos.Resize(nsfunc);
129 totoccupation.Resize(1);
130 tripletorbuu.Resize(nsfunc);
131 tripletorbdd.Resize(nsfunc);
132 singletorb.Resize(nsfunc);
133 unpairedorb.Resize(nsfunc);
134 normalization.Resize(nsfunc);
135 optimize_string.Resize(nsfunc);
136 optimize_total.Resize(nsfunc);
137 noccupied.Resize(nsfunc);
138 optimize_pf=0;
139 ntote_pairs.Resize(nsfunc);
140
141 //orbitals are read in bfwrapper section
142
143 vector <vector <string> > strpairingorbs;
144 vector <string> strpairingorbs_tmp;
145 pos=startpos;
146 while(readsection(words, pos, strpairingorbs_tmp, "PAIRING_ORBITAL")){
147 strpairingorbs.push_back(strpairingorbs_tmp);
148 }
149 if( int(strpairingorbs.size()) < nsfunc)
150 error("Could not find enough PAIRING_ORBITALs");
151 for (int pf=0;pf<nsfunc;pf++){
152 allocate_pairing_orbital(strpairingorbs[pf],pf);
153 }
154
155 // here will start to read the particular pairing functions for each weight
156 if( mpi_info.node==0 )
157 cout<<"Size of pfaffian matrix is "<<npairs(0)+ npairs(1)+ npairs(2)
158 <<"x"<<npairs(0)+ npairs(1)+ npairs(2)<<endl;
159
160 nelectrons.Resize(2);
161
162 nelectrons(0)=sys->nelectrons(0);
163 nelectrons(1)=sys->nelectrons(1);
164
165 tote=nelectrons(0)+nelectrons(1);
166 ndim=3;
167
168 //cout <<"end of pfaffian read"<<endl;
169
170 //molecorb->buildLists(totoccupation);
171
172 }
173
174
175 //----------------------------------------------------------------------
176
177
178 //----------------------------------------------------------------------
179
getOccupation(Array1<Array1<int>> & totaloccupation,int nmoread)180 void Pfaffian_keeper::getOccupation(Array1 <Array1 <int> > & totaloccupation,
181 int nmoread) {
182
183
184 //cout <<"getOccupation "<<endl;
185 nmo=nmoread;
186 Array1 <int> totocc_temp(nmo);
187 totocc_temp=0;
188
189
190 int counter=0;
191 // fill the array for totoccupation
192 for (int mo=0;mo<nmo;mo++){
193 for (int pf=0;pf<nsfunc;pf++){
194 occupation_pos(pf).Resize(occupation(pf).GetSize());
195 for(int orb=0;orb<occupation(pf).GetSize();orb++){
196 if(occupation(pf)(orb)==mo){
197 totocc_temp(mo)=1;
198 }
199 }
200 }
201 if(totocc_temp(mo))
202 counter++;
203 }
204
205
206 totoccupation(0).Resize(counter);
207 totaloccupation.Resize(1);
208 totaloccupation(0).Resize(counter);
209
210 counter=0;
211 for (int mo=0;mo<nmo;mo++){
212 if(totocc_temp(mo)){
213 totoccupation(0)(counter)=mo;
214 totaloccupation(0)(counter)=totoccupation(0)(counter);
215 counter++;
216 }
217 }
218
219 if( mpi_info.node==0 ){
220 cout << "These are all orbitals involved in the Wavefunction: "<<endl;
221 for(int mo=0;mo<totoccupation(0).GetSize();mo++){
222 cout << totoccupation(0)(mo)+1 <<" ";
223 }
224 cout << endl;
225 }
226
227 for (int pf=0;pf<nsfunc;pf++){
228 if( mpi_info.node==0 )
229 cout <<"Orbitals positions in totoccupation for Pfaffian "<<pf<<" : "<<endl;
230 for(int orb=0;orb<occupation(pf).GetSize();orb++){
231 for(int mo=0;mo<totoccupation(0).GetSize();mo++)
232 if(totoccupation(0)(mo)==occupation(pf)(orb))
233 occupation_pos(pf)(orb)=mo;
234 if( mpi_info.node==0 )
235 cout << occupation_pos(pf)(orb)+1<<" ";
236 }
237 if( mpi_info.node==0 )
238 cout <<endl;
239 }
240 }
241
242
243 //----------------------------------------------------------------------
allocate_pairing_orbital(vector<string> & pf_place,int orbint)244 void Pfaffian_keeper::allocate_pairing_orbital(vector <string> & pf_place,
245 int orbint){
246
247 unsigned int startpos=0;
248
249
250 //ORBITALS_IN_PAIRING
251 vector <string> strobs_pairs;
252 unsigned int pos=startpos;
253 if(!readsection(pf_place, pos, strobs_pairs, "ORBITALS_IN_PAIRING")){
254 if( mpi_info.node==0 )
255 cout << "Assuming you want to use all orbitals up to NMO in PFAFFIAN section \n";
256 occupation(orbint).Resize(nmo);
257 for (int i=0;i<nmo;i++){
258 occupation(orbint)(i)=i;
259 }
260 }
261 else {
262 occupation(orbint).Resize(strobs_pairs.size());
263 for (unsigned int i=0;i<strobs_pairs.size();i++){
264 occupation(orbint)(i)=atoi(strobs_pairs[i].c_str())-1;
265 }
266 }
267
268 if( mpi_info.node==0 )
269 cout <<occupation(orbint).GetSize()<<" orbitals involved in pairing for "<<
270 orbint+1<<" pairing orb"<<endl;
271
272 ntote_pairs(orbint)=occupation(orbint).GetSize();
273 if( mpi_info.node==0 )
274 cout << "ntote_pairs("<<orbint<<") : "<<ntote_pairs(orbint)<<endl;
275
276
277 //OPTIMIZE PART
278 vector <string> pf_place2;
279
280 pos=startpos;
281 if(readsection(pf_place, pos, pf_place2, "OPTIMIZE_PF")){
282 optimize_pf=1;
283 optimize_string(orbint).Resize(4);
284 optimize_total(orbint).Resize(4);
285 Pfaffian_optimize_read(pf_place2, orbint, pos);
286 }
287 else {
288 optimize_string(orbint).Resize(0);
289 optimize_total(orbint).Resize(0);
290 }
291
292 vector <string> strtripletorb1,strtripletorb2;
293 vector <string> strsingletorb;
294 vector <string> strunpairedorb;
295 vector <string> strnormalization;
296
297 Array1 <Array1 <doublevar> > unpairedorb_temp;
298
299 normalization(orbint).Resize(3+npairs(2));
300 normalization(orbint)=1.0;
301
302
303 pos=startpos;
304 normalization(orbint).Resize(3+npairs(2));
305
306 if( readsection(pf_place,pos,strnormalization,"ALPHA")){
307 if( mpi_info.node==0 )
308 cout<<"ALPHA section is no longer used"<<endl;
309 }
310
311
312 pos=startpos;
313 tripletorbuu(orbint).Resize(ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2);
314 tripletorbuu(orbint)=0;
315 if(readsection(pf_place, pos, strtripletorb1,"TRIPLET_UU_COEF")){
316 if(int(strtripletorb1.size()) < ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2){
317 if( mpi_info.node==0 )
318 cout <<"Will fill up upto "<< ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2 << " values in TRIPLET_UU_COEF \n";
319 }
320 else if (int(strtripletorb1.size()) > ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2){
321 error("Too many TRIPLET_UU_COEF in PFAFFIAN section");
322 }
323 for(unsigned int i=0; i<strtripletorb1.size(); i++)
324 {
325 tripletorbuu(orbint)(i)=atof(strtripletorb1[i].c_str());
326 }
327 }
328
329 Renormalization(tripletorbuu(orbint), normalization(orbint)(0),coef_eps);
330
331 if (npairs(1)) {
332 pos=startpos;
333 tripletorbdd(orbint).Resize(ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2);
334 tripletorbdd(orbint)=0;
335 if(readsection(pf_place,pos,strtripletorb2,"TRIPLET_DD_COEF")){
336 if(int(strtripletorb2.size()) < ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2){
337 if( mpi_info.node==0 )
338 cout <<"Will fill up upto "<< ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2 << " values in TRIPLET_DD_COEF \n";
339 }
340 else if (int(strtripletorb2.size()) > ntote_pairs(orbint)*(ntote_pairs(orbint)-1)/2){
341 error("Too many TRIPLET_DD_COEF in PFAFFIAN section");
342 }
343 for(unsigned int i=0; i<strtripletorb2.size(); i++)
344 {
345 tripletorbdd(orbint)(i)=atof(strtripletorb2[i].c_str());
346 }
347 }
348 Renormalization(tripletorbdd(orbint), normalization(orbint)(1),coef_eps);
349
350 pos=startpos;
351 singletorb(orbint).Resize(ntote_pairs(orbint)*(ntote_pairs(orbint)+1)/2);
352 singletorb(orbint)=0;
353 int kounter=0;
354 int k=0;
355 for(int i=0;i<ntote_pairs(orbint);i++)
356 for(int j=i;j<ntote_pairs(orbint);j++){
357 if (i==j && k< npairs(1)){
358 singletorb(orbint)(kounter)=1.0;
359 k++;
360 }
361 kounter++;
362 }
363 if(readsection(pf_place,pos,strsingletorb,"SINGLET_COEF")){
364 if(int(strsingletorb.size()) < ntote_pairs(orbint)*(ntote_pairs(orbint)+1)/2)
365 if( mpi_info.node==0 )
366 cout <<"Will fill up upto "<< ntote_pairs(orbint)*(ntote_pairs(orbint)+1)/2 << " values in SINGLET_COEF \n";
367 else if (int(strsingletorb.size()) > ntote_pairs(orbint)*(ntote_pairs(orbint)+1)/2){
368 error("Too many SINGLET_COEF in PFAFFIAN section");
369 }
370 for(unsigned int i=0; i<strsingletorb.size(); i++)
371 singletorb(orbint)(i)=atof(strsingletorb[i].c_str());
372 }
373 Renormalization(singletorb(orbint), normalization(orbint)(2),coef_eps);
374 }
375 else {
376 tripletorbdd(orbint).Resize(0);
377 singletorb(orbint).Resize(0);
378 }
379 if (npairs(2)) {
380 pos=startpos;
381 unpairedorb(orbint).Resize(ntote_pairs(orbint)*npairs(2));
382 unpairedorb(orbint)=0;
383
384 if(!readsection(pf_place,pos,strunpairedorb,"UNPAIRED_COEF")){
385 error("Need UNPAIRED_COEF in PFAFFIAN section");}
386 if(int(strunpairedorb.size()) < ntote_pairs(orbint)*npairs(2))
387 if( mpi_info.node==0 )
388 cout <<"Will fill up upto "<< ntote_pairs(orbint)*npairs(2) << " values in UNPAIRED_COEF \n";
389 else if (int(strunpairedorb.size()) > ntote_pairs(orbint)*npairs(2)){
390 error("Too many UNPAIRED_COEF's in PFAFFIAN section");
391 }
392 unpairedorb_temp.Resize(npairs(2));
393 for(int i=0; i< npairs(2); i++)
394 {
395 unpairedorb_temp(i).Resize(ntote_pairs(orbint));
396 for(int j=0; j<ntote_pairs(orbint); j++)
397 unpairedorb_temp(i)(j)=atof(strunpairedorb[i*ntote_pairs(orbint)+j].c_str());
398 Renormalization(unpairedorb_temp(i), normalization(orbint)(3+i),coef_eps);//we have to normalize for every row
399 for(int j=0; j<ntote_pairs(orbint); j++){
400 unpairedorb(orbint)(i*ntote_pairs(orbint)+j)=unpairedorb_temp(i)(j);
401 }
402 }
403 }
404 else {
405 unpairedorb(orbint).Resize(0);
406 }
407
408
409 }
410
411 //----------------------------------------------------------------------
Pfaffian_optimize_read(vector<string> & pf_place,int pf,unsigned int startpos)412 void Pfaffian_keeper::Pfaffian_optimize_read(vector <string> & pf_place,
413 int pf, unsigned int startpos){
414
415
416 unsigned int pos=startpos;
417
418
419 if(!readvalue(pf_place, pos=0, noccupied(pf), "NOCCUPIED")){
420 noccupied(pf)=npairs(0);
421 }
422 if( mpi_info.node==0 )
423 cout << "NOCCUPIED : " <<noccupied(pf)<<endl;
424
425 optimize_string(pf)(0)="";
426 optimize_string(pf)(1)="";
427 optimize_string(pf)(2)="";
428 optimize_string(pf)(3)="";
429
430
431 optimize_total(pf)(0).Resize(((ntote_pairs(pf)-1)*ntote_pairs(pf))/2);
432 optimize_total(pf)(1).Resize(((ntote_pairs(pf)-1)*ntote_pairs(pf))/2);
433 optimize_total(pf)(2).Resize(((ntote_pairs(pf)+1)*ntote_pairs(pf))/2);
434 optimize_total(pf)(3).Resize(ntote_pairs(pf)*npairs(2));
435
436
437 for (int i=0;i<optimize_total(pf).GetSize();i++)
438 for (int j=0;j<optimize_total(pf)(i).GetSize();j++)
439 optimize_total(pf)(i)(j)=0;
440
441 if( mpi_info.node==0 )
442 cout <<"Pfaffian parameters to optimize: "<<endl;
443 int counter=0;
444 int fullcounter=0;
445
446 //------------Triplets up up------------------------
447 if(haskeyword(pf_place, pos=0, "TRIPLET_UU_DIAG")){
448 counter=fullcounter=0;
449 for (int i=0;i<ntote_pairs(pf);i++){
450 for (int j=i+1;j<ntote_pairs(pf);j++){
451 if((j==i+1)){//&&(i%2==0)){
452 optimize_total(pf)(0)(fullcounter)=1;
453 counter++;
454 }
455 fullcounter++;
456 }
457 }
458 optimize_string(pf)(0)="TRIPLET_UU_DIAG";
459 if( mpi_info.node==0 )
460 cout <<" "<< "TRIPLET_UU_DIAG: "<<counter<<endl;
461 }
462 else if(haskeyword(pf_place, pos=0,"TRIPLET_UU_ALL")){
463 counter=fullcounter=0;
464 for (int i=0;i<ntote_pairs(pf);i++){
465 for (int j=i+1;j<ntote_pairs(pf);j++){
466 //if(i%2==0){
467 optimize_total(pf)(0)(fullcounter)=1;
468 counter++;
469 //}
470 fullcounter++;
471 }
472 }
473 optimize_string(pf)(0)="TRIPLET_UU_ALL";
474 if( mpi_info.node==0 )
475 cout <<" "<< "TRIPLET_UU_ALL: "<<counter<<endl;
476 }
477 else if(haskeyword(pf_place, pos=0, "TRIPLET_UU_VIRTUAL")){
478 counter=fullcounter=0;
479 for (int i=0;i<ntote_pairs(pf);i++){
480 for (int j=i+1;j<ntote_pairs(pf);j++){
481 if((i>=noccupied(pf))){//&&((i-noccupied(pf))%2==0)){
482 optimize_total(pf)(0)(fullcounter)=1;
483 counter++;
484 }
485 fullcounter++;
486 }
487 }
488 optimize_string(pf)(0)="TRIPLET_UU_VIRTUAL";
489 if( mpi_info.node==0 )
490 cout <<" "<< "TRIPLET_UU_VIRTUAL: "<<counter<<endl;
491 }
492 else if(haskeyword(pf_place, pos=0, "TRIPLET_UU_VIRTUAL_DIAG")){
493 counter=fullcounter=0;
494 for (int i=0;i<ntote_pairs(pf);i++){
495 for (int j=i+1;j<ntote_pairs(pf);j++){
496 if((i>=noccupied(pf))&&(j==i+1)){//&&((i-noccupied(pf))%2==0)){
497 optimize_total(pf)(0)(fullcounter)=1;
498 counter++;
499 }
500 fullcounter++;
501 }
502 }
503 optimize_string(pf)(0)="TRIPLET_UU_HF2VIRTUALS";
504 if( mpi_info.node==0 )
505 cout <<" "<< "TRIPLET_UU_VIRTUAL_DIAG: "<<counter<<endl;
506 }
507 else if(haskeyword(pf_place, pos=0, "TRIPLET_UU_HF2VIRTUALS")){
508 counter=fullcounter=0;
509 for (int i=0;i<ntote_pairs(pf);i++){
510 for (int j=i+1;j<ntote_pairs(pf);j++){
511 if((i>= npairs(1))&&(i< npairs(0))&&(j>=npairs(0))){
512 optimize_total(pf)(0)(fullcounter)=1;
513 counter++;
514 }
515 fullcounter++;
516 }
517 }
518 optimize_string(pf)(0)="TRIPLET_UU_HF2VIRTUALS";
519 if( mpi_info.node==0 )
520 cout <<" "<< "TRIPLET_UU_HF2VIRTUALS: "<<counter<<endl;
521 }
522
523
524 //------------Triplets down down------------------------
525 if(haskeyword(pf_place, pos=0, "TRIPLET_DD_DIAG")){
526 counter=fullcounter=0;
527 for (int i=0;i<ntote_pairs(pf);i++){
528 for (int j=i+1;j<ntote_pairs(pf);j++){
529 if((j==i+1)){//&&(i%2==0)){
530 optimize_total(pf)(1)(fullcounter)=1;
531 counter++;
532 }
533 fullcounter++;
534 }
535 }
536 optimize_string(pf)(1)="TRIPLET_DD_DIAG";
537 if( mpi_info.node==0 )
538 cout <<" "<< "TRIPLET_DD_DIAG: "<<counter<<endl;
539 }
540 else if(haskeyword(pf_place, pos=0,"TRIPLET_DD_ALL")){
541 counter=fullcounter=0;
542 for (int i=0;i<ntote_pairs(pf);i++){
543 for (int j=i+1;j<ntote_pairs(pf);j++){
544 //if(i%2==0){
545 optimize_total(pf)(1)(fullcounter)=1;
546 counter++;
547 //}
548 fullcounter++;
549 }
550 }
551 optimize_string(pf)(1)="TRIPLET_DD_ALL";
552 if( mpi_info.node==0 )
553 cout <<" "<< "TRIPLET_DD_ALL: "<<counter<<endl;
554 }
555 else if(haskeyword(pf_place, pos=0, "TRIPLET_DD_VIRTUAL")){
556 counter=fullcounter=0;
557 for (int i=0;i<ntote_pairs(pf);i++){
558 for (int j=i+1;j<ntote_pairs(pf);j++){
559 if((i>=noccupied(pf))&&((i-noccupied(pf))%2==0)){
560 optimize_total(pf)(1)(fullcounter)=1;
561 counter++;
562 }
563 fullcounter++;
564 }
565 }
566 optimize_string(pf)(1)="TRIPLET_DD_VIRTUAL";
567 if( mpi_info.node==0 )
568 cout <<" "<< "TRIPLET_DD_VIRTUAL: "<<counter<<endl;
569 }
570 else if(haskeyword(pf_place, pos=0, "TRIPLET_DD_VIRTUAL_DIAG")){
571 counter=fullcounter=0;
572 for (int i=0;i<ntote_pairs(pf);i++){
573 for (int j=i+1;j<ntote_pairs(pf);j++){
574 if((i>=noccupied(pf))&&(j==i+1)){//&&((i-noccupied(pf))%2==0)){
575 optimize_total(pf)(1)(fullcounter)=1;
576 counter++;
577 }
578 fullcounter++;
579 }
580 }
581 optimize_string(pf)(1)="TRIPLET_DD_VIRTUAL_DIAG";
582 if( mpi_info.node==0 )
583 cout <<" "<< "TRIPLET_DD_VIRTUAL_DIAG: "<<counter<<endl;
584 }
585 //------------Singlets------------------------
586 if(haskeyword(pf_place, pos=0,"SINGLET_DIAG")){
587 counter=fullcounter=0;
588 for (int i=0;i<ntote_pairs(pf);i++){
589 for (int j=i;j<ntote_pairs(pf);j++){
590 if((j==i)){
591 optimize_total(pf)(2)(fullcounter)=1;
592 counter++;
593 }
594 fullcounter++;
595 }
596 }
597 optimize_string(pf)(2)="SINGLET_DIAG";
598 if( mpi_info.node==0 )
599 cout <<" "<< "SINGLET_DIAG: "<<counter<<endl;
600 }
601 else if(haskeyword(pf_place, pos=0,"SINGLET_ALL")){
602 counter=0;
603 for (int i=0;i<ntote_pairs(pf);i++){
604 for (int j=i;j<ntote_pairs(pf);j++){
605 optimize_total(pf)(2)(counter)=1;
606 counter++;
607 }
608 }
609 optimize_string(pf)(2)="SINGLET_ALL";
610 if( mpi_info.node==0 )
611 cout <<" "<< "SINGLET_ALL: "<<counter<<endl;
612 }
613 /*
614 else if(haskeyword(pf_place, pos=0,"SINGLET_DIAG_SIMPLE")){
615 optimize_string(pf)(2)(2)=npairs(1)+1;
616 cout <<" "<< "SINGLET_DIAG_SIMPLE: "<<optimize_string(pf)(2)(2)<<endl;
617 }
618 else if(haskeyword(pf_place, pos=0,"SINGLET_DIAG_XYZ")){
619 optimize_string(pf)(2)(3)=3;
620 cout <<" "<< "SINGLET_DIAG_XYZ: "<<optimize_string(pf)(2)(3)<<endl;
621 }
622 else if(haskeyword(pf_place, pos=0,"SINGLET_XYZ")){
623 optimize_string(pf)(2)(4)=6;
624 cout <<" "<< "SINGLET_XYZ: "<<optimize_string(pf)(2)(4)<<endl;
625 }
626 */
627 else if(haskeyword(pf_place, pos=0, "SINGLET_VIRTUAL")){
628 counter=fullcounter=0;
629 for (int i=0;i<ntote_pairs(pf);i++){
630 for (int j=i;j<ntote_pairs(pf);j++){
631 if((i>=noccupied(pf))){
632 optimize_total(pf)(2)(fullcounter)=1;
633 counter++;
634 }
635 fullcounter++;
636 }
637 }
638 optimize_string(pf)(2)="SINGLET_VIRTUAL";
639 if( mpi_info.node==0 )
640 cout <<" "<< "SINGLET_VIRTUAL: "<<counter<<endl;
641 }
642 else if(haskeyword(pf_place, pos=0, "SINGLET_VIRTUAL_DIAG")){
643 counter=0;
644 for (int i=0;i<ntote_pairs(pf);i++){
645 for (int j=i;j<ntote_pairs(pf);j++){
646 if((i>=noccupied(pf))&&(j==i)){
647 optimize_total(pf)(2)(fullcounter)=1;
648 counter++;
649 }
650 fullcounter++;
651 }
652 }
653 optimize_string(pf)(2)="SINGLET_VIRTUAL_DIAG";
654 if( mpi_info.node==0 )
655 cout <<" "<< "SINGLET_VIRTUAL_DIAG: "<<counter<<endl;
656 }
657 else if(haskeyword(pf_place, pos=0, "SINGLET_SINGLES_AND_DOUBLES")){
658 counter=0;
659 for (int i=0;i<ntote_pairs(pf);i++){
660 for (int j=i;j<ntote_pairs(pf);j++){
661 if((i<noccupied(pf))||(i==j)){
662 optimize_total(pf)(2)(fullcounter)=1;
663 counter++;
664 }
665 fullcounter++;
666 }
667 }
668 optimize_string(pf)(2)="SINGLET_SINGLES_AND_DOUBLES";
669 if( mpi_info.node==0 )
670 cout <<" "<< "SINGLET_SINGLES_AND_DOUBLES: "<<counter<<endl;
671 }
672 else if(haskeyword(pf_place, pos=0, "SINGLET_HF2VIRTUALS")){
673 counter=fullcounter=0;
674 for (int i=0;i<ntote_pairs(pf);i++){
675 for (int j=i;j<ntote_pairs(pf);j++){
676 if((i< npairs(1))&&(j>= npairs(1))){
677 optimize_total(pf)(2)(fullcounter)=1;
678 counter++;
679 }
680 fullcounter++;
681 }
682 }
683 optimize_string(pf)(2)="SINGLET_HF2VIRTUALS";
684 if( mpi_info.node==0 )
685 cout <<" "<< "SINGLET_HF2VIRTUALS: "<<counter<<endl;
686 }
687
688 //------------Unpaired------------------------
689 if(haskeyword(pf_place, pos=0,"UNPAIRED_ALL")){
690 counter=0;
691 for(int i=0;i<ntote_pairs(pf);i++)
692 for(int j=0;j<npairs(2);j++){
693 optimize_total(pf)(3)(counter)=1;
694 counter++;
695 }
696 optimize_string(pf)(3)="UNPAIRED_ALL";
697 if( mpi_info.node==0 )
698 cout <<" "<< "UNPAIRED_ALL: "<<counter<<endl;
699 }
700 //------------Alpha------------------------
701 if(haskeyword(pf_place, pos=0,"ALPHA")){
702 error("ALPHA section of pairing orbital can no longer be optimized");
703 }
704
705 }
706
707
708 //----------------------------------------------------------------------
709
showinfo(ostream & os)710 void Pfaffian_keeper::showinfo(ostream & os ) {
711 string indent=" ";
712
713 os << "Pfaffian " << endl;
714
715 os <<"Number of pairs: ";
716 for (int i=0;i< npairs.GetSize();i++)
717 os <<npairs(i)<<" ";
718 // os << "." << endl;
719 os <<endl;
720
721 os << "Pffafian weights: ";
722 for (int pf=0;pf< npf;pf++)
723 os <<pfwt(pf)<<" ";
724 os << endl << endl;
725
726 os << "Order of pairing orbital functions per pfaffian"<<endl;
727 for(int pf=0;pf<npf;pf++){
728 os << "Pfaffian: " <<pf+1<<endl;
729 for(int i=0;i<order_in_pfaffian(pf).GetDim(0);i++){
730 for(int j=0;j<order_in_pfaffian(pf).GetDim(1);j++){
731 char tmp[10];
732 char tmp2[10];
733 sprintf(tmp, "(%d,%d) ", i+1,j+1);
734 sprintf(tmp2, "(%d)(%d) ", i+1,j+1);
735 if((i<npairs(0)+npairs(1) || j<npairs(0)+npairs(1)) && (i!=j)){
736 if(i<npairs(0)){
737 if(j<npairs(0))
738 os << " Chi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp;
739 else if(j<npairs(0)+npairs(1))
740 os << " Phi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp;
741 else
742 os << " varphi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp2;
743 }
744 else if(i<npairs(0)+npairs(1)){
745 if(j<npairs(0))
746 os << "-Phi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp;
747 else if(j<npairs(0)+npairs(1))
748 os << " Chi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp;
749 else
750 os << " varphi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp2;
751 }
752 else{
753 if(j<npairs(0))
754 os << "-varphi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp2;
755 else
756 os << "-varphi_"<<order_in_pfaffian(pf)(i,j)+1<<tmp2;
757 }
758
759 }
760 else
761 os <<" "<<order_in_pfaffian(pf)(i,j)<<" ";
762 }
763 os <<endl;
764 }
765 os <<endl;
766 }
767
768 for (int pf=0;pf<nsfunc;pf++){
769
770 os <<"Pairing Orbital "<<pf+1<<endl;
771
772 os <<indent<< "Orbitals involved in pairing: \n";
773 os <<indent;
774 for (int i=0;i< occupation(pf).GetSize();i++)
775 os <<occupation(pf)(i)+1<<" ";
776 //os << "." << endl;
777 os << endl<<endl;
778
779
780 os <<indent<< "Triplet up up coeficients: " << endl;
781 int count=0;
782 // os.setf(ios::showpoint|ios::scientific);
783 for(int i=0;i<ntote_pairs(pf)-1;i++){
784 os <<indent;
785 for(int j=0;j<ntote_pairs(pf);j++){
786 if(j>i){
787 os<<tripletorbuu(pf)(count)<< " ";
788 count++;
789 }
790 //else
791 // os << " ";
792 }
793 os <<endl;
794 }
795
796 //os.unsetf(ios::showpoint|ios::scientific);
797 os <<endl;
798
799 if (npairs(1)){
800 os <<indent<< "Triplet down down coeficients: " << endl;
801 count=0;
802 // os.setf(ios::showpoint|ios::scientific);
803 for(int i=0;i<ntote_pairs(pf)-1;i++){
804 os <<indent;
805 for(int j=0;j<ntote_pairs(pf);j++){
806 if(j>i){
807 os<<tripletorbdd(pf)(count)<< " ";
808 count++;
809 }
810 //else
811 // os << " ";
812 }
813 os <<endl;
814 }
815 os <<endl;
816
817
818 count=0;
819 os <<indent<<"Singlet coeficients: " << endl;
820 //os.setf(ios::showpoint|ios::scientific);
821 for(int i=0;i<ntote_pairs(pf);i++){
822 os <<indent;
823 for(int j=0;j<ntote_pairs(pf);j++){
824 if(j>=i){
825 os<<singletorb(pf)(count)<< " ";
826 count++;
827 }
828 // else
829 // os << " ";
830 }
831 os <<endl;
832 }
833 //os.unsetf(ios::showpoint|ios::scientific);
834 os <<endl;
835 }
836
837 if (npairs(2)){
838 count=0;
839 os <<indent<<"Unpaired coeficients: " << endl;
840 //os.setf(ios::showpoint|ios::scientific);
841 for(int i=0;i< npairs(2);i++){
842 os <<indent;
843 for(int j=0;j<ntote_pairs(pf);j++){
844 os<<unpairedorb(pf)(count)<< " ";
845 count++;
846 }
847 os <<endl;
848 }
849 //os.unsetf(ios::showpoint|ios::scientific);
850 os <<endl;
851 }
852 }
853 }
854
855 //----------------------------------------------------------------------
856
writeinput(string & indent,ostream & os)857 void Pfaffian_keeper::writeinput(string & indent, ostream & os ) {
858 os << indent << "NPAIRS { ";
859 for (int i=0;i< npairs.GetSize();i++)
860 os <<npairs(i)<<" ";
861 os << "}" << endl;
862
863 os << indent << "PFWT { ";
864 for (int pf=0;pf< npf;pf++)
865 os <<pfwt(pf)<<" ";
866 os << "}" << endl;
867 if (optimize_pfwt)
868 os << indent << "OPTIMIZE_PFWT" <<endl;
869
870 os << indent << "ORBITAL_ORDER { "<<endl;
871 for(int pf=0;pf<npf;pf++){
872 for(int i=0;i<order_in_pfaffian(pf).GetDim(0);i++){
873 for(int j=i+1;j<order_in_pfaffian(pf).GetDim(1);j++){
874 os << indent << order_in_pfaffian(pf)(i,j)+1 <<" ";
875 }
876 os << indent << endl;
877 }
878 os << indent << endl;
879 }
880 os << indent << "}" << endl;
881
882
883 for(int pf=0;pf<nsfunc;pf++){
884 os << indent << "PAIRING_ORBITAL { "<<endl;
885
886 os << indent << indent << "ORBITALS_IN_PAIRING { ";
887 for (int i=0;i< occupation(pf).GetSize();i++)
888 os <<occupation(pf)(i)+1<<" ";
889 os << "}" << endl;
890
891 if(optimize_pf) {
892 if(optimize_string(pf).GetSize()){
893 os << indent << indent << "OPTIMIZE_PF {" << endl;
894 os << indent << indent << indent << "NOCCUPIED " <<noccupied(pf)<<endl;
895 for(int i=0;i<optimize_string(pf).GetSize();i++){
896 if (optimize_string(pf)(i)!="")
897 os << indent << indent << indent << optimize_string(pf)(i) << endl;
898 }
899 os << indent << indent << "}" << endl;
900 }
901 }
902
903 os << indent<< indent << "TRIPLET_UU_COEF {" << endl;
904 Renormalization(tripletorbuu(pf), normalization(pf)(0),coef_eps);
905 int count=0;
906 // os.setf(ios::showpoint|ios::scientific);
907 for(int i=0;i<ntote_pairs(pf)-1;i++){
908 os << indent << indent;
909 for(int j=0;j<ntote_pairs(pf);j++){
910 if(j>i){
911 os<<tripletorbuu(pf)(count)<< " ";
912 count++;
913 }
914 // else
915 // os << " ";
916 }
917 os <<endl;
918 }
919
920 //os.unsetf(ios::showpoint|ios::scientific);
921 os << indent <<indent << "}"<<endl;
922
923 if (npairs(1)){
924 os << indent << indent << "TRIPLET_DD_COEF {" << endl;
925 Renormalization(tripletorbdd(pf), normalization(pf)(1),coef_eps);
926 count=0;
927 // os.setf(ios::showpoint|ios::scientific);
928 for(int i=0;i<ntote_pairs(pf)-1;i++){
929 os << indent << indent;
930 for(int j=0;j<ntote_pairs(pf);j++){
931 if(j>i){
932 os<<tripletorbdd(pf)(count)<< " ";
933 count++;
934 }
935 // else
936 // os << " ";
937 }
938 os <<endl;
939 }
940
941 //os.unsetf(ios::showpoint|ios::scientific);
942 os << indent << indent << "}"<<endl;
943
944
945
946 count=0;
947 os << indent << indent<< "SINGLET_COEF {" << endl;
948 Renormalization(singletorb(pf), normalization(pf)(2),coef_eps);
949 //os.setf(ios::showpoint|ios::scientific);
950 for(int i=0;i<ntote_pairs(pf);i++){
951 os<< indent << indent;
952 for(int j=0;j<ntote_pairs(pf);j++){
953 if(j>=i){
954 os<<singletorb(pf)(count)<< " ";
955 count++;
956 }
957 // else
958 // os << " ";
959 }
960 os <<endl;
961 }
962 //os.unsetf(ios::showpoint|ios::scientific);
963 os << indent << indent << "}"<<endl;
964 }
965
966 if (npairs(2)){
967 count=0;
968 os << indent << indent << "UNPAIRED_COEF {" << endl;
969 Array1 <Array1 <doublevar> > unpairedorb_temp(npairs(2));
970 for(int i=0; i<npairs(2); i++){
971 unpairedorb_temp(i).Resize(ntote_pairs(pf));
972 for(int j=0; j<ntote_pairs(pf); j++)
973 unpairedorb_temp(i)(j)=unpairedorb(pf)(i*ntote_pairs(pf)+j);
974 Renormalization(unpairedorb_temp(i), normalization(pf)(3+i),coef_eps);//we have to normalize for every row
975 for(int j=0; j<ntote_pairs(pf); j++)
976 unpairedorb(pf)(i*ntote_pairs(pf)+j)=unpairedorb_temp(i)(j);
977 }
978 //os.setf(ios::showpoint|ios::scientific);
979 for(int i=0;i<npairs(2);i++){
980 os<< indent << indent ;
981 for(int j=0;j<ntote_pairs(pf);j++){
982 os<<unpairedorb(pf)(count)<< " ";
983 count++;
984 }
985 os <<endl;
986 }
987 //os.unsetf(ios::showpoint|ios::scientific);
988 os << indent << indent << "}"<<endl;
989 }
990 os << indent << "}"<<endl;
991 }
992 }
993
994
nparms()995 int Pfaffian_keeper::nparms(){
996 if(optimize_pf || optimize_pfwt) {
997
998 int counter=0;
999 for(int pf=0;pf<nsfunc;pf++)
1000 for(int k=0;k<optimize_total(pf).GetSize();k++)
1001 for(int l=0;l<optimize_total(pf)(k).GetSize();l++){
1002 if(optimize_total(pf)(k)(l)){
1003 counter++;
1004 }
1005 }
1006 if (optimize_pfwt){
1007 for(int pf=1; pf< pfwt.GetSize(); pf++) {
1008 counter++;
1009 }
1010 }
1011 return counter;
1012 }
1013 else
1014 return 0;
1015 }
1016
1017 //----------------------------------------------------------------------
1018
getPFCoeff(Array1<doublevar> & parms)1019 void Pfaffian_keeper::getPFCoeff(Array1 <doublevar> & parms){
1020 //cout <<"Start getPFCoeff"<<endl;
1021 int counter=0;
1022 for(int pf=0;pf<nsfunc;pf++)
1023 for(int k=0;k<optimize_total(pf).GetSize();k++)
1024 for(int l=0;l<optimize_total(pf)(k).GetSize();l++){
1025 if(optimize_total(pf)(k)(l)){
1026 switch(k){
1027 case 0:
1028 parms(counter)=tripletorbuu(pf)(l);
1029 break;
1030 case 1:
1031 parms(counter)=tripletorbdd(pf)(l);
1032 break;
1033 case 2:
1034 parms(counter)=singletorb(pf)(l);
1035 break;
1036 case 3:
1037 parms(counter)=unpairedorb(pf)(l);
1038 break;
1039 }
1040 counter++;
1041 }
1042 }
1043 if (optimize_pfwt){
1044 for(int pf=1; pf< pfwt.GetSize(); pf++) {
1045 parms(counter)=pfwt(pf);
1046 counter++;
1047 }
1048 }
1049 /*
1050 if(mpi_info.node==0) {
1051 cout<<"Get parameters values"<<endl;
1052 for (int i=0;i<parms.GetSize();i++)
1053 cout << i+1<<" "<<parms(i)<<endl;
1054 }
1055 */
1056 }
1057 //----------------------------------------------------------------------
1058
setPFCoeff(Array1<doublevar> & parms)1059 void Pfaffian_keeper::setPFCoeff(Array1 <doublevar> & parms){
1060 // cout <<"Start setPFCoeff"<<endl;
1061 int counter=0;
1062 Array1 <Array1 <doublevar> > unpairedorb_temp(npairs(2));
1063
1064 for(int pf=0;pf<nsfunc;pf++){
1065 for(int k=0;k<optimize_total(pf).GetSize();k++){
1066 for(int l=0;l<optimize_total(pf)(k).GetSize();l++){
1067 if(optimize_total(pf)(k)(l)){
1068 switch(k){
1069 case 0:
1070 tripletorbuu(pf)(l)=parms(counter);
1071 break;
1072 case 1:
1073 tripletorbdd(pf)(l)=parms(counter);
1074 break;
1075 case 2:
1076 singletorb(pf)(l)=parms(counter);
1077 break;
1078 case 3:
1079 unpairedorb(pf)(l)=parms(counter);
1080 break;
1081 }
1082 counter++;
1083 }
1084 }
1085 }
1086
1087 Getnormalization(tripletorbuu(pf), normalization(pf)(0));
1088 Getnormalization(tripletorbdd(pf), normalization(pf)(1));
1089 Getnormalization(singletorb(pf), normalization(pf)(2));
1090
1091 //Unpairedorb renormalization
1092 for(int i=0; i<npairs(2); i++){
1093 unpairedorb_temp(i).Resize(ntote_pairs(pf));
1094 for(int j=0; j<ntote_pairs(pf); j++)
1095 unpairedorb_temp(i)(j)=unpairedorb(pf)(i*ntote_pairs(pf)+j);
1096 Getnormalization(unpairedorb_temp(i), normalization(pf)(3+i));//we have to normalize for every row
1097 for(int j=0; j<ntote_pairs(pf); j++)
1098 unpairedorb(pf)(i*ntote_pairs(pf)+j)=unpairedorb_temp(i)(j);
1099 }
1100 }
1101 if (optimize_pfwt){
1102 for(int pf=1; pf< pfwt.GetSize(); pf++) {
1103 pfwt(pf)=parms(counter);
1104 counter++;
1105 }
1106 }
1107
1108 /*
1109 if(mpi_info.node==0) {
1110 cout<<"Set parameters values"<<endl;
1111 for (int i=0;i<parms.GetSize();i++)
1112 cout << i+1<<" "<<parms(i)<<endl;
1113 }
1114 */
1115 }
1116
1117 //######################################################################
1118 //----------------------------------------------------------------------
1119 //######################################################################
1120
1121 /*!
1122 */
read(vector<string> & words,unsigned int & pos,System * sys)1123 void Backflow_pf_wf_data::read(vector <string> & words, unsigned int & pos,
1124 System * sys)
1125 {
1126 unsigned int startpos=pos;
1127 //optimize_backflow=haskeyword(words, pos=startpos,"OPTIMIZE_BACKFLOW");
1128 vector <string> pfwords;
1129 if(!readsection(words, pos=0, pfwords,"PFAFFIAN"))
1130 error("Couldn't find PFAFFIAN section");
1131 vector <string> bfwords;
1132 if(!readsection(words,pos=0, bfwords,"BFLOW"))
1133 error("Couldn't find BACKFLOW section");
1134
1135 bfwrapper.readOrbitals(sys,bfwords);
1136 pfkeeper.read(pfwords,sys);
1137 Array1 <Array1 <int> > totoccupation;
1138 pfkeeper.getOccupation(totoccupation, bfwrapper.nmo());
1139 bfwrapper.init(sys, totoccupation, bfwords);
1140
1141
1142
1143
1144 }
1145
1146 //----------------------------------------------------------------------
1147
supports(wf_support_type support)1148 int Backflow_pf_wf_data::supports(wf_support_type support) {
1149 switch(support) {
1150 case laplacian_update:
1151 return 1;
1152 case density:
1153 return 0;
1154 case parameter_derivatives:
1155 return 0;
1156 default:
1157 return 0;
1158 }
1159 }
1160
1161 //----------------------------------------------------------------------
1162
1163
1164
generateWavefunction(Wavefunction * & wf)1165 void Backflow_pf_wf_data::generateWavefunction(Wavefunction *& wf)
1166 {
1167 //cout <<"generating wf"<<endl;
1168 assert(wf==NULL);
1169 wf=new Backflow_pf_wf;
1170 Backflow_pf_wf * pfaffwf;
1171 recast(wf, pfaffwf);
1172 pfaffwf->init(this);
1173 attachObserver(pfaffwf);
1174 //cout <<"end of generating wf"<<endl;
1175 }
1176
showinfo(ostream & os)1177 int Backflow_pf_wf_data::showinfo(ostream & os)
1178 {
1179 os << "Backflow-Pfaffian" << endl;
1180 pfkeeper.showinfo(os);
1181 bfwrapper.showinfo(os);
1182 return 1;
1183 }
1184
1185 //----------------------------------------------------------------------
1186
writeinput(string & indent,ostream & os)1187 int Backflow_pf_wf_data::writeinput(string & indent, ostream & os)
1188 {
1189
1190
1191 os << indent << "BACKFLOW-PFAFFIAN" << endl;
1192 os << indent << "BFLOW { " << endl;
1193 string indent2=indent+" ";
1194 bfwrapper.writeinput(indent2,os);
1195 os << indent << "}" << endl;
1196 os << indent << "PFAFFIAN { " << endl;
1197 pfkeeper.writeinput(indent2,os);
1198 os << indent << "}" << endl;
1199
1200 return 1;
1201 }
1202
1203 //------------------------------------------------------------------------
getVarParms(Array1<doublevar> & parms)1204 void Backflow_pf_wf_data::getVarParms(Array1 <doublevar> & parms)
1205 {
1206 //cout <<"start getVarParms"<<endl;
1207 // has to complete
1208 // if(pfkeeper.optimize_pf || pfkeeper.optimize_pfwt) {
1209 // getPFCoeff(parms);
1210 //}
1211
1212 Array1 <doublevar> parms_tmp1;
1213 Array1 <doublevar> parms_tmp2;
1214 bfwrapper.jdata.getVarParms(parms_tmp1);
1215 //start: added for ei back-flow
1216 if(bfwrapper.has_electron_ion_bf){
1217 bfwrapper.jgroupdata.getVarParms(parms_tmp2);
1218 }
1219 //end: added for ei back-flow
1220 int sizetotal=bfwrapper.nparms();
1221 int sizejdata=bfwrapper.jdata.nparms();
1222 parms.Resize(sizetotal);
1223 //start: added for ei back-flow
1224 for(int i=0;i<sizetotal;i++){
1225 if(i<sizejdata)
1226 parms(i)=parms_tmp1(i);
1227 else
1228 parms(i)=parms_tmp2(i-sizejdata);
1229 }
1230 //end: added for ei back-flow
1231 //cout <<"done getVarParms"<<endl;
1232 }
1233
1234 //------------------------------------------------------------------------
setVarParms(Array1<doublevar> & parms)1235 void Backflow_pf_wf_data::setVarParms(Array1 <doublevar> & parms)
1236 {
1237 //cout <<"start setVarParms"<<endl;
1238 assert(parms.GetDim(0)==nparms());
1239 //error("need to do parameter optimization");
1240 Array1 <doublevar> parms_tmp1(bfwrapper.jdata.nparms());
1241 //start: added for ei back-flow
1242 Array1 <doublevar> parms_tmp2(bfwrapper.jgroupdata.nparms());
1243 //end: added for ei back-flow
1244
1245 int sizetotal=bfwrapper.nparms();
1246 int sizejdata=bfwrapper.jdata.nparms();
1247 for(int i=0;i<sizetotal;i++){
1248 if(i<sizejdata)
1249 parms_tmp1(i)=parms(i);
1250 else
1251 parms_tmp2(i-sizejdata)=parms(i);
1252 }
1253 bfwrapper.jdata.setVarParms(parms_tmp1);
1254 //start: added for ei back-flow
1255 if(bfwrapper.has_electron_ion_bf){
1256 bfwrapper.jgroupdata.setVarParms(parms_tmp2);
1257 }
1258 //end: added for ei back-flow
1259 int max=wfObserver.size();
1260 //cout << "slatmax " << max << endl;
1261 for(int i=0; i< max; i++) {
1262 wfObserver[i]->notify(all_wf_parms_change, 0);
1263 }
1264 //cout <<"done setVarParms"<<endl;
1265 }
1266 //------------------------------------------------------------------------
1267