1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: complexmomentumbatch.cc
4 // Copyright (C) 2014 Toru Shiozaki
5 //
6 // Author: Ryan D. Reynolds <RyanDReynolds@u.northwestern.edu>
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // This program is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program 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 General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
23 //
24 
25 
26 #include <complex>
27 #include <src/integral/carsphlist.h>
28 #include <src/integral/compos/complexmomentumbatch.h>
29 
30 using namespace std;
31 using namespace bagel;
32 
33 const static CCarSphList carsphlist;
34 
35 
compute()36 void ComplexMomentumBatch::compute() {
37 
38   const CSortList sort_(spherical_);
39 
40   complex<double>* const intermediate_p = stack_->get<complex<double>>(size_block_*3);
41   fill_n(intermediate_p, size_block_*3, 0.0);
42   perform_VRR(intermediate_p);
43 
44   for (int i = 0; i != 3; ++i) {
45     complex<double>* const cdata = data_ + i*size_block_;
46     const complex<double>* const csource = intermediate_p + i*size_block_;
47     complex<double>* const intermediate_c = stack_->get<complex<double>>(cont0_ * cont1_ * asize_intermediate_);
48     fill_n(intermediate_c, cont0_ * cont1_ * asize_intermediate_, 0.0);
49     perform_contraction(asize_intermediate_, csource, prim0_, prim1_, intermediate_c,
50                         basisinfo_[0]->contractions(), basisinfo_[0]->contraction_ranges(), cont0_,
51                         basisinfo_[1]->contractions(), basisinfo_[1]->contraction_ranges(), cont1_);
52 
53     if (basisinfo_[0]->spherical() && basisinfo_[1]->spherical()) {
54       // transform both indices to spherical
55       complex<double>* const intermediate_i = stack_->get<complex<double>>(cont0_ * cont1_ * asize_final_);
56       fill_n(intermediate_i, cont0_ * cont1_ * asize_final_, 0.0);
57       const unsigned int carsph_index = basisinfo_[0]->angular_number() * ANG_HRR_END + basisinfo_[1]->angular_number();
58       const int nloops = cont0_ * cont1_;
59       carsphlist.carsphfunc_call(carsph_index, nloops, intermediate_c, intermediate_i);
60 
61       const unsigned int sort_index = basisinfo_[1]->angular_number() * ANG_HRR_END + basisinfo_[0]->angular_number();
62       sort_.sortfunc_call(sort_index, cdata, intermediate_i, cont1_, cont0_, 1, swap01_);
63       stack_->release(cont0_ * cont1_ * asize_final_, intermediate_i);
64     } else {
65       const unsigned int sort_index = basisinfo_[1]->angular_number() * ANG_HRR_END + basisinfo_[0]->angular_number();
66       sort_.sortfunc_call(sort_index, cdata, intermediate_c, cont1_, cont0_, 1, swap01_);
67     }
68 
69     stack_->release(cont0_ * cont1_ * asize_intermediate_, intermediate_c);
70   }
71   stack_->release(size_block_*3, intermediate_p);
72 
73 }
74 
75 
perform_VRR(complex<double> * intermediate)76 void ComplexMomentumBatch::perform_VRR(complex<double>* intermediate) {
77   const int amax2 = amax1_+1;
78   const int worksize = amax2;
79   complex<double>* worksx = stack_->get<complex<double>>(worksize*worksize);
80   complex<double>* worksy = stack_->get<complex<double>>(worksize*worksize);
81   complex<double>* worksz = stack_->get<complex<double>>(worksize*worksize);
82 
83   complex<double> imag (0.0, 1.0);
84   if (swap01_) imag *= -1.0;
85   const double Bx = basisinfo_[1]->position(0);
86   const double By = basisinfo_[1]->position(1);
87   const double Bz = basisinfo_[1]->position(2);
88   const double A_Bx = basisinfo_[1]->vector_potential(0);
89   const double A_By = basisinfo_[1]->vector_potential(1);
90   const double A_Bz = basisinfo_[1]->vector_potential(2);
91   const double field_x_half = 0.5 * magnetic_field_[0];
92   const double field_y_half = 0.5 * magnetic_field_[1];
93   const double field_z_half = 0.5 * magnetic_field_[2];
94 
95   for (int ii = 0; ii != prim0_*prim1_; ++ii) {
96     // Perform VRR
97     int offset_ii = ii*asize_intermediate_;
98     const double cop = 1.0 / xp_[ii];
99     const double cb = xb_[ii];
100     const complex<double> cxpa = P_[ii*3  ]-basisinfo_[0]->position(0);
101     const complex<double> cypa = P_[ii*3+1]-basisinfo_[0]->position(1);
102     const complex<double> czpa = P_[ii*3+2]-basisinfo_[0]->position(2);
103     complex<double>* current_data = &intermediate[offset_ii];
104 
105     // obtain S(0, 0)
106     worksx[0] = coeffsx_[ii];
107     worksy[0] = coeffsy_[ii];
108     worksz[0] = coeffsz_[ii];
109 
110     // obtain S(1, 0)
111     worksx[1] = cxpa*worksx[0];
112     worksy[1] = cypa*worksy[0];
113     worksz[1] = czpa*worksz[0];
114 
115     for (int i = 2; i != amax2; ++i) {
116       // obtain S(i, 0)
117       worksx[i] = cxpa*worksx[i-1]+0.5*(i-1)*cop*worksx[i-2];
118       worksy[i] = cypa*worksy[i-1]+0.5*(i-1)*cop*worksy[i-2];
119       worksz[i] = czpa*worksz[i-1]+0.5*(i-1)*cop*worksz[i-2];
120     }
121 
122     // peform HRR
123     for (int j = 1; j <= ang1_ + 1; ++j) {
124       for (int i = 0; i != amax2 - j; ++i) {
125         // obtain S(i, j)
126         worksx[j*amax2+i] = AB_[0]*worksx[(j-1)*amax2+i]+worksx[(j-1)*amax2+i+1];
127         worksy[j*amax2+i] = AB_[1]*worksy[(j-1)*amax2+i]+worksy[(j-1)*amax2+i+1];
128         worksz[j*amax2+i] = AB_[2]*worksz[(j-1)*amax2+i]+worksz[(j-1)*amax2+i+1];
129       }
130     }
131 
132     complex<double> Sx, Sy, Sz;    // operator = 1
133     complex<double> ox, oy, oz;    // operator = x
134     complex<double> dx, dy, dz;    // operator = d/dx
135 
136     // now we obtain the output
137     assert((ang0_+1)*(ang0_+2)*(ang1_+1)*(ang1_+2)/4 == asize_intermediate_);
138 
139     int cnt = 0;
140     for (int iz = 0; iz <= ang0_; ++iz) {
141       for (int iy = 0; iy <= ang0_-iz; ++iy) {
142         const int ix = ang0_-iy-iz;
143         if (ix >= 0) {
144           for (int jz = 0; jz <= ang1_; ++jz) {
145             for (int jy = 0; jy <= ang1_-jz; ++jy) {
146               const int jx = ang1_-jy-jz;
147               if (jx >= 0) {
148                 const double jxd = jx;
149                 const double jyd = jy;
150                 const double jzd = jz;
151                 complex<double> minus1x = 0.0;
152                 complex<double> minus1y = 0.0;
153                 complex<double> minus1z = 0.0;
154                 if (jx > 0) minus1x = jxd * worksx[ix+amax2*(jx-1)];
155                 if (jy > 0) minus1y = jyd * worksy[iy+amax2*(jy-1)];
156                 if (jz > 0) minus1z = jzd * worksz[iz+amax2*(jz-1)];
157 
158                 Sx = worksx[ix + amax2 * jx];
159                 Sy = worksy[iy + amax2 * jy];
160                 Sz = worksz[iz + amax2 * jz];
161 
162                 ox = worksx[ix + amax2 * (jx+1)] + Bx * worksx[ix + amax2 * jx];
163                 oy = worksy[iy + amax2 * (jy+1)] + By * worksy[iy + amax2 * jy];
164                 oz = worksz[iz + amax2 * (jz+1)] + Bz * worksz[iz + amax2 * jz];
165 
166                 dx = minus1x - 2.0 * cb * worksx[ix+amax2 * (jx+1)] - imag * A_Bx * worksx[ix + amax2 * jx];
167                 dy = minus1y - 2.0 * cb * worksy[iy+amax2 * (jy+1)] - imag * A_By * worksy[iy + amax2 * jy];
168                 dz = minus1z - 2.0 * cb * worksz[iz+amax2 * (jz+1)] - imag * A_Bz * worksz[iz + amax2 * jz];
169 
170                 current_data[cnt              ] = -imag * dx * Sy * Sz;
171                 current_data[cnt+size_block_  ] = -imag * Sx * dy * Sz;
172                 current_data[cnt+size_block_*2] = -imag * Sx * Sy * dz;
173 
174                 current_data[cnt              ] += field_y_half * Sx * Sy * oz;
175                 current_data[cnt+size_block_  ] += field_z_half * ox * Sy * Sz;
176                 current_data[cnt+size_block_*2] += field_x_half * Sx * oy * Sz;
177                 current_data[cnt              ] -= field_z_half * Sx * oy * Sz;
178                 current_data[cnt+size_block_  ] -= field_x_half * Sx * Sy * oz;
179                 current_data[cnt+size_block_*2] -= field_y_half * ox * Sy * Sz;
180                 ++cnt;
181               }
182             }
183           }
184         }
185       }
186     }
187     assert(cnt == asize_intermediate_);
188 
189   } // end of prim exponent loop
190 
191   stack_->release(worksize*worksize, worksz);
192   stack_->release(worksize*worksize, worksy);
193   stack_->release(worksize*worksize, worksx);
194 }
195 
get_P(const double coord1,const double coord2,const double exp1,const double exp2,const double one12,const int dim,const bool swap)196 complex<double> ComplexMomentumBatch::get_P(const double coord1, const double coord2, const double exp1, const double exp2, const double one12,
197                                                  const int dim, const bool swap) {
198   const double Areal = coord1*exp1;
199   const double Breal = coord2*exp2;
200   const double Aimag = basisinfo_[0]->vector_potential(dim);
201   const double Bimag = basisinfo_[1]->vector_potential(dim);
202   const double imag = 0.5 * (swap ? Bimag - Aimag : Aimag - Bimag );
203   const complex<double> num (Areal + Breal, imag);
204   return num * one12;
205 }
206 
207