1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: reldf.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
26 #include <src/df/reldf.h>
27 #include <src/df/complexdf.h>
28
29 using namespace std;
30 using namespace bagel;
31
RelDF(shared_ptr<const DFDist> df,pair<int,int> cartesian,const vector<int> alpha)32 RelDF::RelDF(shared_ptr<const DFDist> df, pair<int, int> cartesian, const vector<int> alpha) : RelDFBase(cartesian), alpha_(alpha), swap_(false) {
33 if (!dynamic_pointer_cast<const ComplexDFDist>(df)) {
34 dfdata_[0] = df;
35 dfdata_[1] = nullptr;
36 } else {
37 auto cdf = dynamic_pointer_cast<const ComplexDFDist>(df);
38 array<shared_ptr<const DFDist>,2> tmp = cdf->split_real_imag();
39 dfdata_[0] = tmp[0];
40 dfdata_[1] = tmp[1];
41 }
42 set_basis();
43 }
44
RelDF(const RelDF & o,bool coo)45 RelDF::RelDF(const RelDF& o, bool coo) : RelDFBase(o), alpha_(o.alpha_), dfdata_(o.get_data()), swap_(o.swap_) {
46 set_basis();
47 if (coo) {
48 swap_ ^= true;
49 vector<shared_ptr<const SpinorInfo>> newbas;
50 for (auto& i : basis_)
51 newbas.push_back(i->swap());
52 basis_ = newbas;
53 std::swap(cartesian_.first, cartesian_.second);
54 }
55 }
56
57
58 //swap cartesian
swap() const59 shared_ptr<const RelDF> RelDF::swap() const {
60 return make_shared<const RelDF>(*this, true);
61 }
62
63
compute_half_transform(array<shared_ptr<const Matrix>,4> rc,array<shared_ptr<const Matrix>,4> ic) const64 vector<shared_ptr<RelDFHalf>> RelDF::compute_half_transform(array<shared_ptr<const Matrix>,4> rc, array<shared_ptr<const Matrix>,4> ic) const {
65 vector<shared_ptr<RelDFHalf>> out;
66
67 // first make a subset
68 vector<vector<shared_ptr<const SpinorInfo>>> subsets;
69 for (int i = 0; i != 4; ++i) {
70 vector<shared_ptr<const SpinorInfo>> tmp;
71 for (auto& j : basis())
72 if (j->basis(0) == i)
73 tmp.push_back(j);
74 if (!tmp.empty())
75 subsets.push_back(tmp);
76 }
77
78 // transform
79 for (auto& i : subsets)
80 out.push_back(make_shared<RelDFHalf>(shared_from_this(), i, rc, ic));
81 return out;
82 }
83
84
compute_Jop(list<shared_ptr<const RelCDMatrix>> & cd) const85 vector<shared_ptr<ZMatrix>> RelDF::compute_Jop(list<shared_ptr<const RelCDMatrix>>& cd) const {
86
87 vector<shared_ptr<ZVectorB>> sum;
88 for (auto& b : basis_) {
89 sum.push_back(cd.front()->clone());
90 for (auto& i : cd) {
91 if(b->alpha_comp() == i->alpha_comp())
92 *sum.back() += *i;
93 }
94 }
95
96 vector<shared_ptr<ZMatrix>> out;
97
98 // real
99 if (!get_imag()) {
100 for (auto& i : sum) {
101 shared_ptr<const Matrix> rdat = get_real()->compute_Jop_from_cd(i->get_real_part());
102 shared_ptr<const Matrix> idat = get_real()->compute_Jop_from_cd(i->get_imag_part());
103 out.push_back(make_shared<ZMatrix>(*rdat, *idat));
104 }
105
106 // complex (GIAO)
107 } else {
108 // get_real() + get_imag() needed for 3-multiplication algorithm
109 auto dfri = make_shared<DFDist>(get_real());
110 const int n = get_real()->block().size();
111 for (int i=0; i!=n; ++i) {
112 dfri->add_block(get_real()->block(i)->copy());
113 dfri->block(i)->ax_plus_y(1.0, get_imag()->block(i));
114 }
115 for (auto& i : sum) {
116 shared_ptr<const VectorB> cdr = i->get_real_part();
117 shared_ptr<const VectorB> cdi = i->get_imag_part();
118 auto cdri = make_shared<const VectorB> (*cdr + *cdi);
119
120 shared_ptr<Matrix> rdat = get_real()->compute_Jop_from_cd(cdr);
121 shared_ptr<Matrix> tmpdat = get_imag()->compute_Jop_from_cd(cdi);
122 shared_ptr<Matrix> idat = dfri->compute_Jop_from_cd(cdri);
123 *idat -= *rdat;
124 *idat -= *tmpdat;
125 *rdat -= *tmpdat;
126 out.push_back(make_shared<ZMatrix>(*rdat, *idat));
127 }
128 }
129
130 return out;
131 }
132
133
134