1 /*
2
3 Copyright (C) 2007 Michal Bajdich
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 #include "converter.h"
21 #include <iostream>
22 #include <iomanip>
23 #include <cassert>
24 #include <string>
25 #include <vector>
26 #include <fstream>
27 #include <cmath>
28 #include <cstring>
29 #include <cstdlib>
30
31 using namespace std;
32
reorder_orbs(vector<int> & orbs,double & det_weights)33 void reorder_orbs(vector <int> & orbs, double & det_weights){
34 // sort by ascending value
35 int n=orbs.size();
36 int tmp;
37 for(int i=0; i < n; i++) {
38 if(abs(orbs[i])>99)
39 cout <<"Excitation with orbital > 99. GAMESS formated file is not fully usable. Inspect your result! M.B."<<endl;
40 for(int j=0; j < n; j++) {
41 if (orbs[i]>orbs[j]){
42 tmp=orbs[j];
43 orbs[j]=orbs[i];
44 orbs[i]=tmp;
45 det_weights=-det_weights;
46 }
47 }
48 }
49 /*
50 cout <<"reorder_orbs"<<endl;
51 for(int i=0; i < n; i++) {
52 cout <<orbs[i]<<" ";
53 }
54 cout <<endl;
55 */
56
57 }
58
59
normal_reorder_orbs(int & nup,int & ndown,int & ncore,vector<int> & orbs,double & det_weights)60 void normal_reorder_orbs(int & nup, int & ndown, int & ncore, vector <int> & orbs, double & det_weights){
61 // put to normal order
62 int n=orbs.size();
63 int tmp;
64
65 /*
66 cout <<"before normal_reorder_orbs"<<endl;
67 for(int i=0; i < n; i++) {
68 cout <<orbs[i]<<" ";
69 }
70 cout <<endl;
71 */
72
73 for(int i=0; i < n; i++) {
74 if (orbs[i]>0) {
75 if(orbs[i]!=nup-i){
76 for(int j=ncore; j < nup; j++){
77 if(orbs[i]==j+1 && nup-j-1!=i ){
78 tmp=orbs[nup-j-1];
79 orbs[nup-j-1]=orbs[i];
80 orbs[i]=tmp;
81 det_weights=-det_weights;
82 }
83 }
84 }
85 }
86 else{
87 if(-orbs[i]!=(i+1-(nup))){
88 //cout <<"orbs[i] "<<orbs[i]<<"i+1-(nup) "<<i+1-(nup)<<endl;
89 for(int j=ncore; j < ndown; j++){
90 if(-orbs[i]==j+1 && nup-2*ncore+j!=i ){
91 //cout <<" swap pos "<<nup-2*ncore+j<<" with pos "<<i<<endl;
92 tmp=orbs[nup-2*ncore+j];
93 orbs[nup-2*ncore+j]=orbs[i];
94 orbs[i]=tmp;
95 det_weights=-det_weights;
96 }
97 }
98 }
99 }
100 }
101
102 /*
103 cout <<"after normal_reorder_orbs"<<endl;
104 for(int i=0; i < n; i++) {
105 cout <<orbs[i]<<" ";
106 }
107 cout <<endl;
108 */
109
110 }
111
112
spin_separate(vector<int> & vals,vector<int> & spinup,vector<int> & spindown,int & tmp_up,int & tmp_down)113 void spin_separate(vector <int> & vals, vector <int> & spinup, vector <int> & spindown, int & tmp_up, int & tmp_down){
114 int n=vals.size();
115 tmp_up=tmp_down=0;
116 for (int i=0; i < n; i++){
117 if(vals[i]>tmp_up)
118 tmp_up=vals[i];
119 if(vals[i]<tmp_down)
120 tmp_down=vals[i];
121 if(vals[i]>0)
122 spinup.push_back(vals[i]);
123 else
124 spindown.push_back(-vals[i]);
125 }
126 tmp_down=-tmp_down;
127 //cout <<"The largest spin-up orbital is "<<tmp_up<<" and spin-down orbital is "<<-tmp_down<<endl;
128 }
129
130
sort_abs_largest_first(vector<double> & vals,vector<int> & list)131 void sort_abs_largest_first(vector <double> & vals, vector <int> & list){
132 int n=vals.size();
133 list.resize(n);
134 for (int i=0; i < n; i++)
135 list[i] = i;
136
137 for (int i=1; i < n; i++) {
138 const double temp = vals[i];
139 const double abstemp = fabs(vals[i]);
140 int j;
141 for (j=i-1; j>=0 && fabs(vals[j])<abstemp; j--) {
142 vals[j+1] = vals[j];
143 list[j+1] = list[j];
144 }
145 vals[j+1] = temp;
146 list[j+1] = i;
147 }
148 //for (int i=0; i < n; i++)
149 //cout << list[i]<<endl;
150 }
151
152
usage(const char * name)153 void usage(const char * name) {
154 cout << "usage: " << name << " <options> <output> " << endl;
155 cout << "Where options can be: \n";
156 cout << "-csf write CSFs rather than separate determinants\n";
157 cout << "-o desired output file\n";
158 cout << "-wthresh threshold weight value (Default 0.01)\n";
159 cout << "-state write determinants for selected state \n";
160 cout << "-norder orbitals in determinants follow natural order\n";
161 cout << "-print_level (Default 1)\n";
162 exit(1);
163 }
164
main(int argc,char ** argv)165 int main(int argc, char ** argv) {
166 string infilename;
167 string outputname;
168 int computed_states=1;
169 int iroot=1;
170 int nstate=1;
171 double wtresh=0.01;
172 int symmetry=0;
173 int reorder=0;
174 int printout=1;
175
176 for(int i=1; i< argc-1; i++) {
177 if(!strcmp(argv[i], "-wthresh") && argc > i+1) {
178 wtresh=double(atof(argv[++i]));
179 }
180 else if(!strcmp(argv[i], "-o") && argc > i+1) {
181 outputname=argv[++i];
182 }
183 else if(!strcmp(argv[i], "-csf")) {
184 symmetry=1;
185 }
186 else if (!strcmp(argv[i], "-state")) {
187 nstate=atoi(argv[++i]);
188 }
189 else if (!strcmp(argv[i], "-norder")) {
190 reorder=1;
191 }
192 else if(!strcmp(argv[i], "-print_level") && argc > i+1) {
193 printout=int(atoi(argv[++i]));
194 }
195 else {
196 cout << "Didn't understand option " << argv[i] << endl;
197 usage(argv[0]);
198 }
199 }
200
201 if(argc >= 2) {
202 infilename=argv[argc-1];
203 }
204 else { usage(argv[0]); }
205
206 if(outputname=="") {
207 outputname="gamessci2qmc.out";
208 }
209
210 cout << "Starting converter"<<endl;
211 cout << "Using weight treshhold "<<wtresh<<endl;
212
213 ifstream is(infilename.c_str());
214 if(!is) {
215 cout << "Couldn't open" <<infilename<< endl;
216 exit(1);
217 }
218 vector <string> words;
219 string line;
220 string space=" ";
221 int using_guga=0;
222 int using_determinants=0;
223 int csfmax=0; //how many CSFs was there total
224 int number_of_core_orbitals=0;
225 vector <int> nelectrons(2);
226 while(getline(is, line)){
227 words.clear();
228 split(line, space, words);
229 //find GUGA DISTINCT ROW TABLE
230 if(words[0]=="GUGA" && words[1]=="DISTINCT" && words[2]=="ROW" && words[3]=="TABLE"){
231 //cout <<"GUGA CITYP using CSF's"<<endl;
232 using_guga=1;
233 while(getline(is,line)) {
234 words.clear();
235 split(line, space, words);
236 //first read in CORE ALPHA BETA
237 if(words[0]=="NUMBER" && words[2]=="CORE" && words[3]=="MOLECULAR" && words[4]=="ORBITALS"){
238 number_of_core_orbitals=atoi(words[6].c_str());
239 cout <<number_of_core_orbitals<<" core orbitals"<<endl<<endl;
240 }
241 if(words[0]=="NUMBER" && words[2]=="ALPHA" && words[3]=="ELECTRONS"){
242 nelectrons[0]=atoi(words[5].c_str());
243 cout <<nelectrons[0]<<" alpha electrons"<<endl;
244 }
245 if(words[0]=="NUMBER" && words[2]=="BETA" && words[3]=="ELECTRONS"){
246 nelectrons[1]=atoi(words[5].c_str());
247 cout <<nelectrons[1]<<" beta electrons"<<endl;
248 }
249
250
251 //find COMPUTING THE HAMILTONIAN FOR THE ????? CSF-S
252 if(words[0]=="COMPUTING" && words[2]=="HAMILTONIAN" && words[6]=="CSF-S..."){
253 csfmax=atoi(words[5].c_str());
254 cout <<"GUGA CITYP using "<<csfmax<<" CSF's"<<endl;
255 }
256
257 // find NUMBER OF REQUESTED STATES
258 if(words[0]=="DAVIDSON" && words[1]=="METHOD" && words[3]=="DIAGONALIZATION" ){
259 while(getline(is,line)) {
260 words.clear();
261 split(line, space, words);
262 if(words[0]=="NUMBER" && words[2]=="STATES" && words[3]=="REQUESTED"){
263 computed_states=atoi(words[5].c_str());
264 cout << "Number of states computed in CI calculation : "<< computed_states << endl;
265 break;
266 }
267 }
268 }
269
270 // find IROOT value
271 if(words[0]=="PROPERTIES" && words[3]=="COMPUTED" && words[5]=="STATE" && words[6]=="-IROOT-" ){
272 iroot=atoi(words[7].c_str());
273 cout <<"IROOT value in CI calculation : "<< iroot <<endl;
274 break;
275 }
276 }
277 break;
278 }
279
280 //find AMES LABORATORY DETERMINANTAL FULL CI
281 if(words[0]=="AMES" && words[1]=="LABORATORY" && words[2]=="DETERMINANTAL" && words[3]=="FULL" && words[4]=="CI" ){
282 cout <<"Found ALDET CITYP using determinantal CI"<<endl;
283 using_determinants=1;
284 while(getline(is,line)) {
285 words.clear();
286 split(line, space, words);
287 if(words[0]=="NUMBER" && words[3]=="STATES" && words[4]=="REQUESTED" && words[5]=="="){
288 computed_states=atoi(words[6].c_str());
289 cout <<computed_states<<" states computed in CI calculation"<<endl;
290 }
291 if(words[0]=="CI" && words[1]=="PROPERTIES" && words[6]=="ROOT" && words[7]=="NUMBER"){
292 iroot=atoi(words[8].c_str());
293 cout <<"Properties of state "<<iroot<<" were requested in CI"<<endl;
294 break;
295 }
296 }
297 break;
298 }
299 }
300 is.close();
301
302 if (nstate < 1) {
303 cout << "Wrong input for -state . Setting STATE value to default (STATE=IROOT).\n" ;
304 nstate=iroot;
305 }
306 else if (nstate == 0) {
307 cout << "Setting STATE value to default (STATE=IROOT).\n" ;
308 nstate=iroot;
309 }
310 else if (nstate>computed_states) {
311 cout << "STATE value is higher than number of computed states in CI.\n" ;
312 cout << "Setting STATE value to default (STATE=IROOT).\n" ;
313 nstate=iroot;
314 }
315
316
317 if(using_guga){
318 if(symmetry)
319 cout << "Keeping the whole CSFs" <<endl;
320 vector <int> csf_full;
321 vector <int> csf;
322 vector <double> csf_weights;
323 vector <string> csf_occupations_str;
324
325
326 int line_position;
327 is.open(infilename.c_str());
328 while(getline(is, line)) {
329 words.clear();
330 split(line, space, words);
331 if(words[0]=="DAVIDSON" && words[1]=="METHOD"){
332 if(printout>1) cout <<"found DAVIDSON METHOD CI-MATRIX DIAGONALIZATION"<<endl;
333 //read in info on DAVIDSON METHOD CI-MATRIX DIAGONALIZATION
334 while(getline(is, line)) {
335 words.clear();
336 split(line, space, words);
337 if(words[0]=="ITER." && words[2]=="IMPROVED") {
338 while(getline(is, line)) {
339 if(printout>1) cout << line <<endl;
340 if(line=="") break;
341 }
342 }//---done energy info
343 //cout << line<<endl;
344 words.clear();
345 split(line, space, words);
346 //read in contributing CSFs
347 if (words[0]=="STATE"&& words[1]=="#" && words[3]=="ENERGY" && atoi(words[2].c_str())==nstate){
348 cout << "Printing CSFs for state : " << atoi(words[2].c_str()) << endl;
349 int i=0;
350 double smallest_csf_weight=1.0;
351 while(getline(is,line)) {
352 words.clear();
353 split(line, space, words);
354 //cout << "in state : " << line<<endl;
355 if(words.size() >0 and (words[0]=="......" || words[0]=="STATE")) break;
356
357 if(words.size()>1 && atoi(words[0].c_str())>0 && !(words[0]=="---") && !(words[0]=="CSF")){
358 csf.push_back(atoi(words[0].c_str())-1);
359 csf_weights.push_back(atof(words[1].c_str()));
360 csf_occupations_str.push_back(words[2]);
361 if(smallest_csf_weight > abs(csf_weights[i]))
362 smallest_csf_weight=abs(csf_weights[i]);
363 i++;
364 }
365 else if (words.size()==1){
366 csf_occupations_str.back()+=words[0];
367 }
368 }
369 if(csf.size()>0){
370 cout <<"Found "<<csf.size()<<" contributing CSFs with smallest weight of "<<smallest_csf_weight<<endl;
371 if(printout>1)
372 cout <<"MB: to increase printout set PRTTOL smaller than "<<smallest_csf_weight<<" in $GUGDIA"<<endl;
373 }
374 else{
375 cout <<"did not find any CSF COEF OCCUPANCY array in the DAVIDSON METHOD CI-MATRIX DIAGONALIZATION printout"<<endl;
376 exit(1);
377 }
378 }
379 }
380 }//if DAVIDSON METHOD CI-MATRIX DIAGONALIZATION
381 }
382 is.close();
383
384 vector <vector <double> > det_weights(csf.size());
385 vector <vector <string> > det_str(csf.size());
386
387
388
389 int counter_csf=0;
390 is.open(infilename.c_str());
391 while(getline(is, line)) {
392 words.clear();
393 split(line, space, words);
394 //here where all CSFs printout starts
395 if (words[0]=="DETERMINANT" && words[1]=="CONTRIBUTION") {
396 if(printout>1) cout <<"looping to read the relevant CSFs"<<endl;
397 int kk=0;
398 for(int k=0;k<csf.size();k++){
399 if(printout>1) cout <<"looking for CSF with index "<<csf[k]+1<<" .........";
400 while(getline(is,line)) {
401 words.clear();
402 split(line, space, words);
403 if(words[0]=="......" && words[1]=="END" && words[2]=="OF" && words[3]=="-DRT-" && words[4]=="GENERATION"){
404 cout <<"\n found the end of DRT GENERATION"<<endl;
405 if(kk!=csfmax){
406 cout <<"did not found enought CSF'S!, did you run gamess with NPRT=2 in $CIDRT?"<<endl;
407 exit(1);
408 }
409 break;
410 }
411 if(words[0]=="CASE" && words[1]=="VECTOR" && words[2]== "=" && words.size()==4){
412 if(atoi(words[3].c_str()) > csf[k]+1){
413 cout <<"error, the CSF with index "<<csf[k]+1<<" could not be found"<<endl;
414 exit(1);
415 }
416 if(atoi(words[3].c_str())-1==csf[k]){
417 if(printout>1) cout <<"found!"<<endl;
418 counter_csf++;
419 csf_full.push_back(atoi(words[4].c_str())-1);
420
421
422 //read in the determinants in:
423 double sum_of_weights=0.0;
424 int tries=0;
425 int max_tries=1000;
426
427 while(getline(is,line)) {
428 words.clear();
429 split(line, space, words);
430 /*
431 the old way
432 if ( line.size() > 33 ){
433 if (line.substr(13,2)=="C(" ){
434 if(printout>2) cout <<line<<endl;
435 det_weights[k].push_back(atof(line.substr(20,10).c_str()));
436 det_str[k].push_back(line.substr(33));
437 line_position=is.tellg();
438 //if there is a line which is too long this while has to be there
439 while(getline(is,line)) {
440 words.clear();
441 split(line, space, words);
442 if ((words[0]=="C(" ) || (words[0]=="CASE" && words[1]=="VECTOR") || line.size()==0){
443 is.seekg(line_position);
444 break;
445 }
446 else{
447 det_str[k].back() += " ";
448 det_str[k].back() +=line.substr(33);
449 }
450 }
451 }//if
452 }//if
453 */
454 //new way
455 if ( line.size() > 33 ){
456 if (line.substr(13,2)=="C(" ){
457 if(printout>2) cout <<line<<endl;
458 det_weights[k].push_back(atof(line.substr(20,10).c_str()));
459 det_str[k].push_back(line.substr(33));
460 //MB: new way of doing things: if there is a line which is too long this while has to be there
461 int sum_of_string_leghts=det_str[k].back().size();
462 while(sum_of_string_leghts < (nelectrons[0]+nelectrons[1]-2*number_of_core_orbitals)*3){
463 getline(is,line);
464 split(line, space, words);
465 string remainder=line.substr(33);
466 det_str[k].back() += " ";
467 det_str[k].back() +=remainder;
468 sum_of_string_leghts+=remainder.size();
469 }
470 sum_of_weights+=det_weights[k].back()*det_weights[k].back();
471 if(printout>3) cout <<"w: "<<det_weights[k].back()<<" at line number "<<line_position<<" sum of weights: "<<sum_of_weights<<endl;
472 if(abs(sqrt(sum_of_weights)-1.0)<1e-5)
473 break;
474 //M.B.: this was working before, but tellg() was giving negative number for files > 2Gb
475 //if there is a line which is too long this while has to be there
476 /*
477 line_position=is.tellg();
478 while(getline(is,line)) {
479 words.clear();
480 split(line, space, words);
481 tries++;
482 if ((words[0]=="C(" ) || (words[0]=="CASE" && words[1]=="VECTOR") || line.size()==0){
483 is.seekg(line_position);
484 cout <<"going back to line position" <<line_position<<endl;
485 break;
486 }
487 else{
488 if(printout>2) cout <<line<<endl;
489 det_str[k].back() += " ";
490 det_str[k].back() +=line.substr(33);
491 }
492 }//while
493 */
494
495 }
496 if(tries>max_tries){
497 cout<<"could not fully read CSF "<<csf[k]+1<<" ...exiting"<<endl;
498 exit(1);
499 }
500 tries++;
501 }//if(line.size() > 33)
502 }//end of while
503 if(!(abs(sqrt(sum_of_weights)-1.0)<1e-5)){
504 cout<<"could not fully read CSF "<<csf[k]+1<<" ...exiting"<<endl;
505 exit(1);
506 }
507 //need to break in order to find another CSF
508 break;
509 }//end of if CASE VECTOR = csf[k]
510 if(kk>csfmax){
511 cout<<"Number of CSF's is greater than designed limit of "<<csfmax<<"CSF's"<<endl;
512 exit(1);
513 }
514 kk++;
515 }
516 }//while(getline(is,line))
517 }//int k loop
518 }//end of reading all CSFs
519 }
520 is.close();
521
522 if(counter_csf!=csf.size()){
523 cout <<"counter_csf "<<"and "<<csf.size()<<"dont match "<<endl;
524 exit(1);
525 }
526 cout << "done readout"<<endl;
527 //end of read out
528 //storing determinants and reodering orbs
529 det_weights.resize(csf.size());
530 det_str.resize(csf.size());
531 vector < vector < vector <int> > > det(csf.size());
532
533 for(int i=0;i<det_weights.size();i++){
534 if(!det_weights[i].size()){
535 cout<<" CSF with index "<<csf[i]+1<<" was zero determinants in!, exitting"<<endl;
536 exit(1);
537 }
538
539 det[i].resize(det_weights[i].size());
540 for(int j=0;j<det_weights[i].size();j++){
541 for(int k=0; k < det_str[i][j].size(); k=k+3){
542 det[i][j].push_back(atoi(det_str[i][j].substr(k,3).c_str()));
543 }
544
545 reorder_orbs(det[i][j],det_weights[i][j]);
546 if(reorder)
547 normal_reorder_orbs(nelectrons[0],nelectrons[1], number_of_core_orbitals, det[i][j],det_weights[i][j]);
548
549 //if(i>0)
550 // reorder3(det[i][j],det_weights[i][j],det[0][0]);
551 // det_spinup[i].resize(det_weights[i].size());
552 // det_spindown[i].resize(det_weights[i].size());
553 // spin_separate(det[i][j],det_spinup[i][j],det_spindown[i][j] );
554 //cout << "Weight: "<<det_weights[i][j]<<" state: ";
555 //for(int k=0; k < det[i][j].size();k++){
556 // cout << det[i][j][k]<< " ";
557 // }
558 //cout <<endl;
559 }
560 }
561
562 cout << "done storing determinants and reodering orbs"<<endl;
563
564 //storing occupation array
565 vector < vector <int> > csf_occupation(csf.size());
566 for(int i=0;i< csf.size();i++){
567 for(int j=0;j< csf_occupations_str[i].size();j++){
568 csf_occupation[i].push_back(atoi(csf_occupations_str[i].substr(j,1).c_str()));
569 }
570 }
571 cout << "done storing occupation array"<<endl;
572
573 if(symmetry==0){
574 // assign full weight to each determimant
575 for(int i=0;i<csf.size();i++){
576 for(int j=0;j<det_weights[i].size();j++){
577 //if(csf[i]==847){
578 //cout << csf_weights[i]<< " "<<det_weights[csf[i]-1][j]<<endl;
579 //}
580 det_weights[i][j]*=csf_weights[i];
581 }
582 }
583 cout << "done assign full weight to each determimant"<<endl;
584
585 //find the same determimants and add their contributions
586 for(int i=0;i<csf.size();i++){
587 //cout << csf[i]<<" "<< csf_occupations_str[i]<<endl;
588 for(int j=i+1;j<csf.size();j++){
589 if(csf_occupations_str[i]==csf_occupations_str[j]){
590 //cout << "CSF "<<csf[i]<<" and "<<csf[j]<< " are the same"<<endl;
591 for(int k=0;k<det[i].size();k++)
592 for(int l=0;l<det[j].size();l++)
593 if(det[i][k]==det[j][l]){
594 det_weights[i][k]+=det_weights[j][l];
595 det_weights[j][l]=0.0;
596 //cout << "determimant"<<k<<" and "<<l<<" are the same"<<endl;
597 //for(int h=0;h<det[csf[i]-1][k].size();h++)
598 // cout <<det[csf[i]-1][k][h]<<" ";
599 //cout <<endl;
600 //for(int h=0;h<det[csf[j]-1][k].size();h++)
601 // cout <<det[csf[j]-1][l][h]<<" ";
602 //cout <<endl;
603 }
604 }
605 }
606 }
607 cout << "done finding the same determimants and adding their contributions"<<endl;
608 }
609
610
611 //use weight treshold and store for printout
612 vector <double> det_weights_printing;
613 //int ircounter=0;
614 //vector <int> det_weights_bonds;
615 vector < vector <int> > det_printing;
616 vector < vector <int> > det_printing_old;
617 vector <double> csf_weights_printing;
618 vector <int> csf_printing;
619 vector < vector <double> > det_weights_printing_symmetry;
620 for(int i=0;i<csf.size();i++){
621 if(symmetry){
622 if(abs(csf_weights[i])> wtresh){
623 csf_weights_printing.push_back(csf_weights[i]);
624 csf_printing.push_back(i);
625 //det_weights_printing_symmetry.push_back(det_weights[csf[]]);
626 //for(int k=0;k<det[csf[i]].size();k++){
627 // det_weights_printing.push_back(det_weights[csf[i]][k]);
628 //det_printing.push_back(det[csf[i]][k]);
629 //det_weights_bonds.push_back(ircounter);
630 //}
631 //ircounter++;
632 }
633 }
634 else{
635 for(int k=0;k<det[i].size();k++){
636 if(abs(det_weights[i][k])> wtresh){
637 // cout <<"det_weights "<< det_weights[csf[i]][k]<<endl;;
638 det_weights_printing.push_back(det_weights[i][k]);
639 det_printing_old.push_back(det[i][k]);
640 }
641 }
642 }
643 }
644
645 //sorting
646 vector <int> list;
647 if(symmetry){
648 sort_abs_largest_first(csf_weights_printing, list);
649 for(int i=0;i<csf_weights_printing.size();i++){
650 det_weights_printing_symmetry.push_back(det_weights[csf_printing[list[i]]]);
651 for(int k=0;k<det[csf_printing[list[i]]].size();k++){
652 det_weights_printing.push_back(det_weights[csf_printing[list[i]]][k]);
653 det_printing.push_back(det[csf_printing[list[i]]][k]);
654 }
655 }
656 }
657 else{
658 sort_abs_largest_first(det_weights_printing, list);
659 for (int i=0;i<det_weights_printing.size();i++){
660 det_printing.push_back(det_printing_old[list[i]]);
661 }
662
663 }
664
665
666 if(symmetry){
667 cout << "After cutoff applied: found "<<det_weights_printing.size()<<" determinats with weights"<<endl;
668 cout << "and number of independent weights "<<csf_weights_printing.size()<<endl;
669 }
670 else
671 cout << "After cutoff applied: found "<<det_weights_printing.size()<<" unique determinats with weights"<<endl;
672
673
674 vector < vector <int> > det_up(det_weights_printing.size());
675 vector < vector <int> > det_down(det_weights_printing.size());
676 int max_up, max_down, tmp_up,tmp_down;
677 max_up=max_down=0;
678 for (int i=0;i<det_weights_printing.size();i++){
679 spin_separate(det_printing[i], det_up[i], det_down[i], tmp_up, tmp_down);
680 if(tmp_up>max_up)
681 max_up=tmp_up;
682 if(tmp_down>max_down)
683 max_down=tmp_down;
684 }
685 cout <<"the largest spin-up orbital is "<<max_up<<" and spin-down orbital is "<<max_down<<endl;
686
687
688
689
690 //FINAL OUTPUT
691 ofstream output(outputname.c_str());
692 string gap=" ";
693
694 if(symmetry){
695 output <<gap<< "# NCSF = "<<csf_weights_printing.size()<<endl;
696 for (int i=0;i<csf_weights_printing.size();i++){
697 output.precision(7);
698 output.width(10);
699 output <<gap<<"CSF { "<< csf_weights_printing[i] <<" ";
700 for(int j=0;j<det_weights_printing_symmetry[i].size();j++){
701 output.precision(7);
702 output.width(10);
703 output << det_weights_printing_symmetry[i][j]<<" ";
704 }
705 output << "}"<<endl;
706 }
707 }
708 else {
709 output <<gap<< "DETWT {"<<endl;
710 output <<gap<< "# NDET = "<<det_weights_printing.size()<<endl;
711 output <<gap;
712 for (int i=0;i<det_weights_printing.size();i++){
713 output.precision(7);
714 output.width(10);
715 output << det_weights_printing[i] <<" ";
716 if ((i+1)%6==0)
717 output <<endl<<gap;
718 }
719 output << "}"<<endl;
720 }
721 output <<gap<<"STATES {"<<endl<<gap;
722 if(symmetry){
723 int counter=0;
724 for (int i=0;i<csf_weights_printing.size();i++){
725 output <<"# CSF "<<i+1<<": weight: "<<csf_weights_printing[i]<<endl<<gap;
726 for(int j=0;j<det_weights_printing_symmetry[i].size();j++){
727 output <<"# Determinant "<<j+1<<": weight: "<<det_weights_printing_symmetry[i][j]<<endl<<gap;
728 for(int m=0;m<number_of_core_orbitals;m++)
729 output << m+1 <<" ";
730 for(int k=det_up[counter].size();k>0;k--)
731 output << det_up[counter][k-1] <<" ";
732 output <<endl<<gap;
733 for(int m=0;m<number_of_core_orbitals;m++)
734 output << m+1 <<" ";
735 for(int k=0;k<det_down[counter].size();k++)
736 output << det_down[counter][k] <<" ";
737 output <<endl<<gap;
738 counter++;
739 }
740 }
741 }
742 else {
743 for (int i=0;i<det_weights_printing.size();i++){
744 output <<"# Determinant "<<i+1<<": weight: "<<det_weights_printing[i]<<endl<<gap;
745 for(int j=0;j<number_of_core_orbitals;j++)
746 output << j+1 <<" ";
747 for(int k=det_up[i].size();k>0;k--)
748 output << det_up[i][k-1] <<" ";
749 output <<endl<<gap;
750 for(int j=0;j<number_of_core_orbitals;j++)
751 output << j+1 <<" ";
752 for(int k=0;k<det_down[i].size();k++)
753 output << det_down[i][k] <<" ";
754 output <<endl<<gap;
755 }
756 }
757 output << "}"<<endl;
758 output.close();
759 cout <<"done"<<endl;
760 }//end using guga
761 else if(using_determinants){
762 vector <double> det_weights;
763 vector < vector <string> > det_occupations_str(2);
764 int ncorr_orbs;
765 int nacc_orbs;
766 vector <int> nelectrons(2);
767 int orbs_total;
768 int ndet;
769 int cycle=1;
770 if(symmetry)
771 cout << "ignoring option -csf" <<endl;
772 is.open(infilename.c_str());
773 while(getline(is, line) && cycle<2 ) {
774 words.clear();
775 split(line, space, words);
776
777 if(words[0]=="NUMBER" && words[2]=="CORE" && words[3]=="ORBITALS"){
778 ncorr_orbs=atoi(words[5].c_str());
779 cout <<ncorr_orbs<<" core orbitals"<<endl;
780 }
781 if(words[0]=="NUMBER" && words[2]=="ACTIVE" && words[3]=="ORBITALS"){
782 nacc_orbs=atoi(words[5].c_str());
783 cout <<nacc_orbs<<" active orbitals"<<endl;
784 }
785
786 if(words[0]=="NUMBER" && words[2]=="ALPHA" && words[3]=="ELECTRONS"){
787 nelectrons[0]=atoi(words[5].c_str());
788 cout <<nelectrons[0]<<" alpha electrons"<<endl;
789 }
790 if(words[0]=="NUMBER" && words[2]=="BETA" && words[3]=="ELECTRONS"){
791 nelectrons[1]=atoi(words[5].c_str());
792 cout <<nelectrons[1]<<" beta electrons"<<endl;
793 }
794 if(words[0]=="NUMBER" && words[2]=="OCCUPIED" && words[3]=="ORBITALS" && words[4]=="="){
795 orbs_total=atoi(words[5].c_str());
796 cout <<orbs_total<<" total number orbitals"<<endl;
797 }
798
799 if(words[0]=="ITERATION" && words[1]=="ENERGY" && words[2]=="GRADIENT") {
800 while(getline(is, line)) {
801 words.clear();
802 split(line, space, words);
803 cout << line <<endl;
804 if(words[0]=="ALL" && words[1]=="STATES" && words[2]=="CONVERGED.") break;
805 }
806 }//---done energy info
807 //read in actual determinamnts
808
809 if (words[0]=="STATE" && words[2]=="ENERGY=" && atoi(words[1].c_str())==nstate){
810 while(getline(is, line)) {
811 words.clear();
812 split(line, space, words);
813 ndet=0;
814 if(words[0]=="ALPHA" && words[2]=="BETA" && words[4]=="COEFFICIENT") {
815 while(getline(is, line)) {
816 words.clear();
817 split(line, space, words);
818 if(words[1]=="|" && words[3]=="|"){
819 // cout << words[0]<<" "<<words[2]<<" "<<words[4]<<endl;
820 det_weights.push_back(atof(words[4].c_str()));
821 det_occupations_str[0].push_back(words[0]);
822 det_occupations_str[1].push_back(words[2]);
823 ndet++;
824 }
825 //..... DONE WITH DETERMINANT CI COMPUTATION .....
826 //find the end
827 if(words[0]=="....." || line.size()==0) {
828 cout <<"found "<<ndet<<" determinants"<<endl;
829 cycle++;
830 break;
831 }
832 }
833 break;
834 }
835 }
836 }
837 }
838 is.close();
839 cout <<"done readout"<<endl;
840
841 vector < vector < vector <int> > > det_occupations(2);
842
843 for(int spin=0;spin<2;spin++){
844 det_occupations[spin].resize(ndet);
845 for(int det=0;det<ndet;det++){
846 for(int j=0;j<ncorr_orbs;j++){
847 det_occupations[spin][det].push_back(j+1);
848 }
849 for(int j=0;j<nacc_orbs;j++){
850 if(atoi(det_occupations_str[spin][det].substr(j,1).c_str())==1)
851 det_occupations[spin][det].push_back(ncorr_orbs+j+1);
852 }
853 }//ndets
854 }//spin
855
856
857
858
859 //select determinants
860 vector <double> det_weights_print;
861 vector < vector < vector <int> > > det_occupations_print(2);
862 for(int det=0;det<ndet;det++){
863 if(fabs(det_weights[det])>wtresh){
864 det_weights_print.push_back(det_weights[det]);
865 for(int spin=0;spin<2;spin++){
866 det_occupations_print[spin].push_back(det_occupations[spin][det]);
867 }
868 }
869 else{
870 //cout <<" det "<<det<<endl;
871 }
872 }
873 cout <<"done selecting determinants"<<endl;
874 cout <<det_weights_print.size()<<" determinants has weight larger than "<<wtresh<<endl;
875
876 vector <int> list;
877 sort_abs_largest_first(det_weights_print, list);
878
879 //final printout
880 ofstream output(outputname.c_str());
881 string gap=" ";
882 output <<gap<< "DETWT {"<<endl;
883 output <<gap<< "# NDET = "<<det_weights_print.size()<<endl;
884 output <<gap;
885 for (int i=0;i<det_weights_print.size();i++){
886 output.precision(7);
887 output.width(10);
888 output << det_weights_print[i] <<" ";
889 if ((i+1)%6==0)
890 output <<endl<<gap;
891 }
892 output << "}"<<endl;
893 output <<gap<<"STATES {"<<endl<<gap;
894 for (int i=0;i<det_weights_print.size();i++){
895 output <<"# Determinant "<<i+1<<": weight: "<<det_weights_print[i]<<endl<<gap;
896 for(int spin=0;spin<2;spin++){
897 for(int j=0;j<det_occupations_print[spin][list[i]].size();j++)
898 output <<det_occupations_print[spin][list[i]][j]<<" ";
899 output <<endl<<gap;
900 }
901 }
902 output << "}"<<endl;
903
904
905 output.close();
906 cout <<"done"<<endl;
907
908
909
910
911 }
912 else{
913 cout <<"Could not find GUGA nor ALDET type CI";
914 exit(1);
915 }
916 }
917