1 //
2 // obint.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
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 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <util/misc/formio.h>
33
34 #include <math/scmat/block.h>
35 #include <math/scmat/blkiter.h>
36 #include <math/scmat/offset.h>
37
38 #include <chemistry/qc/basis/obint.h>
39 #include <chemistry/qc/basis/basis.h>
40
41 ///////////////////////////////////////////////////////////////////////
42
43 void
44 EfieldDotVectorData::set_position(double*p)
45 {
46 position[0] = p[0];
47 position[1] = p[1];
48 position[2] = p[2];
49 }
50
51 void
52 EfieldDotVectorData::set_vector(double*v)
53 {
54 vector[0] = v[0];
55 vector[1] = v[1];
56 vector[2] = v[2];
57 }
58
59 ///////////////////////////////////////////////////////////////////////
60
61 void
62 DipoleData::set_origin(double*o)
63 {
64 origin[0] = o[0];
65 origin[1] = o[1];
66 origin[2] = o[2];
67 }
68
69 ///////////////////////////////////////////////////////////////////////
70
71 PointChargeData::PointChargeData(int ncharges,
72 const double *const*positions,
73 const double *charges,
74 int copy_data)
75 {
76 ncharges_ = ncharges;
77 if (copy_data) {
78 alloced_positions_ = new double*[ncharges];
79 alloced_charges_ = new double[ncharges];
80 memcpy(alloced_charges_, charges, sizeof(double)*ncharges);
81 double *tmp = new double[ncharges*3];
82 for (int i=0; i<ncharges; i++) {
83 alloced_positions_[i] = tmp;
84 for (int j=0; j<3; j++) {
85 *tmp++ = positions[i][j];
86 }
87 }
88 positions_ = alloced_positions_;
89 charges_ = alloced_charges_;
90 }
91 else {
92 charges_ = charges;
93 alloced_charges_ = 0;
94 alloced_positions_ = 0;
95 charges_ = charges;
96 positions_ = positions;
97 }
98 }
99
100 PointChargeData::~PointChargeData()
101 {
102 if (alloced_positions_) {
103 delete[] alloced_positions_[0];
104 delete[] alloced_positions_;
105 }
106 delete[] alloced_charges_;
107 }
108
109 ///////////////////////////////////////////////////////////////////////
110
111 OneBodyInt::OneBodyInt(Integral* integral,
112 const Ref<GaussianBasisSet>&bs1,
113 const Ref<GaussianBasisSet>&bs2) :
114 integral_(integral),
115 bs1_(bs1), bs2_(bs2)
116 {
117 buffer_ = 0;
118 }
119
120 OneBodyInt::~OneBodyInt()
121 {
122 }
123
124 int
125 OneBodyInt::nbasis() const
126 {
127 return bs1_->nbasis();
128 }
129
130 int
131 OneBodyInt::nbasis1() const
132 {
133 return bs1_->nbasis();
134 }
135
136 int
137 OneBodyInt::nbasis2() const
138 {
139 return bs2_->nbasis();
140 }
141
142 int
143 OneBodyInt::nshell() const
144 {
145 return bs1_->nshell();
146 }
147
148 int
149 OneBodyInt::nshell1() const
150 {
151 return bs1_->nshell();
152 }
153
154 int
155 OneBodyInt::nshell2() const
156 {
157 return bs2_->nshell();
158 }
159
160 Ref<GaussianBasisSet>
161 OneBodyInt::basis()
162 {
163 return bs1_;
164 }
165
166 Ref<GaussianBasisSet>
167 OneBodyInt::basis1()
168 3 0 2 0 0 0 0 0 0 0 {
169 return bs1_;
170 }
171
172 Ref<GaussianBasisSet>
173 OneBodyInt::basis2()
174 1 0 0 0 0 0 0 0 0 0 {
175 return bs2_;
176 }
177
178 const double *
179 OneBodyInt::buffer() const
180 {
181 return buffer_;
182 }
183
184 void
185 OneBodyInt::reinitialize()
186 {
187 }
188
189 ///////////////////////////////////////////////////////////////////////
190
191 ShellPairIter::ShellPairIter()
192 {
193 }
194
195 ShellPairIter::~ShellPairIter()
196 {
197 }
198
199 void
200 ShellPairIter::init(const double * b, int ishell, int jshell,
201 int fi, int fj, int ni, int nj,
202 int red, double scl)
203 1 0 0 0 0 0 0 0 0 0 {
204 1 0 0 0 0 1 0 0 0 0 e12 = ((ishell==jshell) && red);
205
206 ioffset=fi;
207 joffset=fj;
208
209 iend=ni;
210 jend=nj;
211
212 1 0 0 0 0 0 0 0 0 0 buf=b;
213 scale_=scl;
214 }
215
216 ///////////////////////////////////////////////////////////////////////
217
218 OneBodyIntIter::OneBodyIntIter()
219 {
220 }
221
222 OneBodyIntIter::OneBodyIntIter(const Ref<OneBodyInt>& o) :
223 obi(o)
224 {
225 }
226
227 OneBodyIntIter::~OneBodyIntIter()
228 {
229 }
230
231 void
232 OneBodyIntIter::start(int ist, int jst, int ien, int jen)
233 4 0 3 0 1 0 0 0 0 0 {
234 istart=ist;
235 jstart=jst;
236 iend=ien;
237 1 0 0 0 0 0 0 0 0 0 jend=jen;
238
239 icur=istart;
240 jcur=jstart;
241
242 if (!iend) {
243 iend=obi->nshell1();
244 jend=obi->nshell2();
245 }
246
247 1 0 1 0 0 0 0 0 0 0 ij = (icur*(icur+1)>>1) + jcur;
248 }
249
250 static inline int
251 min(int i, int j)
252 0 0 1 0 1 0 0 0 0 0 {
253 return (i<j) ? i : j;
254 }
255
256 void
257 OneBodyIntIter::next()
258 1 0 0 0 0 1 0 0 0 0 {
259 1 0 0 0 0 0 0 0 0 0 int jlast = (redund) ? min(icur,jend-1) : jend-1;
260
261 if (jcur < jlast) {
262 jcur++;
263 ij++;
264 return;
265 }
266
267 jcur=jstart;
268 icur++;
269
270 1 0 0 0 0 0 0 0 0 0 ij = (icur*(icur+1)>>1) + jcur;
271 1 0 0 0 0 0 0 0 0 0 }
272
273 double
274 OneBodyIntIter::scale() const
275 {
276 return 1.0;
277 }
278
279 ShellPairIter&
280 OneBodyIntIter::current_pair()
281 {
282 obi->compute_shell(icur,jcur);
283 spi.init(obi->buffer(), icur, jcur,
284 obi->basis1()->shell_to_function(icur),
285 obi->basis2()->shell_to_function(jcur),
286 obi->basis1()->operator()(icur).nfunction(),
287 obi->basis2()->operator()(jcur).nfunction(),
288 redund, scale()
289 );
290
291 return spi;
292 0 0 0 0 1 0 0 0 0 0 }
293
294 ///////////////////////////////////////////////////////////////////////
295
296 OneBodyIntOp::OneBodyIntOp(const Ref<OneBodyInt>& it)
297 {
298 iter = new OneBodyIntIter(it);
299 }
300
301 OneBodyIntOp::OneBodyIntOp(const Ref<OneBodyIntIter>& it) :
302 iter(it)
303 {
304 }
305
306 OneBodyIntOp::~OneBodyIntOp()
307 {
308 }
309
310 void
311 OneBodyIntOp::process(SCMatrixBlockIter& b)
312 {
313 ExEnv::err() << node0 << indent
314 << "OneBodyIntOp::process: cannot handle generic case\n";
315 abort();
316 }
317
318 void
319 OneBodyIntOp::process_spec_rect(SCMatrixRectBlock* b)
320 {
321 Ref<GaussianBasisSet> bs1 = iter->one_body_int()->basis1();
322 Ref<GaussianBasisSet> bs2 = iter->one_body_int()->basis2();
323
324 // convert basis function indices into shell indices
325 int ishstart = bs1->function_to_shell(b->istart);
326 int jshstart = bs2->function_to_shell(b->jstart);
327
328 int b1end = b->iend;
329 int ishend = (b1end?bs1->function_to_shell(b1end-1) + 1 : 0);
330
331 int b2end = b->jend;
332 int jshend = (b2end?bs2->function_to_shell(b2end-1) + 1 : 0);
333
334 int njdata = b->jend - b->jstart;
335
336 iter->set_redundant(0);
337
338 for (iter->start(ishstart,jshstart,ishend,jshend);
339 iter->ready(); iter->next()) {
340 ShellPairIter& spi = iter->current_pair();
341
342 for (spi.start(); spi.ready(); spi.next()) {
343 int ifn = spi.i();
344 int jfn = spi.j();
345
346 if (ifn < b->istart || ifn >= b->iend ||
347 jfn < b->jstart || jfn >= b->jend)
348 continue;
349
350 int data_index = (ifn - b->istart)*njdata + jfn - b->jstart;
351 b->data[data_index] += spi.val();
352 }
353 }
354 }
355
356 void
357 OneBodyIntOp::process_spec_ltri(SCMatrixLTriBlock* b)
358 {
359 Ref<GaussianBasisSet> bs1 = iter->one_body_int()->basis1();
360
361 // convert basis function indices into shell indices
362 int fnstart = b->start;
363 int fnend = b->end;
364 int shstart = bs1->function_to_shell(fnstart);
365 int shend = (fnend?bs1->function_to_shell(fnend - 1) + 1 : 0);
366
367 iter->set_redundant(1);
368
369 // loop over all needed shells
370 for (iter->start(shstart,shstart,shend,shend); iter->ready(); iter->next()) {
371 ShellPairIter& spi = iter->current_pair();
372
373 // compute a set of shell integrals
374 for (spi.start(); spi.ready(); spi.next()) {
375 int ifn = spi.i();
376 int jfn = spi.j();
377
378 if (ifn < fnstart || ifn >= fnend)
379 continue;
380
381 int ioff = ifn-fnstart;
382 int joff = jfn-fnstart;
383
384 int data_index = i_offset(ioff)+joff;
385
386 b->data[data_index] += spi.val();
387 }
388 }
389 }
390
391 void
392 OneBodyIntOp::process_spec_rectsub(SCMatrixRectSubBlock* b)
393 {
394 Ref<GaussianBasisSet> bs1 = iter->one_body_int()->basis1();
395 Ref<GaussianBasisSet> bs2 = iter->one_body_int()->basis2();
396
397 // convert basis function indices into shell indices
398 int istart = b->istart;
399 int jstart = b->jstart;
400 int iend = b->iend;
401 int jend = b->jend;
402
403 int ishstart = bs1->function_to_shell(istart);
404 int jshstart = bs2->function_to_shell(jstart);
405
406 int ishend = (iend ? bs1->function_to_shell(iend-1) + 1 : 0);
407 int jshend = (jend ? bs2->function_to_shell(jend-1) + 1 : 0);
408
409 int njdata = b->istride;
410
411 iter->set_redundant(0);
412
413 for (iter->start(ishstart,jshstart,ishend,jshend);
414 iter->ready(); iter->next()) {
415 ShellPairIter& spi = iter->current_pair();
416
417 for (spi.start(); spi.ready(); spi.next()) {
418 int ifn = spi.i();
419 int jfn = spi.j();
420
421 if (ifn < istart || ifn >= iend || jfn < jstart || jfn >= jend)
422 continue;
423
424 int data_index = ifn*njdata + jfn;
425 b->data[data_index] += spi.val();
426 }
427 }
428 }
429
430 void
431 OneBodyIntOp::process_spec_ltrisub(SCMatrixLTriSubBlock* b)
432 1 0 0 0 0 0 0 0 0 0 {
433 Ref<GaussianBasisSet> bs1 = iter->one_body_int()->basis1();
434
435 // convert basis function indices into shell indices
436 int istart = b->istart;
437 1 0 1 0 0 0 0 0 0 0 int iend = b->iend;
438
439 int jstart = b->jstart;
440 0 0 0 0 1 0 0 0 0 0 int jend = b->jend;
441
442 int ishstart = bs1->function_to_shell(istart);
443 int jshstart = bs1->function_to_shell(jstart);
444
445 int ishend = (iend ? bs1->function_to_shell(iend-1) + 1 : 0);
446 1 0 0 0 1 0 0 0 0 0 int jshend = (jend ? bs1->function_to_shell(jend-1) + 1 : 0);
447
448 iter->set_redundant(1);
449
450 // loop over all needed shells
451 for (iter->start(ishstart,jshstart,ishend,jshend);
452 iter->ready(); iter->next()) {
453 ShellPairIter& spi = iter->current_pair();
454
455 // compute a set of shell integrals
456 for (spi.start(); spi.ready(); spi.next()) {
457 int ifn = spi.i();
458 int jfn = spi.j();
459
460 0 0 2 1 1 1 0 0 0 0 if (ifn < istart || ifn >= iend || jfn < jstart || jfn >= jend)
461 continue;
462
463 3 0 0 0 1 0 0 0 0 0 int data_index = i_offset(ifn)+jfn;
464 0 0 0 1 0 0 0 0 0 0 b->data[data_index] += spi.val();
465 }
466 }
467 }
468
469 int
470 OneBodyIntOp::has_side_effects()
471 {
472 return 1;
473 }
474
475 ///////////////////////////////////////////////////////////////////////
476
477 OneBody3IntOp::OneBody3IntOp(const Ref<OneBodyInt>& it)
478 {
479 iter = new OneBodyIntIter(it);
480 }
481
482 OneBody3IntOp::OneBody3IntOp(const Ref<OneBodyIntIter>& it) :
483 iter(it)
484 {
485 }
486
487 OneBody3IntOp::~OneBody3IntOp()
488 {
489 }
490
491 void
492 OneBody3IntOp::process(SCMatrixBlockIter&,
493 SCMatrixBlockIter&,
494 SCMatrixBlockIter&)
495 {
496 ExEnv::err() << node0 << indent
497 << "OneBody3IntOp::process(SCMatrixBlockIter&): "
498 << "cannot handle generic case\n";
499 abort();
500 }
501
502 void
503 OneBody3IntOp::process_spec_rect(SCMatrixRectBlock* a,
504 SCMatrixRectBlock* b,
505 SCMatrixRectBlock* c)
506 {
507 Ref<GaussianBasisSet> bs1 = iter->one_body_int()->basis1();
508 Ref<GaussianBasisSet> bs2 = iter->one_body_int()->basis2();
509
510 // convert basis function indices into shell indices
511 int ishstart = bs1->function_to_shell(b->istart);
512 int jshstart = bs2->function_to_shell(b->jstart);
513
514 int ishend = bs1->function_to_shell(b->iend);
515 int jshend = bs2->function_to_shell(b->jend);
516
517 iter->set_redundant(0);
518
519 for (iter->start(ishstart,jshstart,ishend,jshend);
520 iter->ready(); iter->next()) {
521 // compute a set of shell integrals
522 ShellPairIter& spi = iter->current_pair();
523
524 for (spi.start(); spi.ready(); spi.next()) {
525 int ifn = spi.i();
526 int jfn = spi.j();
527
528 if (ifn < b->istart || ifn >= b->iend ||
529 jfn < b->jstart || jfn >= b->jend)
530 continue;
531
532 #if 0
533 for (int i=0; i<nish; i++,ifn++) {
534 if (ifn < b->istart || ifn >= b->iend) {
535 tmp += njsh * 3;
536 }
537 else {
538 int jfn = jfnsave;
539 int data_index = (ifn - b->istart)*njdata + jfn - b->jstart;
540 for (int j=0; j<njsh; j++,jfn++) {
541 if (jfn >= b->jstart && jfn < b->jend) {
542 a->data[data_index] += tmp[0] * scale;
543 b->data[data_index] += tmp[1] * scale;
544 c->data[data_index] += tmp[2] * scale;
545 data_index++;
546 }
547 tmp += 3;
548 }
549 }
550 }
551 #endif
552 }
553 }
554 }
555
556 void
557 OneBody3IntOp::process_spec_ltri(SCMatrixLTriBlock* a,
558 SCMatrixLTriBlock* b,
559 SCMatrixLTriBlock* c)
560 {
561 #if 0
562 Ref<GaussianBasisSet> bs1 = iter->one_body_int()->basis1();
563
564 // convert basis function indices into shell indices
565 int fnstart = b->start;
566 int fnend = b->end;
567 int shstart = bs1->function_to_shell(fnstart);
568 int shend = (fnend?bs1->function_to_shell(fnend - 1) + 1 : 0);
569
570 // loop over all needed shells
571 iter->reset(shstart, shend, 0, 0);
572 for (iter->start_ltri(); iter->ready_ltri(); iter->next_ltri()) {
573 int ish=iter->ishell();
574 int jsh=iter->jshell();
575
576 int nish = bs1->operator[](ish).nfunction();
577 int njsh = bs1->operator[](jsh).nfunction();
578
579 double scale = iter->scale();
580
581 // compute a set of shell integrals
582 compute_shell(ish,jsh,buffer_);
583
584 // take the integrals from buffer and put them into the LTri block
585 double*tmp = buffer_;
586 int ifn = bs1->shell_to_function(ish);
587 int jfnsave = bs1->shell_to_function(jsh);
588 for (int i=0; i<nish; i++,ifn++) {
589 // skip over basis functions that are not needed
590 if (ifn < fnstart || ifn >= fnend) {
591 tmp += njsh * 3;
592 }
593 else {
594 int jfn = jfnsave;
595 int irelfn = ifn - fnstart;
596 int data_index = ((irelfn+1)*irelfn>>1) + jfn - fnstart;
597 for (int j=0; j<njsh; j++,jfn++) {
598 // skip over basis functions that are not needed
599 if (jfn <= ifn && jfn >= fnstart) {
600 a->data[data_index] += tmp[0] * scale;
601 b->data[data_index] += tmp[1] * scale;
602 c->data[data_index] += tmp[2] * scale;
603 data_index++;
604 }
605 tmp += 3;
606 }
607 }
608 }
609 }
610 #endif
611 }
612
613 int
614 OneBody3IntOp::has_side_effects()
615 {
616 return 1;
617 }
618
619 int
620 OneBody3IntOp::has_side_effects_in_arg1()
621 {
622 return 1;
623 }
624
625 int
626 OneBody3IntOp::has_side_effects_in_arg2()
627 {
628 return 1;
629 }
630
631 ///////////////////////////////////////////////////////////////////////
632
633 OneBodyDerivInt::OneBodyDerivInt(Integral *integral,
634 const Ref<GaussianBasisSet>&b) :
635 integral_(integral),
636 bs1(b), bs2(b)
637 {
638 // allocate a buffer
639 int biggest_shell = b->max_nfunction_in_shell();
640 biggest_shell *= biggest_shell * 3;
641
642 if (biggest_shell) {
643 buffer_ = new double[biggest_shell];
644 } else {
645 buffer_ = 0;
646 }
647 }
648
649 OneBodyDerivInt::OneBodyDerivInt(Integral *integral,
650 const Ref<GaussianBasisSet>&b1,
651 const Ref<GaussianBasisSet>&b2) :
652 integral_(integral),
653 bs1(b1), bs2(b2)
654 {
655 buffer_ = 0;
656 }
657
658 OneBodyDerivInt::~OneBodyDerivInt()
659 {
660 }
661
662 int
663 OneBodyDerivInt::nbasis() const
664 {
665 return bs1->nbasis();
666 }
667
668 int
669 OneBodyDerivInt::nbasis1() const
670 {
671 return bs1->nbasis();
672 }
673
674 int
675 OneBodyDerivInt::nbasis2() const
676 {
677 return bs2->nbasis();
678 }
679
680 int
681 OneBodyDerivInt::nshell() const
682 {
683 return bs1->nshell();
684 }
685
686 int
687 OneBodyDerivInt::nshell1() const
688 {
689 return bs1->nshell();
690 }
691
692 int
693 OneBodyDerivInt::nshell2() const
694 {
695 return bs2->nshell();
696 }
697
698 Ref<GaussianBasisSet>
699 OneBodyDerivInt::basis()
700 {
701 return bs1;
702 }
703
704 Ref<GaussianBasisSet>
705 OneBodyDerivInt::basis1()
706 {
707 return bs1;
708 }
709
710 Ref<GaussianBasisSet>
711 OneBodyDerivInt::basis2()
712 {
713 return bs2;
714 }
715
716 const double *
717 OneBodyDerivInt::buffer() const
718 {
719 return buffer_;
720 }
721
722 /////////////////////////////////////////////////////////////////////////////
723
724 // Local Variables:
725 // mode: c++
726 // c-file-style: "ETS"
727 // End:
728