1 /*
2
3 Copyright (C) 2007 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 #include "converter.h"
21 #include "wf_writer.h"
22 #include "basis_writer.h"
23 #include "Pseudo_writer.h"
24 #include <algorithm>
25 #include <fstream>
26 #include <cstring>
27 #include "vecmath.h"
28 #include <string>
29 #include "elements.h"
30 #include <stdlib.h>
31 #include <stdio.h>
32 using namespace std;
33
34 void get_crystal_pseudo(istream & infile,
35 vector <Gaussian_pseudo_writer> & pseudo);
36
37 void get_crystal_latvec(istream & infile,
38 vector< vector < double > > & latvec);
39
40 void get_crystal_atoms(istream & infile,
41 vector < Atom > & atoms);
42
43 void get_crystal_basis(istream & infile,
44 vector <Gaussian_basis_set> & basis);
45
46 // for orbitals with real coefficients
47 void read_crystal_orbital(istream & is,
48 string & fort10file,
49 vector <Atom> & atoms,
50 Slat_wf_writer & slwriter,
51 vector <Gaussian_basis_set> & basis,
52 vector <double> & origin,
53 vector < vector <double> >& latvec,
54 vector < vector < double > > & moCoeff,
55 vector <double> & shift //amount to shift the atoms
56 );
57 void read_kpt_eigenpos(istream & is,
58 vector <string> & rkpoints,
59 vector <long int> & reigen_start,
60 vector <string> & ckpoints,
61 vector <long int> & ceigen_start);
62 // for orbitals with complex coefficients
63 void read_crystal_orbital(istream & is,
64 vector <Atom> & atoms,
65 Slat_wf_writer & slwriter,
66 vector <Gaussian_basis_set> & basis,
67 vector <double> & origin,
68 vector < vector <double> >& latvec,
69 vector < vector < dcomplex > > & moCoeff,
70 vector <double> & shift //amount to shift the atoms
71 );
72
73 void read_crystal_orbital_all(istream & is,
74 string & fort10file,
75 vector <Atom> &atoms,
76 Slat_wf_writer & slwriter,
77 vector <Gaussian_basis_set> & basis,
78 vector <double> & origin,
79 vector < vector <double> >& latvec,
80 vector < vector < vector < double > > > & moCoeff,
81 vector <string> &kpoints,
82 vector < long int > & eigen_start,
83 vector <vector <double> > &kptCoord,
84 vector <double> & shift, //amount to shift the atoms
85 vector < int > shifted,
86 vector < vector <int> > nshift
87 );
88
89 // for orbitals with complex coefficients
90 void read_crystal_orbital_all(istream & is,
91 vector <Atom> &atoms,
92 Slat_wf_writer & slwriter,
93 vector <Gaussian_basis_set> & basis,
94 vector <double> & origin,
95 vector < vector <double> >& latvec,
96 vector < vector < vector < dcomplex > > > & moCoeff,
97 vector <string> &kpoints,
98 vector < long int > & eigen_start,
99 vector < vector <double> > &kptCoord,
100 vector <double> & shift, //amount to shift the atoms
101 vector < int > shifted,
102 vector < vector <int> > nshift
103 );
104
105
106
usage(const char * name)107 void usage(const char * name) {
108 cout << "usage: " << name << " <options> <crystal output> " << endl;
109 cout << "Where options can be: \n\n";
110 cout << "-c Work with complex orbitals. Only real are searched by default.\n\n";
111 cout << "-fort10file Formatted fortran unit 10 from read10 utility.\n";
112 cout << " try to match the vectors in fort.10 to the ones in the output.\n";
113 cout << " It may or may not work correctly..\n\n";
114 cout << "-o Base name for your run case\n";
115 exit(1);
116 }
117
main(int argc,char ** argv)118 int main(int argc, char ** argv) {
119
120 bool had_error=false;
121 bool cmplx=false;
122 string infilename;
123 string outputname, fort10file;
124 vector <double> shift(3);
125 for(int d=0; d< 3; d++) shift[d]=.125;
126
127 for(int i=1; i< argc-1; i++) {
128 if(!strcmp(argv[i], "-o") && argc > i+1) {
129 outputname=argv[++i];
130 }
131 else if(!strcmp(argv[i],"-fort10file")) {
132 if (i+1 < argc) {
133 fort10file=argv[i+1];
134 i++;
135 }
136 else {
137 cout << "-fort10file needs an argument" << endl;
138 exit(1);
139 }
140 }
141 else if(!strcmp(argv[i],"-shift")){
142 if(i+1 < argc) {
143 for(int d=0; d< 3; d++) { shift[d]=atof(argv[i+1]);
144 i++;
145 }
146 }
147 else {
148 cout << "-shiftz needs an argument" << endl;
149 exit(1);
150 }
151 }
152 else if(!strcmp(argv[i], "-noshift")) {
153 for(vector<double>::iterator i=shift.begin();
154 i!= shift.end(); i++) *i=0.0;
155 }
156 else if(!strcmp(argv[i], "-c")) {
157 cmplx=true;
158 }
159
160 else {
161 cout << "Didn't understand option " << argv[i] << endl;
162
163 had_error=true;
164 }
165 }
166
167 if(argc >= 2) {
168 infilename=argv[argc-1];
169 }
170 else { had_error=true; }
171
172 if(outputname == "") {
173 outputname="qwalk";
174 //outputname=infilename;
175 }
176
177 if(had_error) usage(argv[0]);
178
179 if ( cmplx && ( fort10file != "" ) ) {
180 cout << "-fort10file is not supported for complex orbitals (-c)." << endl;
181 exit(1);
182 }
183
184
185 ///////////////////////////////////////////////////////////////
186
187 vector <Atom> atoms;
188 vector <Gaussian_pseudo_writer > pseudo;
189 vector <Gaussian_basis_set > basis;
190
191 Slat_wf_writer slwriter;
192 slwriter.use_global_centers=true;
193 slwriter.write_centers=false;
194 if ( cmplx ) {
195 slwriter.mo_matrix_type="CUTOFF_MO";
196 slwriter.orbtype="CORBITALS";
197 } else {
198 slwriter.mo_matrix_type="CUTOFF_MO";
199 }
200
201 vector < vector < double > > latvec;
202
203
204 ifstream infile(infilename.c_str());
205 if(!infile) {
206 cout << "couldn't open " << infilename << endl;
207 exit(1);
208 }
209 get_crystal_latvec(infile, latvec);
210 infile.close();
211 infile.clear();
212
213 infile.open(infilename.c_str());
214 get_crystal_atoms(infile, atoms);
215 infile.close();
216 infile.clear();
217 if(latvec.size() > 0) {
218 vector <double> atomshift;
219 for(int i=0; i< 3; i++) atomshift.push_back(0);
220 for(int i=0; i< 3; i++) {
221 for(int j=0; j< 3; j++) {
222 atomshift[j]+=shift[i]*latvec[i][j];
223 }
224 }
225 for(int i=0; i < atoms.size(); i++) {
226 atoms[i].pos=atoms[i].pos+atomshift;
227 }
228 }
229 Shifter shiftobj;
230 vector < vector <int> > nshift;
231 nshift.resize(atoms.size());
232 vector <int> shifted(atoms.size());
233 if(latvec.size() > 0) {
234 for (int at = 0; at<atoms.size(); at++)
235 shifted[at]=shiftobj.enforcepbc(atoms[at].pos, latvec, nshift[at]);
236 }
237 infile.open(infilename.c_str());
238 get_crystal_basis(infile, basis);
239 infile.close();
240 infile.clear();
241 infile.open(infilename.c_str());
242 get_crystal_pseudo(infile, pseudo);
243 infile.close();
244 infile.clear();
245
246 int nelectrons=-1;
247 double totspin=-1e8;
248 infile.open(infilename.c_str());
249 string calctype, testword;
250 int nelectrons_up=-1, nelectrons_down=-1;
251 string line, space=" ";
252 vector <string> spl;
253 while(true) {
254 spl.clear();
255 getline(infile,line);
256 split(line,space,spl);
257 if(spl.size() > 2 and spl[1]=="SCF" and spl[2]=="ENDED") break;
258 if(spl.size() > 2 and spl[0]=="TOTAL" and spl[1]=="ATOMIC" and spl[2]=="SPINS") {
259 spl.clear();
260 while(true) {
261 getline(infile,line);
262 if(line[3]=='T') break;
263 split(line,space,spl);
264 }
265 totspin=0.0;
266 for(vector<string>::iterator i=spl.begin(); i!=spl.end(); i++) {
267 totspin+=atof(i->c_str());
268 }
269 }
270 if(spl.size() > 4 and spl[0]=="N." and spl[2]=="ELECTRONS" and spl[4]=="CELL") {
271 nelectrons=atof(spl[5].c_str());
272 }
273 if(spl.size() > 4 and spl[0]=="TYPE" and spl[2]=="CALCULATION") {
274 if(spl[4]=="RESTRICTED") {
275 if(spl[5]=="CLOSED") calctype="RHF";
276 else if(spl[5]=="OPEN") calctype="ROHF";
277 }
278 else if(spl[4]=="UNRESTRICTED") calctype="UHF";
279 else {
280 cout << "Didn't understand calctype " << spl[4] << endl;
281 exit(1);
282 }
283 }
284 }
285 slwriter.calctype=calctype;
286 infile.close();
287
288 if(atoms.size() == 0) {
289 cout << "I couldn't find any atoms!!" << endl;
290 exit(1);
291 }
292
293 //----------------------------------------------------------------------
294 //Done parsing, now organize and link.
295
296 //Figure out which atoms go with which basis sets.
297 int natoms=atoms.size();
298 int nbasis=basis.size();
299 for(int at=0; at<natoms; at++) {
300 int found_basis=0;
301 for(int bas=0; bas< nbasis; bas++) {
302
303 if(atoms[at].name == basis[bas].label) {
304 atoms[at].basis=bas;
305 found_basis=1;
306 break;
307 }
308 }
309 if(!found_basis) {
310 cout << "Couldn't find basis for atom " << atoms[at].name
311 << endl;
312 exit(1);
313 }
314 }
315
316
317 //We need to get the effective charges from the user, since
318 //Crystal isn't nice enough to output them in any consistent way.
319 int testelectrons=0;
320
321 testelectrons=0;
322 for(int at=0; at < natoms; at++) {
323 int isUnique=1;
324 for(unsigned int i=0; i< at; i++) {
325 if(atoms[at].name == atoms[i].name) {
326 isUnique=0;
327 atoms[at].charge=atoms[i].charge;
328 break;
329 }
330 }
331 if(isUnique) {
332 int npseudo=pseudo.size();
333 bool found_in_pseudo=false;
334 for(int j=0; j< npseudo; j++) {
335 if(atoms[at].name==pseudo[j].label) {
336 atoms[at].charge=pseudo[j].effcharge;
337 found_in_pseudo=true;
338 break;
339 }
340 }
341 if(!found_in_pseudo) {
342 cout << "Please enter the effective charge of " << atoms[at].name
343 << endl;
344 cin >> atoms[at].charge;
345 cout << "Thanks" << endl;
346 }
347 }
348 testelectrons+=(int) atoms[at].charge;
349 }
350 if(testelectrons != nelectrons) {
351 cout << "The total number of electrons should be " << nelectrons
352 << ", \nbut the sum due to the effective charge on the ion is "
353 << testelectrons << ". This may be what you intended, but be careful!\n";
354 }
355
356
357
358 //Determine how many electrons are spin up or down if it's not RHF
359 if(calctype=="RHF" ) {
360 if(nelectrons % 2 != 0) {
361 cout << "It seems like you're doing RHF, but there is an odd number of"
362 << " electrons! I'm confused." << endl;
363 exit(1);
364 }
365 slwriter.nup=nelectrons/2;
366 slwriter.ndown=nelectrons/2;
367 }
368 else if(calctype=="ROHF" || calctype=="UHF") {
369 cout << "Detected a ROHF or UHF calculation." << endl;
370 // " What is the spin state?(1=singlet, 2=doublet...)" << endl;
371 cout << "N_up-N_down is found to be " << totspin << endl;
372 int spin=totspin+0.1;
373 if( abs(totspin-spin) > 1e-4) {
374 cout << "WARNING: spin is not close to an integer!" << endl;
375 }
376 slwriter.nup=(nelectrons-spin)/2 + spin;
377 slwriter.ndown=(nelectrons-spin)/2;
378 if(slwriter.nup+slwriter.ndown != nelectrons) {
379 cout << "problem.. nup and ndown are calculated to be "
380 << slwriter.nup << " " << slwriter.ndown
381 << " but they don't add up to be " << nelectrons << endl;
382 exit(1);
383 }
384 }
385 else {
386 cout << "Don't understand calculation type of "
387 << calctype << endl;
388 exit(1);
389 }
390
391 vector <double> origin(3);
392 for(int i=0; i< 3; i++) origin[i]=0;
393
394 vector <vector < vector < double > > > moCoeff;
395 vector < vector < vector < dcomplex > > > CmoCoeff;
396 infile.clear();
397 infile.open(infilename.c_str());
398 vector <string> rkpts;
399 vector <string> ckpts;
400 vector < vector <double> > ckptCoord;
401 vector < vector <double> > rkptCoord;
402 vector <long int > ceigen_start, reigen_start;
403 read_kpt_eigenpos(infile, rkpts, reigen_start, ckpts, ceigen_start);
404 if(cmplx)
405 read_crystal_orbital_all(infile, atoms, slwriter, basis,
406 origin, latvec, CmoCoeff, ckpts, ceigen_start, ckptCoord, shift, shifted, nshift);
407 infile.close();
408 infile.clear();
409 infile.open(infilename.c_str());
410 read_crystal_orbital_all(infile, fort10file, atoms, slwriter, basis,
411 origin, latvec, moCoeff, rkpts, reigen_start, rkptCoord, shift, shifted, nshift);
412 infile.close();
413 natoms=atoms.size();
414
415
416 vector <Center> centers;
417 centers.resize(atoms.size());
418 for(int at=0; at < natoms; at++) {
419 for(int i=0; i< 3; i++) centers[at].pos[i]=atoms[at].pos[i];
420 centers[at].equiv_atom=at;
421 centers[at].name=atoms[at].name;
422 }
423
424 //-------------------------------
425 //print out the qmc input file
426
427 string orboutname=outputname+".orb";
428 slwriter.orbname=orboutname;
429 string basisoutname=outputname+".basis";
430 slwriter.basisname=basisoutname;
431
432 vector <int> nbasis_list;
433 nbasis_list.resize(natoms);
434 for(int at=0; at < natoms; at++) {
435 nbasis_list[at]=basis[atoms[at].basis].nfunc();
436 }
437
438 ofstream basisout(basisoutname.c_str());
439 int nbas=basis.size();
440 for(int bas=0; bas < nbas; bas++) {
441 basisout << "BASIS { \n";
442 basis[bas].print_basis(basisout);
443 basisout << "}\n\n\n";
444 }
445 basisout.close();
446
447
448
449 //--------------------------Jastrow 2 output
450
451 double basis_cutoff;
452 if(latvec.size() >0)
453 basis_cutoff=find_basis_cutoff(latvec);
454 else basis_cutoff=7.5;
455
456 Jastrow2_wf_writer jast2writer;
457 jast2writer.set_atoms(atoms);
458 string jast2outname=outputname+".jast2";
459 ofstream jast2out(jast2outname.c_str());
460 print_std_jastrow2(jast2writer, jast2out, basis_cutoff);
461 jast2out.close();
462
463 // Jastrow 3 output
464 string jast3outname=outputname+".jast3";
465 ofstream jast3out(jast3outname.c_str());
466 vector<string> unique_atoms;
467 find_unique_atoms(atoms, unique_atoms);
468 print_3b_jastrow2(jast3out,unique_atoms,basis_cutoff);
469 jast3out.close();
470
471 //------------------------------------------System output
472
473 double min=1e8;
474 for(vector<Gaussian_basis_set>::iterator bas=basis.begin();
475 bas != basis.end(); bas++) {
476 for(vector <vector<double> >::iterator i=bas->exponents.begin();
477 i!=bas->exponents.end(); i++) {
478 for(vector<double>::iterator j=i->begin(); j!= i->end(); j++) {
479 if(*j < min) min=*j;
480 }
481 }
482 }
483 cout << "minimum exponent " << min << endl;
484 double cutoff_length= sqrt(-log(1e-8)/min);
485 cout << "cutoff length " << cutoff_length << endl;
486
487 cout << "Writing real orbitals ... " << endl;
488 for (int kpt=0; kpt<rkpts.size(); kpt++) {
489 slwriter.orbtype="ORBITALS";
490
491 vector <string> words;
492 string space = " ";
493 split(rkpts[kpt], space, words);
494 stringstream lk;
495 lk << atoi(words[0].c_str())-1;
496
497 string lst = lk.str();
498 string forb=outputname+"_"+lst+".orb";
499 cout << " " << forb + ": ";
500 ofstream orbout(forb.c_str());
501 print_orbitals(orbout, centers, nbasis_list, moCoeff[kpt]);
502 orbout.close();
503 string slateroutname=outputname+"_"+lst+".slater";
504 slwriter.orbname=forb;
505 ofstream slaterout(slateroutname.c_str());
506 slwriter.print_wavefunction(slaterout);
507 slaterout.close();
508
509 string sysoutname=outputname+"_"+lst+".sys";
510 ofstream sysout(sysoutname.c_str());
511 sysout.precision(12);
512 sysout << "SYSTEM { ";
513 if(latvec.size() > 0) sysout << " PERIODIC \n";
514 else sysout << " MOLECULE \n";
515 sysout << " NSPIN { " << slwriter.nup << " "
516 << slwriter.ndown << " } \n";
517
518
519 if(latvec.size() > 0) {
520 sysout << "LATTICEVEC { \n";
521 for(int i=0; i< 3; i++) {
522 for(int j=0; j< 3; j++) {
523 sysout << latvec[i][j] << " ";
524 }
525 sysout << endl;
526 }
527 sysout << " } " << endl;
528 sysout << " origin { " << origin[0] << " "
529 << origin[1] << " " << origin[2] << " } " << endl;
530 sysout << " cutoff_divider "
531 << basis_cutoff*2.0/cutoff_length << endl;
532 sysout << " kpoint { " << rkptCoord[kpt][0]
533 << " " << rkptCoord[kpt][1]
534 << " " << rkptCoord[kpt][2] << " } " << endl;
535 }
536
537 for(int at=0; at <natoms; at++) {
538 atoms[at].print_atom(sysout);
539 }
540 sysout << "}\n\n\n";
541
542
543
544 int npsp=pseudo.size();
545 for(int psp=0; psp < npsp; psp++) {
546 pseudo[psp].print_pseudo(sysout);
547 }
548 sysout.close();
549
550 }
551 // cout << "Now writing complex orbitals" << endl;
552 if(cmplx) {
553 cout << "Writing complex orbitals ... " << endl;
554 for (int kpt=0; kpt<ckpts.size(); kpt++) {
555 slwriter.orbtype="CORBITALS";
556 vector <string> words;
557 string space = " ";
558 split(ckpts[kpt], space, words);
559 // string lst=string(atoi(words[0].c_str())-1);
560 stringstream lk;
561 lk << atoi(words[0].c_str())-1;
562 string lst = lk.str();
563 string forb=outputname+"_"+lst+".orb";
564 cout << " " << forb << ": ";
565 ofstream orbout(forb.c_str());
566 print_orbitals(orbout, centers, nbasis_list, CmoCoeff[kpt]);
567 orbout.close();
568 string slateroutname=outputname+"_"+lst+".slater";
569 slwriter.orbname=forb;
570 ofstream slaterout(slateroutname.c_str());
571 slwriter.print_wavefunction(slaterout);
572 slaterout.close();
573
574 string sysoutname=outputname+"_"+lst+".sys";
575 ofstream sysout(sysoutname.c_str());
576 sysout.precision(12);
577 sysout << "SYSTEM { ";
578 if(latvec.size() > 0) sysout << " PERIODIC \n";
579 else sysout << " MOLECULE \n";
580 sysout << " NSPIN { " << slwriter.nup << " "
581 << slwriter.ndown << " } \n";
582
583
584 if(latvec.size() > 0) {
585 sysout << "LATTICEVEC { \n";
586 for(int i=0; i< 3; i++) {
587 for(int j=0; j< 3; j++) {
588 sysout << latvec[i][j] << " ";
589 }
590 sysout << endl;
591 }
592 sysout << " } " << endl;
593 sysout << " origin { " << origin[0] << " "
594 << origin[1] << " " << origin[2] << " } " << endl;
595 sysout << " cutoff_divider "
596 << basis_cutoff*2.0/cutoff_length << endl;
597 sysout << " kpoint { " << ckptCoord[kpt][0]
598 << " " << ckptCoord[kpt][1]
599 << " " << ckptCoord[kpt][2] << " } " << endl;
600 }
601
602 for(int at=0; at <natoms; at++) {
603 atoms[at].print_atom(sysout);
604 }
605 sysout << "}\n\n\n";
606
607
608
609 int npsp=pseudo.size();
610 for(int psp=0; psp < npsp; psp++) {
611 pseudo[psp].print_pseudo(sysout);
612 }
613 sysout.close();
614 }
615 }
616
617
618 }
619
620 //######################################################################
621
get_crystal_latvec(istream & infile,vector<vector<double>> & latvec)622 void get_crystal_latvec(istream & infile,
623 vector< vector < double > > & latvec) {
624 string testword;
625 while(infile >> testword) {
626 //cheesy way to get the converter to work for a molecular case
627 //I don't know if it would work in general.
628 // if(testword == "MOLECULAR") {
629
630 // vector <double> tmp;
631 // for(int i=0; i< 3; i++) {
632 //tmp.push_back(100);
633 // }
634 // for(int i=0; i< 3; i++) latvec.push_back(tmp);
635 // break;
636 //}
637 if(testword == "DIRECT") {
638 infile >> testword; //LATTICE
639 infile >> testword; //VECTORS
640 infile >> testword; //COMPON.
641 infile >> testword; //(A.U.)
642 if(testword=="(A.U.)") {
643 //cout << "found lattice parms" << endl;
644 infile.ignore(150, '\n');
645 infile >> testword;
646 //cout << testword << endl;
647 infile >> testword;
648 //cout << testword << endl;
649 infile.ignore(150,'\n');
650 latvec.reserve(3);
651 for(int i=0; i< 3; i++) {
652 latvec[i].reserve(3);
653 vector <double> tmp;
654 for(int j=0; j< 3; j++) {
655 double dummy;
656 infile >> dummy;
657 //here we try to avoid having issues with the Ewald summation
658 if(dummy > 100.) {
659 cout << "WARNING: rescaling lattice vector to " << 100 << " from " << dummy
660 << ". This will probably not change your results, but if it does, change EWALD_GMAX in the system input to something higher." << endl;
661 dummy=100.;
662 }
663 tmp.push_back(dummy);
664 }
665 latvec.push_back(tmp);
666 //cout << endl;
667 infile.ignore(150, '\n');
668 }
669 break;
670 }
671 }
672 }
673
674 if(latvec.size() == 0) {
675
676 //cout << "Couldn't find lattice vector in output file." << endl;
677 //exit(1);
678 }
679 }
680
681 //######################################################################
682
erasetail(string s)683 string erasetail(string s) {
684 char ms[] = "1";
685 string ns="";
686 int i=0;
687 while (s[i]!='1' and i < s.length()) {
688 i++;
689 }
690 return s.substr(0, i);
691 }
692
693
get_crystal_atoms(istream & infile,vector<Atom> & atoms)694 void get_crystal_atoms(istream & infile,
695 vector < Atom > & atoms) {
696 string line;
697 string space=" ";
698 Atom temp_atom;
699 double bohr=0.5291772083;
700
701 while(getline(infile, line)) {
702 vector <string> words;
703 split(line, space, words);
704
705 //2003 detection
706 if( words.size() >= 5 &&
707 words[1]=="ATOM" && words[2]=="X(ANGSTROM)") {
708
709
710
711 //cout << "found atom section " << line << endl;
712 infile.ignore(150, '\n'); //ignore line of *'s
713 getline(infile, line);
714 const char endmatch='*';
715 while(true) {
716 vector < string> atomwords;
717 split(line, space, atomwords);
718 if(atomwords.size() > 4) {
719 //cout << "line " << line << endl;
720 temp_atom.name=erasetail(atomwords[2]);
721 temp_atom.pos[0]=atof(atomwords[3].c_str())/bohr;
722 temp_atom.pos[1]=atof(atomwords[4].c_str())/bohr;
723 temp_atom.pos[2]=atof(atomwords[5].c_str())/bohr;
724 atoms.push_back(temp_atom);
725 getline(infile, line);
726 }
727 else break;
728 }
729 break;
730 }
731 //crystal2003 molecules
732 else if(words.size() >=4 && words[0]=="ATOM"
733 && words[1]=="X(ANGSTROM)") {
734 //cout << "found atom section " << line << endl;
735 infile.ignore(150, '\n'); //ignore line of *'s
736 getline(infile, line);
737 const char endmatch='*';
738 while(true) {
739 vector < string> atomwords;
740 split(line, space, atomwords);
741 if(atomwords.size() > 4) {
742 //cout << "line " << line << endl;
743 temp_atom.name=atomwords[3];
744 temp_atom.pos[0]=atof(atomwords[4].c_str())/bohr;
745 temp_atom.pos[1]=atof(atomwords[5].c_str())/bohr;
746 temp_atom.pos[2]=atof(atomwords[6].c_str())/bohr;
747 atoms.push_back(temp_atom);
748 getline(infile, line);
749 }
750 else break;
751 }
752 break;
753 }
754 //crystal98 with COORPRT detection
755 else if(words.size() >= 5 &&
756 words[0]=="ATOM" && words[1]=="X" ) { //1998
757
758 //cout << "found atom section " << line << endl;
759 infile.ignore(150, '\n'); //ignore line of *'s
760 getline(infile, line);
761 const char endmatch='*';
762 while(true) {
763 vector < string> atomwords;
764 split(line, space, atomwords);
765 if(atomwords[0]!="INFORMATION") {
766 //cout << "line " << line << endl;
767 temp_atom.name=atomwords[1];
768 temp_atom.pos[0]=atof(atomwords[2].c_str())/bohr;
769 temp_atom.pos[1]=atof(atomwords[3].c_str())/bohr;
770 temp_atom.pos[2]=atof(atomwords[4].c_str())/bohr;
771 atoms.push_back(temp_atom);
772 getline(infile, line);
773 }
774 else break;
775 }
776 break;
777 }
778
779
780 }
781
782 //We truncate all atom names to two letters, since
783 //crystal does that wrt the basis set & psp
784
785 for(vector<Atom>::iterator at=atoms.begin();
786 at != atoms.end(); at++) {
787 if(at->name.size() > 2)
788 at->name.erase(at->name.begin()+2, at->name.end());
789 }
790
791
792 }
793
794 //##################################################################
795
get_crystal_basis(istream & infile,vector<Gaussian_basis_set> & basis)796 void get_crystal_basis(istream & infile,
797 vector <Gaussian_basis_set> & basis) {
798 string line;
799 string space=" ";
800 vector <string> words;
801 vector < vector < string > > basis_sections;
802 vector <string> basis_labels;
803 vector <string> blank_strvec;
804 while(getline(infile, line)) {
805 words.clear();
806 split(line, space, words);
807 if(words.size() > 4 && words[0]=="ATOM" && words[1] == "X(AU)" &&
808 (words[4]=="NO." or words[4]=="N.")) {
809 infile.ignore(150,'\n');
810
811 cout << "found basis " << line << endl;
812 getline(infile, line);
813 cout << "first line " << line << endl;
814
815 const char endmatch='*';
816 string currname;
817 currname.resize(2);
818 string currtype;
819 currtype.resize(2);
820 int nbasis=-1;
821 while(search_n(line.begin(), line.end(), 5, endmatch) == line.end()
822 && line.size() != 0) {
823
824 if(line[3] != ' ') {
825 currname[0]=line[5];
826 currname[1]=line[6];
827 getline(infile, line);
828 vector <string> this_sec;
829 while(line.size() >= 4 && line[3] == ' ') {
830 this_sec.push_back(line);
831 getline(infile, line);
832 }
833 if(this_sec.size() > 0) {
834 basis_sections.push_back(this_sec);
835 basis_labels.push_back(currname);
836 }
837 }
838 }
839 break;
840 }
841 }
842
843 assert(basis_sections.size() == basis_labels.size());
844
845 //for(unsigned int i=0; i< basis_sections.size(); i++) {
846 // cout << "label " << basis_labels[i] << endl;
847 // for(unsigned int j=0; j< basis_sections[i].size(); j++) {
848 // cout << basis_sections[i][j] << endl;
849 // }
850 //}
851
852 basis.resize(basis_sections.size());
853 for(unsigned int bas=0; bas < basis_sections.size(); bas++) {
854 basis[bas].label=basis_labels[bas];
855 vector < string > indiv_types;
856 vector < vector < string> > indiv_funcs;
857
858 int fnum=-1;
859 string currtype;
860 currtype.resize(2);
861
862 //cout << "parsing the basis" << endl;
863
864 for(unsigned int j=0; j< basis_sections[bas].size(); ) {
865 if(basis_sections[bas][j][36] != ' ') {
866 currtype[0]=basis_sections[bas][j][36];
867 currtype[1]=basis_sections[bas][j][37];
868 indiv_types.push_back(currtype);
869 j++;
870 vector <string> thisfunc;
871 while(j < basis_sections[bas].size() && basis_sections[bas][j][36] == ' ') {
872 thisfunc.push_back(basis_sections[bas][j]);
873 j++;
874 }
875 indiv_funcs.push_back(thisfunc);
876 }
877 }
878
879 assert(indiv_funcs.size()==indiv_types.size());
880
881 //for(unsigned int i=0; i< indiv_funcs.size(); i++) {
882 // cout << "type " << indiv_types[i] << endl;
883 // for(unsigned int j=0; j< indiv_funcs[i].size(); j++) {
884 // cout << "line " << indiv_funcs[i][j] << endl;
885 // }
886 //}
887
888 int nfuncs=0;
889 for(vector<string>::iterator btype=indiv_types.begin();
890 btype != indiv_types.end(); btype++) {
891 if((*btype) == "SP") nfuncs+=2;
892 else nfuncs++;
893 }
894 //cout << nfuncs << " total functions " << endl;
895
896 basis[bas].exponents.resize(nfuncs);
897 basis[bas].coefficients.resize(nfuncs);
898 int currf=0;
899 for(unsigned int f=0; f< indiv_funcs.size(); f++) {
900 if(indiv_types[f]=="SP") {
901 basis[bas].types.push_back("S ");
902 basis[bas].types.push_back("P ");
903 }
904 else if(indiv_types[f]=="D ")
905 basis[bas].types.push_back("5D");
906 else if(indiv_types[f]=="F ")
907 basis[bas].types.push_back("7F_crystal");
908 else basis[bas].types.push_back(indiv_types[f]);
909 //vector <double> tmpexp;
910 //vector < double> tmpcoeff;
911 //vector < double> tmpcoeff2;//For SP functions
912 for(vector <string>::iterator line=indiv_funcs[f].begin();
913 line != indiv_funcs[f].end(); line++) {
914 string exptxt;
915 string::iterator lineb=line->begin();
916 exptxt.assign(lineb+40, lineb+50);
917 double exp=atof(exptxt.c_str());
918 basis[bas].exponents[currf].push_back(exp);
919 //cout << indiv_types[f] << ": "<< exp;
920 string coefftxt, coefftxt2;
921 if(indiv_types[f]=="S ") {
922 coefftxt.assign(lineb+50, lineb+60);
923 }
924 else if(indiv_types[f]=="P ") {
925 coefftxt.assign(lineb+60, lineb+70);
926 }
927 else if(indiv_types[f]=="D ") {
928 coefftxt.assign(lineb+70, lineb+80);
929 }
930 else if(indiv_types[f]=="SP") {
931 coefftxt.assign(lineb+50, lineb+60);
932 coefftxt2.assign(lineb+60, lineb+70);
933 }
934 else if(indiv_types[f]=="F ") {//the position of D/F/G
935 coefftxt.assign(lineb+70, lineb+80);
936 }
937 else {
938 cout << "WARNING!! Don't know type " << indiv_types[f] << endl;
939 }
940 basis[bas].coefficients[currf].push_back(atof(coefftxt.c_str()));
941 // cout << coefftxt.c_str() << endl;
942 //cout << "index: " << currf << endl;
943 //cout << "ORB:" << indiv_types[f] << " "<< atof(coefftxt.c_str())<< endl;
944 if(indiv_types[f]=="SP") {
945 basis[bas].exponents[currf+1].push_back(exp);
946 basis[bas].coefficients[currf+1].push_back(atof(coefftxt2.c_str()));
947 }
948 }
949
950 currf++;
951 if(indiv_types[f]=="SP") currf++;
952
953 }
954 }
955
956 int nbasis=basis.size();
957 //Strip whitespace from basis names and set options
958 for(int bas=0; bas < nbasis; bas++) {
959 if(basis[bas].label[1]==' ') {
960 basis[bas].label.erase(basis[bas].label.end()-1, basis[bas].label.end());
961 }
962 //This isn't necessary in Crystal2009 it seems
963 //basis[bas].options=" NORMTYPE CRYSTAL \n NORENORMALIZE \n";
964
965 //Also strip whitespace from the basis types
966 for(vector <string>::iterator i=basis[bas].types.begin();
967 i!= basis[bas].types.end(); i++) {
968 if( (*i)[1] == ' ') {
969 i->erase(i->begin()+1, i->end());
970 }
971 }
972
973 }
974
975 }
976
977
978 //##################################################################
get_crystal_pseudo(istream & infile,vector<Gaussian_pseudo_writer> & pseudo)979 void get_crystal_pseudo(istream & infile,
980 vector <Gaussian_pseudo_writer> & pseudo) {
981 string testword;
982 while (infile >> testword) {
983 //------------------------------
984 // Pseudopotentials
985 if(testword == "PSEUDOPOTENTIAL") {
986 infile >> testword;
987 if(testword == "INFORMATION") {
988 infile.ignore(125, '\n'); //rest of line
989 infile.ignore(125, '\n'); //Line of **'s
990 infile.ignore(125, '\n'); //empty line
991 string line;
992 getline(infile, line);
993 const char endmatch='*';
994 vector <double> double_blank;
995 vector <int> n_blank;
996 Gaussian_pseudo_writer pseudo_blank;
997 string temp;
998 int currpsp=-1;
999 int currl=-1;
1000 //search loop
1001 vector <string> words;
1002 string space=" ";
1003 //cout << "searching " << line << endl;
1004 while(search_n(line.begin(), line.end(), 20, endmatch) == line.end()) {
1005 words.clear();
1006 split(line,space,words);
1007
1008 //cout << "pre " << line << endl;
1009 if(words[0]=="INFORMATION" or words[0]=="NUCLEAR") break;
1010 if(line.size() > 2 && line[1]=='A' && line[2]=='T') {
1011 //new atom
1012 pseudo.push_back(pseudo_blank);
1013 currpsp++;
1014 temp.assign(line.begin()+15, line.begin()+18);
1015 pseudo[currpsp].atomnum=atoi(temp.c_str())%100;
1016 pseudo[currpsp].label=element_lookup_caps[pseudo[currpsp].atomnum];
1017 temp.assign(line.begin()+34, line.begin()+41);
1018 //pseudo[currpsp].effcharge=int(atof(words[5].c_str()));
1019 //cout << "found " << pseudo[currpsp].label << endl;
1020 cout << "Found " << pseudo[currpsp].label << " : effective charge " << atoi(temp.c_str()) << endl;
1021 pseudo[currpsp].effcharge=int(atoi(temp.c_str()));
1022 currl=-1;
1023 }
1024 else if(words.size() > 1 && words[0]!="TYPE") {
1025 //cout << "LINE: "<< line << " "<< line.size() << line[5] << endl;
1026 if(words[0][0]=='W' || words[0][0] == 'P') {
1027 pseudo[currpsp].exponents.push_back(double_blank);
1028 pseudo[currpsp].coefficients.push_back(double_blank);
1029 pseudo[currpsp].nvalue.push_back(n_blank);
1030 currl++;
1031 words.erase(words.begin());
1032 words.erase(words.begin());
1033 }
1034 if(words.size() > 2) {
1035
1036 string exps, coeffs, ns;
1037 exps.assign(line.begin()+8, line.begin()+22);
1038 coeffs.assign(line.begin()+22, line.begin()+35);
1039 ns.assign(line.begin()+35, line.begin()+39);
1040 if (coeffs.find("*", 0) == 0) {
1041 cout << "****WARNING: Fail to read coefficient (****** occurs) for gaussian function with exponent: "<< exps << endl;
1042 }
1043
1044 // pseudo[currpsp].exponents[currl].push_back(atof(words[0].c_str()));
1045 pseudo[currpsp].exponents[currl].push_back(atof(exps.c_str()));
1046 // pseudo[currpsp].coefficients[currl].push_back(atof(words[1].c_str()));
1047 pseudo[currpsp].coefficients[currl].push_back(atof(coeffs.c_str()));
1048 // pseudo[currpsp].nvalue[currl].push_back(atoi(words[2].c_str()));
1049 pseudo[currpsp].nvalue[currl].push_back(atoi(ns.c_str()));
1050 //cout << "adding " << exps << " " << coeffs << " " << ns << " done " << endl;
1051 if(words.size() > 3) {
1052 exps.assign(line.begin()+39, line.begin()+53);
1053 coeffs.assign(line.begin()+53, line.begin()+66);
1054 ns.assign(line.begin()+66, line.begin()+70);
1055 if (coeffs.find("*", 0) == 0) {
1056 cout << "****WARNING: Fail to read coefficient (****** occurs) for gaussian function with exponent: "<< exps << endl;
1057 }
1058 // pseudo[currpsp].exponents[currl].push_back(atof(words[3].c_str()));
1059 // pseudo[currpsp].coefficients[currl].push_back(atof(words[4].c_str()));
1060 // pseudo[currpsp].nvalue[currl].push_back(atoi(words[5].c_str()));
1061 pseudo[currpsp].exponents[currl].push_back(atof(exps.c_str()));
1062 pseudo[currpsp].coefficients[currl].push_back(atof(coeffs.c_str()));
1063 pseudo[currpsp].nvalue[currl].push_back(atof(ns.c_str()));
1064 //cout << "adding second " << exps << " " << coeffs << " " << ns << endl;
1065
1066 }
1067 }
1068 else {
1069 cout << "WARNING: Pseudopotential output seems corrupted!" << endl;
1070 }
1071
1072 }
1073 getline(infile, line);
1074 }
1075 //----end search loop
1076 }
1077 break;
1078 }
1079 }
1080
1081 // cout << "done " << endl;
1082 int npseud=pseudo.size();
1083 for(int ps=0; ps < npseud; ps++) {
1084 vector <double> tmp=pseudo[ps].exponents[0];
1085 pseudo[ps].exponents.erase(pseudo[ps].exponents.begin());
1086 pseudo[ps].exponents.push_back(tmp);
1087
1088
1089 tmp=pseudo[ps].coefficients[0];
1090 pseudo[ps].coefficients.erase(pseudo[ps].coefficients.begin());
1091 pseudo[ps].coefficients.push_back(tmp);
1092 vector<int> tmp2=pseudo[ps].nvalue[0];
1093 pseudo[ps].nvalue.erase(pseudo[ps].nvalue.begin());
1094 pseudo[ps].nvalue.push_back(tmp2);
1095
1096 // pseudo[ps].exponents.push_back(pseudo[ps].exponents[0]);
1097 // pseudo[ps].nvalue.push_back(pseudo[ps].nvalue[0]);
1098 // pseudo[ps].coefficients.push_back(pseudo[ps].coefficients[0]);
1099 // pseudo[ps].exponents.erase(pseudo[ps].exponents.begin());
1100 // pseudo[ps].coefficients.erase(pseudo[ps].coefficients.begin());
1101 // pseudo[ps].nvalue.erase(pseudo[ps].nvalue.begin());
1102 //cout << "size " << pseudo[ps].exponents.size() << endl;
1103 //pseudo[ps].print_pseudo(cout);
1104 }
1105
1106 }
1107
1108 //######################################################################
1109 /*!
1110 Given a starting position and an input stream, append the MO
1111 coefficients from Crystal to the vector of vectors given.
1112 moCoeff will be in form [mo][coeff].
1113 */
readMO(istream & is,long int start,vector<vector<double>> & moCoeff)1114 int readMO(istream & is, long int start, vector < vector <double> > & moCoeff) {
1115 is.clear();
1116 string line;
1117 is.seekg(start);
1118 // cout << "position " << is.tellg() << endl;
1119 int totmo=moCoeff.size();
1120 string space=" ";
1121 vector <double> emptyVector;
1122 while(getline(is, line) ) {
1123 if(line.size() > 15 && line[5]=='(' && line[15]==')') break;
1124 if(line.size() > 45 && line[35]=='(' && line[45]==')') break;
1125
1126 //check for various words that crystal can end with
1127 if(line[1]=='E' || line[3] == 'T' || line[6]=='B') break;
1128 vector <string> words;
1129 //cout << "start line " << line << endl;
1130 split(line, space, words);
1131 int nmo_this=words.size();
1132 if (words[1]=="NEWK") break;
1133 for(int i=0; i< nmo_this; i++)
1134 moCoeff.push_back(emptyVector);
1135 //cout << "nmo_this " << nmo_this << endl;
1136 //is.ignore(150, '\n'); //empty line
1137 getline(is, line); //sometimes crystal outputs a newline, and sometimes not
1138 if(line.size() < 15) getline(is, line);
1139 while(line.size() > 15 && line[1] !='E') {
1140 words.clear();
1141 split(line, space, words);
1142 //cout << "line " << line << " size " << words.size() << endl;
1143 if(words[0]=="BETA") break;
1144 if(words[1]=="NEWK") break;
1145 if(words.size()!= nmo_this+1) {
1146 cerr << "Problem reading in MOs: expected " << nmo_this << " orbitals, got " << words.size()-1 << endl;
1147 cerr << "Line:" << line << endl;
1148 exit(1);
1149 }
1150 assert(words.size() == nmo_this+1);
1151 for(int i=0; i< nmo_this; i++) {
1152 moCoeff[totmo+i].push_back(atof(words[i+1].c_str()));
1153 }
1154 getline(is, line);
1155 }
1156 totmo+=nmo_this;
1157 }
1158 return totmo;
1159 }
1160
1161 // version for complex orbitals
readMO(istream & is,long int start,vector<vector<dcomplex>> & moCoeff)1162 int readMO(istream & is, long int start,
1163 vector < vector <dcomplex> > & moCoeff) {
1164 is.clear();
1165 string line;
1166 is.seekg(start);
1167 // cout << "position " << is.tellg() << endl;
1168 int totmo=moCoeff.size();
1169 string space=" ";
1170 vector <dcomplex> emptyVector;
1171 while(getline(is, line) ) {
1172 if(line.size() > 15 && line[5]=='(' && line[15]==')') break;
1173 if(line.size() > 45 && line[35]=='(' && line[45]==')') break;
1174 //check for various words that crystal can end with
1175 if(line[1]=='E' || line[3] == 'T' || line[6]=='B') break;
1176 vector <string> words;
1177 //cout << "start line " << line << endl;
1178 split(line, space, words);
1179 if (words[1]=="NEWK") break;
1180 int nmo_this=words.size()/2;
1181 for(int i=0; i< nmo_this; i++)
1182 moCoeff.push_back(emptyVector);
1183 //cout << "nmo_this " << nmo_this << endl;
1184 //is.ignore(150, '\n'); //empty line
1185 getline(is, line); //sometimes crystal outputs a newline, and sometimes not
1186 if(line.size() < 15) getline(is, line);
1187 while(line.size() > 15 && line[1] !='E') {
1188 words.clear();
1189 split(line, space, words);
1190 //cout << "line " << line << " size " << words.size() << endl;
1191 if(words[0]=="BETA") break;
1192 if (words[1]=="NEWK") break;
1193
1194 if(words.size()!= 2*nmo_this+1) {
1195 cerr << "Problem reading in MOs: expected " << nmo_this << " orbitals, got " << words.size()-1 << endl;
1196 cerr << "Line:" << line << endl;
1197 exit(1);
1198 }
1199
1200 for(int i=0; i< nmo_this; i++) {
1201 moCoeff[totmo+i].push_back(
1202 dcomplex( atof(words[2*i+1].c_str()), atof(words[2*i+2].c_str()) )
1203 );
1204 //cout << "imag. part " << atof(words[2*i+2].c_str()) << endl;
1205 }
1206 getline(is, line);
1207 }
1208 totmo+=nmo_this;
1209 }
1210 return totmo;
1211 }
1212
1213 // ---------------------------------------------------------------------------
1214
fort10input(istream & is,string & fort10file,vector<Atom> & atoms,Slat_wf_writer & slwriter,vector<Gaussian_basis_set> & basis,vector<double> & origin,vector<vector<double>> & latvec,vector<vector<double>> & moCoeff,vector<double> & shift)1215 void fort10input(istream & is,
1216 string & fort10file,
1217 vector <Atom> & atoms,
1218 Slat_wf_writer & slwriter,
1219 vector <Gaussian_basis_set> & basis,
1220 vector <double> & origin,
1221 vector < vector <double> > & latvec,
1222 vector < vector <double> > & moCoeff,
1223 vector <double> & shift) {
1224
1225 string calctype=slwriter.calctype;
1226
1227 //reading from the output of the readcrys10.f program to get better CO's.
1228 if(fort10file != "") {
1229 cout << "Reading from the fort10 formatted file " << endl;
1230 ifstream fort10(fort10file.c_str());
1231 if(!fort10) {
1232 cout << "Couldn't open " << fort10file << endl;
1233 exit(1);
1234 }
1235
1236 int nfunctions=moCoeff[0].size();
1237
1238 int nspin=1;
1239 if(calctype=="UHF") {
1240 cout << "fort10file doesn't work for UHF yet.." << endl;
1241 exit(1);
1242 }
1243
1244 string dummy;
1245 int monum;
1246 int nmatch=0;
1247 cout << "Matching the .o orbitals with the .10 orbitals..this may "
1248 << "take some time" << endl;
1249 for(int f=0; f< nfunctions; f++) { //total number of MO's
1250 vector <double> currmo;
1251 currmo.reserve(nfunctions);
1252 fort10.ignore(120, '\n'); //line of =====
1253 fort10.ignore(120, '\n'); //
1254 fort10 >> dummy >> monum;
1255 //cout << "dummy " << dummy << endl;
1256 if(monum != f+1) {
1257 cout << "monum in fort10file doesn't match what it should be: "
1258 << "monum " << monum << " calculated function " << f+1 << endl;
1259 exit(1);
1260 }
1261 int tempfunc;
1262 double tempval;
1263 for(int f2=0; f2 < nfunctions; f2++) {
1264 fort10 >> tempfunc >> tempval;
1265 if(tempfunc != f2+1) {
1266 cout << "functions don't match in fort10file " << tempfunc
1267 << " calculated " << f2+1 << endl;
1268 exit(1);
1269 }
1270 currmo.push_back(tempval);
1271 //now loop through the mo's we've found and find the equivalent.
1272 }
1273 //cout << "currmo size " << currmo.size() << endl;
1274
1275 for(vector < vector <double > >:: iterator mo=moCoeff.begin();
1276 mo != moCoeff.end(); mo++) {
1277 double dot=0, mag_out=0, mag_fort10=0;
1278 vector <double>::iterator curr=currmo.begin();
1279 for(vector<double>::iterator moc=mo->begin();
1280 moc != mo->end(); moc++) {
1281 assert(curr != currmo.end());
1282 dot+=(*moc)*(*curr);
1283 mag_out+=(*moc)*(*moc);
1284 mag_fort10+=(*curr)*(*curr);
1285 curr++;
1286 }
1287 dot /= sqrt(mag_out*mag_fort10);
1288 if(fabs(dot-1) < .01) {
1289 //cout << "match: " << dot << endl;
1290 *mo= currmo;
1291 nmatch++;
1292 }
1293 }
1294 /*
1295 for(int mo=0; mo < totmo; mo++) {
1296 double dot=0;
1297 double mag_out=0, mag_fort10=0;
1298 for(int f=0; f< nfunctions; f++) {
1299 dot+=moCoeff[mo][f]*currmo[f];
1300 mag_out+=moCoeff[mo][f]*moCoeff[mo][f];
1301 mag_fort10+=currmo[f]*currmo[f];
1302 }
1303 dot /= sqrt(mag_out*mag_fort10);
1304
1305 //If they're the same molecular orbital, replace the output file
1306 //one with the read in one.
1307 if(fabs(dot-1) < .01) {
1308 cout << monum << " -> " << mo+1 << " dot " << dot << endl;
1309 for(int f=0; f< nfunctions; f++) {
1310 moCoeff[mo][f]=currmo[f];
1311 }
1312 }
1313
1314 }
1315 */
1316 }
1317 cout << "Matched " << nmatch << " orbitals " << endl;
1318 }
1319 }
1320
1321 //-----------------------------------------------------------------------------
1322
MO_analysis(istream & is,string & fort10file,vector<Atom> & atoms,Slat_wf_writer & slwriter,vector<Gaussian_basis_set> & basis,vector<double> & origin,vector<vector<double>> & latvec,vector<vector<double>> & moCoeff,vector<double> & shift,int totmo,string mo_filename="mo_analysis")1323 void MO_analysis(istream & is,
1324 string & fort10file,
1325 vector <Atom> & atoms,
1326 Slat_wf_writer & slwriter,
1327 vector <Gaussian_basis_set> & basis,
1328 vector <double> & origin,
1329 vector < vector <double> > & latvec,
1330 vector < vector <double> > & moCoeff,
1331 vector <double> & shift,
1332 int totmo, string mo_filename="mo_analysis") {
1333
1334 int natoms=atoms.size();
1335 ofstream an_out(mo_filename.c_str());
1336
1337 const double print_thresh=1e-3;
1338 //const double pi=3.1415926535897932385;
1339 const double pi=3.1415926535;
1340 double snorm=1./sqrt(4.*pi);
1341 double pnorm=sqrt(3.)*snorm; // sqrt(3/4/pi)
1342 vector <double> dnorm;
1343 dnorm.push_back(.5*sqrt(5./(4*pi)));//sqrt(5, 16)
1344 dnorm.push_back(sqrt(15./(4*pi)));
1345 dnorm.push_back(sqrt(15./(4*pi)));
1346 dnorm.push_back(.5*sqrt(15./(4.*pi)));
1347 dnorm.push_back(sqrt(15./(4*pi)));
1348 vector <double> fnorm;
1349 //f orbital normalizations are from <http://winter.group.shef.ac.uk/orbitron/AOs/4f/equations.html>
1350 fnorm.push_back( sqrt( 7./(16.*pi)) );
1351 fnorm.push_back( sqrt(21./(32.*pi)) );
1352 fnorm.push_back( sqrt(21./(32.*pi)) );
1353 fnorm.push_back( sqrt(105./(16.*pi)) );
1354 fnorm.push_back( sqrt(105./(4.*pi)) ); //xyz
1355 fnorm.push_back( sqrt(35./(32.*pi)) );
1356 fnorm.push_back( sqrt(35./(32.*pi)) );
1357
1358 vector <string> dnames(5);
1359 dnames[0]="z2r2";
1360 dnames[1]="xz ";
1361 dnames[2]="yz ";
1362 dnames[3]="x2y2";
1363 dnames[4]="xy ";
1364 vector <string> pnames(3);
1365 pnames[0]="x ";
1366 pnames[1]="y ";
1367 pnames[2]="z ";
1368
1369 vector <string> fnames(7);
1370 fnames[0]="F0 ";
1371 fnames[1]="Fp1 ";
1372 fnames[2]="Fm1 ";
1373 fnames[3]="Fp2 ";
1374 fnames[4]="Fxyz ";
1375 fnames[5]="Fp3 ";
1376 fnames[6]="Fm3 ";
1377 for(int mo=0; mo < totmo; mo++) {
1378 int func=0;
1379 an_out << "\n----------------\n";
1380 an_out << "MO " << mo << endl;
1381 for(int at=0; at < natoms; at++) {
1382 int bas=atoms[at].basis;
1383 int nbasis=basis[bas].types.size();
1384 for(int i=0; i< nbasis; i++) {
1385 if(basis[bas].types[i] == "S") {
1386 moCoeff[mo][func]*=snorm;
1387 if(fabs(moCoeff[mo][func]) > print_thresh) {
1388 an_out << atoms[at].name<< at << " S " << moCoeff[mo][func]
1389 << endl;
1390 }
1391 func++;
1392 }
1393 else if(basis[bas].types[i] == "P") {
1394 for(int j=0; j< 3; j++) {
1395 moCoeff[mo][func]*=pnorm;
1396 if(fabs(moCoeff[mo][func]) > print_thresh) {
1397 an_out << atoms[at].name << at << " " << "P"
1398 << pnames[j] << " " << moCoeff[mo][func]
1399 << endl;
1400 }
1401 func++;
1402 }
1403 }
1404 else if(basis[bas].types[i] == "5D") {
1405 for(int j=0; j< 5; j++) {
1406 moCoeff[mo][func]*=dnorm[j];
1407 if(fabs(moCoeff[mo][func]) > print_thresh) {
1408 an_out << atoms[at].name << at << " " << "D"
1409 << dnames[j] << " " << moCoeff[mo][func]
1410 << endl;
1411 }
1412 func++;
1413 }
1414 }
1415 else if(basis[bas].types[i] == "7F_crystal") {
1416 for(int j=0; j<7; j++) {
1417 moCoeff[mo][func]*=fnorm[j];
1418 if(fabs(moCoeff[mo][func]) > print_thresh) {
1419 an_out << atoms[at].name << at << " " << "F"
1420 << fnames[j] << " " << moCoeff[mo][func]
1421 << endl;
1422 }
1423 func++;
1424 }
1425 }
1426 else {
1427 cout << "Error: unknown basis type in read_crystal_orbital"
1428 << endl;
1429 exit(1);
1430 }
1431 }
1432 }
1433 }
1434
1435 an_out.close();
1436
1437 }
1438
MO_analysis(istream & is,vector<Atom> & atoms,Slat_wf_writer & slwriter,vector<Gaussian_basis_set> & basis,vector<double> & origin,vector<vector<double>> & latvec,vector<vector<dcomplex>> & moCoeff,vector<double> & shift,int totmo,string mo_filename="mo_analysis")1439 void MO_analysis(istream & is,
1440 vector <Atom> & atoms,
1441 Slat_wf_writer & slwriter,
1442 vector <Gaussian_basis_set> & basis,
1443 vector <double> & origin,
1444 vector < vector <double> > & latvec,
1445 vector < vector <dcomplex> > & moCoeff,
1446 vector <double> & shift,
1447 int totmo, string mo_filename="mo_analysis") {
1448
1449 int natoms=atoms.size();
1450 ofstream an_out(mo_filename.c_str());
1451
1452 const double print_thresh=1e-3;
1453 //const double pi=3.1415926535897932385;
1454 const double pi=3.1415926535;
1455 double snorm=1./sqrt(4.*pi);
1456 double pnorm=sqrt(3.)*snorm;
1457 vector <double> dnorm;
1458 dnorm.push_back(.5*sqrt(5./(4*pi)));
1459 dnorm.push_back(sqrt(15./(4*pi)));
1460 dnorm.push_back(sqrt(15./(4*pi)));
1461 dnorm.push_back(.5*sqrt(15./(4.*pi)));
1462 dnorm.push_back(sqrt(15./(4*pi)));
1463 vector <double> fnorm;
1464 //f orbital normalizations are from <http://winter.group.shef.ac.uk/orbitron/AOs/4f/equations.html>
1465 fnorm.push_back( sqrt( 7./(16.*pi)) );
1466 fnorm.push_back( sqrt(21./(32.*pi)) );
1467 fnorm.push_back( sqrt(21./(32.*pi)) );
1468 fnorm.push_back( sqrt(105./(16.*pi)) );
1469 fnorm.push_back( sqrt(105./(4.*pi)) ); //xyz
1470 fnorm.push_back( sqrt(35./(32.*pi)) );
1471 fnorm.push_back( sqrt(35./(32.*pi)) );
1472 vector <string> fnames(7);
1473 fnames[0]="F0 ";
1474 fnames[1]="Fp1 ";
1475 fnames[2]="Fm1 ";
1476 fnames[3]="Fp2 ";
1477 fnames[4]="Fxyz ";
1478 fnames[5]="Fp3 ";
1479 fnames[6]="Fm3 ";
1480
1481 vector <string> dnames(5);
1482 dnames[0]="z2r2";
1483 dnames[1]="xz ";
1484 dnames[2]="yz ";
1485 dnames[3]="x2y2";
1486 dnames[4]="xy ";
1487 vector <string> pnames(3);
1488 pnames[0]="x ";
1489 pnames[1]="y ";
1490 pnames[2]="z ";
1491 for(int mo=0; mo < totmo; mo++) {
1492 int func=0;
1493 an_out << "\n----------------\n";
1494 an_out << "MO " << mo << endl;
1495 for(int at=0; at < natoms; at++) {
1496 int bas=atoms[at].basis;
1497 int nbasis=basis[bas].types.size();
1498 for(int i=0; i< nbasis; i++) {
1499 if(basis[bas].types[i] == "S") {
1500 moCoeff[mo][func]*=snorm;
1501 if(abs(moCoeff[mo][func]) > print_thresh) {
1502 an_out << atoms[at].name<< at << " S " << moCoeff[mo][func]
1503 << endl;
1504 }
1505 func++;
1506 }
1507 else if(basis[bas].types[i] == "P") {
1508 for(int j=0; j< 3; j++) {
1509 moCoeff[mo][func]*=pnorm;
1510 if(abs(moCoeff[mo][func]) > print_thresh) {
1511 an_out << atoms[at].name << at << " " << "P"
1512 << pnames[j] << " " << moCoeff[mo][func]
1513 << endl;
1514 }
1515 func++;
1516 }
1517 }
1518 else if(basis[bas].types[i] == "5D") {
1519 for(int j=0; j< 5; j++) {
1520 moCoeff[mo][func]*=dnorm[j];
1521 if(abs(moCoeff[mo][func]) > print_thresh) {
1522 an_out << atoms[at].name << at << " " << "D"
1523 << dnames[j] << " " << moCoeff[mo][func]
1524 << endl;
1525 }
1526 func++;
1527 }
1528 }
1529 else if(basis[bas].types[i] == "7F_crystal") {
1530 for(int j=0; j< 7; j++) {
1531 moCoeff[mo][func]*=fnorm[j];
1532 if(abs(moCoeff[mo][func]) > print_thresh) {
1533 an_out << atoms[at].name << at << " " << "F"
1534 << fnames[j] << " " << moCoeff[mo][func]
1535 << endl;
1536 }
1537 func++;
1538 }
1539 }
1540 else {
1541 cout << "Error: unknown basis type in read_crystal_orbital"
1542 << endl;
1543 exit(1);
1544 }
1545 }
1546 }
1547 }
1548
1549 an_out.close();
1550
1551 }
1552
1553
1554 // ----------------------------------------------------------------------------
1555
1556 //Reads in all the molecular orbitals that were outputed by Crystal
read_crystal_orbital(istream & is,string & fort10file,vector<Atom> & atoms,Slat_wf_writer & slwriter,vector<Gaussian_basis_set> & basis,vector<double> & origin,vector<vector<double>> & latvec,vector<vector<double>> & moCoeff,vector<double> & shift)1557 void read_crystal_orbital(istream & is,
1558 string & fort10file,
1559 vector <Atom> & atoms,
1560 Slat_wf_writer & slwriter,
1561 vector <Gaussian_basis_set> & basis,
1562 vector <double> & origin,
1563 vector < vector <double> > & latvec,
1564 vector < vector <double> > & moCoeff,
1565 vector <double> & shift) {
1566
1567 assert(atoms.size() > 0);
1568 assert(basis.size() > 0);
1569 string dummy;
1570 string calctype=slwriter.calctype;
1571 //vector <int> nbasis;
1572 //vector <double > emptyVector; //to push onto moCoeff to add a mo.
1573
1574 //Find out how many functions should be there.
1575 int natoms=atoms.size();
1576 int totfunctions=0;
1577 for(int at=0; at < natoms; at++ ) {
1578 int bas=atoms[at].basis;
1579 int nfunc=basis[bas].nfunc();
1580 totfunctions+=nfunc;
1581 }
1582 //cout << "Should be " << totfunctions << " functions " << endl;
1583
1584 string space=" ";
1585 int totmo=0;
1586 vector <long int> eigen_start;
1587 vector <string> kpoints;
1588 while(is >> dummy) {
1589 //cout << "dummy " << dummy << endl;
1590 if(dummy == "FINAL") {
1591 is >> dummy;
1592 if(dummy == "EIGENVECTORS" ) {
1593 is.ignore(125, '\n'); //clear the line with FINAL EIG..
1594 string line;
1595 getline(is, line);
1596 // cout << "line " << line << endl;
1597 while(getline( is,line)) {
1598 if(line.size() > 15 && line[5]=='(' && line[15]==')') {
1599 //cout << line[5] << " " << line[15] << endl;
1600 is.ignore(150, '\n');//two blank lines
1601 is.ignore(150, '\n');
1602 long int pos=is.tellg();
1603 string line2;
1604 getline(is, line2);
1605 vector <string> words;
1606 split(line2, space, words);
1607 if(words[0]!=words[1]) {
1608 //cout << "real k-point line " << line << endl;
1609 kpoints.push_back(line);
1610 eigen_start.push_back(pos);
1611 }
1612 }
1613 }
1614 }
1615 }
1616 }
1617
1618
1619 int nkpts=eigen_start.size();
1620 if(calctype=="UHF") {
1621 assert(nkpts%2==0);
1622 nkpts/=2;
1623 }
1624 cout << "Found " << nkpts << " k-points with real eigenvectors " << endl;
1625 for(int i=0; i< nkpts; i++) {
1626 cout << i << " : " << kpoints[i] << " position " << eigen_start[i]
1627 << endl;
1628 }
1629 int kpt;
1630
1631 if(nkpts > 1) {
1632 cout << "Please choose a point[0-" << nkpts-1 << "]:";
1633 cout.flush();
1634
1635 cin >> kpt;
1636 while(kpt < 0 || kpt >= nkpts) {
1637 cout << "\nout of range..please re-enter:";
1638 cout.flush();
1639 if(!(cin >> kpt)) {
1640 cout << "\nError reading from STDIN" << endl;
1641 exit(1);
1642 }
1643 }
1644 cout << endl;
1645 }
1646 else kpt=0;
1647
1648
1649
1650 is.clear();
1651 totmo=readMO(is, eigen_start[kpt], moCoeff);
1652 if(calctype=="UHF") {
1653 slwriter.spin_dwn_start=moCoeff.size();
1654 totmo=readMO(is, eigen_start[kpt+nkpts], moCoeff);
1655 }
1656 // cout << "nmo's " << moCoeff.size() << endl;
1657
1658 slwriter.kpoint.resize(3);
1659 vector<string> kwords;
1660 kpoints[kpt].erase(kpoints[kpt].find(')'));
1661 split(kpoints[kpt], space, kwords);
1662 assert(kwords.size()>=5);
1663 double max=0;
1664 for(int i=0; i< 3; i++) {
1665 slwriter.kpoint[i]=atoi(kwords[i+2].c_str());
1666 if(slwriter.kpoint[i] > max) max=slwriter.kpoint[i];
1667 }
1668 for(int i=0; i< 3; i++)
1669 if(abs(max) > 1e-5) slwriter.kpoint[i]/=max;
1670
1671
1672 cout << "chosen k-point " << slwriter.kpoint[0] <<" "
1673 << slwriter.kpoint[1] << " " << slwriter.kpoint[2] << endl;
1674
1675
1676 //cout << "Found " << moCoeff.size() << " MO's. and "
1677 //<< moCoeff[0].size() << " functions " << endl;
1678 if(totfunctions != (int) moCoeff[0].size()) {
1679 cout << "The number of basis functions doesn't match between what was"
1680 << "read(" << moCoeff[0].size() << ") and calculated("
1681 << totfunctions << ")" << endl;
1682 exit(1);
1683 }
1684 if(totmo!=moCoeff.size()) {
1685 cout << "in make_orb, totmo is " << totmo
1686 << " and moCoeff.size() is " << moCoeff.size()
1687 << ". They should be equal." << endl;
1688 }
1689
1690
1691 fort10input(is, fort10file, atoms, slwriter, basis, origin,
1692 latvec, moCoeff, shift);
1693
1694 // analysis of band character and NORMALIZATION(!) of coefficients
1695 string mo_filename="xxmo_analysis";
1696 for(int d=0; d< 3; d++) append_number(mo_filename,slwriter.kpoint[d]);
1697 MO_analysis(is, fort10file, atoms, slwriter, basis, origin,
1698 latvec, moCoeff, shift, totmo,mo_filename);
1699
1700
1701 //if our k-point isn't zero, we need to fix any shifts
1702 int f=0;
1703 Shifter shiftobj;
1704 shiftobj.origin=origin;
1705
1706 //------------------------------------
1707 //Shift the atoms away from the edges. Should
1708 //reduce the number of centers outside the cell most of the time.
1709
1710 //vector <double> shift;
1711 //double shift_amount=0;
1712 //shift.push_back(shift_amount); shift.push_back(shift_amount);
1713 //shift.push_back(shift_amount);
1714 if(latvec.size() > 0) {
1715 //Now enforce the pbc's..
1716
1717 for(unsigned int at=0; at< atoms.size(); at++) {
1718 vector <int> nshift;
1719 int bas=atoms[at].basis;
1720 int nfunc=basis[bas].nfunc();
1721 if(shiftobj.enforcepbc(atoms[at].pos, latvec, nshift)) {
1722 //cout << "at " << at << " shifted " << nshift[0] << " " << nshift[1]
1723 // << " " << nshift[2] << endl;
1724 double kdots=0;
1725 for(int d=0; d< 3; d++)
1726 kdots+=slwriter.kpoint[d]*nshift[d];
1727 kdots=cos(pi*kdots);
1728 //cout << "kdots " << kdots << endl;
1729 for(int i=0; i< nfunc; i++) {
1730 for(int mo=0; mo < totmo; mo++) {
1731 moCoeff[mo][f]*=kdots;
1732 }
1733 f++;
1734 }
1735 }
1736 else f+=nfunc;
1737 }
1738 }
1739
1740
1741 }
1742
1743
1744 //Reads in all the molecular orbitals that were outputed by Crystal
1745 //version for complex coefficients
read_crystal_orbital(istream & is,vector<Atom> & atoms,Slat_wf_writer & slwriter,vector<Gaussian_basis_set> & basis,vector<double> & origin,vector<vector<double>> & latvec,vector<vector<dcomplex>> & moCoeff,vector<double> & shift)1746 void read_crystal_orbital(istream & is,
1747 vector <Atom> & atoms,
1748 Slat_wf_writer & slwriter,
1749 vector <Gaussian_basis_set> & basis,
1750 vector <double> & origin,
1751 vector < vector <double> > & latvec,
1752 vector < vector <dcomplex> > & moCoeff,
1753 vector <double> & shift) {
1754
1755 assert(atoms.size() > 0);
1756 assert(basis.size() > 0);
1757 string dummy;
1758 string calctype=slwriter.calctype;
1759 //vector <int> nbasis;
1760
1761 //Find out how many functions should be there.
1762 int natoms=atoms.size();
1763 int totfunctions=0;
1764 for(int at=0; at < natoms; at++ ) {
1765 int bas=atoms[at].basis;
1766 int nfunc=basis[bas].nfunc();
1767 totfunctions+=nfunc;
1768 }
1769 //cout << "Should be " << totfunctions << " functions " << endl;
1770
1771 string space=" ";
1772 int totmo=0;
1773 vector <long int> eigen_start;
1774 vector <string> kpoints;
1775 while(is >> dummy) {
1776 //cout << "dummy " << dummy << endl;
1777 if(dummy == "FINAL") {
1778 is >> dummy;
1779 if(dummy == "EIGENVECTORS" ) {
1780 is.ignore(125, '\n'); //clear the line with FINAL EIG..
1781 string line;
1782 getline(is, line);
1783 // cout << "line " << line << endl;
1784 while(getline( is,line)) {
1785 if(line.size() > 15 && line[5]=='(' && line[15]==')') {
1786 //cout << line[5] << " " << line[15] << endl;
1787 is.ignore(150, '\n');//two blank lines
1788 is.ignore(150, '\n');
1789 long int pos=is.tellg();
1790 string line2;
1791 getline(is, line2);
1792 vector <string> words;
1793 split(line2, space, words);
1794 if(words[0]==words[1]) {
1795 //cout << "complex k-point line " << line << endl;
1796 kpoints.push_back(line);
1797 eigen_start.push_back(pos);
1798 }
1799 }
1800 }
1801 }
1802 }
1803 }
1804
1805
1806 int nkpts=eigen_start.size();
1807 if(calctype=="UHF") {
1808 assert(nkpts%2==0);
1809 nkpts/=2;
1810 }
1811 cout << "Found " << nkpts << " k-points with complex eigenvectors " << endl;
1812 for(int i=0; i< nkpts; i++) {
1813 cout << i << " : " << kpoints[i] << " position " << eigen_start[i]
1814 << endl;
1815 }
1816 int kpt;
1817 if(nkpts > 1) {
1818 cout << "Please choose a point[0-" << nkpts-1 << "]:";
1819 cout.flush();
1820
1821 cin >> kpt;
1822 while(kpt < 0 || kpt >= nkpts) {
1823 cout << "\nout of range..please re-enter:";
1824 cout.flush();
1825 if(!(cin >> kpt)) {
1826 cout << "\nError reading from STDIN" << endl;
1827 exit(1);
1828 }
1829 }
1830 cout << endl;
1831 }
1832 else kpt=0;
1833
1834 is.clear();
1835 is.seekg(1);
1836 string line;
1837 int shrink_fact[3];
1838 while ( getline(is,line) ) {
1839 vector<string> words;
1840 split(line, space, words);
1841 if ( ( words.size()>5 )
1842 && ( words[0]=="SHRINK." )
1843 && ( words[1]=="FACT.(MONKH.)" ) ) {
1844 for (int is = 0; is < 3; is++ ) {
1845 shrink_fact[is]=atoi(words[2+is].c_str());
1846 }
1847 break;
1848 }
1849 }
1850 //cout << "shrinking factor " << shrink_fact << endl;
1851
1852 is.clear();
1853 totmo=readMO(is, eigen_start[kpt], moCoeff);
1854 if(calctype=="UHF") {
1855 slwriter.spin_dwn_start=moCoeff.size();
1856 totmo=readMO(is, eigen_start[kpt+nkpts], moCoeff);
1857 }
1858 cout << "nmo's " << moCoeff.size() << endl;
1859
1860 slwriter.kpoint.resize(3);
1861 vector<string> kwords;
1862 kpoints[kpt].erase(kpoints[kpt].find(')'));
1863 split(kpoints[kpt], space, kwords);
1864 assert(kwords.size()>=5);
1865 for(int i=0; i< 3; i++) {
1866 slwriter.kpoint[i]=atoi(kwords[i+2].c_str());
1867 }
1868 for(int i=0; i< 3; i++)
1869 slwriter.kpoint[i]/=shrink_fact[i]/2.0;
1870
1871
1872 cout << "chosen k-point " << slwriter.kpoint[0] <<" "
1873 << slwriter.kpoint[1] << " " << slwriter.kpoint[2] << endl;
1874
1875
1876 //cout << "Found " << moCoeff.size() << " MO's. and "
1877 //<< moCoeff[0].size() << " functions " << endl;
1878 if(totfunctions != (int) moCoeff[0].size()) {
1879 cout << "The number of basis functions doesn't match between what was"
1880 << "read(" << moCoeff[0].size() << ") and calculated("
1881 << totfunctions << ")" << endl;
1882 exit(1);
1883 }
1884 if(totmo!=moCoeff.size()) {
1885 cout << "in make_orb, totmo is " << totmo
1886 << " and moCoeff.size() is " << moCoeff.size()
1887 << ". They should be equal." << endl;
1888 }
1889
1890 // analysis of band character and NORMALIZATION(!) of coefficients
1891 MO_analysis(is, atoms, slwriter, basis, origin,
1892 latvec, moCoeff, shift, totmo);
1893
1894 //if our k-point isn't zero, we need to fix any shifts
1895 int f=0;
1896 Shifter shiftobj;
1897 shiftobj.origin=origin;
1898
1899 //------------------------------------
1900 //Shift the atoms away from the edges. Should
1901 //reduce the number of centers outside the cell most of the time.
1902
1903 //vector <double> shift;
1904 //double shift_amount=0;
1905 //shift.push_back(shift_amount); shift.push_back(shift_amount);
1906 //shift.push_back(shift_amount);
1907 if(latvec.size() > 0) {
1908 vector <double> atomshift;
1909 for(int i=0; i< 3; i++) atomshift.push_back(0);
1910 for(int i=0; i< 3; i++) {
1911 for(int j=0; j< 3; j++) {
1912 atomshift[j]+=shift[i]*latvec[i][j];
1913 }
1914 }
1915
1916 for(int i=0; i < natoms; i++) {
1917 atoms[i].pos=atoms[i].pos+atomshift;
1918 }
1919 //Now enforce the pbc's..
1920
1921 for(unsigned int at=0; at< atoms.size(); at++) {
1922 vector <int> nshift;
1923 int bas=atoms[at].basis;
1924 int nfunc=basis[bas].nfunc();
1925 if(shiftobj.enforcepbc(atoms[at].pos, latvec, nshift)) {
1926 //cout << "at " << at << " shifted " << nshift[0] << " " << nshift[1]
1927 // << " " << nshift[2] << endl;
1928 dcomplex kdots=0;
1929 for(int d=0; d< 3; d++)
1930 kdots+=slwriter.kpoint[d]*nshift[d];
1931 //kdots=cos(pi*kdots);
1932 kdots=exp(pi*kdots*dcomplex(0.0,1.0));
1933 //cout << "kdots " << kdots << endl;
1934 for(int i=0; i< nfunc; i++) {
1935 for(int mo=0; mo < totmo; mo++) {
1936 moCoeff[mo][f]*=kdots;
1937 }
1938 f++;
1939 }
1940 }
1941 else f+=nfunc;
1942 }
1943 }
1944
1945
1946 }
1947
read_crystal_orbital_all(istream & is,string & fort10file,vector<Atom> & atoms,Slat_wf_writer & slwriter,vector<Gaussian_basis_set> & basis,vector<double> & origin,vector<vector<double>> & latvec,vector<vector<vector<double>>> & moCoeff,vector<string> & kpoints,vector<long int> & eigen_start,vector<vector<double>> & kptCoord,vector<double> & shift,vector<int> shifted,vector<vector<int>> nshift)1948 void read_crystal_orbital_all(istream & is,
1949 string & fort10file,
1950 vector <Atom> &atoms,
1951 Slat_wf_writer & slwriter,
1952 vector <Gaussian_basis_set> & basis,
1953 vector <double> & origin,
1954 vector < vector <double> > & latvec,
1955 vector < vector < vector <double> > > & moCoeff,
1956 vector <string> & kpoints,
1957 vector <long int> & eigen_start,
1958 vector < vector < double > > &kptCoord,
1959 vector <double> & shift,
1960 vector <int> shifted,
1961 vector < vector <int> > nshift) {
1962
1963 assert(atoms.size() > 0);
1964 assert(basis.size() > 0);
1965 string dummy;
1966 string calctype=slwriter.calctype;
1967 //vector <int> nbasis;
1968 //vector <double > emptyVector; //to push onto moCoeff to add a mo.
1969
1970 //Find out how many functions should be there.
1971 int natoms=atoms.size();
1972 int totfunctions=0;
1973 for(int at=0; at < natoms; at++ ) {
1974 int bas=atoms[at].basis;
1975 int nfunc=basis[bas].nfunc();
1976 totfunctions+=nfunc;
1977 }
1978 //cout << "Should be " << totfunctions << " functions " << endl;
1979
1980 string space=" ";
1981 int totmo=0;
1982
1983 int nkpts=eigen_start.size();
1984 if(calctype=="UHF") {
1985 assert(nkpts%2==0);
1986 nkpts/=2;
1987 }
1988 moCoeff.resize(nkpts);
1989 kpoints.resize(nkpts);
1990 kptCoord.resize(nkpts);
1991 cout << "Found " << nkpts << " k-points with real eigenvectors " << endl;
1992 for(int i=0; i< nkpts; i++) {
1993 cout << i << " : " << kpoints[i] << endl;
1994 }
1995 int kpt;
1996
1997 is.clear();
1998 is.seekg(1);
1999 string line;
2000 int shrink_fact[3];
2001 for(int i=0; i<3; i++) shrink_fact[i]=1.0;
2002 while ( getline(is,line) ) {
2003 vector<string> words;
2004 split(line, space, words);
2005 if ( ( words.size()>5 )
2006 && ( words[0]=="SHRINK." )
2007 && ( words[1]=="FACT.(MONKH.)" ) ) {
2008 for (int is=0; is<3;is++)
2009 shrink_fact[is]=atoi(words[2+is].c_str());
2010 break;
2011 }
2012 }
2013 for (kpt=0; kpt<nkpts; kpt++) {
2014 is.clear();
2015 cout << "Reading orbital for " << kpoints[kpt] << "..." << endl;
2016 totmo=readMO(is, eigen_start[kpt], moCoeff[kpt]);
2017 if(calctype=="UHF") {
2018 slwriter.spin_dwn_start=moCoeff[kpt].size();
2019 totmo=readMO(is, eigen_start[kpt+nkpts], moCoeff[kpt]);
2020 }
2021 // cout << "nmo's " << moCoeff[kpt].size() << endl;
2022
2023 slwriter.kpoint.resize(3);
2024 for(int i=0; i<3; i++) slwriter.kpoint[i]=0.0;
2025 vector<string> kwords;
2026 kpoints[kpt].erase(kpoints[kpt].find(')'));
2027 split(kpoints[kpt], space, kwords);
2028 if(kwords.size()<5) {
2029 cout << "Not enough words in " << kpoints[kpt] << endl;
2030 exit(1);
2031 }
2032
2033 double max=0;
2034 for(int i=0; i< 3; i++) {
2035 slwriter.kpoint[i]=atoi(kwords[i+2].c_str());
2036 //if(slwriter.kpoint[i] > max) max=slwriter.kpoint[i];
2037 slwriter.kpoint[i]/=shrink_fact[i]/2.;
2038 }
2039 kptCoord[kpt].resize(3);
2040 for(int i=0; i< 3; i++) {
2041 kptCoord[kpt][i] = slwriter.kpoint[i];
2042 }
2043
2044
2045 //cout << "Found " << moCoeff[kpt].size() << " MO's. and "
2046 //<< moCoeff[kpt][0].size() << " functions " << endl;
2047 if(totfunctions != (int) moCoeff[kpt][0].size()) {
2048 cout << "The number of basis functions doesn't match between what was"
2049 << "read(" << moCoeff[kpt][0].size() << ") and calculated("
2050 << totfunctions << ")" << endl;
2051 exit(1);
2052 }
2053 if(totmo!=moCoeff[kpt].size()) {
2054 cout << "in make_orb, totmo is " << totmo
2055 << " and moCoeff[kpt].size() is " << moCoeff[kpt].size()
2056 << ". They should be equal." << endl;
2057 }
2058
2059
2060 fort10input(is, fort10file, atoms, slwriter, basis, origin,
2061 latvec, moCoeff[kpt], shift);
2062
2063 // analysis of band character and NORMALIZATION(!) of coefficients
2064 string mo_filename="mo_analysis";
2065 for(vector<double>::iterator d=slwriter.kpoint.begin(); d!=slwriter.kpoint.end(); d++)
2066 append_number(mo_filename,*d);
2067 MO_analysis(is, fort10file, atoms, slwriter, basis, origin,
2068 latvec, moCoeff[kpt], shift, totmo,mo_filename);
2069
2070
2071 //if our k-point isn't zero, we need to fix any shifts
2072 int f=0;
2073 Shifter shiftobj;
2074 shiftobj.origin=origin;
2075
2076 //------------------------------------
2077 //Shift the atoms away from the edges. Should
2078 //reduce the number of centers outside the cell most of the time.
2079
2080 //vector <double> shift;
2081 //double shift_amount=0;
2082 //shift.push_back(shift_amount); shift.push_back(shift_amount);
2083 //shift.push_back(shift_amount);
2084 if(latvec.size() > 0) {
2085
2086 //Now enforce the pbc's..
2087
2088 for(unsigned int at=0; at< atoms.size(); at++) {
2089 int bas=atoms[at].basis;
2090 int nfunc=basis[bas].nfunc();
2091 if(shifted[at]) {
2092 //cout << "at " << at << " shifted " << nshift[0] << " " << nshift[1]
2093 // << " " << nshift[2] << endl;
2094 double kdots=0;
2095 for(int d=0; d< 3; d++)
2096 kdots+=slwriter.kpoint[d]*nshift[at][d];
2097 kdots=cos(pi*kdots);
2098 //cout << "kdots " << kdots << endl;
2099 for(int i=0; i< nfunc; i++) {
2100 for(int mo=0; mo < totmo; mo++) {
2101 moCoeff[kpt][mo][f]*=kdots;
2102 }
2103 f++;
2104 }
2105 }
2106 else f+=nfunc;
2107 }
2108 }
2109 }
2110 }
2111
2112
2113
read_crystal_orbital_all(istream & is,vector<Atom> & atoms,Slat_wf_writer & slwriter,vector<Gaussian_basis_set> & basis,vector<double> & origin,vector<vector<double>> & latvec,vector<vector<vector<dcomplex>>> & moCoeff,vector<string> & kpoints,vector<long int> & eigen_start,vector<vector<double>> & kptCoord,vector<double> & shift,vector<int> shifted,vector<vector<int>> nshift)2114 void read_crystal_orbital_all(istream & is,
2115 vector <Atom> &atoms,
2116 Slat_wf_writer & slwriter,
2117 vector <Gaussian_basis_set> & basis,
2118 vector <double> & origin,
2119 vector < vector <double> > & latvec,
2120 vector <vector < vector <dcomplex> > > & moCoeff,
2121 vector <string> & kpoints,
2122 vector <long int> & eigen_start,
2123 vector < vector < double > > &kptCoord,
2124 vector <double> & shift,
2125 vector <int> shifted,
2126 vector < vector <int> > nshift) {
2127
2128 assert(atoms.size() > 0);
2129 assert(basis.size() > 0);
2130 string dummy;
2131 string calctype=slwriter.calctype;
2132 //vector <int> nbasis;
2133 //Find out how many functions should be there.
2134 int natoms=atoms.size();
2135 int totfunctions=0;
2136 for(int at=0; at < natoms; at++ ) {
2137 int bas=atoms[at].basis;
2138 int nfunc=basis[bas].nfunc();
2139 totfunctions+=nfunc;
2140 }
2141 //cout << "Should be " << totfunctions << " functions " << endl;
2142
2143 string space=" ";
2144 int totmo=0;
2145 int nkpts=eigen_start.size();
2146 if(calctype=="UHF") {
2147 assert(nkpts%2==0);
2148 nkpts/=2;
2149 }
2150 kpoints.resize(nkpts);
2151 cout << "Found " << nkpts << " k-points with complex eigenvectors " << endl;
2152 for(int i=0; i< nkpts; i++) {
2153 cout << i << " : " << kpoints[i] << endl;
2154 }
2155 int kpt;
2156 moCoeff.resize(nkpts);
2157 kptCoord.resize(nkpts);
2158
2159 is.clear();
2160 is.seekg(1);
2161 string line;
2162 int shrink_fact[3];
2163 for(int i=0; i<3; i++) shrink_fact[i]=1.0;
2164 while ( getline(is,line) ) {
2165 vector<string> words;
2166 split(line, space, words);
2167 if ( ( words.size()>5 )
2168 && ( words[0]=="SHRINK." )
2169 && ( words[1]=="FACT.(MONKH.)" ) ) {
2170 for (int is = 0; is<3; is++) {
2171 shrink_fact[is]=atoi(words[2+is].c_str());
2172 }
2173 break;
2174 }
2175 }
2176 // cout << "shrinking factor: " << shrink_fact[0] << " " << shrink_fact[1] << " " << shrink_fact[2] << endl;
2177 for (kpt = 0; kpt < nkpts; kpt++) {
2178
2179 //cout << "shrinking factor " << shrink_fact << endl;
2180
2181 is.clear();
2182 cout << "Reading orbital for " << kpoints[kpt] << "..." << endl;
2183 totmo=readMO(is, eigen_start[kpt], moCoeff[kpt]);
2184 if(calctype=="UHF") {
2185 slwriter.spin_dwn_start=moCoeff[kpt].size();
2186 totmo=readMO(is, eigen_start[kpt+nkpts], moCoeff[kpt]);
2187 }
2188 //cout << "nmo's " << moCoeff[kpt].size() << endl;
2189
2190 slwriter.kpoint.resize(3);
2191 for(int i=0;i < 3; i++)
2192 slwriter.kpoint[i]=0.0;
2193
2194 vector<string> kwords;
2195 kpoints[kpt].erase(kpoints[kpt].find(')'));
2196 split(kpoints[kpt], space, kwords);
2197 if(kwords.size()<5) {
2198 cout << "Not enough words in " << kpoints[kpt] << endl;
2199 }
2200 for(int i=0; i< 3; i++) {
2201 slwriter.kpoint[i]=atoi(kwords[i+2].c_str());
2202 }
2203 kptCoord[kpt].resize(3);
2204 for(int i=0; i< 3; i++) {
2205 slwriter.kpoint[i]/=shrink_fact[i]/2.0;
2206 kptCoord[kpt][i] = slwriter.kpoint[i];
2207 }
2208
2209
2210 //cout << "Found " << moCoeff[kpt].size() << " MO's. and "
2211 //<< moCoeff[kpt][0].size() << " functions " << endl;
2212 if(totfunctions != (int) moCoeff[kpt][0].size()) {
2213 cout << "The number of basis functions doesn't match between what was"
2214 << "read(" << moCoeff[kpt][0].size() << ") and calculated("
2215 << totfunctions << ")" << endl;
2216 exit(1);
2217 }
2218 if(totmo!=moCoeff[kpt].size()) {
2219 cout << "in make_orb, totmo is " << totmo
2220 << " and moCoeff[kpt].size() is " << moCoeff[kpt].size()
2221 << ". They should be equal." << endl;
2222 }
2223
2224 // analysis of band character and NORMALIZATION(!) of coefficients
2225 string mo_filename="mo_analysis";
2226 for(vector<double>::iterator d=slwriter.kpoint.begin(); d!=slwriter.kpoint.end(); d++)
2227 append_number(mo_filename,*d);
2228 MO_analysis(is, atoms, slwriter, basis, origin,
2229 latvec, moCoeff[kpt], shift, totmo,mo_filename);
2230
2231 //if our k-point isn't zero, we need to fix any shifts
2232 int f=0;
2233 Shifter shiftobj;
2234 shiftobj.origin=origin;
2235
2236 //------------------------------------
2237 //Shift the atoms away from the edges. Should
2238 //reduce the number of centers outside the cell most of the time.
2239
2240 //vector <double> shift;
2241 //double shift_amount=0;
2242 //shift.push_back(shift_amount); shift.push_back(shift_amount);
2243 //shift.push_back(shift_amount);
2244 if(latvec.size() > 0) {
2245 //Now enforce the pbc's..
2246
2247 for(unsigned int at=0; at< atoms.size(); at++) {
2248
2249 int bas=atoms[at].basis;
2250 int nfunc=basis[bas].nfunc();
2251 if(shifted[at]) {
2252 //cout << "at " << at << " shifted " << nshift[0] << " " << nshift[1]
2253 // << " " << nshift[2] << endl;
2254 dcomplex kdots=0;
2255 for(int d=0; d< 3; d++)
2256 kdots+=slwriter.kpoint[d]*nshift[at][d];
2257 //kdots=cos(pi*kdots);
2258 kdots=exp(pi*kdots*dcomplex(0.0,1.0));
2259 //cout << "kdots " << kdots << endl;
2260 for(int i=0; i< nfunc; i++) {
2261 for(int mo=0; mo < totmo; mo++) {
2262 moCoeff[kpt][mo][f]*=kdots;
2263 }
2264 f++;
2265 }
2266 }
2267 else f+=nfunc;
2268 }
2269 }
2270 }
2271 }
2272
read_kpt_eigenpos(istream & is,vector<string> & rkpoints,vector<long int> & reigen_start,vector<string> & ckpoints,vector<long int> & ceigen_start)2273 void read_kpt_eigenpos(istream & is,
2274 vector <string> & rkpoints,
2275 vector <long int> & reigen_start,
2276 vector <string> & ckpoints,
2277 vector <long int> & ceigen_start) {
2278
2279 string dummy;
2280 string space = " ";
2281 while(is >> dummy) {
2282 //cout << "dummy " << dummy << endl;
2283 if(dummy == "FINAL") {
2284 is >> dummy;
2285 if(dummy == "EIGENVECTORS" ) {
2286 is.ignore(125, '\n'); //clear the line with FINAL EIG..
2287 string line;
2288 getline(is, line);
2289 // cout << "line " << line << endl;
2290 while(getline( is,line)) {
2291 if(line.size() > 15 && line[5]=='(' && line[15]==')') {
2292 //cout << line[5] << " " << line[15] << endl;
2293 is.ignore(150, '\n');//two blank lines
2294 is.ignore(150, '\n');
2295 long int pos=is.tellg();
2296 string line2;
2297 getline(is, line2);
2298 vector <string> words;
2299 split(line2, space, words);
2300 if(words[0]==words[1]) {
2301 //cout << "complex k-point line " << line << endl;
2302 ckpoints.push_back(line);
2303 ceigen_start.push_back(pos);
2304 } else {
2305 rkpoints.push_back(line);
2306 reigen_start.push_back(pos);
2307 }
2308 }
2309 }
2310 }
2311 }
2312 else if(dummy == "NEWK") {
2313 is >> dummy;
2314 if(dummy == "EIGENVECTORS" ) {
2315 is.ignore(125, '\n'); //clear the line with FINAL EIG..
2316 string line;
2317 //getline(is, line);
2318 // cout << "line " << line << endl;
2319 while(getline( is,line)) {
2320 if(line.size() > 45 && line[28]=='K' && line[29]=='='&& line[35]=='(' && line[45]==')') {
2321 //cout << line[5] << " " << line[15] << endl;
2322 is.ignore(150, '\n');//one blank lines
2323 long int pos=is.tellg();
2324 string line2;
2325 getline(is, line2);
2326 vector <string> words1, words;
2327 split(line, space, words1);
2328 string kstr = words1[4] + " ( " + words1[6] + " " + words1[7] + " " + words1[8];
2329 split(line2, space, words);
2330 // cout << words[0] << " " << words[1] << endl;
2331 if(words[0]==words[1]) {
2332 //cout << "complex k-point line " << line << endl;
2333 ckpoints.push_back(kstr);
2334 ceigen_start.push_back(pos);
2335 } else {
2336 rkpoints.push_back(kstr);
2337 reigen_start.push_back(pos);
2338 }
2339 }
2340 }
2341 }
2342 }
2343 }
2344 }
2345
2346 // ===========================================================================
2347