1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: dimer_jop.h
4 // Copyright (C) 2012 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
26
27 #ifndef __BAGEL_DIMER_JOP_H
28 #define __BAGEL_DIMER_JOP_H
29
30 #include <src/ci/fci/mofile.h>
31
32 namespace bagel {
33
34 class DimerJop : public Jop {
35 protected:
36 std::pair<std::shared_ptr<MOFile>, std::shared_ptr<MOFile>> jops_;
37
38 // Array is big enough to store all possible coulomb matrices just for simplicity
39 std::array<std::shared_ptr<const Matrix>, 16> matrices_;
40 std::shared_ptr<const Matrix> cross_mo1e_;
41
42 std::pair<int, int> nact_;
43
44 public:
45 DimerJop(const std::shared_ptr<const Reference> ref, const int nstart, const int nfenceA, const int nfenceB,
46 std::shared_ptr<const Coeff> coeff); // note that in DimerJop, I'm forcing a HZ Jop
47 DimerJop(const int nactA, const int nactB, std::shared_ptr<CSymMatrix> mo1e, std::shared_ptr<Matrix> mo2e);
48
cross_mo1e()49 std::shared_ptr<const Matrix> cross_mo1e() const { return cross_mo1e_; }
50
monomer_jop()51 template<int unit> std::shared_ptr<MOFile> monomer_jop() const { return ( unit == 0 ? jops_.first : jops_.second ); }
52
53 /**
54 \brief Generates a subset of the two-electron integrals \f$\langle ab|cd\rangle = (ac|bd)\f$
55 where each of \f$a,b,c,d\f$ can belong to either monomer \f$A\f$ or \f$B\f$.
56 \details Returns a Matrix(\f$N_A,N_B\f$) where \f$N_A=(\mbox{norb}_A)^d\f$,
57 \f$N_B=(\mbox{norb}_B)^{(4-d)}\f$ and \f$d\f$ is the number of \f$A\f$ operators.
58
59 \tparam A specifies whether the first index belongs to monomer \f$A\f$ (0) or \f$B\f$ (1)
60 \tparam B specifies whether the second index belongs to monomer \f$A\f$ (0) or \f$B\f$ (1)
61 \tparam C specifies whether the third index belongs to monomer \f$A\f$ (0) or \f$B\f$ (1)
62 \tparam D specifies whether the fourth index belongs to monomer \f$A\f$ (0) or \f$B\f$ (1)
63 */
64 template<int A, int B, int C, int D>
65 std::shared_ptr<const Matrix> coulomb_matrix();
66 /// Same as non-const version but will not create the matrix if it doesn't already exist
67 template<int A, int B, int C, int D>
68 std::shared_ptr<const Matrix> coulomb_matrix() const;
69
70 private:
71 void common_init(const int norbA, const int norbB);
72
nact()73 template <int unit> int nact() const { return ( unit == 0 ? nact_.first : nact_.second ); }
active(const int a)74 template <int unit> int active(const int a) const { return a + unit * nact_.first; }
75
76 template <int A, int B, int C, int D>
77 std::pair<int, int> index(const int a, const int b, const int c, const int d) const;
78 };
79
80 template<int A, int B, int C, int D>
coulomb_matrix()81 std::shared_ptr<const Matrix> DimerJop::coulomb_matrix() {
82 // First check to see if it's already stored
83 const int cindex = A + 2*B + 4*C + 8*D;
84 if (matrices_[cindex]) {
85 return matrices_[cindex];
86 }
87 else {
88 const int nactA = nact_.first;
89 const int nactB = nact_.second;
90
91 int ijA = 1;
92 int unitA = 4 - (A + B + C + D);
93 for ( int i = 0; i < unitA; ++i ) ijA *= nactA;
94
95 int ijB = 1;
96 int unitB = A + B + C + D;
97 for ( int i = 0; i < unitB; ++i ) ijB *= nactB;
98
99 auto out = std::make_shared<Matrix>(ijA, ijB);
100
101 for(int d = 0; d < nact<D>(); ++d) {
102 for(int c = 0; c < nact<C>(); ++c) {
103 for(int b = 0; b < nact<B>(); ++b) {
104 for(int a = 0; a < nact<A>(); ++a) {
105 int iA, jB;
106 std::tie(iA, jB) = index<A,B,C,D>(a,b,c,d);
107
108 out->element(iA,jB) = mo2e_hz(active<A>(a), active<B>(b), active<C>(c), active<D>(d));
109 }
110 }
111 }
112 }
113
114 out->localize();
115 matrices_[cindex] = out;
116 return out;
117 }
118 }
119
120 template<int A, int B, int C, int D>
coulomb_matrix()121 std::shared_ptr<const Matrix> DimerJop::coulomb_matrix() const { return matrices_[A + 2*B + 4*C + 8*D]; }
122
index(int a,int b,int c,int d)123 template <int A, int B, int C, int D> std::pair<int, int> DimerJop::index(int a, int b, int c, int d) const {
124 int iA = 0, jB = 0;
125 if (D == 0) iA = d; else jB = d;
126 if (C == 0) iA = c + iA*nact_.first; else jB = c + jB*nact_.second;
127 if (B == 0) iA = b + iA*nact_.first; else jB = b + jB*nact_.second;
128 if (A == 0) iA = a + iA*nact_.first; else jB = a + jB*nact_.second;
129 return {iA,jB};
130 }
131
132 }
133
134 #endif
135