28 #ifndef _chemistry_qc_libint2_boundstimpl_h
29 #define _chemistry_qc_libint2_boundstimpl_h
34 #include <util/misc/scexception.h>
35 #include <chemistry/qc/libint2/int2e.h>
39 struct abssqrt :
public std::unary_function<const T&,T> {
40 T operator()(
const T& x) {
41 if (x < 0)
return std::sqrt(std::negate<T>()(x));
42 else return std::sqrt(x);
49 template <
class Int2e>
50 BoundsLibint2<Int2e>::BoundsLibint2(Integral*integral,
51 const Ref<GaussianBasisSet>& b1,
52 const Ref<GaussianBasisSet>& b2,
53 const Ref<GaussianBasisSet>& b3,
54 const Ref<GaussianBasisSet>& b4,
56 const Ref<IntParams>& params) :
57 nsh2_(b2->nshell()), nsh4_(b4->nshell())
59 equiv_12_34_ = b1->equiv(b3) && b2->equiv(b4);
60 equiv_12_43_ = b1->equiv(b4) && b2->equiv(b3);
61 equiv_1_2_ = b1->equiv(b2);
62 equiv_3_4_ = b3->equiv(b4);
64 Ref<MessageGrp> msg = integral->messagegrp();
65 const int ntasks = msg->n();
66 const int me = msg->me();
68 const int nsh1 = b1->nshell();
69 const int nsh2 = b2->nshell();
70 const int n12 = nsh1*nsh2;
72 libint2::Int2eCreator<Int2e> creator;
73 Ref<Int2e> int12 = creator(integral,b1,b2,b1,b2,storage,params);
74 Q12_.resize(n12);
for(
int i=0; i<n12; ++i) Q12_[i] = (int_bound_t)0;
75 double* buf = int12->buffer(0);
77 for(
int s1=0; s1<nsh1; ++s1) {
78 const int nf1 = b1->shell(s1).nfunction();
79 for(
int s2=0; s2<nsh2; ++s2, ++f12) {
80 const int nf2 = b2->shell(s2).nfunction();
81 const int nf = nf1*nf2*nf1*nf2;
86 int12->compute_quartet(&s1,&s2,&s1,&s2);
89 const double max_elem = *(std::max_element(buf,buf+nf));
90 Q12_[f12] = Log2Bounds::bound_cast(max_elem);
96 msg->sum(&Q12_[0],Q12_.size());
99 if (!equiv_12_34_ && !equiv_12_43_) {
101 const int nsh3 = b3->nshell();
102 const int nsh4 = b4->nshell();
103 const int n34 = nsh3*nsh4;
105 libint2::Int2eCreator<Int2e> creator;
106 Ref<Int2e> int34 = creator(integral,b3,b4,b3,b4,storage,params);
107 Q34_.resize(n34);
for(
int i=0; i<n34; ++i) Q34_[i] = (int_bound_t)0;
108 double* buf = int34->buffer(0);
110 for(
int s3=0; s3<nsh3; ++s3) {
111 const int nf3 = b3->shell(s3).nfunction();
112 for(
int s4=0; s4<nsh4; ++s4, ++f34) {
113 const int nf4 = b4->shell(s4).nfunction();
114 const int nf = nf3*nf4*nf3*nf4;
116 if (f34%ntasks != me)
119 int34->compute_quartet(&s3,&s4,&s3,&s4);
122 const double max_elem = *(std::max_element(buf,buf+nf));
123 Q34_[f34] = Log2Bounds::bound_cast(max_elem);
129 msg->sum(&Q34_[0],Q34_.size());
137 else if (equiv_12_43_) {
140 for(
int s2=0; s2<nsh2; ++s2) {
141 for(
int s1=0; s1<nsh1; ++s1, ++f21) {
142 Q34_[f21] = Q12_[s1*nsh2 + s2];
148 if (debugclass_ > 0) {
149 ExEnv::out0() << indent <<
"BoundsLibint2::Q12 :" << std::endl;
154 for(
int s1=0; s1<nsh1; ++s1) {
155 const int s2max = equiv_1_2_ ? s1 : nsh2-1;
156 for(
int s2=0; s2<=s2max; ++s2) {
157 const int f12 = s1*nsh2 + s2;
159 << int(Q12_[f12]) << std::endl;
164 if (debugclass_ > 0 && not equiv_12_34_ && not equiv_12_43_) {
165 const int nsh3 = b3->nshell();
166 const int nsh4 = b4->nshell();
167 ExEnv::out0() << indent <<
"BoundsLibint2::Q34 :" << std::endl;
172 for(
int s3=0; s3<nsh3; ++s3) {
173 const int s4max = equiv_3_4_ ? s3 : nsh4-1;
174 for(
int s4=0; s4<=s4max; ++s4) {
175 const int f34 = s3*nsh4 + s4;
177 << int(Q34_[f34]) << std::endl;
184 template <
class Int2e>
185 BoundsLibint2<Int2e>::~BoundsLibint2()
189 template <
class Int2e>
193 return Q12_[nsh2_*sh1 + sh2] + Q34_[nsh4_*sh3 + sh4];
197 #endif // header guard