1 //
2 // File: MPQC_IntegralEvaluator4_Impl.cc
3 // Symbol: MPQC.IntegralEvaluator4-v0.2
4 // Symbol Type: class
5 // Babel Version: 0.10.2
6 // Description: Server-side implementation for MPQC.IntegralEvaluator4
7 //
8 // WARNING: Automatically generated; only changes within splicers preserved
9 //
10 // babel-version = 0.10.2
11 //
12 #include "MPQC_IntegralEvaluator4_Impl.hh"
13
14 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4._includes)
15 #include <iostream>
16 #include <sstream>
17 #include <util/class/scexception.h>
18 #include <ccaiter.h>
19
20 using namespace std;
21 using namespace Chemistry::QC::GaussianBasis;
22
23 Ref<GaussianBasisSet> basis_cca_to_sc(Molecular&);
24 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4._includes)
25
26 // user-defined constructor.
_ctor()27 void MPQC::IntegralEvaluator4_impl::_ctor() {
28 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4._ctor)
29 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4._ctor)
30 }
31
32 // user-defined destructor.
_dtor()33 void MPQC::IntegralEvaluator4_impl::_dtor() {
34 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4._dtor)
35 #ifndef INTV3_ORDER
36 if( package_ == "intv3") delete [] temp_buffer_;
37 #endif
38 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4._dtor)
39 }
40
41 // static class initializer.
_load()42 void MPQC::IntegralEvaluator4_impl::_load() {
43 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4._load)
44 // guaranteed to be called at most once before any other method in this class
45 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4._load)
46 }
47
48 // user-defined static methods: (none)
49
50 // user-defined non-static methods:
51 /**
52 * Method: set_integral_package[]
53 */
54 void
set_integral_package(const::std::string & label)55 MPQC::IntegralEvaluator4_impl::set_integral_package (
56 /* in */ const ::std::string& label )
57 throw ()
58 {
59 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4.set_integral_package)
60 package_ = label;
61 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4.set_integral_package)
62 }
63
64 /**
65 * Initialize the evaluator.
66 * @param bs1 Molecular basis on center 1.
67 * @param bs2 Molecular basis on center 2.
68 * @param bs3 Molecular basis on center 3.
69 * @param bs4 Molecular basis on center 4.
70 * @param label String specifying integral type.
71 * @param max_deriv Max derivative to compute.
72 */
73 void
initialize(::Chemistry::QC::GaussianBasis::Molecular bs1,::Chemistry::QC::GaussianBasis::Molecular bs2,::Chemistry::QC::GaussianBasis::Molecular bs3,::Chemistry::QC::GaussianBasis::Molecular bs4,const::std::string & label,int64_t max_deriv)74 MPQC::IntegralEvaluator4_impl::initialize (
75 /* in */ ::Chemistry::QC::GaussianBasis::Molecular bs1,
76 /* in */ ::Chemistry::QC::GaussianBasis::Molecular bs2,
77 /* in */ ::Chemistry::QC::GaussianBasis::Molecular bs3,
78 /* in */ ::Chemistry::QC::GaussianBasis::Molecular bs4,
79 /* in */ const ::std::string& label,
80 /* in */ int64_t max_deriv )
81 throw ()
82 {
83 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4.initialize)
84
85 bufn_ = 0;
86
87 evaluator_label_ = label;
88 int deriv_level = max_deriv;
89
90 bs1_ = basis_cca_to_sc( bs1 );
91 bs2_ = basis_cca_to_sc( bs2 );
92 bs3_ = basis_cca_to_sc( bs3 );
93 bs4_ = basis_cca_to_sc( bs4 );
94
95 max_nshell4_ = bs1_->max_ncartesian_in_shell();
96 max_nshell4_ *= bs2_->max_ncartesian_in_shell();
97 max_nshell4_ *= bs3_->max_ncartesian_in_shell();
98 max_nshell4_ *= bs4_->max_ncartesian_in_shell();
99
100 std::string is_deriv("");
101 if(max_deriv > 0) is_deriv = " derivative";
102 std::cout << " initializing " << package_ << " " << evaluator_label_
103 << is_deriv << " integral evaluator\n";
104 if ( package_ == "intv3" )
105 integral_ = new IntegralV3( bs1_ );
106 #ifdef HAVE_CINTS
107 else if ( package_ == "cints" )
108 integral_ = new IntegralCints( bs1_ );
109 #endif
110 else
111 throw InputError("bad integral package name",
112 __FILE__,__LINE__);
113
114 // a proper solution is required here
115 integral_->set_storage(200000000);
116
117 int error = 0;
118 if(evaluator_label_ == "eri2")
119 switch( deriv_level ) {
120 case 0:
121 { eval_ = integral_->electron_repulsion(); break; }
122 case 1:
123 { deriv_eval_ = integral_->electron_repulsion_deriv();
124 break;
125 }
126 default:
127 ++error;
128 }
129
130 else if(evaluator_label_ == "grt")
131 switch( deriv_level ) {
132 case 0:
133 { eval_ = integral_->grt(); break; }
134 default:
135 ++error;
136 }
137
138 else
139 throw InputError("unsupported integral type",
140 __FILE__,__LINE__);
141
142 if( error )
143 throw InputError("derivative level not supported",
144 __FILE__,__LINE__);
145
146 if( eval_.nonnull() ) {
147 int_type_ = two_body;
148 sc_buffer_ = eval_->buffer();
149 }
150 else if( deriv_eval_.nonnull() ) {
151 int_type_ = two_body_deriv;
152 sc_buffer_ = deriv_eval_->buffer();
153 }
154 else
155 throw ProgrammingError("bad integral evaluator pointer",
156 __FILE__,__LINE__);
157 if( !sc_buffer_ )
158 throw ProgrammingError("buffer not assigned",
159 __FILE__,__LINE__);
160 // get a non-const pointer we can write to
161 buf_ = const_cast<double*>( sc_buffer_ );
162
163 if ( package_ == "intv3" ) {
164 #ifdef INTV3_ORDER
165 std::cout << " using intv3 ordering" << std::endl;
166 #else
167 initialize_reorder_intv3();
168 #endif
169 }
170
171 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4.initialize)
172 }
173
174 /**
175 * Get the buffer pointer.
176 * @return Buffer pointer.
177 */
178 void*
get_buffer()179 MPQC::IntegralEvaluator4_impl::get_buffer ()
180 throw ()
181
182 {
183 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4.get_buffer)
184 return const_cast<double*>( sc_buffer_ );
185 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4.get_buffer)
186 }
187
188 /**
189 * Compute a shell quartet of integrals.
190 * @param shellnum1 Gaussian shell number 1.
191 * @param shellnum2 Gaussian shell number 2.
192 * @param shellnum3 Gaussian shell number 3.
193 * @param shellnum4 Gaussian shell number 4.
194 * @param deriv_level Derivative level.
195 * @param deriv_ctr Derivative center descriptor.
196 */
197 void
compute(int64_t shellnum1,int64_t shellnum2,int64_t shellnum3,int64_t shellnum4,int64_t deriv_level,::Chemistry::QC::GaussianBasis::DerivCenters deriv_ctr)198 MPQC::IntegralEvaluator4_impl::compute (
199 /* in */ int64_t shellnum1,
200 /* in */ int64_t shellnum2,
201 /* in */ int64_t shellnum3,
202 /* in */ int64_t shellnum4,
203 /* in */ int64_t deriv_level,
204 /* in */ ::Chemistry::QC::GaussianBasis::DerivCenters deriv_ctr )
205 throw ()
206 {
207 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4.compute)
208
209 if( int_type_ == two_body ) {
210 eval_->compute_shell( shellnum1, shellnum2,
211 shellnum3, shellnum4);
212 }
213
214 else if( int_type_ == two_body_deriv ) {
215 sc::DerivCenters dc;
216
217 if(deriv_ctr.has_omitted_center() && deriv_ctr.omitted_center() == 0 )
218 dc.add_omitted(0,deriv_ctr.atom(0));
219 else
220 dc.add_center(0,deriv_ctr.atom(0));
221
222 if(deriv_ctr.has_omitted_center() && deriv_ctr.omitted_center() == 1 )
223 dc.add_omitted(1,deriv_ctr.atom(1));
224 else
225 dc.add_center(1,deriv_ctr.atom(1));
226
227 if(deriv_ctr.has_omitted_center() && deriv_ctr.omitted_center() == 2 )
228 dc.add_omitted(2,deriv_ctr.atom(2));
229 else
230 dc.add_center(2,deriv_ctr.atom(2));
231
232 if(deriv_ctr.has_omitted_center() && deriv_ctr.omitted_center() == 3 )
233 dc.add_omitted(3,deriv_ctr.atom(3));
234 else
235 dc.add_center(3,deriv_ctr.atom(3));
236
237 deriv_eval_->compute_shell( shellnum1, shellnum2,
238 shellnum3, shellnum4, dc );
239 }
240 else
241 throw ProgrammingError("bad evaluator type",
242 __FILE__,__LINE__);
243
244 #ifndef INTV3_ORDER
245 if( package_ == "intv3" )
246 reorder_intv3( shellnum1, shellnum2, shellnum3, shellnum4 );
247 #endif
248
249 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4.compute)
250 }
251
252 /**
253 * Compute a shell quartet of integrals and return as a borrowed
254 * sidl array.
255 * @param shellnum1 Gaussian shell number 1.
256 * @param shellnum2 Gaussian shell number 2.
257 * @param shellnum3 Guassian shell number 3.
258 * @param shellnum4 Gaussian shell number 4.
259 * @param deriv_level Derivative level.
260 * @param deriv_ctr Derivative center descriptor.
261 * @return Borrowed sidl array buffer.
262 */
263 ::sidl::array<double>
compute_array(int64_t shellnum1,int64_t shellnum2,int64_t shellnum3,int64_t shellnum4,int64_t deriv_level,::Chemistry::QC::GaussianBasis::DerivCenters deriv_ctr)264 MPQC::IntegralEvaluator4_impl::compute_array (
265 /* in */ int64_t shellnum1,
266 /* in */ int64_t shellnum2,
267 /* in */ int64_t shellnum3,
268 /* in */ int64_t shellnum4,
269 /* in */ int64_t deriv_level,
270 /* in */ ::Chemistry::QC::GaussianBasis::DerivCenters deriv_ctr )
271 throw ()
272 {
273 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4.compute_array)
274
275 compute( shellnum1, shellnum2, shellnum3, shellnum4, deriv_level, deriv_ctr );
276
277 // this creates a proxy SIDL array
278 int lower[1] = {0};
279 int upper[1]; upper[0] = max_nshell4_-1;
280 int stride[1] = {1};
281 sidl_buffer_.borrow( const_cast<double*>(sc_buffer_),
282 1, lower, upper, stride);
283
284 return sidl_buffer_;
285
286 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4.compute_array)
287 }
288
289
290 // DO-NOT-DELETE splicer.begin(MPQC.IntegralEvaluator4._misc)
291
292 void
initialize_reorder_intv3()293 MPQC::IntegralEvaluator4_impl::initialize_reorder_intv3()
294 {
295 if( int_type_ == two_body )
296 temp_buffer_ = new double[max_nshell4_];
297 else if( int_type_ == two_body_deriv )
298 temp_buffer_ = new double[max_nshell4_*3];
299
300 int max12 = max( bs1_->max_angular_momentum(), bs2_->max_angular_momentum() );
301 int max34 = max( bs3_->max_angular_momentum(), bs4_->max_angular_momentum() );
302 int maxam = max( max12, max34 );
303
304 reorder_ = new int*[maxam+1];
305 reorder_[0] = new int[1];
306 reorder_[0][0] = 0;
307 if(maxam==0) return;
308
309 for( int i=1; i<=maxam; ++i) {
310
311 sc::CartesianIter *v3iter = integral_->new_cartesian_iter(i);
312 MPQC::CartesianIterCCA iter(i);
313 MPQC::CartesianIterCCA *ccaiter = &iter;
314 ccaiter->start();
315 int ncf = ccaiter->n();
316
317 reorder_[i] = new int[ncf];
318 v3iter->start();
319 for( int j=0; j<ncf; ++j) {
320 ccaiter->start();
321 for( int k=0; k<ncf; ++k) {
322 if( v3iter->a() == ccaiter->a() &&
323 v3iter->b() == ccaiter->b() &&
324 v3iter->c() == ccaiter->c() ) {
325 reorder_[i][j] = k;
326 k=ncf; //break k loop
327 }
328 else ccaiter->next();
329 }
330 v3iter->next();
331 }
332 }
333
334 }
335
336 void
reorder_intv3(int64_t shellnum1,int64_t shellnum2,int64_t shellnum3,int64_t shellnum4)337 MPQC::IntegralEvaluator4_impl::reorder_intv3(int64_t shellnum1,
338 int64_t shellnum2,
339 int64_t shellnum3,
340 int64_t shellnum4)
341 {
342
343 sc::GaussianShell* s1 = &( bs1_->shell(shellnum1) );
344 sc::GaussianShell* s2 = &( bs2_->shell(shellnum2) );
345 sc::GaussianShell* s3 = &( bs3_->shell(shellnum3) );
346 sc::GaussianShell* s4 = &( bs4_->shell(shellnum4) );
347 int nc1 = s1->ncontraction();
348 int nc2 = s2->ncontraction();
349 int nc3 = s3->ncontraction();
350 int nc4 = s4->ncontraction();
351
352 int reorder_needed=0;
353 for (int i=0; i<nc1; ++i) {
354 if( s1->am(i) == 1) reorder_needed=1;
355 else if( s1->am(i) > 1 && s1->is_cartesian(i) ) reorder_needed=1;
356 }
357 if (!reorder_needed)
358 for (int i=0; i<nc2; ++i) {
359 if( s2->am(i) == 1) reorder_needed=1;
360 else if( s2->am(i) > 1 && s2->is_cartesian(i) ) reorder_needed=1;
361 }
362 if (!reorder_needed)
363 for (int i=0; i<nc3; ++i) {
364 if( s3->am(i) == 1) reorder_needed=1;
365 else if( s3->am(i) > 1 && s3->is_cartesian(i) ) reorder_needed=1;
366 }
367 if (!reorder_needed)
368 for (int i=0; i<nc4; ++i) {
369 if( s4->am(i) == 1) reorder_needed=1;
370 else if( s4->am(i) > 1 && s4->is_cartesian(i) ) reorder_needed=1;
371 }
372 if( !reorder_needed ) return;
373
374 // copy buffer into temp space
375 int nfunc = s1->nfunction() * s2->nfunction() *
376 s3->nfunction() * s4->nfunction();
377 if( int_type_ == two_body_deriv )
378 for( int i=0; i<nfunc*3; ++i)
379 temp_buffer_[i] = sc_buffer_[i];
380 else
381 for( int i=0; i<nfunc; ++i)
382 temp_buffer_[i] = sc_buffer_[i];
383
384 // a derivative buffer is composed of 3 "quartets"
385 int deriv_offset;
386 if( int_type_ == two_body )
387 reorder_quartet( s1, s2, s3, s4, nc1, nc2, nc3, nc4, 0 );
388 else if( int_type_ == two_body_deriv )
389 for( int i=0; i<3; ++i) {
390 deriv_offset = i*nfunc;
391 reorder_quartet( s1, s2, s3, s4, nc1, nc2, nc3, nc4, deriv_offset );
392 }
393
394 }
395
396
397 void
reorder_quartet(sc::GaussianShell * s1,sc::GaussianShell * s2,sc::GaussianShell * s3,sc::GaussianShell * s4,int nc1,int nc2,int nc3,int nc4,int deriv_offset)398 MPQC::IntegralEvaluator4_impl::reorder_quartet( sc::GaussianShell* s1, sc::GaussianShell* s2,
399 sc::GaussianShell* s3, sc::GaussianShell* s4,
400 int nc1, int nc2, int nc3, int nc4,
401 int deriv_offset )
402 {
403
404 int index=deriv_offset, con2_offset=0, con3_offset=0, con4_offset=0, con_offset,
405 local2_offset, local3_offset, local4_offset,
406 c1_base, c2_base, c3_base, c4_base;
407
408 int temp;
409 for( int c4=0; c4<nc4; ++c4 )
410 con4_offset += s4->nfunction(c4);
411
412 temp = 0;
413 con3_offset = con4_offset;
414 for( int c3=0; c3<nc3; ++c3 )
415 temp += s3->nfunction(c3);
416 con3_offset *= temp;
417
418 temp = 0;
419 con2_offset = con3_offset;
420 for( int c2=0; c2<nc2; ++c2 )
421 temp += s2->nfunction(c2);
422 con2_offset *= temp;
423
424 int s1_is_cart, s2_is_cart, s3_is_cart, s4_is_cart,
425 s1_nfunc, s2_nfunc, s3_nfunc, s4_nfunc;
426
427 for( int c1=0; c1<nc1; ++c1 ) {
428
429 c1_base = index;
430 s1_is_cart = s1->is_cartesian(c1);
431 s1_nfunc = s1->nfunction(c1);
432
433 for( int fc1=0; fc1<s1_nfunc; ++fc1 ) {
434
435 if( s1_is_cart )
436 c2_base = c1_base + reorder_[s1->am(c1)][fc1] * con2_offset;
437 else
438 c2_base = c1_base + fc1 * con2_offset;
439
440 local2_offset = 0;
441 for( int c2=0; c2<nc2; ++c2 ) {
442
443 if( c2>0 ) local2_offset += s2->nfunction(c2-1);
444 s2_is_cart = s2->is_cartesian(c2);
445 s2_nfunc = s2->nfunction(c2);
446
447 for( int fc2=0; fc2<s2_nfunc; ++fc2 ) {
448
449 if( s2_is_cart )
450 c3_base = c2_base + (local2_offset + reorder_[s2->am(c2)][fc2])
451 * con3_offset;
452 else
453 c3_base = c2_base + (local2_offset + fc2) * con3_offset;
454
455 local3_offset = 0;
456 for( int c3=0; c3<nc3; ++c3 ) {
457
458 if( c3>0 ) local3_offset += s3->nfunction(c3-1);
459 s3_is_cart = s3->is_cartesian(c3);
460 s3_nfunc = s3->nfunction(c3);
461
462 for( int fc3=0; fc3<s3_nfunc; ++fc3 ) {
463
464 if( s3_is_cart )
465 c4_base = c3_base + (local3_offset + reorder_[s3->am(c3)][fc3])
466 * con4_offset;
467 else
468 c4_base = c3_base + (local3_offset + fc3) * con4_offset;
469
470 local4_offset = 0;
471 for( int c4=0; c4<nc4; ++c4 ) {
472
473 if( c4>0 ) local4_offset += s4->nfunction(c4-1);
474 s4_is_cart = s4->is_cartesian(c4);
475 s4_nfunc = s4->nfunction(c4);
476
477 if( s4_is_cart )
478 for( int fc4=0; fc4<s4_nfunc; ++fc4 ) {
479 buf_[ c4_base + local4_offset + reorder_[s4->am(c4)][fc4] ]
480 = temp_buffer_[index];
481 ++index;
482 }
483 else
484 for( int fc4=0; fc4<s4_nfunc; ++fc4 ) {
485 buf_[ c4_base + local4_offset + fc4 ] = temp_buffer_[index];
486 ++index;
487 }
488 }
489 }
490 }
491 }
492 }
493 }
494 }
495
496 }
497
498 // DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4._misc)
499
500