MPQC  3.0.0-alpha
bounds.timpl.h
1 //
2 // bounds.h
3 //
4 // Copyright (C) 2007 Edward Valeev
5 //
6 // Author: Edward Valeev <evaleev@vt.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifndef _chemistry_qc_libint2_boundstimpl_h
29 #define _chemistry_qc_libint2_boundstimpl_h
30 
31 #include <vector>
32 #include <cmath>
33 #include <algorithm>
34 #include <util/misc/scexception.h>
35 #include <chemistry/qc/libint2/int2e.h>
36 
37 namespace {
38  template <typename T>
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);
43  }
44  };
45 }
46 
47 namespace sc {
48 
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,
55  size_t storage,
56  const Ref<IntParams>& params) :
57  nsh2_(b2->nshell()), nsh4_(b4->nshell())
58  {
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);
63 
64  Ref<MessageGrp> msg = integral->messagegrp();
65  const int ntasks = msg->n();
66  const int me = msg->me();
67 
68  const int nsh1 = b1->nshell();
69  const int nsh2 = b2->nshell();
70  const int n12 = nsh1*nsh2;
71  {
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);
76  int f12 = 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;
82 
83  if (f12%ntasks != me)
84  continue;
85 
86  int12->compute_quartet(&s1,&s2,&s1,&s2);
87 
88  std::transform(buf,buf+nf,buf,abssqrt<double>());
89  const double max_elem = *(std::max_element(buf,buf+nf));
90  Q12_[f12] = Log2Bounds::bound_cast(max_elem);
91  }
92  }
93  }
94  // propagate Q12_ among the nodes
95  {
96  msg->sum(&Q12_[0],Q12_.size());
97  }
98 
99  if (!equiv_12_34_ && !equiv_12_43_) {
100 
101  const int nsh3 = b3->nshell();
102  const int nsh4 = b4->nshell();
103  const int n34 = nsh3*nsh4;
104  {
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);
109  int f34 = 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;
115 
116  if (f34%ntasks != me)
117  continue;
118 
119  int34->compute_quartet(&s3,&s4,&s3,&s4);
120 
121  std::transform(buf,buf+nf,buf,abssqrt<double>());
122  const double max_elem = *(std::max_element(buf,buf+nf));
123  Q34_[f34] = Log2Bounds::bound_cast(max_elem);
124  }
125  }
126  }
127  // propagate Q34_ among the nodes
128  {
129  msg->sum(&Q34_[0],Q34_.size());
130  }
131 
132  }
133  else {
134  if (equiv_12_34_) {
135  Q34_ = Q12_;
136  }
137  else if (equiv_12_43_) {
138  Q34_.resize(n12);
139  int f21 = 0;
140  for(int s2=0; s2<nsh2; ++s2) {
141  for(int s1=0; s1<nsh1; ++s1, ++f21) {
142  Q34_[f21] = Q12_[s1*nsh2 + s2];
143  }
144  }
145  }
146  }
147 
148  if (debugclass_ > 0) {
149  ExEnv::out0() << indent << "BoundsLibint2::Q12 :" << std::endl;
150  ExEnv::out0() << indent << "bs1:" << std::endl;
151  b1->print_brief(ExEnv::out0());
152  ExEnv::out0() << indent << "bs2:" << std::endl;
153  b2->print_brief(ExEnv::out0());
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;
158  ExEnv::out0() << indent << s1 << " " << s2 << " "
159  << int(Q12_[f12]) << std::endl;
160  }
161  }
162  }
163 
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;
168  ExEnv::out0() << indent << "bs3:" << std::endl;
169  b3->print_brief(ExEnv::out0());
170  ExEnv::out0() << indent << "bs4:" << std::endl;
171  b4->print_brief(ExEnv::out0());
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;
176  ExEnv::out0() << indent << s3 << " " << s4 << " "
177  << int(Q34_[f34]) << std::endl;
178  }
179  }
180  }
181 
182  }
183 
184  template <class Int2e>
185  BoundsLibint2<Int2e>::~BoundsLibint2()
186  {
187  }
188 
189  template <class Int2e>
190  int
191  BoundsLibint2<Int2e>::log2_bound(int sh1, int sh2, int sh3, int sh4) const
192  {
193  return Q12_[nsh2_*sh1 + sh2] + Q34_[nsh4_*sh3 + sh4];
194  }
195 }
196 
197 #endif // header guard
198 
199 // Local Variables:
200 // mode: c++
201 // c-file-style: "CLJ"
202 // End:
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc::transform
RefSCMatrix transform(const OrbitalSpace &space2, const OrbitalSpace &space1, const Ref< SCMatrixKit > &kit=SCMatrixKit::default_matrixkit())
transform(s2,s1) returns matrix that transforms s1 to s2.
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::BoundsLibint2::log2_bound
int log2_bound(int sh1, int sh2, int sh3, int sh4) const
computes bound for a given type of integrals
Definition: bounds.timpl.h:191

Generated at Sun Jan 26 2020 23:23:58 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.