1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: ras/rasci_denom.cc
4 // Copyright (C) 2013 Toru Shiozaki
5 //
6 // Author: Shane Parker <shane.parker@u.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/ci/ras/rasci.h>
26 #include <src/ci/ras/denomtask.h>
27
28 using namespace std;
29 using namespace bagel;
30
compute()31 void RAS::DenomTask::compute() {
32 const int nspin = abit_.count() - stringb_->nele();
33 const int nspin2 = nspin*nspin;
34 const int norb = stringb_->norb();
35
36 double* iter = data_;
37 const bitset<nbit__> ia = abit_;
38 for (auto& ib : *stringb_) {
39 const int nopen = (ia^ib).count();
40 const double F = (nopen >> 1) ? (static_cast<double>(nspin2 - nopen)/(nopen*(nopen-1))) : 0.0;
41 *iter = 0.0;
42 for (int i = 0; i < norb; ++i) {
43 const int nia = ia[i];
44 const int nib = ib[i];
45 const int niab = (nia + nib);
46 const int Ni = (nia ^ nib);
47 for (int j = 0; j < i; ++j) {
48 const int nja = ia[j];
49 const int njb = ib[j];
50 const int Nj = nja ^ njb;
51 const int addj = niab * (nja + njb);
52
53 *iter += jop_[j+norb*i] * 2.0 * addj - kop_[j+norb*i] * (F*Ni*Nj + addj);
54 }
55 *iter += h_[i] * niab - kop_[i+norb*i] * 0.5 * (Ni - niab*niab);
56 }
57 ++iter;
58 }
59 }
60
const_denom()61 void RASCI::const_denom() {
62 Timer denom_t;
63 unique_ptr<double[]> h(new double[norb_]);
64 unique_ptr<double[]> jop(new double[norb_*norb_]);
65 unique_ptr<double[]> kop(new double[norb_*norb_]);
66
67 for (int i = 0; i != norb_; ++i) {
68 for (int j = 0; j <= i; ++j) {
69 jop[i*norb_+j] = jop[j*norb_+i] = 0.5*jop_->mo2e_hz(j, i, j, i);
70 kop[i*norb_+j] = kop[j*norb_+i] = 0.5*jop_->mo2e_hz(j, i, i, j);
71 }
72 h[i] = jop_->mo1e(i,i);
73 }
74 denom_t.tick_print("jop, kop");
75
76 denom_ = make_shared<RASCivec>(det());
77
78 size_t tasksize = 0;
79 for (auto& iblock : denom_->blocks()) { if (iblock) tasksize += iblock->lena(); }
80 TaskQueue<RAS::DenomTask> tasks(tasksize);
81
82 for (auto& iblock : denom_->blocks()) {
83 if ( !iblock ) continue;
84 double* iter = iblock->data();
85 for (auto& ia : *iblock->stringsa()) {
86 tasks.emplace_back(iter, ia, iblock->stringsb(), jop.get(), kop.get(), h.get());
87 iter += iblock->lenb();
88 }
89 }
90
91 tasks.compute();
92 denom_t.tick_print("denom");
93 }
94
95
update(shared_ptr<const Coeff> c)96 void RASCI::update(shared_ptr<const Coeff> c) {
97 // iiii file to be created (MO transformation).
98 // now jop_->mo1e() and jop_->mo2e() contains one and two body part of Hamiltonian
99 Timer timer;
100 // Same Jop as used in FCI
101 jop_ = make_shared<Jop>(ref_, ncore_, ncore_+norb_, c, /*store*/false, "HZ");
102
103 // right now full basis is used.
104 cout << " * Integral transformation done. Elapsed time: " << setprecision(2) << timer.tick() << endl << endl;
105
106 const_denom();
107 }
108