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-2015
9 * Copyright (c) 2010-2015, 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 "../checkpoint.h"
18 #include "../mathf.h"
19 #include "../emd/emd_gto.h"
20 #include "../stringutil.h"
21 #include "../settings.h"
22
23 /// Absolute tolerance for normalization of basis functions
24 const double normtol=1e-10;
25
26 /// Absolute tolerance in total energy
27 const double Etol=1e-5;
28 /// Absolute tolerance for density matrix difference
29 const double dPtol=1e-6;
30 /// Absolute tolerance for orbital matrix difference
31 const double dCtol=1e-5;
32 /// Relative tolerance in orbital energies
33 const double dEtol=1e-5;
34
35 /// Tolerance for the number of electrons estimate
36 const double P0TOL=50.0;
37 /// Tolerance for the EMD error estimate
38 const double P2TOL=50.0;
39 /// Orbitals are deemed degenerate if the energy difference is less than
40 const double odegthr=1e-4;
41
42 // Check proper normalization of basis
check_norm(const BasisSet & bas)43 void check_norm(const BasisSet & bas) {
44 size_t Nbf=bas.get_Nbf();
45 arma::mat S=bas.overlap();
46
47 for(size_t i=0;i<Nbf;i++)
48 if(fabs(S(i,i)-1.0)>=normtol) {
49 std::ostringstream oss;
50 ERROR_INFO();
51 fflush(stdout);
52 oss << "Function " << i+1 << " is not normalized: norm is " << S(i,i) << "!.\n";
53 throw std::runtime_error(oss.str());
54 }
55 }
56
fix_signs(arma::mat & C,const arma::mat & Cr,const arma::mat & S)57 void fix_signs(arma::mat & C, const arma::mat & Cr, const arma::mat & S) {
58 for(size_t io=0;io<Cr.n_cols;io++) {
59 // Get dot product
60 double dp=arma::as_scalar(arma::trans(C.col(io))*S*Cr.col(io));
61 // Invert sign?
62 if(dp<0.0)
63 C.col(io)*=-1.0;
64 }
65 }
66
C_diff(const arma::mat & Cr,const arma::mat & S,arma::mat & Cc,const arma::vec & E)67 double C_diff(const arma::mat & Cr, const arma::mat & S, arma::mat & Cc, const arma::vec & E) {
68 size_t io=0;
69 while(io<Cc.n_cols-1) {
70 // Check for degenerate orbitals
71 size_t jo=io;
72 while(jo+1 < Cc.n_cols && fabs(E(io)-E(jo+1))<odegthr)
73 jo++;
74 if(jo!=io) {
75 // Calculate projection of orbitals
76 arma::mat Op(arma::trans(Cc.cols(io,jo))*S*Cr.cols(io,jo));
77 // Rotate orbitals to match original ones
78 Cc.cols(io,jo)=Cc.cols(io,jo)*Op;
79 }
80
81 // Next orbital
82 io=jo+1;
83 }
84
85 // Get norm
86 return rms_norm(Cc-Cr);
87 }
88
E_diff(const arma::mat & Er,const arma::vec & Ec)89 double E_diff(const arma::mat & Er, const arma::vec & Ec) {
90 // Normalization vector
91 arma::vec Enorm(arma::abs(Er));
92 Enorm=arma::max(Enorm,arma::ones(Er.n_elem));
93
94 return arma::max(arma::abs((Er-Ec)/Enorm));
95 }
96
97 Settings settings;
98
main(int argc,char ** argv)99 int main(int argc, char ** argv) {
100 // Get reference directory
101 char * refdir=getenv("ERKALE_REFDIR");
102 if(!refdir)
103 throw std::runtime_error("Need to set reference directory!\n");
104
105 fprintf(stderr,"Reference directory is set to \"%s\".\n",refdir);
106 fflush(stderr);
107
108 if(argc != 2) {
109 std::ostringstream oss;
110 oss << "Usage: " << argv[0] << " chkpt\n";
111 throw std::runtime_error(oss.str());
112 }
113
114 // Current file
115 std::string curfile(argv[1]);
116 if(!file_exists(curfile)) {
117 std::ostringstream oss;
118 oss << "Current checkpoint file \"" << curfile << "\" does not exist!\n";
119 throw std::runtime_error(oss.str());
120 }
121 Checkpoint cur(curfile,false);
122
123 // Reference file
124 std::string reffile=std::string(refdir)+"/"+curfile;
125 if(!file_exists(reffile)) {
126 std::ostringstream oss;
127 oss << "Reference checkpoint file \"" << reffile << "\" does not exist!\n";
128 throw std::runtime_error(oss.str());
129 }
130 Checkpoint ref(reffile,false);
131
132 // Check consistency of run type
133 int rrestr, crestr;
134 ref.read("Restricted",rrestr);
135 cur.read("Restricted",crestr);
136 if(rrestr != crestr)
137 throw std::runtime_error("Run types don't match!\n");
138
139 // Basis sets
140 BasisSet bref, bcur;
141 ref.read(bref);
142 cur.read(bcur);
143
144 // Check that basis sets match
145 if( ! (bref == bcur) )
146 throw std::runtime_error("Basis sets don't match!\n");
147
148 // Check normalization of basis
149 check_norm(bref);
150 check_norm(bcur);
151
152 // Overlap matrix
153 arma::mat S(bcur.overlap());
154
155 if(rrestr) {
156 int Nelr, Nelc;
157 ref.read("Nel",Nelr);
158 cur.read("Nel",Nelc);
159 if(Nelr != Nelc) throw std::runtime_error("Amount of electrons doesn't match!\n");
160
161 // Total energies
162 energy_t Er, Ec;
163 ref.read(Er);
164 cur.read(Ec);
165 double dE=Ec.E-Er.E;
166 printf("Total energy difference %e\n",dE);
167 fflush(stdout);
168 if(fabs(dE) > Etol) throw std::runtime_error("Total energies don't match!\n");
169
170 // Densities
171 arma::cx_mat Pref, Pcur;
172 arma::mat Prefr, Prefi, Pcurr, Pcuri;
173 ref.read("P",Prefr);
174 if(ref.exist("P_im")) {
175 ref.read("P_im",Prefi);
176 Pref=Prefr*COMPLEX1 + Prefi*COMPLEXI;
177 } else
178 Pref=Prefr*COMPLEX1;
179
180 cur.read("P",Pcurr);
181 if(cur.exist("P_im")) {
182 cur.read("P_im",Pcuri);
183 Pcur=Pcurr*COMPLEX1 + Pcuri*COMPLEXI;
184 } else
185 Pcur=Pcurr*COMPLEX1;
186
187 // Check norm
188 double Nelnum=std::real(arma::trace(Pcur*S));
189 double dNel=Nelnum-Nelc;
190 printf("Electron count difference %e\n",dNel);
191 fflush(stdout);
192 if(fabs(dNel)>dPtol) throw std::runtime_error("Norm of density matrix is wrong.\n");
193
194 // Check difference
195 double dP=rms_cnorm(Pcur-Pref);
196 printf("Density matrix difference %e\n",dP);
197 if(dP>dPtol) throw std::runtime_error("Density matrices differ!\n");
198
199 // Orbitals
200 arma::mat Cref, Ccur;
201 ref.read("C",Cref);
202 cur.read("C",Ccur);
203 if(Cref.n_cols != Ccur.n_cols) throw std::runtime_error("Amount of orbitals doesn't match!\n");
204
205 // Fix orbital phase signs
206 fix_signs(Ccur,Cref,S);
207
208 // Orbital energies
209 arma::vec Eref, Ecur;
210 ref.read("E",Eref);
211 cur.read("E",Ecur);
212
213 // Calculate differences of nondegenerate orbitals
214 double dC=C_diff(Cref,S,Ccur,Ecur);
215 printf("Orbital matrix difference %e\n",dC);
216 fflush(stdout);
217 // if(dC>dCtol) throw std::runtime_error("Orbital coefficients differ!\n");
218
219 // Orbital energy differences
220 double dEo=E_diff(Eref,Ecur);
221 printf("Orbital energy difference %e\n",dEo);
222 fflush(stdout);
223 if(dEo > dEtol) throw std::runtime_error("Orbital energies differ!\n");
224
225 // Check EMD
226 GaussianEMDEvaluator eval(bcur,Pcur);
227 EMD emd(&eval, &eval, Nelc, 0, 0);
228 emd.initial_fill(false);
229 emd.find_electrons(false);
230 emd.optimize_moments(false);
231
232 // Get moments
233 arma::mat mom=emd.moments();
234
235 // Compare <p^2> with T and <p^0> with tr(P*S)
236 double p0diff=mom(2,1)-std::real(arma::trace(Pcur*S));
237 double p0err=mom(2,2);
238 printf("<p^0> - N = % e, d<p^0> = %e\n",p0diff,p0err);
239 double p2diff=mom(4,1)-2.0*std::real(arma::trace(Pcur*bcur.kinetic()));
240 double p2err=mom(4,2);
241 printf("<p^2> - T = % e, d<p^2> = %e\n",p2diff,p2err);
242 fflush(stdout);
243 if(fabs(p0diff)>=P0TOL*p0err || fabs(p2diff)>P2TOL*p2err)
244 throw std::runtime_error("EMD failed.\n");
245
246 } else {
247 int Nelar, Nelbr, Nelr;
248 int Nelac, Nelbc, Nelc;
249 ref.read("Nel-a",Nelar);
250 ref.read("Nel-b",Nelbr);
251 ref.read("Nel",Nelr);
252
253 cur.read("Nel-a",Nelac);
254 cur.read("Nel-b",Nelbc);
255 cur.read("Nel",Nelc);
256 if(Nelar != Nelac || Nelbr != Nelbc || Nelr != Nelc)
257 throw std::runtime_error("Amount of electrons doesn't match!\n");
258
259 // Total energies
260 energy_t Er, Ec;
261 ref.read(Er);
262 cur.read(Ec);
263 double dE=Ec.E-Er.E;
264 printf("Total energy difference %e\n",dE);
265 fflush(stdout);
266 if(fabs(dE) > Etol) throw std::runtime_error("Total energies don't match!\n");
267
268 // Densities
269 arma::mat Paref, Pbref, Prefr, Prefi;
270 arma::mat Pacur, Pbcur, Pcurr, Pcuri;
271 arma::cx_mat Pcur, Pref;
272 ref.read("Pa",Paref);
273 ref.read("Pb",Pbref);
274 ref.read("P",Prefr);
275 if(ref.exist("P_im")) {
276 ref.read("P_im",Prefi);
277 Pref=Prefr*COMPLEX1 + Prefi*COMPLEXI;
278 } else
279 Pref=Prefr*COMPLEX1;
280
281 cur.read("Pa",Pacur);
282 cur.read("Pb",Pbcur);
283 cur.read("P",Pcurr);
284 if(cur.exist("P_im")) {
285 cur.read("P_im",Pcuri);
286 Pcur=Pcurr*COMPLEX1 + Pcuri*COMPLEXI;
287 } else
288 Pcur=Pcurr*COMPLEX1;
289
290 // Check norm
291 double Nelanum=arma::trace(Pacur*S);
292 double Nelbnum=arma::trace(Pbcur*S);
293 double Nelnum=std::real(arma::trace(Pcur*S));
294
295 // XRS calculation?
296 bool xrs, xrsspin;
297 std::string xrsmethod;
298 try {
299 cur.read("XRSMethod",xrsmethod);
300 cur.read("XRSSpin",xrsspin);
301 xrs=true;
302 } catch(std::runtime_error &) {
303 xrs=false;
304 }
305
306 // Alpha electron count error
307 double dNela;
308 if(xrs && stricmp(xrsmethod,"TP")==0 && !xrsspin)
309 dNela=Nelanum + 0.5 - Nelac;
310 else if(xrs && stricmp(xrsmethod,"FCH")==0 && !xrsspin)
311 dNela=Nelanum + 1.0 - Nelac;
312 else
313 dNela=Nelanum - Nelac;
314 printf("Alpha electron count difference %e\n",dNela);
315 fflush(stdout);
316 if(fabs(dNela)>dPtol) throw std::runtime_error("Norm of alpha density matrix is wrong.\n");
317
318 // Beta electron count error
319 double dNelb;
320 if(xrs && stricmp(xrsmethod,"TP")==0 && xrsspin)
321 dNelb=Nelbnum + 0.5 - Nelbc;
322 else if(xrs && stricmp(xrsmethod,"FCH")==0 && xrsspin)
323 dNelb=Nelbnum + 1.0 - Nelbc;
324 else
325 dNelb=Nelbnum - Nelbc;
326 printf("Beta electron count difference %e\n",dNela);
327 fflush(stdout);
328 if(fabs(dNelb)>dPtol) throw std::runtime_error("Norm of beta density matrix is wrong.\n");
329
330 double dNel;
331 if(xrs && stricmp(xrsmethod,"TP")==0)
332 dNel=Nelnum + 0.5 - Nelc;
333 else if(xrs && stricmp(xrsmethod,"FCH")==0)
334 dNel=Nelnum + 1.0 - Nelc;
335 else
336 dNel=Nelnum - Nelc;
337 printf("Total electron count difference %e\n",dNela);
338 fflush(stdout);
339 if(fabs(dNel)>dPtol)
340 throw std::runtime_error("Norm of total density matrix is wrong.\n");
341
342 // Check differences
343 double dPa=rms_norm(Pacur-Paref);
344 printf("Alpha density matrix difference %e\n",dPa);
345 fflush(stdout);
346 if(dPa>dPtol) throw std::runtime_error("Alpha density matrices differ!\n");
347 double dPb=rms_norm(Pbcur-Pbref);
348 printf("Beta density matrix difference %e\n",dPb);
349 fflush(stdout);
350 if(dPb>dPtol) throw std::runtime_error("Density matrices differ!\n");
351
352 double dP=rms_cnorm(Pcur-Pref);
353 printf("Total density matrix difference %e\n",dP);
354 fflush(stdout);
355 if(dP >dPtol) throw std::runtime_error("Total density matrices differ!\n");
356
357 // Orbitals
358 arma::mat Caref, Cbref, Cacur, Cbcur;
359 ref.read("Ca",Caref);
360 ref.read("Cb",Cbref);
361 cur.read("Ca",Cacur);
362 cur.read("Cb",Cbcur);
363
364 if(Caref.n_cols != Cacur.n_cols)
365 throw std::runtime_error("Amount of alpha orbitals doesn't match!\n");
366 if(Cbref.n_cols != Cbcur.n_cols)
367 throw std::runtime_error("Amount of beta orbitals doesn't match!\n");
368
369 // Fix orbital phase signs
370 fix_signs(Cacur,Caref,S);
371 fix_signs(Cbcur,Cbref,S);
372
373 // Orbital energies
374 arma::vec Earef, Ebref, Eacur, Ebcur;
375 ref.read("Ea",Earef);
376 ref.read("Eb",Ebref);
377 cur.read("Ea",Eacur);
378 cur.read("Eb",Ebcur);
379
380 // Check differences
381 double dCa=C_diff(Caref,S,Cacur,Eacur);
382 printf("Alpha orbital matrix difference %e\n",dCa);
383 // if(dCa>dCtol) throw std::runtime_error("Alpha orbital coefficients differ!\n");
384 double dCb=C_diff(Cbref,S,Cbcur,Ebcur);
385 printf("Beta orbital matrix difference %e\n",dCb);
386 // if(dCb>dCtol) throw std::runtime_error("Beta orbital coefficients differ!\n");
387
388 double dEoa=E_diff(Earef,Eacur);
389 printf("Alpha orbital energy difference %e\n",dEoa);
390 fflush(stdout);
391 if(dEoa > dEtol) throw std::runtime_error("Alpha orbital energies differ!\n");
392 double dEob=E_diff(Ebref,Ebcur);
393 printf("Beta orbital energy difference %e\n",dEob);
394 fflush(stdout);
395 if(dEob > dEtol) throw std::runtime_error("Beta orbital energies differ!\n");
396
397 // Amount of electrons in density matrix
398 double Pnel(Nelc);
399 if(xrs) {
400 if(stricmp(xrsmethod,"FCH")==0)
401 Pnel-=1.0;
402 else if(stricmp(xrsmethod,"TP")==0)
403 Pnel-=0.5;
404 }
405
406 // Check EMD
407 GaussianEMDEvaluator eval(bcur,Pcur);
408 EMD emd(&eval, &eval, Pnel, 0, 0);
409 emd.initial_fill(false);
410 emd.find_electrons(false);
411 emd.optimize_moments(false);
412 // Get moments
413 arma::mat mom=emd.moments();
414
415 // Compare <p^2> with T and <p^0> with tr(P*S)
416 double p0diff=mom(2,1)-std::real(arma::trace(Pcur*S));
417 double p0err=mom(2,2);
418 printf("<p^0> - N = %e, d<p^0> = %e\n",p0diff,p0err);
419 double p2diff=mom(4,1)-2.0*std::real(arma::trace(Pcur*bcur.kinetic()));
420 double p2err=mom(4,2);
421 printf("<p^2> - T = %e, d<p^2> = %e\n",p2diff,p2err);
422 fflush(stdout);
423 if(fabs(p0diff)>=P0TOL*p0err || fabs(p2diff)>P2TOL*p2err)
424 throw std::runtime_error("EMD failed.\n");
425 }
426 }
427