1 //
2 // BAGEL - Parallel electron correlation program.
3 // Filename: multipolebatch.cc
4 // Copyright (C) 2014 Toru Shiozaki
5 //
6 // Author: Hai-Anh Le <anh@u.northwestern.edu>
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // The BAGEL package 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 3, or (at your option)
14 // any later version.
15 //
16 // The BAGEL package 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 BAGEL package; see COPYING.  If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 
26 
27 #include <src/integral/os/multipolebatch.h>
28 #include <src/integral/carsphlist.h>
29 #include <src/integral/hrrlist.h>
30 #include <src/integral/sortlist.h>
31 
32 using namespace std;
33 using namespace bagel;
34 
35 const static CHRRList hrr;
36 const static CCarSphList carsphlist;
37 
MultipoleBatch(const array<shared_ptr<const Shell>,2> & sh,const array<double,3> centre,const int lmax,shared_ptr<StackMem> stack)38 MultipoleBatch::MultipoleBatch(const array<shared_ptr<const Shell>,2>& sh, const array<double, 3> centre,
39                                const int lmax, shared_ptr<StackMem> stack)
40  : MultipoleBatch_base(sh, centre, lmax, stack) {
41 
42  const double integral_thresh = PRIM_SCREEN_THRESH;
43  allocate_arrays(prim0_ * prim1_);
44 
45  compute_ss(integral_thresh);
46 }
47 
compute()48 void MultipoleBatch::compute() {
49 
50   // First, do (ix, iy, iz; l; m) = (n|O_lm|0) using VRR
51 
52   complex<double>* intermediate_p = stack_->get<complex<double>>(size_alloc_);
53   const complex<double> imag = swap01_ ? complex<double>(0.0, -1.0) : complex<double>(0.0, 1.0);
54 
55   const int nmul = num_multipoles_;
56   vector<complex<double>*> bkup(nmul);
57   for (int imul = 0; imul != nmul; ++imul)
58     bkup[imul] = intermediate_p + imul * size_block_;
59 
60   for (int iprim = 0; iprim != prim0_ * prim1_; ++iprim) {
61 
62     vector<complex<double>> workz(amax1_ * nmul); // (0, 0, iz, l, m)
63 
64     // (0, 0, 0, l, m)
65     for (int imul = 0; imul != nmul; ++imul)
66       workz[imul] = multipole_[imul * prim0_ * prim1_ + iprim];
67 
68     int iang = 0;
69     for (int iz = 0; iz <= amax_; ++iz) {
70       // (0, 0, iz, l, m)
71       if (iz > 0) {
72         int imul = 0;
73         for (int l = 0; l <= lmax_; ++l) {
74           for (int m = 0; m <= 2 * l; ++m, ++imul) {
75             assert(l * l + m == imul);
76             workz[iz * nmul + imul] = -2.0 * xb_[iprim] * AB_[2] * workz[(iz-1)*nmul+imul];
77             if (iz > 1)
78               workz[iz * nmul + imul] += (iz-1.0) * workz[(iz-2)*nmul+imul];
79             if (l > 0 && abs(m-l) < l) {
80               const int i1 = (l-1)*(l-1) + m - 1;
81               workz[iz * nmul + imul] += workz[(iz-1)*nmul+i1];
82             }
83             workz[iz * nmul + imul] *= 0.5 / xp_[iprim];
84           }
85         }
86       }
87 
88       vector<complex<double>> worky((amax1_ - iz) * nmul);
89       for (int imul = 0; imul != nmul; ++imul)
90         worky[imul] = workz[iz * nmul + imul];
91 
92       for (int iy = 0; iy <= amax_ - iz; ++iy) {
93         // (0, iy, iz, l, m)
94         for (int b = 1; b <= iy; ++b) {
95           int imul = 0;
96           for (int l = 0; l <= lmax_; ++l) {
97             for (int m = 0; m <= 2 * l; ++m, ++imul) {
98               worky[b * nmul + imul] = -2.0 * xb_[iprim] * AB_[1] * worky[(b-1)*nmul+imul];
99               if (b > 1)
100                 worky[b * nmul + imul] += (b-1.0) * worky[(b-2)*nmul+imul];
101               if (l > 0) {
102                 const int i1 = (l-1)*(l-1) + m - 1;
103                 if (abs(m - l - 1) < l)
104                   worky[b * nmul + imul] -= imag * 0.5 * worky[(b-1)*nmul+i1-1];
105                 if (abs(m - l + 1) < l)
106                   worky[b * nmul + imul] -= imag * 0.5 * worky[(b-1)*nmul+i1+1];
107               }
108               worky[b * nmul + imul] *= 0.5 / xp_[iprim];
109             }
110           }
111         }
112 
113         vector<complex<double>> workx((amax1_ - iy - iz) * nmul);
114         for (int imul = 0; imul != nmul; ++imul)
115           workx[imul] = worky[iy * nmul + imul];
116 
117         for (int ix = max(0, amin_ - iy - iz); ix <= amax_ - iy - iz; ++ix, ++iang) {
118           // (ix, iy, iz, l, m)
119           for (int a = 1; a <= ix; ++a) {
120             int imul = 0;
121             for (int l = 0; l <= lmax_; ++l) {
122               for (int m = 0; m <= 2 * l; ++m, ++imul) {
123                 workx[a * nmul + imul] = -2.0 * xb_[iprim] * AB_[0] * workx[(a-1)*nmul+imul];
124                 if (a > 1)
125                   workx[a * nmul + imul] += (a-1.0) * workx[(a-2)*nmul+imul];
126                 if (l > 0) {
127                   const int i1 = (l-1)*(l-1) + m - 1;
128                   if (abs(m - l - 1) < l)
129                     workx[a * nmul + imul] -= 0.5 * workx[(a-1)*nmul+i1-1];
130                   if (abs(m - l + 1) < l)
131                     workx[a * nmul + imul] += 0.5 * workx[(a-1)*nmul+i1+1];
132                 }
133                 workx[a * nmul + imul] *= 0.5 / xp_[iprim];
134               }
135             }
136           }
137 
138           int imul = 0;
139           for (int l = 0; l <= lmax_; ++l) {
140             for (int m = 0; m <= 2 * l; ++m, ++imul) {
141               complex<double>* current_data = bkup[imul] + iprim * asize_;
142               const int pos = amapping_[ix + amax1_ * (iy + amax1_ * iz)];
143               if (!swap01_) {
144                 current_data[pos] = workx[ix * nmul + imul];
145               } else {
146                 //current_data[pos] = workx[ix * nmul + imul];
147                 current_data[pos] = conj(workx[ix * nmul + imul]);
148               }
149             }
150           }
151         } // ix
152       } //iy
153     } //iz
154     assert(iang <= asize_);
155 
156   } // end loop over primitives - done (n|O|0)
157 
158   const CSortList sort(spherical_);
159   complex<double>* data_start = intermediate_p;
160   complex<double>* data_final = data_;
161   fill_n(data_, size_alloc_, 0.0);
162 
163 
164   for (int imul = 0; imul != nmul; ++imul, data_start += size_block_, data_final += size_block_) {
165 
166     complex<double>* const intermediate_c = stack_->get<complex<double>>(size_block_);
167     perform_contraction(asize_, data_start, prim0_, prim1_, intermediate_c,
168                         basisinfo_[0]->contractions(), basisinfo_[0]->contraction_ranges(), cont0_,
169                         basisinfo_[1]->contractions(), basisinfo_[1]->contraction_ranges(), cont1_);
170 
171     complex<double>* const intermediate_fi = stack_->get<complex<double>>(size_block_);
172 
173     // now get (a|O_lm|b) using HRR
174     if (basisinfo_[1]->angular_number() != 0) {
175       const int hrr_index = basisinfo_[0]->angular_number() * ANG_HRR_END + basisinfo_[1]->angular_number();
176       hrr.hrrfunc_call(hrr_index, cont0_ * cont1_, intermediate_c, AB_, intermediate_fi);
177     } else {
178       copy_n(intermediate_c, size_block_, intermediate_fi);
179     }
180 
181     if (spherical_) {
182       complex<double>* const intermediate_i = stack_->get<complex<double>>(size_block_);
183       const unsigned int carsph_index = basisinfo_[0]->angular_number() * ANG_HRR_END + basisinfo_[1]->angular_number();
184       carsphlist.carsphfunc_call(carsph_index, cont0_ * cont1_, intermediate_fi, intermediate_i);
185 
186       const unsigned int sort_index = basisinfo_[1]->angular_number() * ANG_HRR_END + basisinfo_[0]->angular_number();
187       sort.sortfunc_call(sort_index, data_final, intermediate_i, cont1_, cont0_, 1, swap01_);
188       stack_->release(size_block_, intermediate_i);
189     } else {
190       const unsigned int sort_index = basisinfo_[1]->angular_number() * ANG_HRR_END + basisinfo_[0]->angular_number();
191       sort.sortfunc_call(sort_index, data_final, intermediate_fi, cont1_, cont0_, 1, swap01_);
192     }
193 
194     stack_->release(size_block_, intermediate_fi);
195     stack_->release(size_block_, intermediate_c);
196   } // end loop over multipoles
197 
198   // release resources
199   stack_->release(size_alloc_, intermediate_p);
200 }
201