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-2011
9 * Copyright (c) 2010-2011, 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 <cstdio>
18 #include "optimize_completeness.h"
19 #include "../basislibrary.h"
20 #include "../global.h"
21 #include "../settings.h"
22
23 #ifdef _OPENMP
24 #include <omp.h>
25 #endif
26
27 #ifdef SVNRELEASE
28 #include "../version.h"
29 #endif
30
31 Settings settings;
32
main_guarded(int argc,char ** argv)33 int main_guarded(int argc, char **argv) {
34 #ifdef _OPENMP
35 printf("ERKALE - Completeness optimization from Hel. OpenMP version, running on %i cores.\n",omp_get_max_threads());
36 #else
37 printf("ERKALE - Completeness optimization from Hel. Serial version.\n");
38 #endif
39 print_copyright();
40 print_license();
41 #ifdef SVNRELEASE
42 printf("At svn revision %s.\n\n",SVNREVISION);
43 #endif
44 print_hostname();
45
46 settings.add_int("am","angular momentum of shell to optimize for",0);
47 settings.add_int("n","moment to optimize for: 1 for maximal area, 2 for minimal rms deviation",1);
48 settings.add_double("min","lower limit of exponent range in log10",-2,true);
49 settings.add_double("max","upper limit of exponent range in log10",6,true);
50 settings.add_double("tol","Optimize and add functions until tolerance is achieved",0.0);
51 settings.add_int("nfunc","Fixed number of functions to optimize",0);
52 settings.add_int("nfull","Number of functions at each side to fully optimize",4);
53 settings.add_bool("coulomb","Use Coulomb metric? (Use only for RI basis sets)",false);
54 settings.add_double("LinDepThresh","Basis set linear dependence threshold",1e-5);
55 settings.add_string("Output","Output file to use","optimized.gbs");
56
57 if(argc!=2) {
58 settings.print();
59 printf("Usage: %s runfile\n",argv[0]);
60 return 1;
61 }
62 settings.parse(std::string(argv[1]),true);
63 settings.print();
64
65 // Get parameters
66 int am=settings.get_int("am");
67 int n=settings.get_int("n");
68 double min=settings.get_double("min");
69 double max=settings.get_double("max");
70
71 // Form optimized set of primitives
72 arma::vec exps;
73
74 // Did we get a tolerance, or a number of functions?
75 double tol=settings.get_double("tol");
76 int nfunc=settings.get_int("nfunc");
77 int nfull=settings.get_int("nfull");
78 bool coulomb=settings.get_bool("coulomb");
79 // The Coulomb metric is equivalent to the normal metric with am-1
80 if(coulomb)
81 am--;
82
83 double tau;
84
85 if(tol != 0.0 && nfunc != 0)
86 throw std::logic_error("Can't specify both wanted tolerance and number of functions!\n");
87 if(tol == 0.0 && nfunc == 0)
88 throw std::logic_error("Neither wanted tolerance or number of functions was given!\n");
89
90 if(tol!=0.0) {
91 exps=get_exponents(am,min,max,tol,n,true,nfull);
92 } else {
93 // Number of functions given.
94 exps=optimize_completeness(am,min,max,nfunc,n,true,&tau,nfull);
95 }
96
97 // Sort into ascending order
98 std::sort(exps.begin(),exps.end());
99
100 // Return the original value if Coulomb metric was used
101 if(coulomb)
102 am++;
103
104 // Create a basis set out of it. Print exponents in descending order
105 ElementBasisSet el("El");
106 for(size_t i=exps.size()-1;i<exps.size();i--) {
107 // Create shell of functions
108 FunctionShell tmp(am);
109 tmp.add_exponent(1.0,exps[i]);
110 // and add it to the basis set
111 el.add_function(tmp);
112 }
113
114 BasisSetLibrary baslib;
115 baslib.add_element(el);
116 baslib.save_gaussian94(settings.get_string("Output"));
117
118 printf("\nCompleteness-optimized basis set saved to %s.\n",settings.get_string("Output").c_str());
119
120 return 0;
121 }
122
main(int argc,char ** argv)123 int main(int argc, char **argv) {
124 try {
125 return main_guarded(argc, argv);
126 } catch (const std::exception &e) {
127 std::cerr << "error: " << e.what() << std::endl;
128 return 1;
129 }
130 }
131