1 //
2 // fock.cc
3 //
4 // Copyright (C) 2004 Edward Valeev
5 //
6 // Author: Edward Valeev <edward.valeev@chemistry.gatech.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit 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 Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <stdexcept>
33 #include <chemistry/qc/mbptr12/r12int_eval.h>
34
35 using namespace std;
36 using namespace sc;
37
38 RefSCMatrix
fock_(const Ref<MOIndexSpace> & occ_space,const Ref<MOIndexSpace> & bra_space,const Ref<MOIndexSpace> & ket_space,double scale_J,double scale_K)39 R12IntEval::fock_(const Ref<MOIndexSpace>& occ_space, const Ref<MOIndexSpace>& bra_space,
40 const Ref<MOIndexSpace>& ket_space, double scale_J, double scale_K)
41 {
42 const Ref<GaussianBasisSet> bs1 = bra_space->basis();
43 const Ref<GaussianBasisSet> bs2 = ket_space->basis();
44 const bool bs1_eq_bs2 = (bs1 == bs2);
45 int nshell1 = bs1->nshell();
46 int nshell2 = bs2->nshell();
47
48 RefSCMatrix vec1t = bra_space->coefs().t();
49 RefSCMatrix vec2 = ket_space->coefs();
50
51 Ref<Integral> localints = r12info_->integral()->clone();
52 localints->set_basis(bs1,bs2);
53
54 Ref<OneBodyInt> h_ints = localints->hcore();
55
56 // form AO moment matrices
57 RefSCDimension aodim1 = vec1t.coldim();
58 RefSCDimension aodim2 = vec2.rowdim();
59 Ref<SCMatrixKit> aokit = bs1->so_matrixkit();
60 RefSCMatrix h(aodim1, aodim2, aokit);
61 h.assign(0.0);
62
63 for(int sh1=0; sh1<nshell1; sh1++) {
64 int bf1_offset = bs1->shell_to_function(sh1);
65 int nbf1 = bs1->shell(sh1).nfunction();
66
67 int sh2max;
68 if (bs1_eq_bs2)
69 sh2max = sh1;
70 else
71 sh2max = nshell2-1;
72
73 for(int sh2=0; sh2<=sh2max; sh2++) {
74 int bf2_offset = bs2->shell_to_function(sh2);
75 int nbf2 = bs2->shell(sh2).nfunction();
76
77 h_ints->compute_shell(sh1,sh2);
78 const double *hintsptr = h_ints->buffer();
79
80 int bf1_index = bf1_offset;
81 for(int bf1=0; bf1<nbf1; bf1++, bf1_index++, hintsptr+=nbf2) {
82 int bf2_index = bf2_offset;
83 const double *ptr = hintsptr;
84 int bf2max;
85 if (bs1_eq_bs2 && sh1 == sh2)
86 bf2max = bf1;
87 else
88 bf2max = nbf2-1;
89 for(int bf2=0; bf2<=bf2max; bf2++, bf2_index++, ptr++) {
90
91 h.set_element(bf1_index, bf2_index, *ptr);
92
93 }
94 }
95 }
96 }
97
98 // Symmetrize matrices, if necessary
99 if (bs1_eq_bs2) {
100 const int nbasis = bs1->nbasis();
101 for(int bf1=0; bf1<nbasis; bf1++)
102 for(int bf2=0; bf2<=bf1; bf2++) {
103 h(bf2,bf1) = h(bf1,bf2);
104 }
105 }
106
107 // finally, transform
108 RefSCMatrix F = vec1t * h * vec2;
109
110 // add coulomb and exchange parts
111 if (scale_J != 0.0) {
112 RefSCMatrix J = coulomb_(occ_space,bra_space,ket_space);
113 J.scale(2.0*scale_J); F.accumulate(J); J = 0;
114 }
115 if (scale_K != 0.0) {
116 RefSCMatrix K = exchange_(occ_space,bra_space,ket_space);
117 K.scale(-1.0*scale_K); F.accumulate(K); K = 0;
118 }
119
120 // and clean up a bit
121 h_ints = 0;
122 h = 0;
123
124 if (debug_ > 1) {
125 F.print("Fock matrix");
126 }
127
128 return F;
129 }
130
131 ///////////////////////////////////////////////////////////////
132
133 // Local Variables:
134 // mode: c++
135 // c-file-style: "CLJ"
136 // End:
137