1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * HF/DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2019
9 * Copyright (c) 2010-2019, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16
17 #include "../basislibrary.h"
18 #include "../stringutil.h"
19 #include "../eriworker.h"
20 #include "../settings.h"
21 #include "pivoted_cholesky_basis.h"
22 #include "../completeness/completeness_profile.h"
23 #ifdef SVNRELEASE
24 #include "../version.h"
25 #endif
26
27 std::string cmds[]={"augdiffuse", "augsteep", "cholesky", "choleskybasis", "completeness", "composition", "daug", "decontract", "densityfit", "dump", "dumpdec", "genbas", "gendecbas", "merge", "norm", "orth", "overlap", "Porth", "prodset", "save", "savecfour", "savedalton", "savemolpro", "sort", "taug"};
28
29
help()30 void help() {
31 printf("Valid commands:\n");
32 for(size_t i=0;i<sizeof(cmds)/sizeof(cmds[0]);i++)
33 printf("\t%s\n",cmds[i].c_str());
34 }
35
36 Settings settings;
37
main_guarded(int argc,char ** argv)38 int main_guarded(int argc, char **argv) {
39 printf("ERKALE - Basis set tools from Hel.\n");
40 print_copyright();
41 print_license();
42 #ifdef SVNRELEASE
43 printf("At svn revision %s.\n\n",SVNREVISION);
44 #endif
45 print_hostname();
46
47 if(argc<3) {
48 printf("Usage: %s input.gbs command\n\n",argv[0]);
49 help();
50 return 0;
51 }
52
53 // Get filename
54 std::string filein(argv[1]);
55 // Load input
56 BasisSetLibrary bas;
57 bas.load_basis(filein);
58
59 // Get command
60 std::string cmd(argv[2]);
61 // and determine what to do.
62 if(stricmp(cmd,"augdiffuse")==0) {
63 // Augment basis set
64
65 if(argc!=5) {
66 printf("\nUsage: %s input.gbs %s nexp output.gbs\n",argv[0],tolower(cmd).c_str());
67 return 1;
68 }
69
70 int naug=atoi(argv[3]);
71 std::string fileout(argv[4]);
72
73 bas.augment_diffuse(naug);
74 bas.save_gaussian94(fileout);
75
76 } else if(stricmp(cmd,"augsteep")==0) {
77 // Augment basis set
78
79 if(argc!=5) {
80 printf("\nUsage: %s input.gbs %s nexp output.gbs\n",argv[0],tolower(cmd).c_str());
81 return 1;
82 }
83
84 int naug=atoi(argv[3]);
85 std::string fileout(argv[4]);
86
87 bas.augment_steep(naug);
88 bas.save_gaussian94(fileout);
89
90 } else if(stricmp(cmd,"cholesky")==0) {
91 // Form Cholesky fitting basis set
92
93 if(argc!=7) {
94 printf("\nUsage: %s input.gbs cholesky thr maxam cholthr output.gbs\n",argv[0]);
95 return 1;
96 }
97
98 double thr(atof(argv[3]));
99 int maxam(atoi(argv[4]));
100 double ovlthr(atof(argv[5]));
101 std::string outfile(argv[6]);
102
103 if(maxam>=LIBINT_MAX_AM) {
104 printf("Setting maxam = %i because limitations in used version of LIBINT.\n",LIBINT_MAX_AM-1);
105 maxam=LIBINT_MAX_AM-1;
106 }
107
108 init_libint_base();
109 BasisSetLibrary ret=bas.cholesky_set(thr,maxam,ovlthr);
110 ret.save_gaussian94(outfile);
111
112 } else if(stricmp(cmd,"choleskybasis")==0) {
113 if(argc!=7) {
114 printf("\nUsage: %s input.gbs choleskybasis system.xyz thr uselm output.gbs\n",argv[0]);
115 return 1;
116 }
117
118 std::vector<atom_t> atoms=load_xyz(argv[3],false);
119 double thr(atof(argv[4]));
120 int uselm(atof(argv[5]));
121 std::string outfile(argv[6]);
122 settings.add_scf_settings();
123 settings.set_bool("UseLM",uselm);
124
125 BasisSetLibrary ret=pivoted_cholesky_basis(atoms,bas,thr);
126 ret.save_gaussian94(outfile);
127
128 } else if(stricmp(cmd,"completeness")==0) {
129 // Print completeness profile.
130
131 if(argc!=5 && argc!=6) {
132 printf("\nUsage: %s input.gbs completeness element output.dat (coulomb)\n",argv[0]);
133 return 1;
134 }
135
136 std::string el(argv[3]);
137 std::string fileout(argv[4]);
138 bool coulomb=false;
139 if(argc==6)
140 coulomb=atoi(argv[5]);
141
142 // Get wanted element from basis
143 ElementBasisSet elbas=bas.get_element(el);
144
145 // Compute completeness profile
146 compprof_t prof=compute_completeness(elbas,-10.0,15.0,3001,coulomb);
147
148 // Print profile in output file
149 FILE *out=fopen(fileout.c_str(),"w");
150 for(size_t i=0;i<prof.lga.size();i++) {
151 // Value of scanning exponent
152 fprintf(out,"%13e",prof.lga[i]);
153 // Print completeness of shells
154 for(size_t j=0;j<prof.shells.size();j++)
155 fprintf(out,"\t%13e",prof.shells[j].Y[i]);
156 fprintf(out,"\n");
157 }
158 fclose(out);
159
160 } else if(stricmp(cmd,"composition")==0) {
161 // Determine composition of basis set.
162
163 if(argc!=3 && argc!=4) {
164 printf("\nUsage: %s input.gbs composition (El)\n",argv[0]);
165 return 1;
166 }
167
168 // Elemental basis sets
169 std::vector<ElementBasisSet> elbases;
170
171 if(argc==4)
172 elbases.push_back(bas.get_element(argv[3]));
173 else
174 elbases=bas.get_elements();
175
176 printf("\n");
177 printf("el at# [npr|nbf] [primitive|contracted(?)]\n");
178 printf("-------------------------------------------\n");
179
180 // Loop over elements
181 for(size_t iel=0;iel<elbases.size();iel++) {
182 // Get the basis set
183 ElementBasisSet elbas=elbases[iel];
184
185 // Decontracted basis
186 ElementBasisSet eldec(elbas);
187 eldec.decontract();
188
189 // Get the shells
190 std::vector<FunctionShell> sh=elbas.get_shells();
191 std::vector<FunctionShell> decsh=eldec.get_shells();
192
193 // Count the shells
194 arma::imat Nsh(max_am,2);
195 Nsh.zeros();
196 for(size_t ish=0;ish<decsh.size();ish++)
197 Nsh(decsh[ish].get_am(),0)++;
198 for(size_t ish=0;ish<sh.size();ish++)
199 Nsh(sh[ish].get_am(),1)++;
200
201 // Determine if basis set is contracted and the amount of
202 // functions
203 bool contr=false;
204 size_t nbf=0;
205 size_t nprim=0;
206 for(int am=0;am<max_am;am++) {
207 // Number of primitives
208 nprim+=Nsh(am,0)*(2*am+1);
209 // Number of contracted functions
210 nbf+=Nsh(am,1)*(2*am+1);
211 }
212 if(nbf!=nprim)
213 contr=true;
214
215 // Print composition
216 printf("%-2s %3i ",elbas.get_symbol().c_str(),(int) elbas.get_number());
217 if(contr) {
218 // Print amount of functions
219 char cmp[20];
220 sprintf(cmp,"[%i|%i]",(int) nprim,(int) nbf);
221 printf("%10s [",cmp);
222
223 // Print primitives
224 for(int am=0;am<max_am;am++)
225 if(Nsh(am,0)>0)
226 printf("%i%c",(int) Nsh(am,0),tolower(shell_types[am]));
227 // Print contractions
228 printf("|");
229 for(int am=0;am<max_am;am++)
230 if(Nsh(am,0)!=Nsh(am,1))
231 printf("%i%c",(int) Nsh(am,1),tolower(shell_types[am]));
232 printf("]\n");
233 } else {
234 printf("%10i ",(int) nbf);
235 for(int am=0;am<max_am;am++)
236 if(Nsh(am,0)>0)
237 printf("%i%c",(int) Nsh(am,0),tolower(shell_types[am]));
238 printf("\n");
239 }
240 }
241 } else if(stricmp(cmd,"daug")==0 || stricmp(cmd,"taug")==0) {
242 // Augment basis set
243
244 if(argc!=4) {
245 printf("\nUsage: %s input.gbs %s output.gbs\n",argv[0],tolower(cmd).c_str());
246 return 1;
247 }
248
249 int naug;
250 if(stricmp(cmd,"daug")==0)
251 naug=1;
252 else
253 naug=2;
254
255 std::string fileout(argv[3]);
256 bas.augment_diffuse(naug);
257 bas.save_gaussian94(fileout);
258 } else if(stricmp(cmd,"decontract")==0) {
259 // Decontract basis set.
260
261 if(argc!=4) {
262 printf("\nUsage: %s input.gbs decontract output.gbs\n",argv[0]);
263 return 1;
264 }
265
266 std::string fileout(argv[3]);
267 bas.decontract();
268 bas.save_gaussian94(fileout);
269
270 } else if(stricmp(cmd,"densityfit")==0) {
271 // Generate density fitted set
272
273 if(argc!=6) {
274 printf("\nUsage: %s input.gbs densityfit lval fsam output.gbs\n",argv[0]);
275 return 1;
276 }
277
278 int lval(atoi(argv[3]));
279 double fsam(atof(argv[4]));
280 std::string fileout(argv[5]);
281 BasisSetLibrary dfit(bas.density_fitting(lval,fsam));
282 dfit.save_gaussian94(fileout);
283
284 } else if(stricmp(cmd,"dump")==0) {
285 // Dump wanted element.
286
287 if(argc!=5 && argc!=6) {
288 printf("\nUsage: %s input.gbs dump element output.gbs (number)\n",argv[0]);
289 return 1;
290 }
291
292 std::string el(argv[3]);
293 std::string fileout(argv[4]);
294
295 int no=0;
296 if(argc==6)
297 no=atoi(argv[5]);
298
299 // Save output
300 BasisSetLibrary elbas;
301 elbas.add_element(bas.get_element(el,no));
302 elbas.save_gaussian94(fileout);
303
304 } else if(stricmp(cmd,"dumpdec")==0) {
305 // Dump wanted element in decontracted form.
306
307 if(argc!=5 && argc!=6) {
308 printf("\nUsage: %s input.gbs dumpdec element output.gbs (number)\n",argv[0]);
309 return 1;
310 }
311
312 std::string el(argv[3]);
313 std::string fileout(argv[4]);
314
315 int no=0;
316 if(argc==6)
317 no=atoi(argv[5]);
318
319 // Save output
320 BasisSetLibrary elbas;
321 bas.decontract();
322 elbas.add_element(bas.get_element(el,no));
323 elbas.save_gaussian94(fileout);
324
325 } else if(stricmp(cmd,"genbas")==0) {
326 // Generate basis set for xyz file
327
328 if(argc!=5) {
329 printf("\nUsage: %s input.gbs genbas system.xyz output.gbs\n",argv[0]);
330 return 1;
331 }
332
333 // Load atoms from xyz file
334 std::vector<atom_t> atoms=load_xyz(argv[3],false);
335 // Output file
336 std::string fileout(argv[4]);
337 // Save output
338 BasisSetLibrary elbas;
339
340 // Collect elements
341 std::vector<ElementBasisSet> els=bas.get_elements();
342 // Loop over atoms in system
343 for(size_t iat=0;iat<atoms.size();iat++) {
344 bool found=false;
345
346 // First, check if there is a special basis for the atom.
347 for(size_t iel=0;iel<els.size();iel++)
348 if(stricmp(atoms[iat].el,els[iel].get_symbol())==0 && atoms[iat].num == els[iel].get_number()) {
349 // Yes, add it.
350 elbas.add_element(els[iel]);
351 found=true;
352 break;
353 }
354
355 // Otherwise, check if a general basis is already in the basis
356 if(!found) {
357 std::vector<ElementBasisSet> added=elbas.get_elements();
358 for(size_t j=0;j<added.size();j++)
359 if(added[j].get_number()==0 && stricmp(atoms[iat].el,added[j].get_symbol())==0)
360 found=true;
361 }
362
363 // If general basis not found, add it.
364 if(!found) {
365 for(size_t iel=0;iel<els.size();iel++)
366 if(stricmp(atoms[iat].el,els[iel].get_symbol())==0 && els[iel].get_number()==0) {
367 // Yes, add it.
368 elbas.add_element(els[iel]);
369 found=true;
370 break;
371 }
372 }
373
374 if(!found) {
375 std::ostringstream oss;
376 oss << "Basis set for element " << atoms[iat].el << " does not exist in " << filein << "!\n";
377 throw std::runtime_error(oss.str());
378 }
379 }
380 elbas.sort();
381 elbas.save_gaussian94(fileout);
382
383 } else if(stricmp(cmd,"gendecbas")==0) {
384 // Generate decontracted basis set for xyz file
385
386 if(argc!=5) {
387 printf("\nUsage: %s input.gbs gendecbas system.xyz output.gbs\n",argv[0]);
388 return 1;
389 }
390
391 // Load atoms from xyz file
392 std::vector<atom_t> atoms=load_xyz(argv[3],false);
393 // Output file
394 std::string fileout(argv[4]);
395 // Save output
396 BasisSetLibrary elbas;
397
398 // Collect elements
399 std::vector<ElementBasisSet> els=bas.get_elements();
400 // Loop over atoms in system
401 for(size_t iat=0;iat<atoms.size();iat++) {
402 bool found=false;
403
404 // First, check if there is a special basis for the atom.
405 for(size_t iel=0;iel<els.size();iel++)
406 if(stricmp(atoms[iat].el,els[iel].get_symbol())==0 && atoms[iat].num == els[iel].get_number()) {
407 // Yes, add it.
408 elbas.add_element(els[iel]);
409 found=true;
410 break;
411 }
412
413 // Otherwise, check if a general basis is already in the basis
414 if(!found) {
415 std::vector<ElementBasisSet> added=elbas.get_elements();
416 for(size_t j=0;j<added.size();j++)
417 if(added[j].get_number()==0 && stricmp(atoms[iat].el,added[j].get_symbol())==0)
418 found=true;
419 }
420
421 // If general basis not found, add it.
422 if(!found) {
423 for(size_t iel=0;iel<els.size();iel++)
424 if(stricmp(atoms[iat].el,els[iel].get_symbol())==0 && els[iel].get_number()==0) {
425 // Yes, add it.
426 elbas.add_element(els[iel]);
427 found=true;
428 break;
429 }
430 }
431
432 if(!found) {
433 std::ostringstream oss;
434 oss << "Basis set for element " << atoms[iat].el << " does not exist in " << filein << "!\n";
435 throw std::runtime_error(oss.str());
436 }
437 }
438 elbas.decontract();
439 elbas.sort();
440 elbas.save_gaussian94(fileout);
441
442 } else if(stricmp(cmd,"merge")==0) {
443 // Merge functions with too big overlap
444
445 if(argc!=5) {
446 printf("\nUsage: %s input.gbs merge cutoff output.gbs\n",argv[0]);
447 return 1;
448 }
449
450 // Cutoff value
451 double cutoff=atof(argv[3]);
452 bas.merge(cutoff);
453 bas.save_gaussian94(argv[4]);
454
455 } else if(stricmp(cmd,"norm")==0) {
456 // Normalize basis
457
458 if(argc!=4) {
459 printf("\nUsage: %s input.gbs norm output.gbs\n",argv[0]);
460 return 1;
461 }
462
463 std::string fileout=argv[3];
464 bas.normalize();
465 bas.save_gaussian94(fileout);
466
467 } else if(stricmp(cmd,"orth")==0) {
468 // Orthogonalize basis
469
470 if(argc!=4) {
471 printf("\nUsage: %s input.gbs orth output.gbs\n",argv[0]);
472 return 1;
473 }
474
475 std::string fileout=argv[3];
476 bas.orthonormalize();
477 bas.save_gaussian94(fileout);
478
479 } else if(stricmp(cmd,"overlap")==0) {
480 // Primitive overlap
481
482 if(argc!=4) {
483 printf("\nUsage: %s input.gbs overlap element\n",argv[0]);
484 return 1;
485 }
486
487 // Get element basis set
488 ElementBasisSet elbas=bas.get_element(argv[3]);
489 elbas.decontract();
490
491 // Loop over angular momentum
492 for(int am=0;am<=elbas.get_max_am();am++) {
493 // Get primitives
494 arma::vec exps;
495 arma::mat contr;
496 elbas.get_primitives(exps,contr,am);
497
498 // Compute overlap matrix
499 arma::mat S=overlap(exps,exps,am);
500
501 // Print out overlap
502 printf("*** %c shell ***\n",shell_types[am]);
503 exps.t().print("Exponents");
504 printf("\n");
505
506 S.print("Overlap");
507 printf("\n");
508 }
509
510 } else if(stricmp(cmd,"Porth")==0) {
511 // P-orthogonalize basis
512
513 if(argc!=6) {
514 printf("\nUsage: %s input.gbs Porth cutoff Cortho output.gbs\n",argv[0]);
515 return 1;
516 }
517
518 double cutoff=atof(argv[3]);
519 double Cortho=atof(argv[4]);
520 std::string fileout=argv[5];
521 bas.P_orthogonalize(cutoff,Cortho);
522 bas.save_gaussian94(fileout);
523
524 } else if(stricmp(cmd,"prodset")==0) {
525 // Generate product set
526
527 if(argc!=6) {
528 printf("\nUsage: %s input.gbs prodset lval fsam output.gbs\n",argv[0]);
529 return 1;
530 }
531
532 int lval(atoi(argv[3]));
533 double fsam(atof(argv[4]));
534 std::string fileout(argv[5]);
535 BasisSetLibrary dfit(bas.product_set(lval,fsam));
536 dfit.save_gaussian94(fileout);
537
538 } else if(stricmp(cmd,"save")==0) {
539 // Save basis
540
541 if(argc!=4) {
542 printf("\nUsage: %s input.gbs save output.gbs\n",argv[0]);
543 return 1;
544 }
545
546 std::string fileout=argv[3];
547 bas.save_gaussian94(fileout);
548
549 } else if(stricmp(cmd,"savecfour")==0) {
550 // Save basis in CFOUR format
551
552 if(argc!=5) {
553 printf("\nUsage: %s input.gbs savecfour name basis.cfour\n",argv[0]);
554 return 1;
555 }
556
557 std::string fileout=argv[3];
558 std::string name=argv[4];
559 bas.save_cfour(name,fileout);
560
561 } else if(stricmp(cmd,"savedalton")==0) {
562 // Save basis in Dalton format
563
564 if(argc!=4) {
565 printf("\nUsage: %s input.gbs savedalton output.dal\n",argv[0]);
566 return 1;
567 }
568
569 std::string fileout=argv[3];
570 bas.save_dalton(fileout);
571
572 } else if(stricmp(cmd,"savemolpro")==0) {
573 // Save basis in Molpro format
574
575 if(argc!=4) {
576 printf("\nUsage: %s input.gbs savemolpro output.mol\n",argv[0]);
577 return 1;
578 }
579
580 std::string fileout=argv[3];
581 bas.save_molpro(fileout);
582
583 } else if(stricmp(cmd,"sort")==0) {
584 // Sort basis set
585
586 if(argc!=4) {
587 printf("\nUsage: %s input.gbs sort output.gbs\n",argv[0]);
588 return 1;
589 }
590
591 std::string fileout=argv[3];
592 bas.sort();
593 bas.save_gaussian94(fileout);
594 } else {
595 printf("\nInvalid command.\n");
596
597 help();
598 }
599
600 return 0;
601 }
602
main(int argc,char ** argv)603 int main(int argc, char **argv) {
604 try {
605 return main_guarded(argc, argv);
606 } catch (const std::exception &e) {
607 std::cerr << "error: " << e.what() << std::endl;
608 return 1;
609 }
610 }
611