1 /*
2
3 Copyright (C) 2007 Jindrich Kolorenc
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
21
22 #include <iostream>
23 #include <fstream>
24 #include <sstream>
25 #include <iomanip> // formating output stream
26 #include <cmath> // pow, abs
27 #include <limits> // to figure out precision of double
28 #include <cstdlib>
29 #include <cstring>
30
31 #include "indexx.h" // Quicksort from Numerical Recipes
32
33 //#define PRINT_VECS 1 // to print k-vectors of basis functions
34 //#define PRINT_SHELLS 1 // to print (closed) shells
35 #define TETRAHEDRON 1 // to use tetrahedron method to generate k-points
36
37 using namespace std;
38
39 void write_basis(string oname, string twist,
40 double* kx, double* ky, double* kz, int* idx,
41 int basissize, int outprec);
42 void write_slater(string oname, string twist, int* Ne);
43 void write_orb(string oname, string twist, int nmo);
44 void write_sys(string oname, string twist, int* Ne, double L,
45 double jKx, double jKy, double jKz, double weight_tw,
46 int outprec);
47 void write_jast2(string oname, double L);
48
49 extern "C" {
50 // soubroutine from TETPACK library
51 void setk01_(int* na, double* a, double* ptk, int* nptk, int*idef,
52 int* ntet, int* nkmax, int* ntmax);
53 }
54
main(int argc,char ** argv)55 int main(int argc, char ** argv) {
56 // {{{ .
57
58 const int ikmax=8; // upper limit for planewave search
59 const double pi=3.1415926535897932385;
60
61 // find out how many decimal places has double and set precision
62 // for printing k-vectors and cell size
63 typedef numeric_limits< double > dl;
64 int outprec=dl::digits10;
65 //cout << "double has " << outprec << " digits" << endl;
66
67 // base for output files is given on the command line
68 string oname;
69 bool oflag=false;
70 for(int i=1; i< argc; i++) {
71 if(!strcmp(argv[i], "-o") && argc > i+1) {
72 oname=argv[++i];
73 oflag=true;
74 }
75 else if(argv[i][0]=='-') {
76 cerr << "Unknown option: " << argv[i] << endl;
77 }
78 }
79
80 if ( !oflag ) {
81 cout << "Usage: heg2qmcC -o <output base> [-bcs]" << endl;
82 cout << endl;
83 cout << "Takes N_up, N_down and rs from standard input and ";
84 cout << "writes" << endl;
85 cout << " <output base>.i.sys" << endl;
86 cout << " <output base>.i.slater" << endl;
87 cout << " <output base>.i.basis" << endl;
88 cout << " <output base>.i.orb" << endl;
89 cout << " <output base>.centers" << endl;
90 cout << " <output base>.jast2" << endl;
91 cout << "where i indexes k-points (twists of boundary conditions) and ";
92 cout << "runs from 0 to 3;" << endl;
93 cout << "0 is the Gamma-point." << endl;
94 cout << endl;
95 return 1;
96 }
97
98 // figure out the parameters we need (read from standard input)
99 int Ne[2];
100 int Nup, Ndown;
101 cout << "Number of spin-up electrons" << endl;
102 cin >> Nup;
103 Ne[0]=Nup;
104 cout << "Number of spin-down electrons" << endl;
105 cin >> Ndown;
106 Ne[1]=Ndown;
107 int nmo; // number of molecular orbitals
108 if ( Ne[0] > Ne[1] ) { nmo=Ne[0]; } else { nmo=Ne[1]; }
109
110 cout << "rs (in atomic units, i.e., in bohrs)" << endl;
111 double rs;
112 cin >> rs;
113 #ifdef TETRAHEDRON
114 cout << "K-mesh discretization parameter (tetrahedron method)" << endl;
115 #else
116 cout << "Number of K-points in each direction" << endl;
117 #endif
118 int NK;
119 cin >> NK;
120 cout << endl;
121 double L=rs*pow(4.0/3*(Ne[0]+Ne[1])*pi,1.0/3);
122 double k0=2*pi/L;
123
124 // for (NK=2; NK<40; NK++) {
125 // cout << "---------------------------------------------------------"
126 // << "-----------------------" << endl;
127 // cout << "NK=" << NK << endl;
128
129 double Ekin=0; // Kinetic energy for constructed slater wavefunction
130 double weight=0;
131
132 int NKpt;
133
134 #ifdef TETRAHEDRON
135
136 int ntet;
137 int nkmax=(NK+1)*(NK+2)*(NK+3)/6;
138 int ntmax=NK*NK*NK;
139 int idef[ntmax*5];
140 double Kpt[nkmax*4];
141 setk01_(&NK,&L,Kpt,&NKpt,idef,&ntet,&nkmax,&ntmax);
142 cout << endl;
143
144 #else // TETRAHEDRON
145
146 NKpt=NK*NK*NK;
147 double Kpt[NKpt*4];
148 int Kcount=0;
149 for (int iKx=0; iKx<NK; iKx++) {
150 for (int iKy=0; iKy<NK; iKy++) {
151 for (int iKz=0; iKz<NK; iKz++) {
152 Kpt[Kcount*4+0]=iKx*k0/NK;
153 Kpt[Kcount*4+1]=iKy*k0/NK;
154 Kpt[Kcount*4+2]=iKz*k0/NK;
155 Kpt[Kcount*4+3]=1.0;
156 Kcount++;
157 }
158 }
159 }
160
161 #endif // TETRAHEDRON
162
163 cout << "NKpt=" << NKpt << endl;
164
165 for (int iK=0; iK<NKpt; iK++) {
166 stringstream twist;
167 twist << iK;
168 double Kx=Kpt[iK*4+0];
169 double Ky=Kpt[iK*4+1];
170 double Kz=Kpt[iK*4+2];
171 double weight_tw=Kpt[iK*4+3];
172
173 double Ekin_tw=0;
174
175 // now we construct all k-vectors in a cube and sort them with
176 // respect to their magnitude; then we occuppy as many of them as
177 // needed
178
179 int ikmax3=(2*ikmax-1)*(2*ikmax-1)*(2*ikmax-1);
180 double kx[ikmax3], ky[ikmax3], kz[ikmax3], k2[ikmax3];
181 int idx[ikmax3], shells[ikmax3];
182
183 int seq=0;
184
185 for ( int i=-ikmax+1; i<ikmax; i++) {
186 for ( int j=-ikmax+1; j<ikmax; j++) {
187 for ( int k=-ikmax+1; k<ikmax; k++) {
188 kx[seq]=Kx+i*k0;
189 ky[seq]=Ky+j*k0;
190 kz[seq]=Kz+k*k0;
191 k2[seq]=kx[seq]*kx[seq]+ky[seq]*ky[seq]+kz[seq]*kz[seq];
192 seq++;
193 }
194 }
195 }
196
197 if ( seq != ikmax3 ) {
198 cerr << "Something's wrong: got incorrect number of k-vectors." << endl;
199 cerr << seq << " instead of " << ikmax3 << endl;
200 exit(1);
201 }
202
203 // order according to length of vectors (i.e., kinetic energy)
204 // LKW note..this could be done better using objects and the STL sort,
205 //but I'll leave it for now.
206 // JK note..be my guest
207 indexx(k2,idx,ikmax3);
208
209 // find closed shells
210 int shell=0; // counter of closed shells
211 int inshell=0; // counter of states in between consecutive shells
212 for ( seq=0; seq < ikmax3; seq++ ) {
213
214 // we do not have complete shells for vectors from corners of the
215 // "test" cube, throw these vectors away (the -2 might be too conservative
216 // but there is the K-point shift we have to consider
217 if ( k2[idx[seq]] > (ikmax-2)*(ikmax-2)*k0*k0 ) break;
218
219 #ifdef PRINT_VECS
220 cout << setw(4) << seq
221 << setw(4) << kx[idx[seq]]
222 << setw(4) << ky[idx[seq]]
223 << setw(4) << kz[idx[seq]]
224 << setw(5) << k2[idx[seq]] << endl;
225 #endif
226
227 inshell++;
228
229 if ( seq+1 < ikmax3 && k2[idx[seq+1]] > k2[idx[seq]] ) {
230 if ( shell > 0 ) {
231 shells[shell]=inshell+shells[shell-1];
232 } else {
233 shells[0]=inshell;
234 }
235 #ifdef PRINT_VECS
236 cout << "----- " << setw(2) << shell+1 << " ("
237 << setw(4) << shells[shell] << ") -----" << endl;
238 #endif
239 inshell=0;
240 shell++;
241 }
242 } // for ( seq=0; seq < ikmax3; seq++ )
243 int nshells=shell;
244
245 #ifdef PRINT_SHELLS
246 // print out closed shells
247 cout << "first " << nshells << " closed shells:" << endl;
248 const int lcut=12;
249 if ( nshells < lcut ) {
250 for (int shell=0; shell < nshells; shell++) {
251 cout << setw(6) << shells[shell];
252 }
253 cout << endl;
254 } else {
255 shell=0;
256 for (int i=0; i < nshells/lcut+1; i++) {
257 int jmax=nshells-shell;
258 if ( jmax > lcut ) jmax=lcut;
259 for(int j=0; j < jmax; j++) {
260 cout << setw(6) << shells[shell];
261 shell++;
262 }
263 cout << endl;
264 }
265 }
266 cout << endl;
267 #endif // PRINT_SHELLS
268
269
270 // test whether user requested more states than we have found
271 if ( Ne[0] > shells[nshells-1] || Ne[1] > shells[nshells-1] ) {
272 cerr << "Maximum number of electrons of one kind is currently "
273 << shells[nshells-1] << endl;
274 cerr << "To go higher, change hardcoded limit of basis size (ikmax)."
275 << endl;
276 exit(1);
277 }
278
279 int basissize=nmo;
280
281 // write basis set file
282 write_basis(oname, twist.str(), kx, ky, kz, idx, basissize, outprec);
283
284 // write slater determinant
285 write_slater(oname, twist.str(), Ne);
286
287 // write orb file
288 write_orb(oname, twist.str(), nmo);
289
290 // write sys file
291 write_sys(oname, twist.str(), Ne, L, 2*Kx/k0, 2*Ky/k0, 2*Kz/k0,
292 weight_tw, outprec);
293
294 // calculate the (kinetic) energy of non-interacting particles
295 // to have an idea of this part of finite-size errors
296 for (int i=0; i<2; i++) {
297 for (int j=0; j<Ne[i]; j++) {
298 double dEkin=k2[idx[j]]/2;
299 //cout << dEkin << endl;
300 Ekin_tw+=dEkin;
301 }
302 }
303
304 // kinetic energy for a given twist
305 /*
306 cout << Kx*NK*L/pi << " " << Ky*NK*L/pi << " " << Kz*NK*L/pi
307 << " (weight "
308 << setiosflags(ios::fixed) << weight_tw << resetiosflags(ios::fixed)
309 << "), ";
310 cout << "kinetic energy: "
311 << Ekin_tw << " ("
312 << Ekin_tw/(Ne[0]+Ne[1]) << " per particle)" << endl;
313 */
314 //cout << endl << "---" << endl << endl;
315 Ekin+=Ekin_tw*weight_tw;
316 weight+=weight_tw;
317 //cout << weight_tw << " " << weight <<endl;;
318
319 Ekin_tw=0;
320
321 }
322
323 // write the file with definition of a single center
324 string centername=oname+".centers";
325 ofstream centerout(centername.c_str());
326 centerout << "1" << endl;
327 centerout << "origin 0.0 0.0 0.0" << endl;
328 centerout.close();
329
330 // write simple two-body Jastrow
331 write_jast2(oname, L);
332
333 // compare our approximate kinetic energy with exact one
334 double Ekin_ex=0;
335 for (int i=0; i<2; i++) {
336 double kF=pow(6*pi*pi*Ne[i],1.0/3)/L;
337 Ekin_ex+=L*L*L*pow(kF,5.0)/(20*pi*pi);
338 }
339 Ekin=Ekin/weight;
340 cout << setprecision(8);
341 cout << endl << "Total kinetic energy (hartree):" << endl;
342 cout << " approximate: "
343 << setw(12) << Ekin << " ("
344 << setw(12) << Ekin/(Ne[0]+Ne[1])
345 << " per particle)" << endl;
346 cout << " exact : "
347 << setw(12) << Ekin_ex << " ("
348 << setw(12) << Ekin_ex/(Ne[0]+Ne[1])
349 << " per particle)" << endl;
350 cout << " error : "
351 << setw(12) << 100*abs((Ekin-Ekin_ex)/Ekin) << "%"
352 << endl << endl;
353
354 Ekin=0.0;
355
356 // } for (NK=2; NK<40; NK++)
357
358 cout << "---" << endl << endl;
359 cout << "BTW, magnification factor in *.slater might need to be lowered for "
360 << "large" << endl << "number of particles (cca 200+ per spin channel)."
361 << endl << endl;
362
363 return 0;
364
365 // }}}
366 }
367
368
write_basis(string oname,string twist,double * kx,double * ky,double * kz,int * idx,int basissize,int outprec)369 void write_basis(string oname, string twist,
370 double* kx, double* ky, double* kz, int* idx,
371 int basissize, int outprec) {
372 // {{{ .
373
374 string basisname=oname+"."+twist+".basis";
375 ofstream basisout(basisname.c_str());
376 basisout << " CBASIS {" << endl;
377 basisout << " origin" << endl;
378 basisout << " CPLANEWAVE" << endl;
379 basisout << " GVECTOR {" << endl;
380
381 for (int i=0; i < basissize; i++) {
382 basisout << " ";
383 basisout << setprecision(outprec) << setw(outprec+6) << kx[idx[i]]
384 << setprecision(outprec) << setw(outprec+6) << ky[idx[i]]
385 << setprecision(outprec) << setw(outprec+6) << kz[idx[i]]
386 << endl;
387 }
388
389 basisout << " }" << endl;
390 basisout << " }" << endl;
391 basisout.close();
392
393 // }}}
394 }
395
write_slater(string oname,string twist,int * Ne)396 void write_slater(string oname, string twist, int* Ne) {
397 // {{{ .
398
399 int nmo_slater;
400 if ( Ne[0] > Ne[1] ) { nmo_slater=Ne[0]; } else { nmo_slater=Ne[1]; }
401
402 string slatername=oname+"."+twist+".slater";
403 ofstream slaterout(slatername.c_str());
404 slaterout << " SLATER" << endl;
405 slaterout << " CORBITALS {" << endl;
406 slaterout << " CBASFUNC_MO" << endl;
407 slaterout << " MAGNIFY 0.5" << endl;
408 slaterout << " INCLUDE " << oname;
409 slaterout << "." << twist;
410 slaterout << ".basis" << endl;
411 slaterout << " CENTERS { READ " << oname << ".centers }" << endl;
412 slaterout << " NMO " << nmo_slater << endl;
413 slaterout << " ORBFILE " << oname;
414 slaterout << "." << twist;
415 slaterout << ".orb" << endl;
416 slaterout << " }" << endl;
417 slaterout << " DETWT { 1.0 }" << endl;
418 slaterout << " STATES {" << endl;
419
420 for(int k=0; k < 2; k++) {
421 if ( k==0 ) {
422 slaterout << " # spin-ups" << endl;
423 } else {
424 slaterout << " # spin-downs" << endl;
425 }
426 if ( Ne[k] < 10 ) {
427 slaterout << " " ;
428 for (int orb=0; orb < Ne[k]; orb++) {
429 slaterout << setw(6) << orb+1;
430 }
431 slaterout << endl;
432 } else {
433 int orb=0;
434 for (int i=0; i < Ne[k]/10+1; i++) {
435 int jmax=Ne[k]-orb;
436 if ( jmax > 10 ) jmax=10;
437 if ( jmax > 0 ) slaterout << " " ;
438 for(int j=0; j < jmax; j++) {
439 orb++;
440 slaterout << setw(6) << orb;
441 }
442 if ( jmax > 0 ) slaterout << endl;
443 }
444 }
445 }
446
447 slaterout << " }" << endl;
448 slaterout.close();
449
450 // }}}
451 }
452
453
write_orb(string oname,string twist,int nmo)454 void write_orb(string oname, string twist, int nmo) {
455 // {{{ .
456 string orbname=oname+"."+twist+".orb";
457 ofstream orbout(orbname.c_str());
458 int istart=1;
459 int istop=nmo+1;
460 for( int i=istart; i < istop; i++) {
461 orbout << i << " " << 1 << endl;
462 }
463 orbout.close();
464
465 // }}}
466 }
467
write_sys(string oname,string twist,int * Ne,double L,double jKx,double jKy,double jKz,double weight_tw,int outprec)468 void write_sys(string oname, string twist, int* Ne, double L,
469 double jKx, double jKy, double jKz, double weight_tw,
470 int outprec) {
471 // {{{ .
472 string sysname=oname+"."+twist+".sys";
473 ofstream sysout(sysname.c_str());
474 sysout << "# weight of this k-point (twist): " << weight_tw << endl;
475 sysout << "SYSTEM { HEG" << endl;
476 sysout << " NSPIN { " << Ne[0] << " " << Ne[1] << " }" << endl;
477 // There are a few places in HEG system (notably calculation of e-e
478 // distances) that are optimized for orthogonal simulation cell
479 // and general lattice vectors would lead to incorrect answers.
480 sysout << " boxsize {" << endl;
481 sysout << setprecision(outprec) << setw(outprec+6) << L
482 << setprecision(outprec) << setw(outprec+6) << L
483 << setprecision(outprec) << setw(outprec+6) << L << endl;
484 sysout << " }" << endl;
485 sysout << " origin { 0.0 0.0 0.0 }" << endl;
486 // all kpoint related information should be in the basis itself
487 // <-- this is probably NOT TRUE - phase change across boundary
488 sysout << " kpoint { " << jKx << " " << jKy << " " << jKz << " }" << endl;
489 sysout << " interaction { truncCoul }" << endl;
490 sysout << "# interaction { coulomb } " << endl;
491 sysout << "}" << endl;
492 sysout.close();
493 // }}}
494 }
495
write_jast2(string oname,double L)496 void write_jast2(string oname, double L) {
497 // {{{ we could certainly use RPA to obtain better initial guess,
498 // but right now I just need something
499 string jastname=oname+".jast2";
500 ofstream jastout(jastname.c_str());
501 double cutoff=L/2.001; // will work as long as we use the simple cubic cell
502 jastout << "JASTROW2" << endl;
503 jastout << " GROUP {" << endl;
504 jastout << " OPTIMIZEBASIS" << endl;
505 jastout << " TWOBODY_SPIN {" << endl;
506 jastout << " FREEZE" << endl;
507 jastout << " LIKE_COEFFICIENTS { 0.25 0.00 }" << endl;
508 jastout << " UNLIKE_COEFFICIENTS { 0.00 0.50 }" << endl;
509 jastout << " }" << endl;
510 jastout << " EEBASIS {" << endl;
511 jastout << " EE" << endl;
512 jastout << " CUTOFF_CUSP" << endl;
513 jastout << " GAMMA 0.1" << endl; // maybe zero here
514 jastout << " CUSP 1" << endl;
515 jastout << " CUTOFF " << cutoff << endl;
516 jastout << " }" << endl;
517 jastout << " EEBASIS {" << endl;
518 jastout << " EE" << endl;
519 jastout << " CUTOFF_CUSP" << endl;
520 jastout << " GAMMA 0.5" << endl;
521 jastout << " CUSP 1" << endl;
522 jastout << " CUTOFF " << cutoff << endl;
523 jastout << " }" << endl;
524 jastout << " }" << endl;
525 jastout << " GROUP {" << endl;
526 jastout << " OPTIMIZEBASIS" << endl;
527 jastout << " TWOBODY {" << endl;
528 jastout << " COEFFICIENTS { 0.0 }" << endl;
529 jastout << " }" << endl;
530 jastout << " EEBASIS {" << endl;
531 jastout << " EE" << endl;
532 jastout << " POLYPADE" << endl;
533 jastout << " RCUT " << cutoff << endl;
534 jastout << " BETA0 -0.5" << endl;
535 jastout << " NFUNC 1" << endl;
536 jastout << " }" << endl;
537 jastout << " }" << endl;
538 jastout.close();
539 // }}}
540 }
541
542
543 // Local variables:
544 // folded-file: t
545