clhftmpl.h Source File

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