tchftmpl.h
namespace sc {
class LocalTCContribution {
private:
double * const gmata;
double * const gmatb;
double * const kmata;
double * const kmatb;
double * const pmata;
double * const pmatb;
double * const opmata;
double * const opmatb;
public:
LocalTCContribution(double *ga, double *pa, double *gb, double *pb,
double *ka, double *opa, double *kb, double *opb) :
gmata(ga), gmatb(gb), kmata(ka), kmatb(kb),
pmata(pa), pmatb(pb), opmata(opa), opmatb(opb) {}
~LocalTCContribution() {}
void set_bound(double,double) {}
inline void cont1(int ij, int kl, double val) {
gmata[ij] += val*pmata[kl];
gmata[kl] += val*pmata[ij];
gmatb[ij] += val*pmatb[kl];
gmatb[kl] += val*pmatb[ij];
}
inline void cont2(int ij, int kl, double val) {
val *= 0.25;
gmata[ij] -= val*pmata[kl];
gmata[kl] -= val*pmata[ij];
gmatb[ij] -= val*pmatb[kl];
gmatb[kl] -= val*pmatb[ij];
kmata[ij] += val*opmata[kl];
kmata[kl] += val*opmata[ij];
kmatb[ij] += val*opmatb[kl];
kmatb[kl] += val*opmatb[ij];
}
inline void cont3(int ij, int kl, double val) {
val *= 0.5;
gmata[ij] -= val*pmata[kl];
gmata[kl] -= val*pmata[ij];
gmatb[ij] -= val*pmatb[kl];
gmatb[kl] -= val*pmatb[ij];
kmata[ij] += val*opmata[kl];
kmata[kl] += val*opmata[ij];
kmatb[ij] += val*opmatb[kl];
kmatb[kl] += val*opmatb[ij];
}
inline void cont4(int ij, int kl, double val) {
gmata[ij] += 0.75*val*pmata[kl];
gmata[kl] += 0.75*val*pmata[ij];
gmatb[ij] += 0.75*val*pmatb[kl];
gmatb[kl] += 0.75*val*pmatb[ij];
kmata[ij] += 0.25*val*opmata[kl];
kmata[kl] += 0.25*val*opmata[ij];
kmatb[ij] += 0.25*val*opmatb[kl];
kmatb[kl] += 0.25*val*opmatb[ij];
}
inline void cont5(int ij, int kl, double val) {
val *= 0.5;
gmata[ij] += val*pmata[kl];
gmata[kl] += val*pmata[ij];
gmatb[ij] += val*pmatb[kl];
gmatb[kl] += val*pmatb[ij];
kmata[ij] += val*opmata[kl];
kmata[kl] += val*opmata[ij];
kmatb[ij] += val*opmatb[kl];
kmatb[kl] += val*opmatb[ij];
}
};
class LocalTCEnergyContribution {
private:
double * const pmata;
double * const pmatb;
double * const opmata;
double * const opmatb;
public:
double eca;
double exa;
double ecb;
double exb;
double ecab;
double exab;
LocalTCEnergyContribution(double *pa, double *pb, double *opa, double *opb) :
pmata(pa), pmatb(pb), opmata(opa), opmatb(opb) {
exa=eca=0;
exb=ecb=0;
exab=ecab=0;
}
~LocalTCEnergyContribution() {}
void set_bound(double,double) {};
inline void cont1(int ij, int kl, double val) {
eca += val*pmata[ij]*pmata[kl];
ecb += val*pmatb[ij]*pmatb[kl];
ecab += val*(pmata[ij]*pmatb[kl]-pmatb[ij]*pmata[kl]);
}
inline void cont2(int ij, int kl, double val) {
exa -= 0.25*val*pmata[ij]*pmata[kl];
exb -= 0.25*val*pmatb[ij]*pmatb[kl];
exab -= 0.25*val*(pmata[ij]*pmatb[kl]-pmatb[ij]*pmata[kl]);
}
inline void cont3(int ij, int kl, double val) {
exa -= 0.5*val*pmata[ij]*pmata[kl];
exb -= 0.5*val*pmatb[ij]*pmatb[kl];
exab -= 0.5*val*(pmata[ij]*pmatb[kl]-pmatb[ij]*pmata[kl]);
}
inline void cont4(int ij, int kl, double val) {
eca += val*pmata[ij]*pmata[kl];
ecb += val*pmatb[ij]*pmatb[kl];
ecab += val*(pmata[ij]*pmatb[kl]-pmatb[ij]*pmata[kl]);
exa -= 0.25*val*pmata[ij]*pmata[kl];
exb -= 0.25*val*pmatb[ij]*pmatb[kl];
exab -= 0.25*val*(pmata[ij]*pmatb[kl]-pmatb[ij]*pmata[kl]);
}
inline void cont5(int ij, int kl, double val) {
eca += val*pmata[ij]*pmata[kl];
ecb += val*pmatb[ij]*pmatb[kl];
ecab += val*(pmata[ij]*pmatb[kl]-pmatb[ij]*pmata[kl]);
exa -= 0.5*val*pmata[ij]*pmata[kl];
exb -= 0.5*val*pmatb[ij]*pmatb[kl];
exab -= 0.5*val*(pmata[ij]*pmatb[kl]-pmatb[ij]*pmata[kl]);
}
};
class LocalTCGradContribution {
private:
double * const pmat;
double * const pmata;
double * const pmatb;
double c1sq;
double c2sq;
double c1c2;
public:
LocalTCGradContribution(double *p, double *pa, double *pb,
double c1, double c2) :
pmat(p), pmata(pa), pmatb(pb)
{
c1sq = c1*c1;
c2sq = c2*c2;
c1c2 = c1*c2;
}
~LocalTCGradContribution() {}
inline double cont1(int ij, int kl) {
return pmat[ij]*pmat[kl] +
c1sq*(pmata[ij]*pmat[kl] + pmat[ij]*pmata[kl]) +
c2sq*(pmatb[ij]*pmat[kl] + pmat[ij]*pmatb[kl]) +
0.5*c1sq*pmata[ij]*pmata[kl] +
0.5*c2sq*pmatb[ij]*pmatb[kl];
}
inline double cont2(int ij, int kl) {
return pmat[ij]*pmat[kl] +
c1sq*(pmata[ij]*pmat[kl] + pmat[ij]*pmata[kl]) +
c2sq*(pmatb[ij]*pmat[kl] + pmat[ij]*pmatb[kl]) -
c1c2*(pmata[ij]*pmatb[kl] + pmatb[ij]*pmata[kl]);
}
};
}