1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: asd_ras_sigma.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/util/prim_op.h>
26 #include <src/asd/asd_ras.h>
27 #include <src/ci/ras/form_sigma.h>
28
29 using namespace std;
30 using namespace bagel;
31
form_sigma(shared_ptr<const RASDvec> ccvec,shared_ptr<const MOFile> jop) const32 shared_ptr<RASDvec> ASD_RAS::form_sigma(shared_ptr<const RASDvec> ccvec, shared_ptr<const MOFile> jop) const {
33 FormSigmaRAS form;
34 vector<int> conv(ccvec->ij(), static_cast<int>(false));
35 return form(ccvec, jop, conv);
36 }
37
form_sigma_1e(shared_ptr<const RASDvec> ccvec,const double * modata) const38 shared_ptr<RASDvec> ASD_RAS::form_sigma_1e(shared_ptr<const RASDvec> ccvec, const double* modata) const {
39 FormSigmaRAS form;
40 const int norb = ccvec->det()->norb();
41 auto mo1e = make_shared<Matrix>(norb, norb);
42 copy_n(modata, norb*norb, mo1e->data());
43 return form(ccvec, mo1e, nullptr, vector<int>(ccvec->ij(), static_cast<int>(false)));
44 }
45
46
compute_rdm12_monomer(shared_ptr<const RASDvec> civec,const int i) const47 tuple<shared_ptr<RDM<1>>,shared_ptr<RDM<2>>> ASD_RAS::compute_rdm12_monomer(shared_ptr<const RASDvec> civec, const int i) const {
48 shared_ptr<const RASCivec> cbra = civec->data(i);
49 const int norb = cbra->det()->norb();
50
51 auto dbra = make_shared<RASDvec>(cbra->det(), norb*norb);
52 dbra->zero();
53 sigma_2a(cbra, dbra);
54
55 return compute_rdm12_last_step(cbra, dbra);
56 }
57
sigma_2a(shared_ptr<const RASCivec> cc,shared_ptr<RASDvec> d) const58 void ASD_RAS::sigma_2a(shared_ptr<const RASCivec> cc, shared_ptr<RASDvec> d) const {
59 shared_ptr<const RASDeterminants> det = cc->det();
60
61 for (auto& block : cc->blocks()) {
62 if (block) {
63 const size_t a_offset = block->stringsa()->offset();
64 const size_t b_offset = block->stringsb()->offset();
65
66 for (size_t ia = 0, ab = 0; ia < block->stringsa()->size(); ++ia) {
67 const auto abit = block->string_bits_a(ia);
68 for (size_t jb = 0; jb < block->stringsb()->size(); ++jb, ++ab) {
69 const auto bbit = block->string_bits_b(jb);
70 const double coef = block->element(ab);
71
72 if (fabs(coef) < numerical_zero__) continue;
73 for (auto& phi : det->phib(jb + b_offset)) {
74 assert(phi.target == jb + b_offset);
75 const auto sbit = det->string_bits_b(phi.source);
76 if(det->allowed(abit,sbit))
77 d->data(phi.ij)->element(sbit,abit) += static_cast<double>(phi.sign) * coef; // i^+ j => (j,i)
78 }
79 for (auto& phi : det->phia(ia + a_offset)) {
80 assert(phi.target == ia + a_offset);
81 const auto sbit = det->string_bits_a(phi.source);
82 if(det->allowed(sbit,bbit))
83 d->data(phi.ij)->element(bbit,sbit) += static_cast<double>(phi.sign) * coef;
84 }
85
86 }
87 }
88 }
89 }
90
91 }
92
compute_rdm12_last_step(shared_ptr<const RASCivec> cibra,shared_ptr<const RASDvec> dbra) const93 tuple<shared_ptr<RDM<1>>, shared_ptr<RDM<2>>> ASD_RAS::compute_rdm12_last_step(shared_ptr<const RASCivec> cibra, shared_ptr<const RASDvec> dbra) const {
94 const int norb = cibra->det()->norb();
95 const int nri = dbra->data(0)->size();
96 auto cimat = make_shared<Matrix>(nri,norb*norb);
97 for (int ij = 0; ij != norb*norb; ++ij)
98 copy_n(dbra->data(ij)->data(), nri, cimat->element_ptr(0,ij));
99
100 // 1RDM
101 auto rdm1 = make_shared<RDM<1>>(norb);
102 dgemv_("T", nri, norb*norb, 1.0, cimat->element_ptr(0,0), nri, cibra->data(), 1, 0.0, rdm1->data(), 1);
103
104 // 2RDM
105 auto rdm2 = make_shared<RDM<2>>(norb);
106 {
107 auto rdmt = rdm2->clone();
108 dgemm_("T", "N", norb*norb, norb*norb, nri, 1.0, cimat->element_ptr(0,0), nri, cimat->element_ptr(0,0), nri, 0.0, rdmt->data(), norb*norb);
109
110 unique_ptr<double[]> buf(new double[norb*norb]);
111 for (int i = 0; i != norb; ++i)
112 for (int k = 0; k != norb; ++k) {
113 copy_n(rdmt->element_ptr(0,0,k,i), norb*norb, buf.get());
114 blas::transpose(buf.get(), norb, norb, rdmt->element_ptr(0,0,k,i));
115 }
116
117 sort_indices<2,3,0,1, 0,1, 1,1>(rdmt->data(), rdm2->data(), norb, norb, norb, norb);
118
119 // put in diagonal into 2RDM
120 // Gamma{i+ k+ l j} = Gamma{i+ j k+ l} - delta_jk Gamma{i+ l}
121 for (int i = 0; i != norb; ++i)
122 for (int k = 0; k != norb; ++k)
123 for (int j = 0; j != norb; ++j)
124 rdm2->element(j,k,k,i) -= rdm1->element(j,i);
125
126 //RDM2 symmetrize (out-of-excitation free parts are copied)
127 for (int i = 0, ij = 0; i != norb; ++i)
128 for (int j = 0; j != norb; ++j, ++ij)
129 for (int k = 0, kl = 0; k != norb; ++k)
130 for (int l = 0; l != norb; ++l, ++kl)
131 if (kl > ij) rdm2->element(i,j,k,l) = rdm2->element(k,l,i,j);
132 }
133
134 return tie(rdm1, rdm2);
135 }
136
contract_I(shared_ptr<const RASDvec> A,shared_ptr<Matrix> adiabats,int ioff,int nstA,int nstB,int kst) const137 shared_ptr<RASDvec> ASD_RAS::contract_I(shared_ptr<const RASDvec> A, shared_ptr<Matrix> adiabats, int ioff, int nstA, int nstB, int kst) const {
138 auto out = make_shared<RASDvec>(A->det(), nstB);
139 for (int ij = 0; ij != nstB; ++ij) {
140 out->data(ij)->zero();
141 }
142
143 for (int j = 0; j != nstB; ++j) {
144 for (int i = 0; i != nstA; ++i) {
145 const int ij = i + (j*nstA);
146 double u_ij = adiabats->element(ioff+ij,kst);
147
148 out->data(j)->ax_plus_y(u_ij, A->data(i));
149
150 }
151 }
152 return out;
153 }
154
contract_J(shared_ptr<const RASDvec> B,shared_ptr<Matrix> adiabats,int ioff,int nstA,int nstB,int kst) const155 shared_ptr<RASDvec> ASD_RAS::contract_J(shared_ptr<const RASDvec> B, shared_ptr<Matrix> adiabats, int ioff, int nstA, int nstB, int kst) const {
156 auto out = make_shared<RASDvec>(B->det(), nstA);
157 for (int ij = 0; ij != nstA; ++ij) {
158 out->data(ij)->zero();
159 }
160
161 for (int i = 0; i != nstA; ++i) {
162 for (int j = 0; j != nstB; ++j) {
163 const int ij = i + (j*nstA);
164 double u_ij = adiabats->element(ioff+ij,kst);
165
166 out->data(i)->ax_plus_y(u_ij, B->data(j));
167
168 }
169 }
170 return out;
171 }
172