1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: asd/asd_compute_rdm.hpp
4 // Copyright (C) 2014 Toru Shiozaki
5 //
6 // Author: Inkoo Kim <inkoo.kim@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 #ifdef ASD_HEADERS
26
27 #ifndef BAGEL_ASD_COMPUTE_RDM_H
28 #define BAGEL_ASD_COMPUTE_RDM_H
29
30 template <class VecType>
compute_rdm12()31 void ASD<VecType>::compute_rdm12() {
32 Timer rdmtime;
33 std::cout << std::endl << " ===== ASD RDM Computation ==== " << std::endl;
34 const int norbA = dimer_->active_refs().first->nact();
35 const int norbB = dimer_->active_refs().second->nact();
36
37 if (rdm1_av_ == nullptr && nstates_ > 1) {
38 rdm1_av_ = std::make_shared<RDM<1>>(norbA+norbB);
39 rdm2_av_ = std::make_shared<RDM<2>>(norbA+norbB);
40 }
41 if (nstates_ > 1) {
42 rdm1_av_->zero();
43 rdm2_av_->zero();
44 }
45
46 compute_rdm12_dimer(); //allocation takes place
47 std::cout << " o Dimer RDM - " << std::setw(9) << std::fixed << std::setprecision(2) << rdmtime.tick() << std::endl;
48
49 compute_rdm12_monomer();
50 std::cout << " o Monomer RDM - " << std::setw(9) << std::fixed << std::setprecision(2) << rdmtime.tick() << std::endl;
51
52 //State-average RDM
53 if (nstates_ != 1) {
54 for (int i = 0; i != nstates_; ++i) {
55 rdm1_av_->ax_plus_y(weight_[i], rdm1_[i]);
56 rdm2_av_->ax_plus_y(weight_[i], rdm2_[i]);
57 }
58 } else {
59 rdm1_av_ = rdm1_[0];
60 rdm2_av_ = rdm2_[0];
61 }
62
63 //print information
64 if (print_info_) {
65 for (int i = 0; i != nstates_; ++i) {
66 print_rdm_info(rdm1_[i], rdm2_[i], i);
67 print_energy_info(rdm1_[i], rdm2_[i], i);
68 }
69 }
70
71 std::cout << " ============================== " << std::endl;
72 }
73
74
75 template <class VecType>
compute_rdm12_monomer()76 void ASD<VecType>::compute_rdm12_monomer() {
77 const int nactA = dimer_->active_refs().first->nact();
78 const int nactB = dimer_->active_refs().second->nact();
79
80 std::vector<std::shared_ptr<RDM<1>>> rdm1A; rdm1A.resize(nstates_);
81 std::vector<std::shared_ptr<RDM<2>>> rdm2A; rdm2A.resize(nstates_);
82 std::vector<std::shared_ptr<RDM<1>>> rdm1B; rdm1B.resize(nstates_);
83 std::vector<std::shared_ptr<RDM<2>>> rdm2B; rdm2B.resize(nstates_);
84
85 //allocate first (filled with zero)
86 for (int kst = 0; kst != nstates_; ++kst) {
87 rdm1A[kst] = std::make_shared<RDM<1>>(nactA);
88 rdm2A[kst] = std::make_shared<RDM<2>>(nactA);
89 rdm1B[kst] = std::make_shared<RDM<1>>(nactB);
90 rdm2B[kst] = std::make_shared<RDM<2>>(nactB);
91 }
92
93 const int n = subspaces_.size() * nstates_;
94 const int numproc = mpi__->size();
95 const int myid = mpi__->rank();
96
97 std::vector<std::pair<int,int>> proc;
98 for (int i = 0; i != subspaces_.size(); ++i)
99 for (int j = 0; j != nstates_; ++j)
100 proc.emplace_back(i,j);
101 assert(proc.size() == n);
102
103 int mystart, myend;
104 if (n >= numproc) {
105 mystart = (n / numproc) * myid + ((n % numproc) < myid ? (n % numproc) : myid);
106 myend = mystart + (n / numproc) + ((n % numproc) > myid);
107 }
108 else {
109 mystart = myid < n ? myid : -1;
110 myend = myid < n ? (myid + 1) : 0;
111 }
112
113 //diagonal dimer subspace
114 for (int job = mystart; job != myend; ++job) {
115
116 if (job < 0) break; //MPI workers without allocated jobs, fill zero and wait for reduction
117
118 const int isub = proc[job].first;
119 const int kst = proc[job].second;
120
121 auto& AB = subspaces_[isub];
122 std::shared_ptr<const VecType> A = AB.template ci<0>(); //CASDvec, RASDvec
123 std::shared_ptr<const VecType> B = AB.template ci<1>();
124 const int ioff = AB.offset();
125 const int nstA = A->ij(); //XDvec size
126 const int nstB = B->ij();
127 assert(nstA == AB.nstatesA());
128 assert(nstB == AB.nstatesB());
129
130 std::shared_ptr<RDM<1>> rdm1;
131 std::shared_ptr<RDM<2>> rdm2;
132
133 {//MonomerA i.e. delta_JJ'
134 auto dket = contract_I(A, adiabats_, ioff, nstA, nstB, kst);
135 for (int j = 0; j != nstB; ++j) {
136 std::tie(rdm1, rdm2) = compute_rdm12_monomer(dket, j);
137 rdm1A[kst]->ax_plus_y(1.0, rdm1);
138 rdm2A[kst]->ax_plus_y(1.0, rdm2);
139 }
140 }
141
142 {//MonomerB i.e. delta_II'
143 auto dket = contract_J(B, adiabats_, ioff, nstA, nstB, kst);
144 for (int i = 0; i != nstA; ++i) {
145 std::tie(rdm1, rdm2) = compute_rdm12_monomer(dket, i);
146 rdm1B[kst]->ax_plus_y(1.0, rdm1);
147 rdm2B[kst]->ax_plus_y(1.0, rdm2);
148 }
149 }
150 } //subspaces
151
152 #ifdef HAVE_MPI_H
153 for (int kst = 0; kst != nstates_; ++kst) {
154 mpi__->allreduce(rdm1A[kst]->data(), rdm1A[kst]->size());
155 mpi__->allreduce(rdm2A[kst]->data(), rdm2A[kst]->size());
156 mpi__->allreduce(rdm1B[kst]->data(), rdm1B[kst]->size());
157 mpi__->allreduce(rdm2B[kst]->data(), rdm2B[kst]->size());
158 }
159 #endif
160
161 for (int istate = 0; istate != nstates_; ++istate) {
162 auto rdm1 = std::make_shared<RDM<1>>(nactA+nactB);
163 {//A
164 auto low = {0,0};
165 auto up = {nactA,nactA};
166 auto outv = make_rwview(rdm1->range().slice(low,up), rdm1->storage());
167 copy(rdm1A[istate]->begin(), rdm1A[istate]->end(), outv.begin());
168 }
169 {//B
170 auto low = {nactA,nactA};
171 auto up = {nactA+nactB,nactA+nactB};
172 auto outv = make_rwview(rdm1->range().slice(low,up), rdm1->storage());
173 copy(rdm1B[istate]->begin(), rdm1B[istate]->end(), outv.begin());
174 }
175 auto rdm2 = std::make_shared<RDM<2>>(nactA+nactB);
176 {//A
177 auto low = {0,0,0,0};
178 auto up = {nactA,nactA,nactA,nactA};
179 auto outv = make_rwview(rdm2->range().slice(low,up), rdm2->storage());
180 copy(rdm2A[istate]->begin(), rdm2A[istate]->end(), outv.begin());
181 }
182 {//B
183 auto low = {nactA,nactA,nactA,nactA};
184 auto up = {nactA+nactB,nactA+nactB,nactA+nactB,nactA+nactB};
185 auto outv = make_rwview(rdm2->range().slice(low,up), rdm2->storage());
186 copy(rdm2B[istate]->begin(), rdm2B[istate]->end(), outv.begin());
187 }
188
189 //update rdms
190 *rdm1_[istate] += *rdm1;
191 *rdm2_[istate] += *rdm2;
192 }
193
194 }
195 #endif
196
197 #endif
198