source: ThirdParty/mpqc_open/src/lib/chemistry/qc/cints/cintstest.cc@ e9f307

Candidate_v1.7.0 stable
Last change on this file since e9f307 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: 44.3 KB
Line 
1//
2// cintstest.cc
3//
4// Copyright (C) 2001 Edward Valeev
5//
6// Author: Edward Valeev <edward.valeev@chemistry.gatech.edu>
7// Maintainer: EV
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#include <stdlib.h>
29#include <string.h>
30
31#include <util/misc/formio.h>
32#include <util/misc/regtime.h>
33#include <util/keyval/keyval.h>
34#include <util/group/message.h>
35#include <util/group/pregtime.h>
36#include <chemistry/qc/basis/integral.h>
37#include <chemistry/qc/intv3/int1e.h>
38#include <chemistry/qc/intv3/int2e.h>
39#include <chemistry/qc/intv3/intv3.h>
40#include <chemistry/qc/intv3/cartitv3.h>
41#include <chemistry/qc/cints/cints.h>
42#include <chemistry/qc/cints/int2e.h>
43#include <chemistry/qc/cints/cartit.h>
44
45
46using namespace std;
47using namespace sc;
48
49#define CINTS
50
51void compare_1e_cints_vs_v3(Ref<OneBodyInt>& obcints, Ref<OneBodyInt>& obv3);
52void compare_2e_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
53void compare_2e_puream_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
54void compare_2e_bufsum_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
55void compare_2e_unique_bufsum_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
56void print_grt_ints(Ref<TwoBodyInt>& tbcints);
57void compare_2e_permute(Ref<Integral>& cints);
58void test_int_shell_1e(const Ref<KeyVal>&, const Ref<Int1eV3> &int1ev3,
59 void (Int1eV3::*int_shell_1e)(int,int),
60 int permute);
61void test_3_center(const Ref<KeyVal>&, const Ref<Int2eV3> &);
62void test_4_center(const Ref<KeyVal>& keyval, const Ref<Int2eV3> &int2ev3);
63void test_4der_center(const Ref<KeyVal>&, const Ref<Int2eV3> &int2ev3);
64
65#define maxint 9
66
67void
68testint(const Ref<OneBodyInt>& in)
69{
70 if (in.null()) {
71 cout << "null integral generator" << endl;
72 abort();
73 }
74}
75
76void
77testint(const Ref<OneBodyDerivInt>& in)
78{
79 if (in.null()) {
80 cout << "null integral generator" << endl;
81 abort();
82 }
83}
84
85void
86testint(const Ref<TwoBodyInt>& in)
87{
88 if (in.null()) {
89 cout << "null integral generator" << endl;
90 abort();
91 }
92}
93
94void
95testint(const Ref<TwoBodyDerivInt>& in)
96{
97 if (in.null()) {
98 cout << "null integral generator" << endl;
99 abort();
100 }
101}
102
103/*void
104do_bounds_stats(const Ref<KeyVal>& keyval,
105 const Ref<Int2eV3> &int2ev3)
106{
107 int i,j;
108 int nshell = int2ev3->basis()->nshell();
109 int eps = -10;
110 int *nonzero = new int[nshell];
111 for (i=0; i<nshell; i++) {
112 if (i==0) nonzero[i] = 0;
113 else nonzero[i] = nonzero[i-1];
114 for (j=0; j<=i; j++) {
115 if (int2ev3->erep_4bound(i,j,-1,-1) > eps) {
116 nonzero[i]++;
117 }
118 }
119 }
120 for (i=0; i<nshell; i++) {
121 int natom = 1 + int2ev3->basis()->shell_to_center(i);
122 int npq = (i*(i+1))/2;
123 cout<<scprintf("nsh=%2d nat=%2d npq=%4d npq>eps=%4d npq>eps/nsh=%9.4f /nat=%9.4f",
124 i, natom, npq, nonzero[i], double(nonzero[i])/i,
125 double(nonzero[i])/natom)
126 << endl;
127 }
128 delete[] nonzero;
129}
130*/
131
132int main(int argc, char **argv)
133{
134 int ii, i,j,k,l,m,n;
135
136 Ref<MessageGrp> msg = MessageGrp::initial_messagegrp(argc,argv);
137 if (msg.null()) msg = new ProcMessageGrp();
138 MessageGrp::set_default_messagegrp(msg);
139
140 Ref<RegionTimer> tim = new ParallelRegionTimer(msg,"cintstest", 1, 1);
141
142 char *infile = new char[strlen(SRCDIR)+strlen("/cintstest.in")+1];
143 sprintf(infile,SRCDIR "/cintstest.in");
144 if (argc == 2) {
145 infile = argv[1];
146 }
147
148 Ref<KeyVal> pkv(new ParsedKeyVal(infile));
149 Ref<KeyVal> tkeyval(new PrefixKeyVal(":test", pkv));
150
151 Ref<GaussianBasisSet> basis = require_dynamic_cast<GaussianBasisSet*>(
152 tkeyval->describedclassvalue("basis").pointer(),"main\n");
153 Ref<Molecule> mol = basis->molecule();
154
155 int tproc = tkeyval->intvalue("test_processor");
156 if (tproc >= msg->n()) tproc = 0;
157 int me = msg->me();
158
159 if (me == tproc) cout << "testing on processor " << tproc << endl;
160
161 int storage = tkeyval->intvalue("storage");
162 cout << "storage = " << storage << endl;
163 /* Ref<Integral> intgrlv3 = new IntegralV3(basis,basis,basis,basis);
164 Ref<Int1eV3> int1ev3 = new Int1eV3(intgrlv3.pointer(),basis,basis,1);
165 Ref<Int2eV3> int2ev3 = new Int2eV3(intgrlv3.pointer(),basis,basis,basis,basis,
166 1, storage);
167
168
169 int permute = tkeyval->booleanvalue("permute");
170 tim->enter("overlap");
171 if (me == tproc && tkeyval->booleanvalue("overlap")) {
172 cout << scprintf("testing overlap:\n");
173 test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::overlap, permute);
174 }
175 tim->change("kinetic");
176 if (me == tproc && tkeyval->booleanvalue("kinetic")) {
177 cout << scprintf("testing kinetic:\n");
178 test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::kinetic, permute);
179 }
180 tim->change("hcore");
181 if (me == tproc && tkeyval->booleanvalue("hcore")) {
182 cout << scprintf("testing hcore:\n");
183 test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::hcore, permute);
184 }
185 tim->change("nuclear");
186 if (me == tproc && tkeyval->booleanvalue("nuclear")) {
187 cout << scprintf("testing nuclear:\n");
188 test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::nuclear, permute);
189 }
190 tim->change("3 center");
191 if (me == tproc && tkeyval->booleanvalue("3")) {
192 test_3_center(tkeyval, int2ev3);
193 }
194 tim->change("4 center");
195 if (me == tproc && tkeyval->booleanvalue("4")) {
196 test_4_center(tkeyval, int2ev3);
197 }
198 tim->change("4 center der");
199 if (me == tproc && tkeyval->booleanvalue("4der")) {
200 test_4der_center(tkeyval, int2ev3);
201 }
202 tim->change("bound stats");
203 if (me == tproc && tkeyval->booleanvalue("boundstats")) {
204 do_bounds_stats(tkeyval, int2ev3);
205 }
206
207 tim->change("IntegralV3");*/
208
209 tim->enter("Integral");
210 Ref<Integral> integral = new IntegralV3(basis);
211#ifdef CINTS
212 Ref<Integral> integralcints = new IntegralCints(basis);
213#endif
214
215 Ref<OneBodyInt> overlapv3 = integral->overlap();
216 Ref<OneBodyInt> kineticv3 = integral->kinetic();
217 Ref<OneBodyInt> nuclearv3 = integral->nuclear();
218 Ref<OneBodyInt> hcorev3 = integral->hcore();
219
220#ifdef CINTS
221 Ref<OneBodyInt> overlapcints = integralcints->overlap();
222 testint(overlapcints);
223 Ref<OneBodyInt> kineticcints = integralcints->kinetic();
224 testint(kineticcints);
225 Ref<OneBodyInt> nuclearcints = integralcints->nuclear();
226 testint(nuclearcints);
227 Ref<OneBodyInt> hcorecints = integralcints->hcore();
228 testint(hcorecints);
229#endif
230
231 Ref<TwoBodyInt> erepv3 = integral->electron_repulsion();
232
233#ifdef CINTS
234 int storage_needed = integralcints->storage_required_eri(basis);
235 cout << scprintf("Need %d bytes to create EriCints\n",storage_needed);
236 Ref<TwoBodyInt> erepcints = integralcints->electron_repulsion();
237 testint(erepcints);
238 storage_needed = integralcints->storage_required_grt(basis);
239 cout << scprintf("Need %d bytes to create GRTCints\n",storage_needed);
240 Ref<TwoBodyInt> grtcints = integralcints->grt();
241 testint(grtcints);
242#endif
243 tim->exit();
244
245 // Test iterators
246 /* CartesianIterCints citer(3);
247 cout << "Cartesian f-shell:" << endl;
248 for ( citer.start(); int(citer) ; citer.next() )
249 cout << "nx = " << citer.a() << " ny = " << citer.b() << " nz = " << citer.c() << endl;
250 RedundantCartesianIterCints rciter(3);
251 cout << "Redundant Cartesian f-shell:" << endl;
252 for ( rciter.start(); int(rciter) ; rciter.next() )
253 cout << "nx = " << rciter.a() << " ny = " << rciter.b() << " nz = " << rciter.c() << endl;
254 */
255
256 //cout << "Testing Cints' overlap integrals against IntV3's" << endl;
257 // compare_1e_cints_vs_v3(overlapcints,overlapv3);
258 //cout << "Testing Cints' kinetic energy integrals against IntV3's" << endl;
259 //compare_1e_cints_vs_v3(kineticcints,kineticv3);
260 //cout << "Testing Cints' nuclear attraction integrals against IntV3's" << endl;
261 //compare_1e_cints_vs_v3(nuclearcints,nuclearv3);
262 //cout << "Testing Cints' core hamiltonian integrals against IntV3's" << endl;
263 //compare_1e_cints_vs_v3(hcorecints,hcorev3);
264
265 // compare_2e_permute(integralcints);
266
267 cout << "Testing Cints' ERIs against IntV3's" << endl;
268 compare_2e_cints_vs_v3(erepcints,erepv3);
269 //compare_2e_puream_cints_vs_v3(erepcints,erepv3);
270 cout << "Testing Cints' ERIs (from GRTCints) against IntV3's" << endl;
271 compare_2e_cints_vs_v3(grtcints,erepv3);
272
273#ifdef CINTS
274 cout << "Testing sums of Cints' ERIs against IntV3's" << endl;
275 compare_2e_bufsum_cints_vs_v3(erepcints,erepv3);
276 cout << "Testing sums of Cints' ERIs (from GRTCints) against IntV3's" << endl;
277 compare_2e_bufsum_cints_vs_v3(grtcints,erepv3);
278
279 cout << "Testing sums of unique Cints' ERIs against IntV3's" << endl;
280 erepcints->set_redundant(0);
281 erepv3->set_redundant(0);
282 grtcints->set_redundant(0);
283 compare_2e_unique_bufsum_cints_vs_v3(erepcints,erepv3);
284 cout << "Testing sums of unique Cints' ERIs (from GRTCints) against IntV3's" << endl;
285 compare_2e_unique_bufsum_cints_vs_v3(grtcints,erepv3);
286
287 cout << "Printing GRT integrals" << endl;
288 print_grt_ints(grtcints);
289#endif
290
291 // tim->print();
292 return 0;
293}
294
295void
296compare_1e_cints_vs_v3(Ref<OneBodyInt>& obcints, Ref<OneBodyInt>& obv3)
297{
298 Ref<GaussianBasisSet> basis = obcints->basis();
299 for (int sh1=4; sh1<basis->nshell(); sh1++)
300 for (int sh2=0; sh2<basis->nshell(); sh2++) {
301 int nbf2 = basis->shell(sh2).nfunction();
302 obv3->compute_shell(sh1,sh2);
303 obcints->compute_shell(sh1,sh2);
304 const double *buffercints = obcints->buffer();
305 const double *bufferv3 = obv3->buffer();
306
307 int bf1_offset = 0;
308 for (int gc1=0; gc1<basis->shell(sh1).ncontraction(); gc1++) {
309 int am1 = basis->shell(sh1).am(gc1);
310 CartesianIterCints citer1(am1);
311 CartesianIterV3 iter1(am1);
312 for ( citer1.start(); int(citer1) ; citer1.next() ) {
313 int bf1cints = bf1_offset + citer1.bfn();
314 int bf1v3;
315 for( iter1.start(); int(iter1) ; iter1.next() ) {
316 if (iter1.a() == citer1.a() &&
317 iter1.b() == citer1.b() &&
318 iter1.c() == citer1.c()) {
319 bf1v3 = bf1_offset + iter1.bfn();
320 break;
321 }
322 }
323
324 int bf2_offset = 0;
325 for (int gc2=0; gc2<basis->shell(sh2).ncontraction(); gc2++) {
326 int am2 = basis->shell(sh2).am(gc2);
327 CartesianIterCints citer2(am2);
328 CartesianIterV3 iter2(am2);
329
330 for ( citer2.start(); int(citer2) ; citer2.next() ) {
331 int bf2cints = bf2_offset + citer2.bfn();
332 int bf2v3;
333 for( iter2.start(); int(iter2) ; iter2.next() ) {
334 if (iter2.a() == citer2.a() &&
335 iter2.b() == citer2.b() &&
336 iter2.c() == citer2.c()) {
337 bf2v3 = bf2_offset + iter2.bfn();
338 break;
339 }
340 }
341
342 double valuecints = buffercints[bf1cints*nbf2 + bf2cints];
343 double valuev3 = bufferv3[bf1v3*nbf2 + bf2v3];
344 if (fabs(valuecints-valuev3) > 1E-13) {
345 cout << scprintf("Discrepancy in OEInt(sh1 = %d, sh2 = %d)\n",sh1,sh2);
346 cout << scprintf("bf1 = %d bf2 = %d OEIntegral(cints) = %20.15lf\n",bf1cints,bf2cints,valuecints);
347 cout << scprintf("bf1 = %d bf2 = %d OEIntegral(V3) = %20.15lf\n\n",bf1v3,bf2v3,valuev3);
348 }
349 }
350 bf2_offset += basis->shell(sh2).nfunction(gc2);
351 }
352 }
353 bf1_offset += basis->shell(sh1).nfunction(gc1);
354 }
355 }
356}
357
358void
359compare_2e_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3)
360{
361 const double *buffercints = tbcints->buffer();
362 const double *bufferv3 = tbv3->buffer();
363
364 Ref<GaussianBasisSet> basis = tbcints->basis();
365 for (int sh1=0; sh1<basis->nshell(); sh1++)
366 for (int sh2=0; sh2<basis->nshell(); sh2++)
367 for (int sh3=0; sh3<basis->nshell(); sh3++)
368 for (int sh4=0; sh4<basis->nshell(); sh4++)
369 {
370 //sh1=0;sh2=0;sh3=8;sh4=3;
371 //cout << scprintf("Computing TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4);
372 tbv3->compute_shell(sh1,sh2,sh3,sh4);
373 tbcints->compute_shell(sh1,sh2,sh3,sh4);
374
375 int nbf2 = basis->shell(sh2).nfunction();
376 int nbf3 = basis->shell(sh3).nfunction();
377 int nbf4 = basis->shell(sh4).nfunction();
378
379 int bf1_offset = 0;
380 for(int gc1=0;gc1<basis->shell(sh1).ncontraction(); gc1++) {
381 int am1 = basis->shell(sh1).am(gc1);
382 CartesianIterCints citer1(am1);
383 CartesianIterV3 iter1(am1);
384 for ( citer1.start(); int(citer1) ; citer1.next() ) {
385 int bf1cints = citer1.bfn();
386 int bf1v3;
387 for( iter1.start(); int(iter1) ; iter1.next() ) {
388 if (iter1.a() == citer1.a() &&
389 iter1.b() == citer1.b() &&
390 iter1.c() == citer1.c()) {
391 bf1v3 = iter1.bfn();
392 break;
393 }
394 }
395 bf1cints += bf1_offset;
396 bf1v3 += bf1_offset;
397
398 int bf2_offset = 0;
399 for(int gc2=0;gc2<basis->shell(sh2).ncontraction(); gc2++) {
400 int am2 = basis->shell(sh2).am(gc2);
401 CartesianIterCints citer2(am2);
402 CartesianIterV3 iter2(am2);
403 for ( citer2.start(); int(citer2) ; citer2.next() ) {
404 int bf2cints = citer2.bfn();
405 int bf2v3;
406 for( iter2.start(); int(iter2) ; iter2.next() ) {
407 if (iter2.a() == citer2.a() &&
408 iter2.b() == citer2.b() &&
409 iter2.c() == citer2.c()) {
410 bf2v3 = iter2.bfn();
411 break;
412 }
413 }
414 bf2cints += bf2_offset;
415 bf2v3 += bf2_offset;
416
417 int bf3_offset = 0;
418 for(int gc3=0;gc3<basis->shell(sh3).ncontraction(); gc3++) {
419 int am3 = basis->shell(sh3).am(gc3);
420 CartesianIterCints citer3(am3);
421 CartesianIterV3 iter3(am3);
422 for ( citer3.start(); int(citer3) ; citer3.next() ) {
423 int bf3cints = citer3.bfn();
424 int bf3v3;
425 for( iter3.start(); int(iter3) ; iter3.next() ) {
426 if (iter3.a() == citer3.a() &&
427 iter3.b() == citer3.b() &&
428 iter3.c() == citer3.c()) {
429 bf3v3 = iter3.bfn();
430 break;
431 }
432 }
433 bf3cints += bf3_offset;
434 bf3v3 += bf3_offset;
435
436 int bf4_offset = 0;
437 for(int gc4=0;gc4<basis->shell(sh4).ncontraction(); gc4++) {
438 int am4 = basis->shell(sh4).am(gc4);
439 CartesianIterCints citer4(am4);
440 CartesianIterV3 iter4(am4);
441 for ( citer4.start(); int(citer4) ; citer4.next() ) {
442 int bf4cints = citer4.bfn();
443 int bf4v3;
444 for( iter4.start(); int(iter4) ; iter4.next() ) {
445 if (iter4.a() == citer4.a() &&
446 iter4.b() == citer4.b() &&
447 iter4.c() == citer4.c()) {
448 bf4v3 = iter4.bfn();
449 break;
450 }
451 }
452 bf4cints += bf4_offset;
453 bf4v3 += bf4_offset;
454
455 double valuecints = buffercints[((bf1cints*nbf2 + bf2cints)*nbf3 + bf3cints)*nbf4 + bf4cints];
456 double valuev3 = bufferv3[((bf1v3*nbf2 + bf2v3)*nbf3 + bf3v3)*nbf4 + bf4v3];
457 if (fabs(valuecints-valuev3) > 1E-12) {
458 cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4);
459 cout << scprintf("bf1 = %d bf2 = %d bf3 = %d bf4 = %d TEIntegral(cints) = %20.15lf\n",
460 bf1cints,bf2cints,bf3cints,bf4cints,valuecints);
461 cout << scprintf("bf1 = %d bf2 = %d bf3 = %d bf4 = %d TEIntegral(V3) = %20.15lf\n\n",
462 bf1v3,bf2v3,bf3v3,bf4v3,valuev3);
463 }
464 }
465 bf4_offset+=basis->shell(sh4).nfunction(gc4);
466 }
467 }
468 bf3_offset+=basis->shell(sh3).nfunction(gc3);
469 }
470 }
471 bf2_offset+=basis->shell(sh2).nfunction(gc2);
472 }
473 }
474 bf1_offset+=basis->shell(sh1).nfunction(gc1);
475 }
476 //return;
477 }
478}
479
480
481void
482compare_2e_puream_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3)
483{
484 const double *buffercints = tbcints->buffer();
485 const double *bufferv3 = tbv3->buffer();
486
487 Ref<GaussianBasisSet> basis = tbcints->basis();
488 for (int sh1=0; sh1<basis->nshell(); sh1++)
489 for (int sh2=0; sh2<basis->nshell(); sh2++)
490 for (int sh3=0; sh3<basis->nshell(); sh3++)
491 for (int sh4=0; sh4<basis->nshell(); sh4++)
492 {
493 // sh1 = 0; sh2 = 0; sh3 = 6; sh4 = 13;
494 /* if ( !((basis->shell(sh1).has_pure() || basis->shell(sh1).max_am()==0) &&
495 (basis->shell(sh2).has_pure() || basis->shell(sh2).max_am()==0) &&
496 (basis->shell(sh3).has_pure() || basis->shell(sh3).max_am()==0) &&
497 (basis->shell(sh4).has_pure() || basis->shell(sh4).max_am()==0)) )
498 continue;*/
499 tbcints->compute_shell(sh1,sh2,sh3,sh4);
500 tbv3->compute_shell(sh1,sh2,sh3,sh4);
501
502 int nbf1 = basis->shell(sh1).nfunction();
503 int nbf2 = basis->shell(sh2).nfunction();
504 int nbf3 = basis->shell(sh3).nfunction();
505 int nbf4 = basis->shell(sh4).nfunction();
506
507 cout << scprintf("Computing TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4);
508 cout << scprintf("size = %d\n",nbf1*nbf2*nbf3*nbf4);
509 for(int ijkl=0;ijkl<nbf1*nbf2*nbf3*nbf4;ijkl++) {
510 double valuecints = buffercints[ijkl];
511 double valuev3 = bufferv3[ijkl];
512 if (fabs(valuecints-valuev3) > 1E-13) {
513 cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4);
514 cout << scprintf("1234 = %d TEIntegral(cints) = %20.15lf\n",
515 ijkl,valuecints);
516 cout << scprintf("TEIntegral(V3) = %20.15lf\n\n",
517 valuev3);
518 }
519 }
520 // return;
521 }
522}
523
524
525void
526compare_2e_bufsum_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3)
527{
528 Ref<GaussianBasisSet> basis = tbcints->basis();
529 const double *buffercints = tbcints->buffer();
530 const double *bufferv3 = tbv3->buffer();
531
532 for (int sh1=0; sh1<basis->nshell(); sh1++)
533 for (int sh2=0; sh2<basis->nshell(); sh2++)
534 for (int sh3=0; sh3<basis->nshell(); sh3++)
535 for (int sh4=0; sh4<basis->nshell(); sh4++)
536 {
537 // sh1=12;sh2=12;sh3=12;sh4=12;
538 tbcints->compute_shell(sh1,sh2,sh3,sh4);
539 tbv3->compute_shell(sh1,sh2,sh3,sh4);
540
541 int nbf1 = basis->shell(sh1).nfunction();
542 int nbf2 = basis->shell(sh2).nfunction();
543 int nbf3 = basis->shell(sh3).nfunction();
544 int nbf4 = basis->shell(sh4).nfunction();
545
546 double sum_cints = 0.0;
547 double sum_v3 = 0.0;
548
549 int index = 0;
550 for (int i=0; i<nbf1; i++) {
551 for (int j=0; j<nbf2; j++) {
552 for (int k=0; k<nbf3; k++) {
553 for (int l=0; l<nbf4; l++) {
554 sum_cints += buffercints[index];
555 sum_v3 += bufferv3[index];
556 /* cout << scprintf("index = %d TEIntegral(cints) = %20.15lf\n",
557 index,buffercints[index]);
558 cout << scprintf("index = %d TEIntegral(V3) = %20.15lf\n\n",
559 index,bufferv3[index]);
560 */
561 index++;
562 }
563 }
564 }
565 }
566
567 if (fabs(sum_cints-sum_v3) > 1E-12) {
568 cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4);
569 cout << scprintf("TEIntegralSum(cints) = %20.15lf\n",
570 sum_cints);
571 cout << scprintf("TEIntegralSum(V3) = %20.15lf\n\n",
572 sum_v3);
573 }
574 //return;
575 }
576}
577
578void
579compare_2e_unique_bufsum_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3)
580{
581 Ref<GaussianBasisSet> basis = tbcints->basis();
582 const double *buffercints = tbcints->buffer();
583 const double *bufferv3 = tbv3->buffer();
584
585 for (int sh1=0; sh1<basis->nshell(); sh1++)
586 for (int sh2=0; sh2<=sh1; sh2++)
587 for (int sh3=0; sh3<=sh1; sh3++)
588 for (int sh4=0; sh4 <= ((sh1==sh3) ? sh2 : sh3) ; sh4++)
589 {
590 tbcints->compute_shell(sh1,sh2,sh3,sh4);
591 tbv3->compute_shell(sh1,sh2,sh3,sh4);
592
593 int nbf1 = basis->shell(sh1).nfunction();
594 int nbf2 = basis->shell(sh2).nfunction();
595 int nbf3 = basis->shell(sh3).nfunction();
596 int nbf4 = basis->shell(sh4).nfunction();
597
598 double sum_cints = 0.0;
599 double sum_v3 = 0.0;
600
601 int e12 = (sh1 == sh2) ? 1 : 0;
602 int e34 = (sh3 == sh4) ? 1 : 0;
603 int e13e24 = (((sh1 == sh3)&&(sh2==sh4))||((sh1==sh4)&&(sh2==sh3))) ? 1 : 0;
604
605 int index = 0;
606 for (int i=0; i<nbf1; i++) {
607 int jmax = e12 ? i : nbf2-1;
608 for (int j=0; j<=jmax; j++) {
609 int kmax = e13e24 ? i : nbf3-1;
610 for (int k=0; k<=kmax; k++) {
611 int lmax = e34 ? ( (e13e24&&(i==k)) ? j : k) : ( (e13e24&&(i==k)) ? j : nbf4-1);
612 for (int l=0; l<=lmax; l++) {
613 sum_cints += buffercints[index];
614 sum_v3 += bufferv3[index];
615 /* cout << scprintf("index = %d TEIntegral(cints) = %20.15lf\n",
616 index,buffercints[index]);
617 cout << scprintf("index = %d TEIntegral(V3) = %20.15lf\n\n",
618 index,bufferv3[index]);
619 */
620 index++;
621 }
622 }
623 }
624 }
625
626 if (fabs(sum_cints-sum_v3) > 1E-12) {
627 cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4);
628 cout << scprintf("TEIntegralSum(cints) = %20.15lf\n",
629 sum_cints);
630 cout << scprintf("TEIntegralSum(V3) = %20.15lf\n\n",
631 sum_v3);
632 }
633 }
634}
635
636void
637print_grt_ints(Ref<TwoBodyInt>& tbcints)
638{
639 Ref<GaussianBasisSet> basis = tbcints->basis();
640 const double *buffer[4];
641 buffer[0] = tbcints->buffer(TwoBodyInt::eri);
642 buffer[1] = tbcints->buffer(TwoBodyInt::r12);
643 buffer[2] = tbcints->buffer(TwoBodyInt::r12t1);
644 buffer[3] = tbcints->buffer(TwoBodyInt::r12t2);
645 char teout_filename[] = "teout0.dat";
646 FILE *teout[4];
647
648 for(int te_type=0;te_type<4;te_type++) {
649 teout_filename[5] = te_type + '0';
650 teout[te_type] = fopen(teout_filename,"w");
651 }
652
653 for (int ush1=0; ush1<basis->nshell(); ush1++)
654 for (int ush2=0; ush2<=ush1; ush2++)
655 for (int ush3=0; ush3<=ush2; ush3++)
656 for (int ush4=0; ush4 <=ush3 ; ush4++)
657 {
658 int S1[3], S2[3], S3[3], S4[4];
659 int num = 1;
660 S1[0] = ush1;
661 S2[0] = ush2;
662 S3[0] = ush3;
663 S4[0] = ush4;
664 if (ush1==ush2 && ush1==ush3 || ush2==ush3 && ush2==ush4)
665 num=1;
666 else if (ush1==ush3 || ush2==ush4) {
667 num = 2;
668 S1[1] = ush1;
669 S2[1] = ush3;
670 S3[1] = ush2;
671 S4[1] = ush4;
672 }
673 else if (ush2==ush3) {
674 num = 2;
675 S1[1] = ush1;
676 S2[1] = ush4;
677 S3[1] = ush2;
678 S4[1] = ush3;
679 }
680 else if (ush1==ush2 || ush3==ush4) {
681 num = 2;
682 S1[1] = ush1;
683 S2[1] = ush3;
684 S3[1] = ush2;
685 S4[1] = ush4;
686 }
687 else {
688 num = 3;
689 S1[1] = ush1;
690 S2[1] = ush3;
691 S3[1] = ush2;
692 S4[1] = ush4;
693 S1[2] = ush1;
694 S2[2] = ush4;
695 S3[2] = ush2;
696 S4[2] = ush3;
697 }
698
699 for(int uq=0;uq<num;uq++) {
700 int sh1 = S1[uq];
701 int sh2 = S2[uq];
702 int sh3 = S3[uq];
703 int sh4 = S4[uq];
704
705 tbcints->compute_shell(sh1,sh2,sh3,sh4);
706
707 int nbf1 = basis->shell(sh1).nfunction();
708 int nbf2 = basis->shell(sh2).nfunction();
709 int nbf3 = basis->shell(sh3).nfunction();
710 int nbf4 = basis->shell(sh4).nfunction();
711
712 int e12 = (sh1 == sh2) ? 1 : 0;
713 int e34 = (sh3 == sh4) ? 1 : 0;
714 int e13e24 = (((sh1 == sh3)&&(sh2==sh4))||((sh1==sh4)&&(sh2==sh3))) ? 1 : 0;
715
716 int index = 0;
717 for (int i=0; i<nbf1; i++) {
718 int jmax = e12 ? i : nbf2-1;
719 for (int j=0; j<=jmax; j++) {
720 int kmax = e13e24 ? i : nbf3-1;
721 for (int k=0; k<=kmax; k++) {
722 int lmax = e34 ? ( (e13e24&&(i==k)) ? j : k) : ( (e13e24&&(i==k)) ? j : nbf4-1);
723 for (int l=0; l<=lmax; l++) {
724
725 double integral = buffer[0][index];
726 if (fabs(integral) > 1E-15) {
727 fprintf(teout[0], "%5d%5d%5d%5d%20.10lf\n",
728 basis->shell_to_function(sh1) + i+1,
729 basis->shell_to_function(sh2) + j+1,
730 basis->shell_to_function(sh3) + k+1,
731 basis->shell_to_function(sh4) + l+1,
732 integral);
733 }
734 integral = buffer[1][index];
735 if (fabs(integral) > 1E-15) {
736 fprintf(teout[1], "%5d%5d%5d%5d%20.10lf\n",
737 basis->shell_to_function(sh1) + i+1,
738 basis->shell_to_function(sh2) + j+1,
739 basis->shell_to_function(sh3) + k+1,
740 basis->shell_to_function(sh4) + l+1,
741 integral);
742 }
743 integral = buffer[2][index];
744 if (fabs(integral) > 1E-15) {
745 fprintf(teout[2], "%5d%5d%5d%5d%20.10lf\n",
746 basis->shell_to_function(sh1) + i+1,
747 basis->shell_to_function(sh2) + j+1,
748 basis->shell_to_function(sh3) + k+1,
749 basis->shell_to_function(sh4) + l+1,
750 integral);
751 }
752 integral = buffer[3][index];
753 if (fabs(integral) > 1E-15) {
754 fprintf(teout[3], "%5d%5d%5d%5d%20.10lf\n",
755 basis->shell_to_function(sh1) + i+1,
756 basis->shell_to_function(sh2) + j+1,
757 basis->shell_to_function(sh3) + k+1,
758 basis->shell_to_function(sh4) + l+1,
759 integral);
760 }
761 index++;
762 }
763 }
764 }
765 }
766 }
767 }
768 for(int te_type=0;te_type<4;te_type++)
769 fclose(teout[te_type]);
770}
771
772void
773compare_2e_permute(Ref<Integral>& cints)
774{
775 Ref<TwoBodyInt> tb1 = cints->electron_repulsion();
776 Ref<TwoBodyInt> tb2 = cints->electron_repulsion();
777 Ref<GaussianBasisSet> basis = tb1->basis();
778 const double *buffer1 = tb1->buffer();
779 const double *buffer2 = tb2->buffer();
780
781 int sh1 = 0;
782 int sh2 = 0;
783 int sh3 = 4;
784 int sh4 = 0;
785
786 tb1->compute_shell(sh1,sh2,sh3,sh4);
787 tb2->compute_shell(sh1,sh2,sh4,sh3);
788
789 int nbf1 = basis->shell(sh1).nfunction();
790 int nbf2 = basis->shell(sh2).nfunction();
791 int nbf3 = basis->shell(sh3).nfunction();
792 int nbf4 = basis->shell(sh4).nfunction();
793
794 for(int index = 0; index<nbf1*nbf2*nbf3*nbf4; index++)
795 if (fabs(buffer1[index]-buffer2[index]) > 1E-13)
796 {
797 cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4);
798 cout << scprintf("TEIntegral(cints1) = %20.15lf\n",buffer1[index]);
799 cout << scprintf("TEIntegral(cints2) = %20.15lf\n\n",buffer2[index]);
800 }
801}
802
803
804void
805do_shell_test_1e(const Ref<Int1eV3> &int1ev3,
806 void (Int1eV3::*int_shell_1e)(int,int),
807 int permute, int i, int j, int na, int nb,
808 double *buf, double *pbuf)
809{
810 int ii = 0;
811 int a;
812 double *buffer = int1ev3->buffer();
813 (int1ev3->*int_shell_1e)(i, j);
814 for (a=0; a<na*nb; a++) {
815 buf[a] = buffer[a];
816 }
817 (int1ev3->*int_shell_1e)(j, i);
818 for (a=0; a<na*nb; a++) {
819 pbuf[a] = buffer[a];
820 }
821 for (a=0; a<na; a++) {
822 for (int b=0; b<nb; b++) {
823 if (fabs(buf[ii] - pbuf[a + na*b]) > 1.0e-13) {
824 cout << scprintf("----- 1e perm failed:"
825 "<%d %d|%d %d>:"
826 " %18.14f != %18.14f "
827 "<%d %d|%d %d>\n",
828 i, a, j, b,
829 buf[ii],
830 pbuf[a + na*b],
831 j, b, i, a);
832 }
833 if (fabs(buf[ii]) > 1.0e-15) {
834 cout << scprintf(" <(%d %d)|(%d %d)> = %15.11f\n",
835 i,a,j,b, buf[ii]);
836 }
837 ii++;
838 }
839 }
840}
841
842void
843test_int_shell_1e(const Ref<KeyVal>& keyval, const Ref<Int1eV3> &int1ev3,
844 void (Int1eV3::*int_shell_1e)(int,int),
845 int permute)
846{
847 int flags = 0;
848 Ref<GaussianBasisSet> basis = int1ev3->basis();
849 int maxfunc = basis->max_nfunction_in_shell();
850 int size = maxfunc * maxfunc;
851 double *buf = new double[size];
852 double *pbuf = new double[size];
853 int nshell = int1ev3->basis()->nshell();
854
855 for (int i=0; i<nshell; i++) {
856 int na = basis->shell(i).nfunction();
857 for (int j=0; j<nshell; j++) {
858 int nb = basis->shell(j).nfunction();
859 do_shell_test_1e(int1ev3, int_shell_1e, permute,
860 i, j, na, nb, buf, pbuf);
861
862 }
863 }
864
865 delete[] buf;
866 delete[] pbuf;
867}
868
869void
870test_3_center(const Ref<KeyVal>& keyval, const Ref<Int2eV3> &int2ev3)
871{
872 int ii, i,j,k,l,m,n;
873
874 int2ev3->set_redundant(1);
875 int2ev3->set_permute(0);
876 double *buffer = int2ev3->buffer();
877 int nshell = int2ev3->basis()->nshell();
878
879 for (i=0; i<nshell; i++) {
880 for (j=0; j<nshell; j++) {
881 int sh[2], sizes[2];
882 sh[0] = i;
883 sh[1] = j;
884 int2ev3->erep_2center(sh,sizes);
885 ii = 0;
886 for (k=0; k<sizes[0]; k++) {
887 for (l=0; l<sizes[1]; l++) {
888 if (fabs(buffer[ii])>1.0e-15)
889 cout << scprintf(" ((%d %d)|(%d %d)) = %15.11f\n",
890 sh[0],k,sh[1],l, buffer[ii]);
891 ii++;
892 }
893 }
894 }
895 }
896
897 for (i=0; i<nshell; i++) {
898 for (j=0; j<nshell; j++) {
899 for (m=0; m<nshell; m++) {
900 int sh[3], sizes[3];
901 sh[0] = i;
902 sh[1] = j;
903 sh[2] = m;
904 int2ev3->erep_3center(sh,sizes);
905 ii = 0;
906 for (k=0; k<sizes[0]; k++) {
907 for (l=0; l<sizes[1]; l++) {
908 for (n=0; n<sizes[2]; n++) {
909 if (fabs(buffer[ii])>1.0e-15)
910 cout << scprintf(
911 " ((%d %d)|(%d %d)(%d %d)) = %15.11f\n",
912 sh[0],k,sh[1],l,sh[2],n, buffer[ii]);
913 ii++;
914 }
915 }
916 }
917 }
918 }
919 }
920
921}
922
923void
924init_shell_perm(const Ref<Int2eV3> &int2ev3, double *integrals,
925 double buff[maxint][maxint][maxint][maxint],
926 int sh[4], int sizes[4])
927{
928 int i, j, k, l;
929 int oldp = int2ev3->permute();
930 int2ev3->set_permute(0);
931 int2ev3->erep(sh, sizes);
932 int2ev3->set_permute(oldp);
933 for (i=0; i<sizes[0]; i++) {
934 for (j=0; j<sizes[1]; j++) {
935 for (k=0; k<sizes[2]; k++) {
936 for (l=0; l<sizes[3]; l++) {
937 buff[i][j][k][l] = *integrals++;
938 }
939 }
940 }
941 }
942}
943
944void
945check_shell_perm(const Ref<Int2eV3> &int2ev3, double *integrals,
946 double buff[maxint][maxint][maxint][maxint],
947 int sh[4], int sizes[4], int p0, int p1, int p2, int p3)
948{
949 int ip[4], p[4];
950 int psizes[4];
951 int psh[4];
952 int index = 0;
953 int i[4];
954 p[0] = p0;
955 p[1] = p1;
956 p[2] = p2;
957 p[3] = p3;
958 ip[p0] = 0;
959 ip[p1] = 1;
960 ip[p2] = 2;
961 ip[p3] = 3;
962 psh[0] = sh[p0];
963 psh[1] = sh[p1];
964 psh[2] = sh[p2];
965 psh[3] = sh[p3];
966 int oldp = int2ev3->permute();
967 int2ev3->set_permute(0);
968 int2ev3->erep(psh, psizes);
969 int2ev3->set_permute(oldp);
970 for (i[0]=0; i[0]<psizes[0]; i[0]++) {
971 for (i[1]=0; i[1]<psizes[1]; i[1]++) {
972 for (i[2]=0; i[2]<psizes[2]; i[2]++) {
973 for (i[3]=0; i[3]<psizes[3]; i[3]++) {
974 if (fabs(buff[i[ip[0]]][i[ip[1]]][i[ip[2]]][i[ip[3]]]
975 - integrals[index]) > 1.0e-13) {
976 cout << scprintf("perm %d %d %d %d failed:"
977 "((%d %d)(%d %d)|(%d %d)(%d %d)):"
978 " %18.14f != %18.14f "
979 "((%d %d)(%d %d)|(%d %d)(%d %d))\n",
980 p0, p1, p2, p3,
981 sh[0],i[0], sh[1],i[1], sh[2],i[2], sh[3],i[3],
982 buff[i[ip[0]]][i[ip[1]]][i[ip[2]]][i[ip[3]]],
983 integrals[index],
984 psh[0],i[p[0]], psh[1],i[p[1]],
985 psh[2],i[p[2]], psh[3],i[p[3]]);
986 }
987 index++;
988 }
989 }
990 }
991 }
992}
993
994void
995do_shell_quartet_test(const Ref<Int2eV3> &int2ev3,
996 int print, int printbounds, int bounds, int permute,
997 const Ref<KeyVal>& keyval,
998 int i, int j, int k, int l)
999{
1000 int sh[4], sizes[4];
1001 int ibuf;
1002 int ii, jj, kk, ll;
1003 sh[0] = i;
1004 sh[1] = j;
1005 sh[2] = k;
1006 sh[3] = l;
1007 double maxintegral, integralbound;
1008 int boundijkl;
1009 if (bounds) {
1010 integralbound
1011 = int2ev3->logbound_to_bound(
1012 (boundijkl = int2ev3->erep_4bound(i,j,k,l))
1013 );
1014 }
1015 double *buffer = int2ev3->buffer();
1016 int2ev3->erep(sh,sizes);
1017 ibuf = 0;
1018 maxintegral = 0.0;
1019 for (ii=0; ii<sizes[0]; ii++) {
1020 for (jj=0; jj<sizes[1]; jj++) {
1021 for (kk=0; kk<sizes[2]; kk++) {
1022 for (ll=0; ll<sizes[3]; ll++) {
1023 double absint = fabs(buffer[ibuf]);
1024 if (absint > maxintegral) {
1025 maxintegral = absint;
1026 }
1027 if (bounds && absint > integralbound) {
1028 cout << scprintf("((%d %d)(%d %d)|(%d %d)(%d %d)) = %15.11f, "
1029 "bound = %15.11f\n",
1030 sh[0], ii, sh[1], jj, sh[2], kk, sh[3], ll,
1031 buffer[ibuf], integralbound);
1032 abort();
1033 }
1034 if (print && (absint > 1.0e-9
1035 || (bounds && integralbound > 1.0e-9))) {
1036 cout << scprintf(" ((%d %d)(%d %d)|(%d %d)(%d %d))"
1037 " = %15.11f",
1038 sh[0],ii,
1039 sh[1],jj,
1040 sh[2],kk,
1041 sh[3],ll,
1042 buffer[ibuf]);
1043 if (bounds) {
1044 cout << scprintf(" (%2d%% of bound)",
1045 (int)(100*(absint/integralbound)));
1046 }
1047 cout << scprintf("\n");
1048 }
1049 ibuf++;
1050 }
1051 }
1052 }
1053 }
1054
1055 if (permute) {
1056 double buff1[maxint][maxint][maxint][maxint];
1057 sh[0] = i;
1058 sh[1] = j;
1059 sh[2] = k;
1060 sh[3] = l;
1061 init_shell_perm(int2ev3, buffer, buff1, sh, sizes);
1062 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 0, 1, 2, 3);
1063 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 1, 0, 2, 3);
1064 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 0, 1, 3, 2);
1065 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 1, 0, 3, 2);
1066 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 2, 3, 0, 1);
1067 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 2, 3, 1, 0);
1068 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 3, 2, 0, 1);
1069 check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 3, 2, 1, 0);
1070 }
1071
1072 if (bounds) {
1073 int boundij = int2ev3->erep_4bound(i,j,-1,-1);
1074 int boundkl = int2ev3->erep_4bound(-1,-1,k,l);
1075 int badbound = 0;
1076 if (boundij < boundijkl || boundkl < boundijkl) {
1077 badbound = 1;
1078 }
1079 if (badbound || printbounds) {
1080 cout << scprintf("max(%d,%d,%d,%d)=%7.4f, bnd=%7.4f, "
1081 "bnd(%d,%d,*,*)=%7.4f, bnd(*,*,%d,%d)=%7.4f\n",
1082 i, j, k, l, maxintegral, integralbound,
1083 i,j, int2ev3->logbound_to_bound(boundij),
1084 k,l, int2ev3->logbound_to_bound(boundkl));
1085 }
1086 if (badbound) {
1087 cout << scprintf("ERROR: bad bound\n");
1088 abort();
1089 }
1090 }
1091}
1092
1093void
1094do_4_center_test(const Ref<Int2eV3> &int2ev3, int print, int printbounds,
1095 int bounds, int permute,
1096 const Ref<KeyVal>& keyval)
1097{
1098 int ii,jj,kk,ll, i,j,k,l, ibuf;
1099 int nshell = int2ev3->basis()->nshell();
1100 int unique = keyval->booleanvalue("unique");
1101 int timestats = keyval->booleanvalue("timestats");
1102 Ref<RegionTimer> timer = new RegionTimer();
1103
1104 if (!timestats) {
1105 for (i=0; i<nshell; i++) {
1106 int jmax = nshell - 1;
1107 if (unique) jmax = i;
1108 for (j=0; j<=jmax; j++) {
1109 int kmax = nshell - 1;
1110 if (unique) kmax = i;
1111 for (k=0; k<=kmax; k++) {
1112 int lmax = nshell - 1;
1113 if (unique) {
1114 if (k==i) lmax = j;
1115 else lmax = k;
1116 }
1117 for (l=0; l<=lmax; l++) {
1118 do_shell_quartet_test(int2ev3, print, printbounds,
1119 bounds, permute,
1120 keyval, i, j, k, l);
1121 }
1122 }
1123 }
1124 }
1125 }
1126 if (timestats && nshell) {
1127 unsigned short seed = 1234;
1128 seed48(&seed);
1129 const int nsample = 5000;
1130 const int ntrials = 50;
1131 double times[ntrials];
1132 for (i=0; i<ntrials; i++) {
1133 double t1 = timer->get_cpu_time();
1134 for (j=0; j<nsample; j++) {
1135 // pick an integral at random
1136 int ish = int(drand48()*nshell);
1137 int jsh = int(drand48()*ish);
1138 int ksh = int(drand48()*ish);
1139 int lsh;
1140 if (ish==ksh) lsh = int(drand48()*jsh);
1141 else lsh = int(drand48()*ksh);
1142 int sh[4], sizes[4];
1143 if (ish >= nshell) ish = nshell-1;
1144 if (jsh >= nshell) jsh = nshell-1;
1145 if (ksh >= nshell) ksh = nshell-1;
1146 if (lsh >= nshell) lsh = nshell-1;
1147 sh[0] = ish;
1148 sh[1] = jsh;
1149 sh[2] = ksh;
1150 sh[3] = lsh;
1151 int2ev3->erep(sh,sizes);
1152 }
1153 double t2 = timer->get_cpu_time();
1154 times[i] = t2-t1;
1155 }
1156 double ave = 0.0;
1157 for (i=0; i<ntrials; i++) {
1158 ave += times[i];
1159 }
1160 ave /= ntrials;
1161 double sigma2 = 0.0;
1162 for (i=0; i<ntrials; i++) {
1163 double diff = times[i] - ave;
1164 sigma2 += diff*diff;
1165 }
1166 double sigma = sqrt(sigma2/ntrials);
1167 // adjust sigma and ave from the trial average results to
1168 // the integral results
1169 ave /= nsample;
1170 sigma /= sqrt(double(nsample));
1171 cout << scprintf(" ave = %10.8f sigma = %10.8f (microsecs)\n"
1172 " sigma/ave = %10.4f",
1173 ave*1e6, sigma*1e6,
1174 sigma/ave)
1175 << endl;
1176 }
1177}
1178
1179void
1180test_4_center(const Ref<KeyVal>& keyval, const Ref<Int2eV3> &int2ev3)
1181{
1182 int i;
1183
1184 cout << scprintf("4 center test:\n");
1185 cout << scprintf(" on entry int2ev3 used %d bytes\n", int2ev3->used_storage());
1186
1187 int2ev3->set_permute(0);
1188 int2ev3->set_redundant(1);
1189
1190 int storage = keyval->intvalue("storage") - int2ev3->used_storage();
1191 if (storage < 0) storage = 0;
1192 if (keyval->booleanvalue("store_integrals")) storage = 0;
1193 int niter = keyval->intvalue("niter");
1194 int print = keyval->booleanvalue("print");
1195 int bounds = keyval->booleanvalue("bounds");
1196 int permute = keyval->booleanvalue("permute");
1197 int printbounds = keyval->booleanvalue("printbounds");
1198
1199 cout << scprintf(" storage = %d\n", storage);
1200 cout << scprintf(" niter = %d\n", niter);
1201 cout << scprintf(" print = %d\n", print);
1202 cout << scprintf(" bounds = %d\n", bounds);
1203 cout << scprintf(" permute = %d\n", permute);
1204 cout << scprintf("printbounds = %d\n", printbounds);
1205
1206 if (bounds) int2ev3->init_bounds();
1207
1208 int2ev3->init_storage(storage);
1209
1210 for (i=0; i<niter; i++) {
1211 do_4_center_test(int2ev3, print, printbounds, bounds, permute, keyval);
1212 }
1213
1214 if (keyval->count("quartet") == 4) {
1215 do_shell_quartet_test(int2ev3, print, printbounds, bounds, permute,
1216 keyval,
1217 keyval->intvalue("quartet", 0),
1218 keyval->intvalue("quartet", 1),
1219 keyval->intvalue("quartet", 2),
1220 keyval->intvalue("quartet", 3));
1221 }
1222
1223 int2ev3->done_storage();
1224 int2ev3->done_bounds();
1225}
1226
1227void
1228do_shell_quartet_der_test(const Ref<Int2eV3> &int2ev3,
1229 double* buffer, int print, int printbounds,
1230 int bounds, int permute,
1231 const Ref<KeyVal>& keyval,
1232 int i, int j, int k, int l)
1233{
1234 int ii,jj,kk,ll, ibuf, ider, xyz;
1235 der_centersv3_t dercenters;
1236
1237 int sh[4], sizes[4];
1238 sh[0] = i;
1239 sh[1] = j;
1240 sh[2] = k;
1241 sh[3] = l;
1242 double maxintegral = 0.0, integralbound;
1243 int boundijkl;
1244 if (bounds) {
1245 integralbound
1246 = int2ev3->logbound_to_bound(
1247 (boundijkl = int2ev3->erep_4bound_1der(i,j,k,l))
1248 );
1249 }
1250 int2ev3->erep_all1der(sh,sizes,&dercenters);
1251 ibuf = 0;
1252 for (ider=0; ider<dercenters.n; ider++) {
1253 for (xyz=0; xyz<3; xyz++) {
1254 for (ii=0; ii<sizes[0]; ii++) {
1255 for (jj=0; jj<sizes[1]; jj++) {
1256 for (kk=0; kk<sizes[2]; kk++) {
1257 for (ll=0; ll<sizes[3]; ll++) {
1258 double absint = fabs(buffer[ibuf]);
1259 if (absint > maxintegral) {
1260 maxintegral = absint;
1261 }
1262 if (bounds && absint > integralbound) {
1263 cout << scprintf("((%d %d)(%d %d)|(%d %d)(%d %d))"
1264 " = %15.11f, bound = %15.11f\n",
1265 sh[0], ii, sh[1], jj,
1266 sh[2], kk, sh[3], ll,
1267 buffer[ibuf], integralbound);
1268 abort();
1269 }
1270 if (print && absint > 1.0e-15) {
1271 cout << scprintf(" ((%d %d)(%d %d)"
1272 "|(%d %d)(%d %d))(%d %d)"
1273 " = %15.11f\n",
1274 sh[0],ii,
1275 sh[1],jj,
1276 sh[2],kk,
1277 sh[3],ll,
1278 dercenters.num[ider], xyz,
1279 buffer[ibuf]
1280 );
1281 }
1282 ibuf++;
1283 }
1284 }
1285 }
1286 }
1287 }
1288 }
1289
1290 if (bounds) {
1291 int boundij = int2ev3->erep_4bound_1der(i,j,-1,-1);
1292 int boundkl = int2ev3->erep_4bound_1der(-1,-1,k,l);
1293 int badbound = 0;
1294 if (boundij < boundijkl || boundkl < boundijkl) {
1295 badbound = 1;
1296 }
1297 if (badbound || printbounds) {
1298 cout << scprintf("max(%d,%d,%d,%d)=%7.4f, bnd=%7.4f, "
1299 "bnd(%d,%d,*,*)=%8.4f, bnd(*,*,%d,%d)=%8.4f\n",
1300 i, j, k, l, maxintegral, integralbound,
1301 i,j, int2ev3->logbound_to_bound(boundij),
1302 k,l, int2ev3->logbound_to_bound(boundkl));
1303 }
1304 if (badbound) {
1305 cout << scprintf("ERROR: bad bound\n");
1306 abort();
1307 }
1308 }
1309}
1310
1311void
1312do_test_4der_center(const Ref<Int2eV3> &int2ev3,
1313 double* buffer, int print, int printbounds,
1314 int bounds, int permute,
1315 const Ref<KeyVal>& keyval)
1316{
1317 int i,j,k,l;
1318 int nshell = int2ev3->basis()->nshell();
1319 for (i=0; i<nshell; i++) {
1320 for (j=0; j<nshell; j++) {
1321 for (k=0; k<nshell; k++) {
1322 for (l=0; l<nshell; l++) {
1323 do_shell_quartet_der_test(int2ev3, buffer,
1324 print, printbounds,
1325 bounds, permute,
1326 keyval,
1327 i, j, k, l);
1328 }
1329 }
1330 }
1331 }
1332}
1333
1334void
1335test_4der_center(const Ref<KeyVal>& keyval, const Ref<Int2eV3> &int2ev3)
1336{
1337 int i;
1338
1339 int2ev3->set_permute(0);
1340 int2ev3->set_redundant(1);
1341 double *buffer = int2ev3->buffer();
1342
1343 int niter = keyval->intvalue("niter");
1344 int print = keyval->booleanvalue("print");
1345 int bounds = keyval->booleanvalue("bounds");
1346 int printbounds = keyval->booleanvalue("printbounds");
1347 int permute = keyval->booleanvalue("permute");
1348
1349 cout << scprintf("4 center derivative test:\n");
1350 cout << scprintf(" niter = %d\n", niter);
1351 cout << scprintf(" print = %d\n", print);
1352 cout << scprintf(" bounds = %d\n", bounds);
1353 cout << scprintf("printbounds = %d\n", printbounds);
1354 cout << scprintf(" permute = %d\n", permute);
1355
1356 if (bounds) int2ev3->init_bounds_1der();
1357
1358 for (i=0; i<niter; i++) {
1359 do_test_4der_center(int2ev3, buffer,
1360 print, printbounds, bounds, permute, keyval);
1361 }
1362
1363 if (keyval->count("quartet") == 4) {
1364 do_shell_quartet_der_test(int2ev3, buffer, print, printbounds,
1365 bounds, permute,
1366 keyval,
1367 keyval->intvalue("quartet", 0),
1368 keyval->intvalue("quartet", 1),
1369 keyval->intvalue("quartet", 2),
1370 keyval->intvalue("quartet", 3));
1371 }
1372
1373 if (bounds) int2ev3->done_bounds_1der();
1374}
1375
1376/////////////////////////////////////////////////////////////////////////////
1377
1378// Local Variables:
1379// mode: c++
1380// c-file-style: "CLJ"
1381// End:
Note: See TracBrowser for help on using the repository browser.