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 "basislibrary.h"
18 #include "basis.h"
19 #include "checkpoint.h"
20 #include "dftgrid.h"
21 #include "elements.h"
22 #include "find_molecules.h"
23 #include "guess.h"
24 #include "linalg.h"
25 #include "mathf.h"
26 #include "xyzutils.h"
27 #include "properties.h"
28 #include "sap.h"
29 #include "scf.h"
30 #include "settings.h"
31 #include "stringutil.h"
32 #include "timer.h"
33
34 // Needed for libint init
35 #include "eriworker.h"
36
37 #include <armadillo>
38 #include <cstdio>
39 #include <cstdlib>
40 #include <cfloat>
41
42 #ifdef _OPENMP
43 #include <omp.h>
44 #endif
45
46 #ifdef SVNRELEASE
47 #include "version.h"
48 #endif
49
print_header()50 void print_header() {
51 #ifdef _OPENMP
52 printf("ERKALE - HF/DFT from Hel, OpenMP version, running on %i cores.\n",omp_get_max_threads());
53 #else
54 printf("ERKALE - HF/DFT from Hel, serial version.\n");
55 #endif
56 print_copyright();
57 print_license();
58 #ifdef SVNRELEASE
59 printf("At svn revision %s.\n\n",SVNREVISION);
60 #endif
61 print_hostname();
62 }
63
diag(arma::vec & E,arma::mat & C,const arma::mat & H,const arma::mat & Sinvh)64 void diag(arma::vec & E, arma::mat & C, const arma::mat & H, const arma::mat & Sinvh) {
65 arma::eig_sym(E,C,Sinvh.t()*H*Sinvh);
66 C=Sinvh*C;
67
68 //E.print("Eigenvalues");
69 }
70
71 Settings settings;
72
main_guarded(int argc,char ** argv)73 int main_guarded(int argc, char **argv) {
74 print_header();
75
76 if(argc!=2) {
77 printf("Usage: $ %s runfile\n",argv[0]);
78 return 0;
79 }
80
81 // Initialize libint
82 init_libint_base();
83
84 Timer t;
85 t.print_time();
86
87 // Parse settings
88 settings.add_scf_settings();
89 settings.add_string("LoadChk","File to load old results from","");
90 settings.add_string("SaveChk","File to save guess to","");
91 settings.add_double("LinDepThr","Linear dependency threshold",1e-5);
92 settings.add_bool("FON","Fermi occupation numbers",false);
93 settings.add_string("FONscan","Fermi occupation scan","");
94 settings.add_double("T","Temperature in K",1000);
95 settings.add_double("K","GWH approximation scaling factor",1.75);
96 settings.add_int("nrad","Number of radial shells for SAP",99);
97 settings.add_int("lmax","Angular rule for SAP (defaults to l=41 i.e. 590 points)",41);
98 settings.parse(std::string(argv[1]),true);
99 settings.print();
100
101 // SAP quadrature rule
102 int nrad(settings.get_int("nrad"));
103 int lmax(settings.get_int("lmax"));
104 // GWH constant
105 double K(settings.get_double("K"));
106
107 // Basis set
108 BasisSet basis;
109 // Get checkpoint file
110 std::string chkf(settings.get_string("LoadChk"));
111 std::string savef(settings.get_string("SaveChk"));
112 Checkpoint chk(chkf,false);
113 chk.read(basis);
114
115 int restr;
116 chk.read("Restricted",restr);
117
118 // Read in density matrix
119 arma::mat Pa, Pb;
120 if(restr) {
121 chk.read("P",Pa);
122 Pa/=2.0;
123 Pb=Pa;
124 } else {
125 chk.read("Pa",Pa);
126 chk.read("Pb",Pb);
127 }
128
129 // Number of electrons
130 int Nela;
131 chk.read("Nel-a",Nela);
132 int Nelb;
133 chk.read("Nel-b",Nelb);
134
135 // Overlap matrix
136 arma::mat S(basis.overlap());
137 // and its half-inverse
138 double linthr(settings.get_double("LinDepThr"));
139 arma::mat Sinvh(CanonicalOrth(S,linthr));
140
141 // Guess energy and orbitals
142 arma::vec E;
143 arma::mat C;
144 arma::mat Pag, Pbg;
145 // Form core Hamiltonian
146 arma::mat V(basis.nuclear());
147 arma::mat T(basis.kinetic());
148 arma::mat Hcore(T+V);
149
150 // Form density matrix from guess orbitals?
151 bool formP=true;
152
153 // Compute guess density matrix
154 std::string guess(settings.get_string("Guess"));
155 bool FON(settings.get_bool("FON"));
156 std::string FONscan(settings.get_string("FONscan"));
157 double temp(settings.get_double("T"));
158 if(stricmp(guess,"core")==0) {
159 // Diagonalize it
160 diag(E,C,Hcore,Sinvh);
161
162 } else if(stricmp(guess,"gwh")==0) {
163
164 // Form GWH matrix
165 arma::mat Hgwh(Hcore);
166 for(size_t i=0;i<Hcore.n_rows;i++) {
167 Hgwh(i,i)=Hcore(i,i);
168 for(size_t j=0;j<Hcore.n_cols;j++) {
169 Hgwh(i,j)=0.5*K*S(i,j)*(Hcore(i,i)+Hcore(j,j));
170 Hgwh(j,i)=Hgwh(i,j);
171 }
172 }
173 // and diagonalize it
174 diag(E,C,Hgwh,Sinvh);
175
176 } else if(stricmp(guess,"sap")==0) {
177 DFTGrid grid(&basis);
178 bool grad=false;
179 bool tau=false;
180 bool lapl=false;
181 bool strict=false;
182 bool nl=false;
183 grid.construct(nrad,lmax,grad,tau,lapl,strict,nl);
184
185 // Get SAP potential
186 arma::mat Vsap(grid.eval_SAP());
187 // Approximate Hamiltonian is
188 diag(E,C,Hcore+Vsap,Sinvh);
189
190 } else if(stricmp(guess,"sad")==0 || stricmp(guess,"no")==0) {
191 // Get SAD guess
192 Pag=sad_guess(basis)/2.0;
193 Pbg=Pag;
194
195 if(stricmp(guess,"no")==0) {
196 // Build Fock operator from SAD density matrix. Get natural orbitals
197 arma::vec occ;
198 form_NOs(Pag,S,C,occ);
199 } else {
200 if((Nela+Nelb)!=basis.Ztot()) {
201 printf("Warning: SAD density doesn't integrate to wanted number of electrons.\n");
202 }
203 if(Nela!=Nelb) {
204 printf("Warning: SAD density doesn't correspond to wanted spin state.\n");
205 }
206 formP=false;
207 }
208
209 } else if(stricmp(guess,"gsap")==0) {
210 // Get SAP guess
211 diag(E,C,Hcore+sap_guess(basis),Sinvh);
212
213 } else if(stricmp(guess,"huckel")==0) {
214 // Get Huckel guess
215 diag(E,C,huckel_guess(basis,K),Sinvh);
216
217 } else if(stricmp(guess,"minsap")==0) {
218 DFTGrid grid(&basis);
219 bool grad=false;
220 bool tau=false;
221 bool lapl=false;
222 bool strict=false;
223 bool nl=false;
224 grid.construct(nrad,lmax,grad,tau,lapl,strict,nl);
225
226 // Get SAP potential
227 arma::mat Vsap(grid.eval_SAP());
228 // SAP Hamiltonian
229 arma::mat Hsap(Hcore+Vsap);
230
231 // Get minimal basis
232 arma::mat Cmin(minimal_basis_projection(basis));
233
234 // Convert to minimal basis
235 Hsap=Cmin.t()*Hsap*Cmin;
236 // Scale off-diagonal elements
237 for(size_t j=0;j<Hsap.n_cols;j++)
238 for(size_t i=0;i<j;i++) {
239 Hsap(i,j) *= 0.5*K;
240 Hsap(j,i) *= 0.5*K;
241 }
242 // Go back to AO basis
243 Hsap=Cmin*Hsap*Cmin.t();
244
245 // Diagonalize
246 diag(E,C,Hsap,Sinvh);
247
248 } else {
249 throw std::logic_error("Unsupported guess!\n");
250 }
251
252 if(formP) {
253 arma::vec occa(C.n_cols), occb(C.n_cols);
254 if(FON) {
255 if(!formP)
256 throw std::logic_error("Can't do Fermi occupations without orbitals!\n");
257 if(Nela)
258 printf("Alpha gap is %e i.e. % 6.3f eV\n",(E(Nela)-E(Nela-1)),27.2114*(E(Nela)-E(Nela-1)));
259 if(Nelb && Nela!=Nelb)
260 printf("Beta gap is %e i.e. % 6.3f eV\n",(E(Nelb)-E(Nelb-1)),27.2114*(E(Nelb)-E(Nelb-1)));
261 occa=FermiON(E,Nela,temp);
262 occb=FermiON(E,Nelb,temp);
263
264 if(Nela)
265 printf("Alpha NOON gap is %e\n",(occa(Nela-1)-occa(Nela)));
266 if(Nelb && Nela!=Nelb)
267 printf("Beta NOON gap is %e\n",(occb(Nelb-1)-occb(Nelb)));
268 } else {
269 occa.zeros();
270 occb.zeros();
271 occa.subvec(0,Nela-1).ones();
272 if(Nelb)
273 occb.subvec(0,Nelb-1).ones();
274 }
275 //E.print("Orbital energies");
276 //occa.print("Alpha occupation");
277 //occb.print("Beta occupation");
278 Pag=C*arma::diagmat(occa)*C.t();
279 if(occb.n_elem)
280 Pbg=C*arma::diagmat(occb)*C.t();
281 }
282
283 // Calculate the projection
284 double aproj(arma::trace(Pa*S*Pag*S));
285 double bproj(arma::trace(Pb*S*Pbg*S));
286 if(std::abs(aproj-bproj)>=sqrt(DBL_EPSILON)) {
287 printf("Alpha projection of guess onto SCF density is %e i.e. %5.2f %%\n",aproj,aproj/Nela*100.0);
288 printf("Beta projection of guess onto SCF density is %e i.e. %5.2f %%\n",bproj,bproj/Nelb*100.0);
289 }
290 printf("Projection of guess onto SCF density is %e i.e. %5.2f %%\n",aproj+bproj,(aproj+bproj)/(Nela+Nelb)*100.0);
291
292 // Save result?
293 if(savef.size()) {
294 std::string cmd("cp " + chkf + " " + savef);
295 (void) system(cmd.c_str());
296
297 Checkpoint save(savef,true,false);
298 save.write("P",Pag+Pbg);
299 if(restr) {
300 save.write("C",C);
301 } else {
302 save.write("Ca",C);
303 save.write("Cb",C);
304 save.write("Pa",Pag);
305 save.write("Pb",Pbg);
306 }
307 }
308
309 if(FONscan.size()) {
310 if(!formP)
311 throw std::logic_error("Can't do Fermi temperature scan without orbitals!\n");
312
313 // Scan FONs
314 arma::vec Ts(arma::linspace<arma::vec>(0,20000,101));
315
316 // Data
317 arma::mat fon(Ts.n_elem,2);
318 fon.col(0)=Ts;
319 for(size_t i=0;i<Ts.n_elem;i++) {
320 arma::vec occa(FermiON(E,Nela,Ts(i)));
321 arma::vec occb(FermiON(E,Nelb,Ts(i)));
322 Pag=C*arma::diagmat(occa)*C.t();
323 Pbg=C*arma::diagmat(occb)*C.t();
324 aproj=arma::trace(Pa*S*Pag*S);
325 bproj=arma::trace(Pb*S*Pbg*S);
326 fon(i,1)=(aproj+bproj)/(Nela+Nelb)*100.0;
327 }
328 fon.save(FONscan,arma::raw_ascii);
329 printf("FON scan saved in %s\n",FONscan.c_str());
330 }
331
332 printf("\nRunning program took %s.\n",t.elapsed().c_str());
333
334 return 0;
335 }
336
main(int argc,char ** argv)337 int main(int argc, char **argv) {
338 try {
339 return main_guarded(argc, argv);
340 } catch (const std::exception &e) {
341 std::cerr << "error: " << e.what() << std::endl;
342 return 1;
343 }
344 }
345