source: ThirdParty/mpqc_open/src/lib/chemistry/cca/MPQC_IntegralEvaluator4_Impl.cc@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 14.5 KB
Line 
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
20using namespace std;
21using namespace Chemistry::QC::GaussianBasis;
22
23Ref<GaussianBasisSet> basis_cca_to_sc(Molecular&);
24// DO-NOT-DELETE splicer.end(MPQC.IntegralEvaluator4._includes)
25
26// user-defined constructor.
27void 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.
33void 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.
42void 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 */
54void
55MPQC::IntegralEvaluator4_impl::set_integral_package (
56 /* in */ const ::std::string& label )
57throw ()
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 */
73void
74MPQC::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 )
81throw ()
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 */
178void*
179MPQC::IntegralEvaluator4_impl::get_buffer ()
180throw ()
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 */
197void
198MPQC::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 )
205throw ()
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>
264MPQC::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 )
271throw ()
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
292void
293MPQC::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
336void
337MPQC::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
397void
398MPQC::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
Note: See TracBrowser for help on using the repository browser.