1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * 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
18
19 #include <cfloat>
20 #include "diis.h"
21 #include "lbfgs.h"
22 #include "linalg.h"
23 #include "mathf.h"
24 #include "stringutil.h"
25
26 // Maximum allowed absolute weight for a Fock matrix
27 #define MAXWEIGHT 10.0
28 // Trigger cooloff if energy rises more than
29 #define COOLTHR 0.1
30
operator <(const diis_pol_entry_t & lhs,const diis_pol_entry_t & rhs)31 bool operator<(const diis_pol_entry_t & lhs, const diis_pol_entry_t & rhs) {
32 return lhs.E < rhs.E;
33 }
34
operator <(const diis_unpol_entry_t & lhs,const diis_unpol_entry_t & rhs)35 bool operator<(const diis_unpol_entry_t & lhs, const diis_unpol_entry_t & rhs) {
36 return lhs.E < rhs.E;
37 }
38
DIIS(const arma::mat & S_,const arma::mat & Sinvh_,bool usediis_,double diiseps_,double diisthr_,bool useadiis_,bool verbose_,size_t imax_)39 DIIS::DIIS(const arma::mat & S_, const arma::mat & Sinvh_, bool usediis_, double diiseps_, double diisthr_, bool useadiis_, bool verbose_, size_t imax_) {
40 S=S_;
41 Sinvh=Sinvh_;
42 usediis=usediis_;
43 useadiis=useadiis_;
44 verbose=verbose_;
45 imax=imax_;
46
47 // Start mixing in DIIS weight when error is
48 diiseps=diiseps_;
49 // and switch to full DIIS when error is
50 diisthr=diisthr_;
51
52 // No cooloff
53 cooloff=0;
54 }
55
rDIIS(const arma::mat & S_,const arma::mat & Sinvh_,bool usediis_,double diiseps_,double diisthr_,bool useadiis_,bool verbose_,size_t imax_)56 rDIIS::rDIIS(const arma::mat & S_, const arma::mat & Sinvh_, bool usediis_, double diiseps_, double diisthr_, bool useadiis_, bool verbose_, size_t imax_) : DIIS(S_,Sinvh_,usediis_,diiseps_,diisthr_,useadiis_,verbose_,imax_) {
57 }
58
uDIIS(const arma::mat & S_,const arma::mat & Sinvh_,bool combine_,bool usediis_,double diiseps_,double diisthr_,bool useadiis_,bool verbose_,size_t imax_)59 uDIIS::uDIIS(const arma::mat & S_, const arma::mat & Sinvh_, bool combine_, bool usediis_, double diiseps_, double diisthr_, bool useadiis_, bool verbose_, size_t imax_) : DIIS(S_,Sinvh_,usediis_,diiseps_,diisthr_,useadiis_,verbose_,imax_), combine(combine_) {
60 }
61
~DIIS()62 DIIS::~DIIS() {
63 }
64
~rDIIS()65 rDIIS::~rDIIS() {
66 }
67
~uDIIS()68 uDIIS::~uDIIS() {
69 }
70
clear()71 void rDIIS::clear() {
72 stack.clear();
73 }
74
clear()75 void uDIIS::clear() {
76 stack.clear();
77 }
78
erase_last()79 void rDIIS::erase_last() {
80 stack.erase(stack.begin());
81 }
82
erase_last()83 void uDIIS::erase_last() {
84 stack.erase(stack.begin());
85 }
86
update(const arma::mat & F,const arma::mat & P,double E,double & error)87 void rDIIS::update(const arma::mat & F, const arma::mat & P, double E, double & error) {
88 // New entry
89 diis_unpol_entry_t hlp;
90 hlp.F=F;
91 hlp.P=P;
92 hlp.E=E;
93
94 // Compute error matrix
95 arma::mat errmat(F*P*S);
96 // FPS - SPF
97 errmat-=arma::trans(errmat);
98 // and transform it to the orthonormal basis (1982 paper, page 557)
99 errmat=arma::trans(Sinvh)*errmat*Sinvh;
100 // and store it
101 hlp.err=arma::vectorise(errmat);
102
103 // DIIS error is
104 error=arma::max(arma::max(arma::abs(errmat)));
105
106 // Is stack full?
107 if(stack.size()==imax) {
108 erase_last();
109 }
110 // Add to stack
111 stack.push_back(hlp);
112
113 // Update ADIIS helpers
114 PiF_update();
115 }
116
PiF_update()117 void rDIIS::PiF_update() {
118 const arma::mat & Fn=stack[stack.size()-1].F;
119 const arma::mat & Pn=stack[stack.size()-1].P;
120
121 // Update matrices
122 PiF.zeros(stack.size());
123 for(size_t i=0;i<stack.size();i++)
124 PiF(i)=arma::trace((stack[i].P-Pn)*Fn);
125
126 PiFj.zeros(stack.size(),stack.size());
127 for(size_t i=0;i<stack.size();i++)
128 for(size_t j=0;j<stack.size();j++)
129 PiFj(i,j)=arma::trace((stack[i].P-Pn)*(stack[j].F-Fn));
130 }
131
update(const arma::mat & Fa,const arma::mat & Fb,const arma::mat & Pa,const arma::mat & Pb,double E,double & error)132 void uDIIS::update(const arma::mat & Fa, const arma::mat & Fb, const arma::mat & Pa, const arma::mat & Pb, double E, double & error) {
133 // New entry
134 diis_pol_entry_t hlp;
135 hlp.Fa=Fa;
136 hlp.Fb=Fb;
137 hlp.Pa=Pa;
138 hlp.Pb=Pb;
139 hlp.E=E;
140
141 // Compute error matrices
142 arma::mat errmata(Fa*Pa*S);
143 arma::mat errmatb(Fb*Pb*S);
144 // FPS - SPF
145 errmata-=arma::trans(errmata);
146 errmatb-=arma::trans(errmatb);
147 // and transform them to the orthonormal basis (1982 paper, page 557)
148 errmata=arma::trans(Sinvh)*errmata*Sinvh;
149 errmatb=arma::trans(Sinvh)*errmatb*Sinvh;
150 // and store it
151 if(combine) {
152 hlp.err=arma::vectorise(errmata+errmatb);
153 } else {
154 hlp.err.zeros(errmata.n_elem+errmatb.n_elem);
155 hlp.err.subvec(0,errmata.n_elem-1)=arma::vectorise(errmata);
156 hlp.err.subvec(errmata.n_elem,hlp.err.n_elem-1)=arma::vectorise(errmatb);
157 }
158
159 // DIIS error is
160 error=arma::max(arma::abs(hlp.err));
161
162 // Is stack full?
163 if(stack.size()==imax) {
164 erase_last();
165 }
166 // Add to stack
167 stack.push_back(hlp);
168
169 // Update ADIIS helpers
170 PiF_update();
171 }
172
PiF_update()173 void uDIIS::PiF_update() {
174 const arma::mat & Fan=stack[stack.size()-1].Fa;
175 const arma::mat & Fbn=stack[stack.size()-1].Fb;
176 const arma::mat & Pan=stack[stack.size()-1].Pa;
177 const arma::mat & Pbn=stack[stack.size()-1].Pb;
178
179 // Update matrices
180 PiF.zeros(stack.size());
181 for(size_t i=0;i<stack.size();i++)
182 PiF(i)=arma::trace((stack[i].Pa-Pan)*Fan) + arma::trace((stack[i].Pb-Pbn)*Fbn);
183
184 PiFj.zeros(stack.size(),stack.size());
185 for(size_t i=0;i<stack.size();i++)
186 for(size_t j=0;j<stack.size();j++)
187 PiFj(i,j)=arma::trace((stack[i].Pa-Pan)*(stack[j].Fa-Fan))+arma::trace((stack[i].Pb-Pbn)*(stack[j].Fb-Fbn));
188 }
189
get_energies() const190 arma::vec rDIIS::get_energies() const {
191 arma::vec E(stack.size());
192 for(size_t i=0;i<stack.size();i++)
193 E(i)=stack[i].E;
194 return E;
195 }
196
get_diis_error() const197 arma::mat rDIIS::get_diis_error() const {
198 arma::mat err(stack[0].err.n_elem,stack.size());
199 for(size_t i=0;i<stack.size();i++)
200 err.col(i)=stack[i].err;
201 return err;
202 }
203
get_energies() const204 arma::vec uDIIS::get_energies() const {
205 arma::vec E(stack.size());
206 for(size_t i=0;i<stack.size();i++)
207 E(i)=stack[i].E;
208 return E;
209 }
get_diis_error() const210 arma::mat uDIIS::get_diis_error() const {
211 arma::mat err(stack[0].err.n_elem,stack.size());
212 for(size_t i=0;i<stack.size();i++)
213 err.col(i)=stack[i].err;
214 return err;
215 }
216
get_w_wrk()217 arma::vec DIIS::get_w_wrk() {
218 // DIIS error
219 arma::mat de=get_diis_error();
220 double err=arma::max(arma::abs(de.col(de.n_cols-1)));
221
222 // Weight
223 arma::vec w;
224
225 if(useadiis && !usediis) {
226 w=get_w_adiis();
227 if(verbose) {
228 printf("ADIIS weights\n");
229 print_mat(w.t(),"% .2e ");
230 }
231 } else if(!useadiis && usediis) {
232 // Use DIIS only if error is smaller than threshold
233 if(err>diisthr)
234 throw std::runtime_error("DIIS error too large for only DIIS to converge wave function.\n");
235
236 w=get_w_diis();
237
238 if(verbose) {
239 printf("DIIS weights\n");
240 print_mat(w.t(),"% .2e ");
241 }
242 } else if(useadiis && usediis) {
243 // Sliding scale: DIIS weight
244 double diisw=std::max(std::min(1.0 - (err-diisthr)/(diiseps-diisthr), 1.0), 0.0);
245 // ADIIS weght
246 double adiisw=1.0-diisw;
247
248 // Determine cooloff
249 if(cooloff>0) {
250 diisw=0.0;
251 cooloff--;
252 } else {
253 // Check if energy has increased
254 arma::vec E=get_energies();
255 if(E.n_elem>1 && E(E.n_elem-1)-E(E.n_elem-2) > COOLTHR) {
256 cooloff=2;
257 diisw=0.0;
258 }
259 }
260
261 w.zeros(de.n_cols);
262
263 // DIIS and ADIIS weights
264 arma::vec wd, wa;
265 if(diisw!=0.0) {
266 wd=get_w_diis();
267 w+=diisw*wd;
268 }
269 if(adiisw!=0.0) {
270 wa=get_w_adiis();
271 w+=adiisw*wa;
272 }
273
274 if(verbose) {
275 if(adiisw!=0.0) {
276 printf("ADIIS weights\n");
277 print_mat(wa.t(),"% .2e ");
278 }
279 if(diisw!=0.0) {
280 printf("CDIIS weights\n");
281 print_mat(wd.t(),"% .2e ");
282 }
283 if(adiisw!=0.0 && diisw!=0.0) {
284 printf(" DIIS weights\n");
285 print_mat(w.t(),"% .2e ");
286 }
287 }
288
289 } else
290 throw std::runtime_error("Nor DIIS or ADIIS has been turned on.\n");
291
292 return w;
293 }
294
get_w()295 arma::vec DIIS::get_w() {
296 arma::vec sol;
297 while(true) {
298 sol=get_w_wrk();
299 if(std::abs(sol(sol.n_elem-1))<=sqrt(DBL_EPSILON)) {
300 if(sol.n_elem > 2) {
301 if(verbose) printf("Weight on last matrix too small, dropping the oldest matrix.\n");
302 erase_last();
303 PiF_update();
304 } else {
305 if(verbose) printf("Weight on last matrix still too small; switching to damping.\n");
306 sol.zeros();
307
308 double damp=0.1;
309 sol(0)=1-damp;
310 sol(1)=damp;
311 break;
312 }
313 } else
314 break;
315 }
316
317 return sol;
318 }
319
get_w_diis() const320 arma::vec DIIS::get_w_diis() const {
321 arma::mat errs=get_diis_error();
322 return get_w_diis_wrk(errs);
323 }
324
get_w_diis_wrk(const arma::mat & errs) const325 arma::vec DIIS::get_w_diis_wrk(const arma::mat & errs) const {
326 // Size of LA problem
327 int N=(int) errs.n_cols;
328
329 // Array holding the errors
330 arma::mat B(N,N);
331 B.zeros();
332 // Compute errors
333 for(int i=0;i<N;i++)
334 for(int j=0;j<N;j++) {
335 B(i,j)=arma::dot(errs.col(i),errs.col(j));
336 }
337
338 /*
339 The C1-DIIS method is equivalent to solving the group of linear
340 equations
341 B w = lambda 1 (1)
342
343 where B is the error matrix, w are the DIIS weights, lambda is the
344 Lagrange multiplier that guarantees that the weights sum to unity,
345 and 1 stands for a unit vector (1 1 ... 1)^T.
346
347 By rescaling the weights as w -> w/lambda, equation (1) is
348 reverted to the form
349 B w = 1 (2)
350
351 which can easily be solved using SVD techniques.
352
353 Finally, the weights are renormalized to satisfy
354 \sum_i w_i = 1
355 which takes care of the Lagrange multipliers.
356 */
357
358 // Right-hand side of equation is
359 arma::vec rh(N);
360 rh.ones();
361
362 // Singular value decomposition
363 arma::mat U, V;
364 arma::vec sval;
365 if(!arma::svd(U,sval,V,B,"std")) {
366 throw std::logic_error("SVD failed in DIIS.\n");
367 }
368 //sval.print("Singular values");
369
370 // Form solution vector
371 arma::vec sol(N);
372 sol.zeros();
373 for(int i=0;i<N;i++) {
374 #if 0
375 // Perform Tikhonov regularization for the singular
376 // eigenvalues. This doesn't appear to work as well as just
377 // cutting out the spectrum.
378 double s(sval(i));
379 double t=1e-4;
380 double invs=s/(s*s + t*t);
381 sol += arma::dot(U.col(i),rh)*invs * V.col(i);
382 #else
383 // Just cut out the problematic part. Weirdly, it appears that any
384 // kind of screening of the singular eigenvalues results in poorer
385 // convergence behavior(!)
386 if(sval(i) != 0.0)
387 sol += arma::dot(U.col(i),rh)/sval(i) * V.col(i);
388 #endif
389 }
390
391 // Sanity check
392 if(arma::sum(sol)==0.0)
393 sol.ones();
394
395 // Normalize solution
396 //printf("Sum of weights is %e\n",arma::sum(sol));
397 sol/=arma::sum(sol);
398
399 return sol;
400 }
401
solve_F(arma::mat & F)402 void rDIIS::solve_F(arma::mat & F) {
403 arma::vec sol(get_w());
404 // Form weighted Fock matrix
405 F.zeros();
406 for(size_t i=0;i<stack.size();i++)
407 F+=sol(i)*stack[i].F;
408 }
409
solve_F(arma::mat & Fa,arma::mat & Fb)410 void uDIIS::solve_F(arma::mat & Fa, arma::mat & Fb) {
411 arma::vec sol(get_w());
412
413 // Form weighted Fock matrix
414 Fa.zeros();
415 Fb.zeros();
416 for(size_t i=0;i<stack.size();i++) {
417 Fa+=sol(i)*stack[i].Fa;
418 Fb+=sol(i)*stack[i].Fb;
419 }
420 }
421
solve_P(arma::mat & P)422 void rDIIS::solve_P(arma::mat & P) {
423 arma::vec sol(get_w());
424
425 // Form weighted density matrix
426 P.zeros();
427 for(size_t i=0;i<stack.size();i++)
428 P+=sol(i)*stack[i].P;
429 }
430
solve_P(arma::mat & Pa,arma::mat & Pb)431 void uDIIS::solve_P(arma::mat & Pa, arma::mat & Pb) {
432 arma::vec sol(get_w());
433
434 // Form weighted density matrix
435 Pa.zeros();
436 Pb.zeros();
437 for(size_t i=0;i<stack.size();i++) {
438 Pa+=sol(i)*stack[i].Pa;
439 Pb+=sol(i)*stack[i].Pb;
440 }
441 }
442
find_minE(const std::vector<std::pair<double,double>> & steps,double & Emin,size_t & imin)443 static void find_minE(const std::vector< std::pair<double,double> > & steps, double & Emin, size_t & imin) {
444 Emin=steps[0].second;
445 imin=0;
446 for(size_t i=1;i<steps.size();i++)
447 if(steps[i].second < Emin) {
448 Emin=steps[i].second;
449 imin=i;
450 }
451 }
452
compute_c(const arma::vec & x)453 static arma::vec compute_c(const arma::vec & x) {
454 // Compute contraction coefficients
455 return x%x/arma::dot(x,x);
456 }
457
compute_jac(const arma::vec & x)458 static arma::mat compute_jac(const arma::vec & x) {
459 // Compute jacobian of transformation: jac(i,j) = dc_i / dx_j
460
461 // Compute coefficients
462 arma::vec c(compute_c(x));
463 double xnorm=arma::dot(x,x);
464
465 arma::mat jac(c.n_elem,c.n_elem);
466 for(size_t i=0;i<c.n_elem;i++) {
467 double ci=c(i);
468 double xi=x(i);
469
470 for(size_t j=0;j<c.n_elem;j++) {
471 double xj=x(j);
472
473 jac(i,j)=-ci*2.0*xj/xnorm;
474 }
475
476 // Extra term on diagonal
477 jac(i,i)+=2.0*xi/xnorm;
478 }
479
480 return jac;
481 }
482
get_w_adiis() const483 arma::vec DIIS::get_w_adiis() const {
484 // Number of parameters
485 size_t N=PiF.n_elem;
486
487 if(N==1) {
488 // Trivial case.
489 arma::vec ret(1);
490 ret.ones();
491 return ret;
492 }
493
494 // Starting point: equal weights on all matrices
495 arma::vec x=arma::ones<arma::vec>(N)/N;
496
497 // BFGS accelerator
498 LBFGS bfgs;
499
500 // Step size
501 double steplen=0.01, fac=2.0;
502
503 for(size_t iiter=0;iiter<1000;iiter++) {
504 // Get gradient
505 //double E(get_E_adiis(x));
506 arma::vec g(get_dEdx_adiis(x));
507 if(arma::norm(g,2)<=1e-7) {
508 break;
509 }
510
511 // Search direction
512 bfgs.update(x,g);
513 arma::vec sd(-bfgs.solve());
514
515 // Do a line search on the search direction
516 std::vector< std::pair<double, double> > steps;
517 // First, we try a fraction of the current step length
518 {
519 std::pair<double, double> p;
520 p.first=steplen/fac;
521 p.second=get_E_adiis(x+sd*p.first);
522 steps.push_back(p);
523 }
524 // Next, we try the current step length
525 {
526 std::pair<double, double> p;
527 p.first=steplen;
528 p.second=get_E_adiis(x+sd*p.first);
529 steps.push_back(p);
530 }
531
532 // Minimum energy and index
533 double Emin;
534 size_t imin;
535
536 while(true) {
537 // Sort the steps in length
538 std::sort(steps.begin(),steps.end());
539
540 // Find the minimum energy
541 find_minE(steps,Emin,imin);
542
543 // Where is the minimum?
544 if(imin==0 || imin==steps.size()-1) {
545 // Need smaller step
546 std::pair<double,double> p;
547 if(imin==0) {
548 p.first=steps[imin].first/fac;
549 if(steps[imin].first<DBL_EPSILON)
550 break;
551 } else {
552 p.first=steps[imin].first*fac;
553 }
554 p.second=get_E_adiis(x+sd*p.first);
555 steps.push_back(p);
556 } else {
557 // Optimum is somewhere in the middle
558 break;
559 }
560 }
561
562 if((imin!=0) && (imin!=steps.size()-1)) {
563 // Interpolate: A b = y
564 arma::mat A(3,3);
565 arma::vec y(3);
566 for(size_t i=0;i<3;i++) {
567 A(i,0)=1.0;
568 A(i,1)=steps[imin+i-1].first;
569 A(i,2)=std::pow(A(i,1),2);
570
571 y(i)=steps[imin+i-1].second;
572 }
573
574 arma::mat b;
575 if(arma::solve(b,A,y) && b(2)>sqrt(DBL_EPSILON)) {
576 // Success in solution and parabola gives minimum.
577
578 // The minimum of the parabola is at
579 double x0=-b(1)/(2*b(2));
580
581 // Is this an interpolation?
582 if(A(0,1) < x0 && x0 < A(2,1)) {
583 // Do the calculation with the interpolated step
584 std::pair<double,double> p;
585 p.first=x0;
586 p.second=get_E_adiis(x+sd*p.first);
587 steps.push_back(p);
588
589 // Find the minimum energy
590 find_minE(steps,Emin,imin);
591 }
592 }
593 }
594
595 if(steps[imin].first<DBL_EPSILON)
596 break;
597
598 // Switch to the minimum geometry
599 x+=steps[imin].first*sd;
600 // Store optimal step length
601 steplen=steps[imin].first;
602
603 //printf("Step %i: energy decreased by %e, gradient norm %e\n",(int) iiter+1,steps[imin].second-E,arma::norm(g,2)); fflush(stdout);
604 }
605
606 // Calculate weights
607 return compute_c(x);
608 }
609
get_E_adiis(const arma::vec & x) const610 double DIIS::get_E_adiis(const arma::vec & x) const {
611 // Consistency check
612 if(x.n_elem != PiF.n_elem) {
613 throw std::domain_error("Incorrect number of parameters.\n");
614 }
615
616 arma::vec c(compute_c(x));
617
618 // Compute energy
619 double Eval=0.0;
620 Eval+=2.0*arma::dot(c,PiF);
621 Eval+=arma::as_scalar(arma::trans(c)*PiFj*c);
622
623 return Eval;
624 }
625
get_dEdx_adiis(const arma::vec & x) const626 arma::vec DIIS::get_dEdx_adiis(const arma::vec & x) const {
627 // Compute contraction coefficients
628 arma::vec c(compute_c(x));
629
630 // Compute derivative of energy
631 arma::vec dEdc=2.0*PiF + PiFj*c + arma::trans(PiFj)*c;
632
633 // Compute jacobian of transformation: jac(i,j) = dc_i / dx_j
634 arma::mat jac(compute_jac(x));
635
636 // Finally, compute dEdx by plugging in Jacobian of transformation
637 // dE/dx_i = dc_j/dx_i dE/dc_j
638 return arma::trans(jac)*dEdc;
639 }
640