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 <iomanip>
22 #include <cstdlib>
23 #include <algorithm>
24
25 using namespace std;
26
split(string & text,string & separators,vector<string> & words)27 void split
28 (string & text, string & separators, vector<string> & words)
29 {
30 int n = text.length();
31 int start, stop;
32
33 start = text.find_first_not_of(separators);
34 while ((start >= 0) && (start < n)) {
35 stop = text.find_first_of(separators, start);
36 if ((stop < 0) || (stop > n)) stop = n;
37 words.push_back(text.substr(start, stop - start));
38 start = text.find_first_not_of(separators, stop+1);
39 }
40 }
41
append_number(string & str,int num)42 void append_number(string & str, int num)
43 {
44 char strbuff[40];
45 sprintf(strbuff, "%d", num);
46 str+=strbuff;
47 }
48
49
50
make_uniform(vector<double> & r,vector<double> & vals,vector<double> & ur,vector<double> & uvals,double spacing,double cutoff)51 void make_uniform(vector <double> & r,
52 vector <double> & vals,
53 vector <double> & ur,
54 vector <double> & uvals,
55 double spacing,
56 double cutoff) {
57 if(cutoff < 0) {
58 cutoff=r[r.size()-1]; //mostly for convenience..we work in terms of npts below
59 //this can be set quite high, since QWalk will find where to cut it off anyway
60 //(and fairly agressively)
61 }
62 int npts=int(cutoff/spacing);
63
64 int interval=0;
65 for(int i=0; i< npts; i++) {
66 double rad=i*spacing;
67 double val=0.0;
68 if(i==0) {
69 val=vals[0];
70 }
71 else {
72 while(r[interval+1] < rad && interval < r.size()-1) interval++;
73 double x=(rad-r[interval])/(r[interval+1]-r[interval]);
74 //cout << "rad " << rad << " x " << x << " r1 " << r[interval]
75 // << " r2 " << r[interval+1] << endl;
76 val=(1-x)*vals[interval]+x*vals[interval+1];
77 }
78 ur.push_back(rad);
79 uvals.push_back(val);
80 }
81
82 }
83
84
85 //--------------------------------------------------------
86
87
88
print_orbitals(ostream & os,std::vector<Center> & centers,std::vector<int> & nbasis,std::vector<std::vector<double>> & moCoeff)89 void print_orbitals(ostream & os,
90 std::vector <Center> & centers,
91 std::vector <int> & nbasis,
92 std::vector < std::vector < double > > & moCoeff) {
93
94 cout << "writing orb file. This could take a few seconds." << endl;
95 os.precision(12);
96 const double threshold=1e-7;
97 int ncenters=centers.size();
98 int totcoeff=0;
99 int totmo=moCoeff.size();
100
101 vector <int> totfunc_start;
102 int natoms=nbasis.size();
103 int totbasis=0;
104
105 for(int at=0; at< natoms; at++) {
106 totfunc_start.push_back(totbasis);
107 totbasis+=nbasis[at];
108 }
109
110 for(int mo=0; mo < totmo; mo++) {
111 for(int cent=0; cent < ncenters; cent++) {
112 int at=centers[cent].equiv_atom;
113 for(int ao=0; ao < nbasis[at]; ao++) {
114 //if(fabs(moCoeff[mo][totfunc_start[at]+ao]) > threshold) {
115 os << setw(10) << mo+1
116 << setw(10) << ao+1
117 << setw(10) << cent+1
118 << setw(10) << totcoeff+1 << endl;
119 totcoeff++;
120 // }
121 }
122 }
123 }
124
125 totcoeff=0;
126 os << "COEFFICIENTS\n";
127 for(int mo=0; mo < totmo; mo++) {
128 for(int cent=0; cent < ncenters; cent++) {
129 int at=centers[cent].equiv_atom;
130 for(int ao=0; ao < nbasis[at]; ao++) {
131 //if(fabs(moCoeff[mo][totfunc_start[at]+ao]) > threshold) {
132 os << moCoeff[mo][totfunc_start[at]+ao] << " ";
133 totcoeff++;
134 if(totcoeff%5==0) os << endl;
135 // }
136 }
137 }
138 }
139
140 }
141
142
print_orbitals(ostream & os,std::vector<Center> & centers,std::vector<int> & nbasis,std::vector<std::vector<dcomplex>> & moCoeff)143 void print_orbitals(ostream & os,
144 std::vector <Center> & centers,
145 std::vector <int> & nbasis,
146 std::vector < std::vector < dcomplex > > & moCoeff) {
147
148 cout << "writing orb file. This could take a few seconds." << endl;
149 os.precision(12);
150 const double threshold=1e-7;
151 int ncenters=centers.size();
152 int totcoeff=0;
153 int totmo=moCoeff.size();
154
155 vector <int> totfunc_start;
156 int natoms=nbasis.size();
157 int totbasis=0;
158
159 for(int at=0; at< natoms; at++) {
160 totfunc_start.push_back(totbasis);
161 totbasis+=nbasis[at];
162 }
163
164 for(int mo=0; mo < totmo; mo++) {
165 for(int cent=0; cent < ncenters; cent++) {
166 int at=centers[cent].equiv_atom;
167 for(int ao=0; ao < nbasis[at]; ao++) {
168 //if(fabs(moCoeff[mo][totfunc_start[at]+ao]) > threshold) {
169 os << setw(10) << mo+1
170 << setw(10) << ao+1
171 << setw(10) << cent+1
172 << setw(10) << totcoeff+1 << endl;
173 totcoeff++;
174 // }
175 }
176 }
177 }
178
179 totcoeff=0;
180 os << "COEFFICIENTS\n";
181 for(int mo=0; mo < totmo; mo++) {
182 for(int cent=0; cent < ncenters; cent++) {
183 int at=centers[cent].equiv_atom;
184 for(int ao=0; ao < nbasis[at]; ao++) {
185 //if(fabs(moCoeff[mo][totfunc_start[at]+ao]) > threshold) {
186 os << moCoeff[mo][totfunc_start[at]+ao] << " ";
187 totcoeff++;
188 if(totcoeff%5==0) os << endl;
189 // }
190 }
191 }
192 }
193
194 }
195
196
197 //###########################################################################
find_unique_atoms(const vector<Atom> & atoms,vector<string> & unique_atoms)198 void find_unique_atoms(const vector<Atom> & atoms, vector<string> & unique_atoms) {
199 unique_atoms.clear();
200 for(vector<Atom>::const_iterator at=atoms.begin(); at != atoms.end();
201 at++) {
202 if(find(unique_atoms.begin(), unique_atoms.end(),at->name)==unique_atoms.end()){
203 unique_atoms.push_back(at->name);
204 }
205 }
206 }
207
208 //######################################################################
209
210
print_centers(ostream & os,vector<Center> & centers)211 void print_centers(ostream & os, vector <Center> & centers) {
212 int ncenters=centers.size();
213 os.precision(12);
214 os << ncenters << endl;
215 for(int cent=0; cent < ncenters; cent++) {
216 os << centers[cent].name
217 << " ";
218 for(int i=0; i< 3; i++) os << centers[cent].pos[i] << " ";
219 os << endl;
220 }
221
222 }
223
224
225 //######################################################################
print_vmc_section(ostream & os,string & outputname,double eref)226 void print_vmc_section( ostream & os, string & outputname, double eref) {
227 os << "METHOD { VMC } \n";
228 }
229
230
print_opt_section(ostream & os,string & outputname,double eref)231 void print_opt_section( ostream & os, string & outputname, double eref) {
232 os << "METHOD { OPTIMIZE } \n";
233 }
234
print_dmc_section(ostream & os,string & outputname,double eref)235 void print_dmc_section( ostream & os, string & outputname, double eref) {
236 os << "METHOD { DMC timestep .01 } \n";
237 }
238
239
240
241
242 //######################################################################
243 #include "vecmath.h"
244
find_basis_cutoff(std::vector<std::vector<double>> & latvec)245 double find_basis_cutoff(std::vector <std::vector <double> > & latvec) {
246 double cutoff_divider=2.000001;
247 assert(latvec.size()==3);
248 assert(latvec[0].size()==3);
249 assert(latvec[1].size()==3);
250 assert(latvec[2].size()==3);
251
252
253 //vector <double> origin;
254 vector <double> cross01, cross12, cross02;
255 for(int i=0; i< 3; i++) {
256 //origin.push_back(0);
257 cross01.push_back(0);
258 cross12.push_back(0);
259 cross02.push_back(0);
260 }
261
262
263 //TODO-Check the cutoff determination..
264 cross01=cross(latvec[0], latvec[1]);
265 cross12=cross(latvec[1], latvec[2]);
266 cross02=cross(latvec[0], latvec[2]);
267
268 double height0=fabs(projection(latvec[0], cross12));
269 double height1=fabs(projection(latvec[1], cross02));
270 double height2=fabs(projection(latvec[2], cross01));
271 return min(min(height0, height1), height2)/cutoff_divider;
272 }
273
274 //--------------------------------------------------------------------
275
find_centers(vector<double> & origin,vector<vector<double>> & latvec,vector<Atom> & atoms,vector<Center> & centers)276 double find_centers(vector <double> & origin,
277 vector <vector <double> > & latvec,
278 vector <Atom> & atoms,
279 vector <Center> & centers ) {
280 const double cutoff_divider=1.000001;
281
282 assert(origin.size() == 3);
283 assert(latvec.size()==3);
284 assert(latvec[0].size()==3);
285 assert(latvec[1].size()==3);
286 assert(latvec[2].size()==3);
287
288
289 //vector <double> origin;
290 vector <double> cross01, cross12, cross02;
291 for(int i=0; i< 3; i++) {
292 //origin.push_back(0);
293 cross01.push_back(0);
294 cross12.push_back(0);
295 cross02.push_back(0);
296 }
297
298
299 //TODO-Check the cutoff determination..
300 cross01=cross(latvec[0], latvec[1]);
301 cross12=cross(latvec[1], latvec[2]);
302 cross02=cross(latvec[0], latvec[2]);
303
304 double height0=fabs(projection(latvec[0], cross12));
305 double height1=fabs(projection(latvec[1], cross02));
306 double height2=fabs(projection(latvec[2], cross01));
307 double cutoff=min(min(height0, height1), height2)/cutoff_divider;
308
309 double basis_cutoff=cutoff;
310
311 //cout << "heights " << height0 << " " << height1 << " " << height2 << endl;
312 //cout << "cutoff " << cutoff << endl;
313
314 vector <vector <double > > cutoff_piped;
315 for(int i=0; i< 3; i++) {
316 vector <double> tmp;
317 double norm=sqrt(dot(latvec[i], latvec[i]));
318 for(int j=0; j< 3; j++) {
319 origin[j]-=cutoff*latvec[i][j]/norm;
320 tmp.push_back(latvec[i][j]*(1+2*cutoff/norm));
321 }
322 cutoff_piped.push_back(tmp);
323 }
324
325 //the origin and cutoff_piped now define the parallelpiped that
326 //we want to find centers in.
327
328 //cout << "origin";
329 //for(int i=0; i< 3; i++) cout << origin[i] << " ";
330 //cout << endl;
331
332 //cout << "box " << endl;
333 //for(int i=0; i< 3; i++) {
334 // for(int j=0; j< 3; j++) {
335 // cout << cutoff_piped[i][j] << " ";
336 // }
337 // cout << endl;
338 //}
339
340
341 const int nsearch=1;//number of adjacent cells to search
342 int natoms=atoms.size();
343 vector <double > pos;
344 for(int i=0; i< 3; i++) pos.push_back(0);
345 vector <double> effpos; //position shifted by origin
346 for(int i=0; i< 3; i++) effpos.push_back(0);
347 vector <double> temp; //temporary vector
348 for(int i=0; i< 3; i++) temp.push_back(0);
349 Center tmpcenter;
350
351 int total=0;
352
353 int ncorner=0;
354 int nedge=0;
355 int nside=0;
356
357
358 for(int ii=-nsearch; ii < nsearch+1; ii++) {
359 for(int jj=-nsearch; jj < nsearch+1; jj++) {
360 for(int kk=-nsearch; kk < nsearch+1; kk++) {
361 for(int at=0; at < natoms; at++) {
362 for(int i=0; i< 3; i++) {
363 pos[i]=atoms[at].pos[i]
364 +(latvec[0][i])*ii
365 +(latvec[1][i])*jj
366 +(latvec[2][i])*kk;
367 }
368 effpos=pos-origin;
369 int use=0;
370 if(is_enclosed(effpos, cutoff_piped) ) {
371
372 //if all three are non-zero, we're in a corner
373 if(ii*jj*kk !=0) {
374 vector <double> corner;
375 for(int i=0; i< 3; i++) corner.push_back(0);
376
377 //find the closest corner
378 for(int i=0; i< 3; i++) {
379 corner[i]+=(ii+1)*latvec[0][i]/2;
380 corner[i]+=(jj+1)*latvec[1][i]/2;
381 corner[i]+=(kk+1)*latvec[2][i]/2;
382 }
383 if( distance_vec(pos, corner) < cutoff ) {
384 use=1;
385 ncorner++;
386 }
387 }
388
389 //if we've moved two indices at the same time,
390 //we're near the edge of the cell
391 else if(ii*jj != 0) { //a and b move
392 double aproj=projection(pos, latvec[0]);
393 double len=length_vec(latvec[0]);
394 if(aproj > len) {
395 aproj-=len;
396 }
397 double bproj=projection(pos, latvec[1]);
398 len=length_vec(latvec[1]);
399 if(bproj > len) {
400 bproj -= len;
401 }
402 if(sqrt(aproj*aproj+bproj*bproj) < cutoff) {
403 use=1;
404 nedge++;
405 }
406 else {
407 use=0;
408 }
409 }
410 else if(jj*kk != 0) { //b and c move
411 double aproj=projection(pos, latvec[1]);
412 double len=length_vec(latvec[1]);
413 if(aproj > len) {
414 aproj-=len;
415 }
416 double bproj=projection(pos, latvec[2]);
417 len=length_vec(latvec[2]);
418 if(bproj > len) {
419 bproj -= len;
420 }
421 if(sqrt(aproj*aproj+bproj*bproj) < cutoff) {
422 use=1;
423 nedge++;
424 }
425 else {
426 use=0;
427 }
428 }
429 else if(ii*kk != 0) { //a and c move
430 double aproj=projection(pos, latvec[0]);
431 double len=length_vec(latvec[0]);
432 if(aproj > len) {
433 aproj-=len;
434 }
435 double bproj=projection(pos, latvec[2]);
436 len=length_vec(latvec[2]);
437 if(bproj > len) {
438 bproj -= len;
439 }
440 if(sqrt(aproj*aproj+bproj*bproj) < cutoff) {
441 use=1;
442 nedge++;
443 }
444 else {
445 use=0;
446 }
447 }
448 else {
449 nside++;
450 use=1;
451 }
452 }
453
454 if(use ) {
455 tmpcenter.pos=pos;
456 tmpcenter.name=atoms[at].name;
457 tmpcenter.equiv_atom=at;
458 tmpcenter.basis=atoms[at].basis;
459 centers.push_back(tmpcenter);
460 total++;
461 //cout << atoms[at].name
462 // << " enclosed " << pos[0] << " " << pos[1]
463 // << " " << pos[2] << " " << at+1 << endl;
464 }
465 }
466 }
467 }
468 }
469
470 cout << "nedge " << nedge << " nside " << nside << " ncorner " << ncorner << endl;
471 cout << "total of " << total << " centers " << endl;
472 return basis_cutoff;
473 }
474
475
476 //#####################################################################
477
478
compare_mo(vector<vector<double>> & oldMOCoeff,vector<vector<double>> & moCoeff,vector<int> & compare_list)479 bool compare_mo(vector <vector < double> > & oldMOCoeff,
480 vector <vector < double> > & moCoeff,
481 vector <int> & compare_list ) {
482
483
484 int nfunctions=oldMOCoeff[0].size();
485
486 int ncompare=compare_list.size();
487 if(moCoeff.size() < ncompare) {
488 cout << "the punch file doesn't have enough molecular orbitals to compare." << endl;
489 exit(1);
490 }
491 if(oldMOCoeff.size() < ncompare) {
492 cout << "the old punch file doesn't have enough MO's to compare." << endl;
493 }
494
495 if(moCoeff[0].size() != nfunctions) {
496 cout << "Number of functions don't match between new and old mo\n";
497 exit(1);
498 }
499 vector <int> unresolved_mos;
500
501 //First check to see if the mo's are in the same place
502 //(most should be)
503 for(int i=0; i< ncompare; i++) {
504 int mo=compare_list[i];
505 double dot=0, mag_old=0, mag_new=0;
506 for(int f=0; f< nfunctions; f++) {
507 dot+=moCoeff[mo][f]*oldMOCoeff[mo][f];
508 mag_old+=oldMOCoeff[mo][f]*oldMOCoeff[mo][f];
509 mag_new+=moCoeff[mo][f]*moCoeff[mo][f];
510 }
511 dot /= sqrt(mag_old*mag_new);
512 dot =fabs(dot);
513 cout << "mo " << mo << " dot " << dot << endl;
514 if(fabs(dot-1) > .01) {
515 unresolved_mos.push_back(mo);
516 }
517 }
518
519 int nunresolved=unresolved_mos.size();
520 for(int i=0; i< nunresolved; i++) {
521 cout << "not matched: " << unresolved_mos[i] << endl;
522 }
523
524
525 bool are_same=true;
526 //See if any just swapped..
527 for(int i=0; i< nunresolved; i++) {
528 int mo1=unresolved_mos[i];
529 bool resolved_swapping=false;
530
531 for(int j=0; j< nunresolved; j++) {
532 int mo2=unresolved_mos[j];
533
534 double dot=0, mag_old=0, mag_new=0;
535 for(int f=0; f< nfunctions; f++) {
536 dot+=moCoeff[mo1][f]*oldMOCoeff[mo2][f];
537 mag_old+=oldMOCoeff[mo2][f]*oldMOCoeff[mo2][f];
538 mag_new+=moCoeff[mo1][f]*moCoeff[mo1][f];
539 }
540 dot /= sqrt(mag_old*mag_new);
541 dot=fabs(dot);
542 if(fabs(dot-1) < .01) {
543 cout << "switched orbital: mo " << mo2 << " went to " << mo1
544 << " dot product " << dot
545 << endl;
546 resolved_swapping=true;
547 }
548 }
549
550 if(!resolved_swapping) {
551 cout << "Unresolvable change in mo " << mo1 << endl;
552 are_same=false;
553 }
554
555 }
556
557 return are_same;
558 }
559
560
561
562 //----------------------------------------------------------------------
563