MPQC 2.3.1

- Main Page
- Related Pages
- Classes
- Files

- File List

clhftmpl.h
1
2 #include <scconfig.h>
3 #include <util/misc/exenv.h>
4
5 #undef SCF_CHECK_BOUNDS
6
7 #ifdef SCF_CHECK_BOUNDS
8 #define CHECK(ival,pval,ij,kl,con) check(ival,pval,ij,kl,con)
9 #else
10 #define CHECK(ival,pval,ij,kl,con)
11 #endif
12
13 namespace sc {
14
15 class LocalCLHFContribution {
16 private:
17 double * const restrictxx gmat;
18 double * const restrictxx pmat;
19
20 double ibound_;
21 double pbound_;
22
23 public:
24 LocalCLHFContribution(double *g, double *p) : gmat(g), pmat(p) {}
25 ~LocalCLHFContribution() {}
26
27 void set_bound(double i, double p) { ibound_ = i; pbound_ = p; }
28 void check(double ival, double pval, int ij, int kl, const char *contrib)
29 {
30 int bad = 0;
31 if ( 1.000001 * ibound_ < (ival > 0 ? ival : -ival)) {
32 ExEnv::errn() << "BAD INTEGRAL BOUND" << std::endl;
33 ExEnv::errn() << " bound = " << ibound_ << std::endl;
34 ExEnv::errn() << " value = " << ival << std::endl;
35 bad = 1;
36 }
37 if ( 1.000001 * pbound_ < (pval > 0 ? pval : -pval)) {
38 ExEnv::errn() << "BAD DENSITY BOUND" << std::endl;
39 ExEnv::errn() << " bound = " << pbound_ << std::endl;
40 ExEnv::errn() << " value = " << pval << std::endl;
41 bad = 1;
42 }
43 if (bad) {
44 ExEnv::errn() << " ij = " << ij << std::endl;
45 ExEnv::errn() << " kl = " << kl << std::endl;
46 ExEnv::errn() << " cont = " << contrib << std::endl;
47 abort();
48 }
49 }
50
51 inline void cont1(int ij, int kl, double val) {
52 gmat[ij] += val*pmat[kl]; CHECK(val,pmat[kl],ij,kl,"cont1a");
53 gmat[kl] += val*pmat[ij]; CHECK(val,pmat[ij],ij,kl,"cont1b");
54 }
55
56 inline void cont2(int ij, int kl, double val) {
57 val *= -0.25;
58 gmat[ij] += val*pmat[kl]; CHECK(4*val,0.25*pmat[kl],ij,kl,"cont2a");
59 gmat[kl] += val*pmat[ij]; CHECK(4*val,0.25*pmat[ij],ij,kl,"cont2b");
60 }
61
62 inline void cont3(int ij, int kl, double val) {
63 val *= -0.5;
64 gmat[ij] += val*pmat[kl]; CHECK(2*val,0.5*pmat[kl],ij,kl,"cont3a");
65 gmat[kl] += val*pmat[ij]; CHECK(2*val,0.5*pmat[ij],ij,kl,"cont3b");
66 }
67
68 inline void cont4(int ij, int kl, double val) {
69 val *= 0.75;
70 gmat[ij] += val*pmat[kl]; CHECK(4./3.*val,0.75*pmat[kl],ij,kl,"cont4a");
71 gmat[kl] += val*pmat[ij]; CHECK(4./3.*val,0.75*pmat[ij],ij,kl,"cont4b");
72 }
73
74 inline void cont5(int ij, int kl, double val) {
75 val *= 0.5;
76 gmat[ij] += val*pmat[kl]; CHECK(2*val,0.5*pmat[kl],ij,kl,"cont5a");
77 gmat[kl] += val*pmat[ij]; CHECK(2*val,0.5*pmat[ij],ij,kl,"cont5b");
78 }
79 };
80
81 class LocalCLHFEnergyContribution {
82 private:
83 double * const pmat;
84
85 public:
86 double ec;
87 double ex;
88
89 void set_bound(double,double) {}
90
91 LocalCLHFEnergyContribution(double *p) : pmat(p) {
92 ec=ex=0;
93 }
94 ~LocalCLHFEnergyContribution() {}
95
96 inline void cont1(int ij, int kl, double val) {
97 ec += val*pmat[ij]*pmat[kl];
98 }
99
100 inline void cont2(int ij, int kl, double val) {
101 ex -= 0.25*val*pmat[ij]*pmat[kl];
102 }
103
104 inline void cont3(int ij, int kl, double val) {
105 ex -= 0.5*val*pmat[ij]*pmat[kl];
106 }
107
108 inline void cont4(int ij, int kl, double val) {
109 ec += val*pmat[ij]*pmat[kl];
110 ex -= 0.25*val*pmat[ij]*pmat[kl];
111 }
112
113 inline void cont5(int ij, int kl, double val) {
114 ec += val*pmat[ij]*pmat[kl];
115 ex -= 0.5*val*pmat[ij]*pmat[kl];
116 }
117 };
118
119 class LocalCLHFGradContribution {
120 private:
121 double * const pmat;
122
123 public:
124 LocalCLHFGradContribution(double *p) : pmat(p) {}
125 ~LocalCLHFGradContribution() {}
126
127 inline double cont1(int ij, int kl) {
128 return pmat[ij]*pmat[kl];
129 }
130
131 inline double cont2(int ij, int kl) {
132 return pmat[ij]*pmat[kl];
133 }
134 };
135
136 }