1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: cphf.cc
4 // Copyright (C) 2012 Toru Shiozaki
5 //
6 // Author: Toru Shiozaki <shiozaki@northwestern.edu>
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // This program is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program. If not, see <http://www.gnu.org/licenses/>.
23 //
24
25 #include <src/grad/cphf.h>
26 #include <src/util/math/linearRM.h>
27
28 using namespace std;
29 using namespace bagel;
30
CPHF(const shared_ptr<const Matrix> grad,const VectorB & eig,const shared_ptr<const DFHalfDist> h,const shared_ptr<const Reference> r)31 CPHF::CPHF(const shared_ptr<const Matrix> grad, const VectorB& eig, const shared_ptr<const DFHalfDist> h,
32 const shared_ptr<const Reference> r)
33 : grad_(grad), eig_(eig), halfjj_(h), ref_(r), geom_(r->geom()) {
34
35 }
36
37
solve(const double zthresh,const int zmaxiter)38 shared_ptr<Matrix> CPHF::solve(const double zthresh, const int zmaxiter) {
39
40 LinearRM<Matrix> solver(zmaxiter, grad_);
41
42 const size_t nmobasis = ref_->coeff()->mdim();
43 const size_t nocca = ref_->nocc();
44 const size_t nvirt = nmobasis - nocca;
45
46 const MatView ocoeff = ref_->coeff()->slice(0, nocca);
47 const MatView vcoeff = ref_->coeff()->slice(nocca, nmobasis);
48
49 auto t = make_shared<Matrix>(nmobasis, nmobasis);
50 for (int i = 0; i != nocca; ++i)
51 for (int a = nocca; a != nvirt+nocca; ++a)
52 t->element(a,i) = grad_->element(a,i) / (eig_(a)-eig_(i));
53 t->scale(1.0/t->norm());
54
55 cout << " === Z-vector iteration ===" << endl << endl;
56
57 Timer timer;
58 for (int iter = 0; iter != zmaxiter; ++iter) {
59
60 auto sigma = make_shared<Matrix>(nmobasis, nmobasis);
61 // one electron part
62 for (int i = 0; i != nocca; ++i)
63 for (int a = nocca; a != nocca+nvirt; ++a)
64 (*sigma)(a,i) = (eig_(a)-eig_(i)) * t->element(a,i);
65
66 // J part
67 shared_ptr<const Matrix> tvo = t->get_submatrix(nocca, 0, nvirt, nocca);
68 auto pbmao = make_shared<Matrix>(ocoeff ^ (vcoeff * *tvo));
69 pbmao->symmetrize();
70 Matrix jri = *geom_->df()->compute_Jop(pbmao) * ocoeff;
71 Matrix jai = (vcoeff % jri) * 4.0;
72
73 // K part
74 // halfjj is an half transformed DF integral with J^{-1}_{DE}, given by the constructor
75 shared_ptr<const Matrix> kir = halfjj_->compute_Kop_1occ(pbmao, -2.0);
76 Matrix kia = *kir * vcoeff;
77
78 for (int i = 0; i != nocca; ++i)
79 for (int a = 0; a != nvirt; ++a)
80 (*sigma)(a+nocca,i) += jai(a,i) + kia(i,a);
81
82 t = solver.compute_residual(t, sigma);
83
84 cout << setw(7) << iter << " " << setw(20) << setprecision(14) << t->rms() << setw(15) << setprecision(2) << timer.tick() << endl;
85 if (t->rms() < zthresh) break;
86
87 for (int i = 0; i != nocca; ++i)
88 for (int a = nocca; a != nvirt+nocca; ++a)
89 t->element(a,i) /= (eig_(a)-eig_(i));
90 t->scale(1.0/t->norm());
91 }
92
93 cout << endl;
94 t = solver.civec();
95 t->fill_upper();
96 return t;
97
98 }
99