2 #include <mpqc_config.h>
3 #include <util/misc/exenv.h>
5 #undef SCF_CHECK_BOUNDS
7 #ifdef SCF_CHECK_BOUNDS
8 #define CHECK(ival,pval,ij,kl,con) check(ival,pval,ij,kl,con)
10 #define CHECK(ival,pval,ij,kl,con)
17 double *
const RESTRICT gmat;
18 double *
const RESTRICT pmat;
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)
31 if ( 1.000001 * ibound_ < (ival > 0 ? ival : -ival)) {
33 ExEnv::errn() <<
" bound = " << ibound_ << std::endl;
37 if ( 1.000001 * pbound_ < (pval > 0 ? pval : -pval)) {
39 ExEnv::errn() <<
" bound = " << pbound_ << std::endl;
46 ExEnv::errn() <<
" cont = " << contrib << std::endl;
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");
56 inline void cont2(
int ij,
int kl,
double val) {
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");
62 inline void cont3(
int ij,
int kl,
double val) {
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");
68 inline void cont4(
int ij,
int kl,
double val) {
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");
74 inline void cont5(
int ij,
int kl,
double val) {
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");
89 void set_bound(
double,
double) {}
96 inline void cont1(
int ij,
int kl,
double val) {
97 ec += val*pmat[ij]*pmat[kl];
100 inline void cont2(
int ij,
int kl,
double val) {
101 ex -= 0.25*val*pmat[ij]*pmat[kl];
104 inline void cont3(
int ij,
int kl,
double val) {
105 ex -= 0.5*val*pmat[ij]*pmat[kl];
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];
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];
127 inline double cont1(
int ij,
int kl) {
128 return pmat[ij]*pmat[kl];
131 inline double cont2(
int ij,
int kl) {
132 return pmat[ij]*pmat[kl];