00001 00002 #include <scconfig.h> 00003 #include <util/misc/exenv.h> 00004 00005 #undef SCF_CHECK_BOUNDS 00006 00007 #ifdef SCF_CHECK_BOUNDS 00008 #define CHECK(ival,pval,ij,kl,con) check(ival,pval,ij,kl,con) 00009 #else 00010 #define CHECK(ival,pval,ij,kl,con) 00011 #endif 00012 00013 namespace sc { 00014 00015 class LocalCLHFContribution { 00016 private: 00017 double * const restrictxx gmat; 00018 double * const restrictxx pmat; 00019 00020 double ibound_; 00021 double pbound_; 00022 00023 public: 00024 LocalCLHFContribution(double *g, double *p) : gmat(g), pmat(p) {} 00025 ~LocalCLHFContribution() {} 00026 00027 void set_bound(double i, double p) { ibound_ = i; pbound_ = p; } 00028 void check(double ival, double pval, int ij, int kl, const char *contrib) 00029 { 00030 int bad = 0; 00031 if ( 1.000001 * ibound_ < (ival > 0 ? ival : -ival)) { 00032 ExEnv::errn() << "BAD INTEGRAL BOUND" << std::endl; 00033 ExEnv::errn() << " bound = " << ibound_ << std::endl; 00034 ExEnv::errn() << " value = " << ival << std::endl; 00035 bad = 1; 00036 } 00037 if ( 1.000001 * pbound_ < (pval > 0 ? pval : -pval)) { 00038 ExEnv::errn() << "BAD DENSITY BOUND" << std::endl; 00039 ExEnv::errn() << " bound = " << pbound_ << std::endl; 00040 ExEnv::errn() << " value = " << pval << std::endl; 00041 bad = 1; 00042 } 00043 if (bad) { 00044 ExEnv::errn() << " ij = " << ij << std::endl; 00045 ExEnv::errn() << " kl = " << kl << std::endl; 00046 ExEnv::errn() << " cont = " << contrib << std::endl; 00047 abort(); 00048 } 00049 } 00050 00051 inline void cont1(int ij, int kl, double val) { 00052 gmat[ij] += val*pmat[kl]; CHECK(val,pmat[kl],ij,kl,"cont1a"); 00053 gmat[kl] += val*pmat[ij]; CHECK(val,pmat[ij],ij,kl,"cont1b"); 00054 } 00055 00056 inline void cont2(int ij, int kl, double val) { 00057 val *= -0.25; 00058 gmat[ij] += val*pmat[kl]; CHECK(4*val,0.25*pmat[kl],ij,kl,"cont2a"); 00059 gmat[kl] += val*pmat[ij]; CHECK(4*val,0.25*pmat[ij],ij,kl,"cont2b"); 00060 } 00061 00062 inline void cont3(int ij, int kl, double val) { 00063 val *= -0.5; 00064 gmat[ij] += val*pmat[kl]; CHECK(2*val,0.5*pmat[kl],ij,kl,"cont3a"); 00065 gmat[kl] += val*pmat[ij]; CHECK(2*val,0.5*pmat[ij],ij,kl,"cont3b"); 00066 } 00067 00068 inline void cont4(int ij, int kl, double val) { 00069 val *= 0.75; 00070 gmat[ij] += val*pmat[kl]; CHECK(4./3.*val,0.75*pmat[kl],ij,kl,"cont4a"); 00071 gmat[kl] += val*pmat[ij]; CHECK(4./3.*val,0.75*pmat[ij],ij,kl,"cont4b"); 00072 } 00073 00074 inline void cont5(int ij, int kl, double val) { 00075 val *= 0.5; 00076 gmat[ij] += val*pmat[kl]; CHECK(2*val,0.5*pmat[kl],ij,kl,"cont5a"); 00077 gmat[kl] += val*pmat[ij]; CHECK(2*val,0.5*pmat[ij],ij,kl,"cont5b"); 00078 } 00079 }; 00080 00081 class LocalCLHFEnergyContribution { 00082 private: 00083 double * const pmat; 00084 00085 public: 00086 double ec; 00087 double ex; 00088 00089 void set_bound(double,double) {} 00090 00091 LocalCLHFEnergyContribution(double *p) : pmat(p) { 00092 ec=ex=0; 00093 } 00094 ~LocalCLHFEnergyContribution() {} 00095 00096 inline void cont1(int ij, int kl, double val) { 00097 ec += val*pmat[ij]*pmat[kl]; 00098 } 00099 00100 inline void cont2(int ij, int kl, double val) { 00101 ex -= 0.25*val*pmat[ij]*pmat[kl]; 00102 } 00103 00104 inline void cont3(int ij, int kl, double val) { 00105 ex -= 0.5*val*pmat[ij]*pmat[kl]; 00106 } 00107 00108 inline void cont4(int ij, int kl, double val) { 00109 ec += val*pmat[ij]*pmat[kl]; 00110 ex -= 0.25*val*pmat[ij]*pmat[kl]; 00111 } 00112 00113 inline void cont5(int ij, int kl, double val) { 00114 ec += val*pmat[ij]*pmat[kl]; 00115 ex -= 0.5*val*pmat[ij]*pmat[kl]; 00116 } 00117 }; 00118 00119 class LocalCLHFGradContribution { 00120 private: 00121 double * const pmat; 00122 00123 public: 00124 LocalCLHFGradContribution(double *p) : pmat(p) {} 00125 ~LocalCLHFGradContribution() {} 00126 00127 inline double cont1(int ij, int kl) { 00128 return pmat[ij]*pmat[kl]; 00129 } 00130 00131 inline double cont2(int ij, int kl) { 00132 return pmat[ij]*pmat[kl]; 00133 } 00134 }; 00135 00136 }