// // cintstest.cc // // Copyright (C) 2001 Edward Valeev // // Author: Edward Valeev // Maintainer: EV // // This file is part of the SC Toolkit. // // The SC Toolkit is free software; you can redistribute it and/or modify // it under the terms of the GNU Library General Public License as published by // the Free Software Foundation; either version 2, or (at your option) // any later version. // // The SC Toolkit is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Library General Public License for more details. // // You should have received a copy of the GNU Library General Public License // along with the SC Toolkit; see the file COPYING.LIB. If not, write to // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. // // The U.S. Government is granted a limited license as per AL 91-7. // #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; using namespace sc; #define CINTS void compare_1e_cints_vs_v3(Ref& obcints, Ref& obv3); void compare_2e_cints_vs_v3(Ref& tbcints, Ref& tbv3); void compare_2e_puream_cints_vs_v3(Ref& tbcints, Ref& tbv3); void compare_2e_bufsum_cints_vs_v3(Ref& tbcints, Ref& tbv3); void compare_2e_unique_bufsum_cints_vs_v3(Ref& tbcints, Ref& tbv3); void print_grt_ints(Ref& tbcints); void compare_2e_permute(Ref& cints); void test_int_shell_1e(const Ref&, const Ref &int1ev3, void (Int1eV3::*int_shell_1e)(int,int), int permute); void test_3_center(const Ref&, const Ref &); void test_4_center(const Ref& keyval, const Ref &int2ev3); void test_4der_center(const Ref&, const Ref &int2ev3); #define maxint 9 void testint(const Ref& in) { if (in.null()) { cout << "null integral generator" << endl; abort(); } } void testint(const Ref& in) { if (in.null()) { cout << "null integral generator" << endl; abort(); } } void testint(const Ref& in) { if (in.null()) { cout << "null integral generator" << endl; abort(); } } void testint(const Ref& in) { if (in.null()) { cout << "null integral generator" << endl; abort(); } } /*void do_bounds_stats(const Ref& keyval, const Ref &int2ev3) { int i,j; int nshell = int2ev3->basis()->nshell(); int eps = -10; int *nonzero = new int[nshell]; for (i=0; ierep_4bound(i,j,-1,-1) > eps) { nonzero[i]++; } } } for (i=0; ibasis()->shell_to_center(i); int npq = (i*(i+1))/2; cout<eps=%4d npq>eps/nsh=%9.4f /nat=%9.4f", i, natom, npq, nonzero[i], double(nonzero[i])/i, double(nonzero[i])/natom) << endl; } delete[] nonzero; } */ int main(int argc, char **argv) { int ii, i,j,k,l,m,n; Ref msg = MessageGrp::initial_messagegrp(argc,argv); if (msg.null()) msg = new ProcMessageGrp(); MessageGrp::set_default_messagegrp(msg); Ref tim = new ParallelRegionTimer(msg,"cintstest", 1, 1); char *infile = new char[strlen(SRCDIR)+strlen("/cintstest.in")+1]; sprintf(infile,SRCDIR "/cintstest.in"); if (argc == 2) { infile = argv[1]; } Ref pkv(new ParsedKeyVal(infile)); Ref tkeyval(new PrefixKeyVal(":test", pkv)); Ref basis = require_dynamic_cast( tkeyval->describedclassvalue("basis").pointer(),"main\n"); Ref mol = basis->molecule(); int tproc = tkeyval->intvalue("test_processor"); if (tproc >= msg->n()) tproc = 0; int me = msg->me(); if (me == tproc) cout << "testing on processor " << tproc << endl; int storage = tkeyval->intvalue("storage"); cout << "storage = " << storage << endl; /* Ref intgrlv3 = new IntegralV3(basis,basis,basis,basis); Ref int1ev3 = new Int1eV3(intgrlv3.pointer(),basis,basis,1); Ref int2ev3 = new Int2eV3(intgrlv3.pointer(),basis,basis,basis,basis, 1, storage); int permute = tkeyval->booleanvalue("permute"); tim->enter("overlap"); if (me == tproc && tkeyval->booleanvalue("overlap")) { cout << scprintf("testing overlap:\n"); test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::overlap, permute); } tim->change("kinetic"); if (me == tproc && tkeyval->booleanvalue("kinetic")) { cout << scprintf("testing kinetic:\n"); test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::kinetic, permute); } tim->change("hcore"); if (me == tproc && tkeyval->booleanvalue("hcore")) { cout << scprintf("testing hcore:\n"); test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::hcore, permute); } tim->change("nuclear"); if (me == tproc && tkeyval->booleanvalue("nuclear")) { cout << scprintf("testing nuclear:\n"); test_int_shell_1e(tkeyval, int1ev3, &Int1eV3::nuclear, permute); } tim->change("3 center"); if (me == tproc && tkeyval->booleanvalue("3")) { test_3_center(tkeyval, int2ev3); } tim->change("4 center"); if (me == tproc && tkeyval->booleanvalue("4")) { test_4_center(tkeyval, int2ev3); } tim->change("4 center der"); if (me == tproc && tkeyval->booleanvalue("4der")) { test_4der_center(tkeyval, int2ev3); } tim->change("bound stats"); if (me == tproc && tkeyval->booleanvalue("boundstats")) { do_bounds_stats(tkeyval, int2ev3); } tim->change("IntegralV3");*/ tim->enter("Integral"); Ref integral = new IntegralV3(basis); #ifdef CINTS Ref integralcints = new IntegralCints(basis); #endif Ref overlapv3 = integral->overlap(); Ref kineticv3 = integral->kinetic(); Ref nuclearv3 = integral->nuclear(); Ref hcorev3 = integral->hcore(); #ifdef CINTS Ref overlapcints = integralcints->overlap(); testint(overlapcints); Ref kineticcints = integralcints->kinetic(); testint(kineticcints); Ref nuclearcints = integralcints->nuclear(); testint(nuclearcints); Ref hcorecints = integralcints->hcore(); testint(hcorecints); #endif Ref erepv3 = integral->electron_repulsion(); #ifdef CINTS int storage_needed = integralcints->storage_required_eri(basis); cout << scprintf("Need %d bytes to create EriCints\n",storage_needed); Ref erepcints = integralcints->electron_repulsion(); testint(erepcints); storage_needed = integralcints->storage_required_grt(basis); cout << scprintf("Need %d bytes to create GRTCints\n",storage_needed); Ref grtcints = integralcints->grt(); testint(grtcints); #endif tim->exit(); // Test iterators /* CartesianIterCints citer(3); cout << "Cartesian f-shell:" << endl; for ( citer.start(); int(citer) ; citer.next() ) cout << "nx = " << citer.a() << " ny = " << citer.b() << " nz = " << citer.c() << endl; RedundantCartesianIterCints rciter(3); cout << "Redundant Cartesian f-shell:" << endl; for ( rciter.start(); int(rciter) ; rciter.next() ) cout << "nx = " << rciter.a() << " ny = " << rciter.b() << " nz = " << rciter.c() << endl; */ //cout << "Testing Cints' overlap integrals against IntV3's" << endl; // compare_1e_cints_vs_v3(overlapcints,overlapv3); //cout << "Testing Cints' kinetic energy integrals against IntV3's" << endl; //compare_1e_cints_vs_v3(kineticcints,kineticv3); //cout << "Testing Cints' nuclear attraction integrals against IntV3's" << endl; //compare_1e_cints_vs_v3(nuclearcints,nuclearv3); //cout << "Testing Cints' core hamiltonian integrals against IntV3's" << endl; //compare_1e_cints_vs_v3(hcorecints,hcorev3); // compare_2e_permute(integralcints); cout << "Testing Cints' ERIs against IntV3's" << endl; compare_2e_cints_vs_v3(erepcints,erepv3); //compare_2e_puream_cints_vs_v3(erepcints,erepv3); cout << "Testing Cints' ERIs (from GRTCints) against IntV3's" << endl; compare_2e_cints_vs_v3(grtcints,erepv3); #ifdef CINTS cout << "Testing sums of Cints' ERIs against IntV3's" << endl; compare_2e_bufsum_cints_vs_v3(erepcints,erepv3); cout << "Testing sums of Cints' ERIs (from GRTCints) against IntV3's" << endl; compare_2e_bufsum_cints_vs_v3(grtcints,erepv3); cout << "Testing sums of unique Cints' ERIs against IntV3's" << endl; erepcints->set_redundant(0); erepv3->set_redundant(0); grtcints->set_redundant(0); compare_2e_unique_bufsum_cints_vs_v3(erepcints,erepv3); cout << "Testing sums of unique Cints' ERIs (from GRTCints) against IntV3's" << endl; compare_2e_unique_bufsum_cints_vs_v3(grtcints,erepv3); cout << "Printing GRT integrals" << endl; print_grt_ints(grtcints); #endif // tim->print(); return 0; } void compare_1e_cints_vs_v3(Ref& obcints, Ref& obv3) { Ref basis = obcints->basis(); for (int sh1=4; sh1nshell(); sh1++) for (int sh2=0; sh2nshell(); sh2++) { int nbf2 = basis->shell(sh2).nfunction(); obv3->compute_shell(sh1,sh2); obcints->compute_shell(sh1,sh2); const double *buffercints = obcints->buffer(); const double *bufferv3 = obv3->buffer(); int bf1_offset = 0; for (int gc1=0; gc1shell(sh1).ncontraction(); gc1++) { int am1 = basis->shell(sh1).am(gc1); CartesianIterCints citer1(am1); CartesianIterV3 iter1(am1); for ( citer1.start(); int(citer1) ; citer1.next() ) { int bf1cints = bf1_offset + citer1.bfn(); int bf1v3; for( iter1.start(); int(iter1) ; iter1.next() ) { if (iter1.a() == citer1.a() && iter1.b() == citer1.b() && iter1.c() == citer1.c()) { bf1v3 = bf1_offset + iter1.bfn(); break; } } int bf2_offset = 0; for (int gc2=0; gc2shell(sh2).ncontraction(); gc2++) { int am2 = basis->shell(sh2).am(gc2); CartesianIterCints citer2(am2); CartesianIterV3 iter2(am2); for ( citer2.start(); int(citer2) ; citer2.next() ) { int bf2cints = bf2_offset + citer2.bfn(); int bf2v3; for( iter2.start(); int(iter2) ; iter2.next() ) { if (iter2.a() == citer2.a() && iter2.b() == citer2.b() && iter2.c() == citer2.c()) { bf2v3 = bf2_offset + iter2.bfn(); break; } } double valuecints = buffercints[bf1cints*nbf2 + bf2cints]; double valuev3 = bufferv3[bf1v3*nbf2 + bf2v3]; if (fabs(valuecints-valuev3) > 1E-13) { cout << scprintf("Discrepancy in OEInt(sh1 = %d, sh2 = %d)\n",sh1,sh2); cout << scprintf("bf1 = %d bf2 = %d OEIntegral(cints) = %20.15lf\n",bf1cints,bf2cints,valuecints); cout << scprintf("bf1 = %d bf2 = %d OEIntegral(V3) = %20.15lf\n\n",bf1v3,bf2v3,valuev3); } } bf2_offset += basis->shell(sh2).nfunction(gc2); } } bf1_offset += basis->shell(sh1).nfunction(gc1); } } } void compare_2e_cints_vs_v3(Ref& tbcints, Ref& tbv3) { const double *buffercints = tbcints->buffer(); const double *bufferv3 = tbv3->buffer(); Ref basis = tbcints->basis(); for (int sh1=0; sh1nshell(); sh1++) for (int sh2=0; sh2nshell(); sh2++) for (int sh3=0; sh3nshell(); sh3++) for (int sh4=0; sh4nshell(); sh4++) { //sh1=0;sh2=0;sh3=8;sh4=3; //cout << scprintf("Computing TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4); tbv3->compute_shell(sh1,sh2,sh3,sh4); tbcints->compute_shell(sh1,sh2,sh3,sh4); int nbf2 = basis->shell(sh2).nfunction(); int nbf3 = basis->shell(sh3).nfunction(); int nbf4 = basis->shell(sh4).nfunction(); int bf1_offset = 0; for(int gc1=0;gc1shell(sh1).ncontraction(); gc1++) { int am1 = basis->shell(sh1).am(gc1); CartesianIterCints citer1(am1); CartesianIterV3 iter1(am1); for ( citer1.start(); int(citer1) ; citer1.next() ) { int bf1cints = citer1.bfn(); int bf1v3; for( iter1.start(); int(iter1) ; iter1.next() ) { if (iter1.a() == citer1.a() && iter1.b() == citer1.b() && iter1.c() == citer1.c()) { bf1v3 = iter1.bfn(); break; } } bf1cints += bf1_offset; bf1v3 += bf1_offset; int bf2_offset = 0; for(int gc2=0;gc2shell(sh2).ncontraction(); gc2++) { int am2 = basis->shell(sh2).am(gc2); CartesianIterCints citer2(am2); CartesianIterV3 iter2(am2); for ( citer2.start(); int(citer2) ; citer2.next() ) { int bf2cints = citer2.bfn(); int bf2v3; for( iter2.start(); int(iter2) ; iter2.next() ) { if (iter2.a() == citer2.a() && iter2.b() == citer2.b() && iter2.c() == citer2.c()) { bf2v3 = iter2.bfn(); break; } } bf2cints += bf2_offset; bf2v3 += bf2_offset; int bf3_offset = 0; for(int gc3=0;gc3shell(sh3).ncontraction(); gc3++) { int am3 = basis->shell(sh3).am(gc3); CartesianIterCints citer3(am3); CartesianIterV3 iter3(am3); for ( citer3.start(); int(citer3) ; citer3.next() ) { int bf3cints = citer3.bfn(); int bf3v3; for( iter3.start(); int(iter3) ; iter3.next() ) { if (iter3.a() == citer3.a() && iter3.b() == citer3.b() && iter3.c() == citer3.c()) { bf3v3 = iter3.bfn(); break; } } bf3cints += bf3_offset; bf3v3 += bf3_offset; int bf4_offset = 0; for(int gc4=0;gc4shell(sh4).ncontraction(); gc4++) { int am4 = basis->shell(sh4).am(gc4); CartesianIterCints citer4(am4); CartesianIterV3 iter4(am4); for ( citer4.start(); int(citer4) ; citer4.next() ) { int bf4cints = citer4.bfn(); int bf4v3; for( iter4.start(); int(iter4) ; iter4.next() ) { if (iter4.a() == citer4.a() && iter4.b() == citer4.b() && iter4.c() == citer4.c()) { bf4v3 = iter4.bfn(); break; } } bf4cints += bf4_offset; bf4v3 += bf4_offset; double valuecints = buffercints[((bf1cints*nbf2 + bf2cints)*nbf3 + bf3cints)*nbf4 + bf4cints]; double valuev3 = bufferv3[((bf1v3*nbf2 + bf2v3)*nbf3 + bf3v3)*nbf4 + bf4v3]; if (fabs(valuecints-valuev3) > 1E-12) { cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4); cout << scprintf("bf1 = %d bf2 = %d bf3 = %d bf4 = %d TEIntegral(cints) = %20.15lf\n", bf1cints,bf2cints,bf3cints,bf4cints,valuecints); cout << scprintf("bf1 = %d bf2 = %d bf3 = %d bf4 = %d TEIntegral(V3) = %20.15lf\n\n", bf1v3,bf2v3,bf3v3,bf4v3,valuev3); } } bf4_offset+=basis->shell(sh4).nfunction(gc4); } } bf3_offset+=basis->shell(sh3).nfunction(gc3); } } bf2_offset+=basis->shell(sh2).nfunction(gc2); } } bf1_offset+=basis->shell(sh1).nfunction(gc1); } //return; } } void compare_2e_puream_cints_vs_v3(Ref& tbcints, Ref& tbv3) { const double *buffercints = tbcints->buffer(); const double *bufferv3 = tbv3->buffer(); Ref basis = tbcints->basis(); for (int sh1=0; sh1nshell(); sh1++) for (int sh2=0; sh2nshell(); sh2++) for (int sh3=0; sh3nshell(); sh3++) for (int sh4=0; sh4nshell(); sh4++) { // sh1 = 0; sh2 = 0; sh3 = 6; sh4 = 13; /* if ( !((basis->shell(sh1).has_pure() || basis->shell(sh1).max_am()==0) && (basis->shell(sh2).has_pure() || basis->shell(sh2).max_am()==0) && (basis->shell(sh3).has_pure() || basis->shell(sh3).max_am()==0) && (basis->shell(sh4).has_pure() || basis->shell(sh4).max_am()==0)) ) continue;*/ tbcints->compute_shell(sh1,sh2,sh3,sh4); tbv3->compute_shell(sh1,sh2,sh3,sh4); int nbf1 = basis->shell(sh1).nfunction(); int nbf2 = basis->shell(sh2).nfunction(); int nbf3 = basis->shell(sh3).nfunction(); int nbf4 = basis->shell(sh4).nfunction(); cout << scprintf("Computing TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4); cout << scprintf("size = %d\n",nbf1*nbf2*nbf3*nbf4); for(int ijkl=0;ijkl 1E-13) { cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4); cout << scprintf("1234 = %d TEIntegral(cints) = %20.15lf\n", ijkl,valuecints); cout << scprintf("TEIntegral(V3) = %20.15lf\n\n", valuev3); } } // return; } } void compare_2e_bufsum_cints_vs_v3(Ref& tbcints, Ref& tbv3) { Ref basis = tbcints->basis(); const double *buffercints = tbcints->buffer(); const double *bufferv3 = tbv3->buffer(); for (int sh1=0; sh1nshell(); sh1++) for (int sh2=0; sh2nshell(); sh2++) for (int sh3=0; sh3nshell(); sh3++) for (int sh4=0; sh4nshell(); sh4++) { // sh1=12;sh2=12;sh3=12;sh4=12; tbcints->compute_shell(sh1,sh2,sh3,sh4); tbv3->compute_shell(sh1,sh2,sh3,sh4); int nbf1 = basis->shell(sh1).nfunction(); int nbf2 = basis->shell(sh2).nfunction(); int nbf3 = basis->shell(sh3).nfunction(); int nbf4 = basis->shell(sh4).nfunction(); double sum_cints = 0.0; double sum_v3 = 0.0; int index = 0; for (int i=0; i 1E-12) { cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4); cout << scprintf("TEIntegralSum(cints) = %20.15lf\n", sum_cints); cout << scprintf("TEIntegralSum(V3) = %20.15lf\n\n", sum_v3); } //return; } } void compare_2e_unique_bufsum_cints_vs_v3(Ref& tbcints, Ref& tbv3) { Ref basis = tbcints->basis(); const double *buffercints = tbcints->buffer(); const double *bufferv3 = tbv3->buffer(); for (int sh1=0; sh1nshell(); sh1++) for (int sh2=0; sh2<=sh1; sh2++) for (int sh3=0; sh3<=sh1; sh3++) for (int sh4=0; sh4 <= ((sh1==sh3) ? sh2 : sh3) ; sh4++) { tbcints->compute_shell(sh1,sh2,sh3,sh4); tbv3->compute_shell(sh1,sh2,sh3,sh4); int nbf1 = basis->shell(sh1).nfunction(); int nbf2 = basis->shell(sh2).nfunction(); int nbf3 = basis->shell(sh3).nfunction(); int nbf4 = basis->shell(sh4).nfunction(); double sum_cints = 0.0; double sum_v3 = 0.0; int e12 = (sh1 == sh2) ? 1 : 0; int e34 = (sh3 == sh4) ? 1 : 0; int e13e24 = (((sh1 == sh3)&&(sh2==sh4))||((sh1==sh4)&&(sh2==sh3))) ? 1 : 0; int index = 0; for (int i=0; i 1E-12) { cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4); cout << scprintf("TEIntegralSum(cints) = %20.15lf\n", sum_cints); cout << scprintf("TEIntegralSum(V3) = %20.15lf\n\n", sum_v3); } } } void print_grt_ints(Ref& tbcints) { Ref basis = tbcints->basis(); const double *buffer[4]; buffer[0] = tbcints->buffer(TwoBodyInt::eri); buffer[1] = tbcints->buffer(TwoBodyInt::r12); buffer[2] = tbcints->buffer(TwoBodyInt::r12t1); buffer[3] = tbcints->buffer(TwoBodyInt::r12t2); char teout_filename[] = "teout0.dat"; FILE *teout[4]; for(int te_type=0;te_type<4;te_type++) { teout_filename[5] = te_type + '0'; teout[te_type] = fopen(teout_filename,"w"); } for (int ush1=0; ush1nshell(); ush1++) for (int ush2=0; ush2<=ush1; ush2++) for (int ush3=0; ush3<=ush2; ush3++) for (int ush4=0; ush4 <=ush3 ; ush4++) { int S1[3], S2[3], S3[3], S4[4]; int num = 1; S1[0] = ush1; S2[0] = ush2; S3[0] = ush3; S4[0] = ush4; if (ush1==ush2 && ush1==ush3 || ush2==ush3 && ush2==ush4) num=1; else if (ush1==ush3 || ush2==ush4) { num = 2; S1[1] = ush1; S2[1] = ush3; S3[1] = ush2; S4[1] = ush4; } else if (ush2==ush3) { num = 2; S1[1] = ush1; S2[1] = ush4; S3[1] = ush2; S4[1] = ush3; } else if (ush1==ush2 || ush3==ush4) { num = 2; S1[1] = ush1; S2[1] = ush3; S3[1] = ush2; S4[1] = ush4; } else { num = 3; S1[1] = ush1; S2[1] = ush3; S3[1] = ush2; S4[1] = ush4; S1[2] = ush1; S2[2] = ush4; S3[2] = ush2; S4[2] = ush3; } for(int uq=0;uqcompute_shell(sh1,sh2,sh3,sh4); int nbf1 = basis->shell(sh1).nfunction(); int nbf2 = basis->shell(sh2).nfunction(); int nbf3 = basis->shell(sh3).nfunction(); int nbf4 = basis->shell(sh4).nfunction(); int e12 = (sh1 == sh2) ? 1 : 0; int e34 = (sh3 == sh4) ? 1 : 0; int e13e24 = (((sh1 == sh3)&&(sh2==sh4))||((sh1==sh4)&&(sh2==sh3))) ? 1 : 0; int index = 0; for (int i=0; i 1E-15) { fprintf(teout[0], "%5d%5d%5d%5d%20.10lf\n", basis->shell_to_function(sh1) + i+1, basis->shell_to_function(sh2) + j+1, basis->shell_to_function(sh3) + k+1, basis->shell_to_function(sh4) + l+1, integral); } integral = buffer[1][index]; if (fabs(integral) > 1E-15) { fprintf(teout[1], "%5d%5d%5d%5d%20.10lf\n", basis->shell_to_function(sh1) + i+1, basis->shell_to_function(sh2) + j+1, basis->shell_to_function(sh3) + k+1, basis->shell_to_function(sh4) + l+1, integral); } integral = buffer[2][index]; if (fabs(integral) > 1E-15) { fprintf(teout[2], "%5d%5d%5d%5d%20.10lf\n", basis->shell_to_function(sh1) + i+1, basis->shell_to_function(sh2) + j+1, basis->shell_to_function(sh3) + k+1, basis->shell_to_function(sh4) + l+1, integral); } integral = buffer[3][index]; if (fabs(integral) > 1E-15) { fprintf(teout[3], "%5d%5d%5d%5d%20.10lf\n", basis->shell_to_function(sh1) + i+1, basis->shell_to_function(sh2) + j+1, basis->shell_to_function(sh3) + k+1, basis->shell_to_function(sh4) + l+1, integral); } index++; } } } } } } for(int te_type=0;te_type<4;te_type++) fclose(teout[te_type]); } void compare_2e_permute(Ref& cints) { Ref tb1 = cints->electron_repulsion(); Ref tb2 = cints->electron_repulsion(); Ref basis = tb1->basis(); const double *buffer1 = tb1->buffer(); const double *buffer2 = tb2->buffer(); int sh1 = 0; int sh2 = 0; int sh3 = 4; int sh4 = 0; tb1->compute_shell(sh1,sh2,sh3,sh4); tb2->compute_shell(sh1,sh2,sh4,sh3); int nbf1 = basis->shell(sh1).nfunction(); int nbf2 = basis->shell(sh2).nfunction(); int nbf3 = basis->shell(sh3).nfunction(); int nbf4 = basis->shell(sh4).nfunction(); for(int index = 0; index 1E-13) { cout << scprintf("Discrepancy in TEInt(sh1 = %d, sh2 = %d, sh3 = %d, sh4 = %d)\n",sh1,sh2,sh3,sh4); cout << scprintf("TEIntegral(cints1) = %20.15lf\n",buffer1[index]); cout << scprintf("TEIntegral(cints2) = %20.15lf\n\n",buffer2[index]); } } void do_shell_test_1e(const Ref &int1ev3, void (Int1eV3::*int_shell_1e)(int,int), int permute, int i, int j, int na, int nb, double *buf, double *pbuf) { int ii = 0; int a; double *buffer = int1ev3->buffer(); (int1ev3->*int_shell_1e)(i, j); for (a=0; a*int_shell_1e)(j, i); for (a=0; a 1.0e-13) { cout << scprintf("----- 1e perm failed:" "<%d %d|%d %d>:" " %18.14f != %18.14f " "<%d %d|%d %d>\n", i, a, j, b, buf[ii], pbuf[a + na*b], j, b, i, a); } if (fabs(buf[ii]) > 1.0e-15) { cout << scprintf(" <(%d %d)|(%d %d)> = %15.11f\n", i,a,j,b, buf[ii]); } ii++; } } } void test_int_shell_1e(const Ref& keyval, const Ref &int1ev3, void (Int1eV3::*int_shell_1e)(int,int), int permute) { int flags = 0; Ref basis = int1ev3->basis(); int maxfunc = basis->max_nfunction_in_shell(); int size = maxfunc * maxfunc; double *buf = new double[size]; double *pbuf = new double[size]; int nshell = int1ev3->basis()->nshell(); for (int i=0; ishell(i).nfunction(); for (int j=0; jshell(j).nfunction(); do_shell_test_1e(int1ev3, int_shell_1e, permute, i, j, na, nb, buf, pbuf); } } delete[] buf; delete[] pbuf; } void test_3_center(const Ref& keyval, const Ref &int2ev3) { int ii, i,j,k,l,m,n; int2ev3->set_redundant(1); int2ev3->set_permute(0); double *buffer = int2ev3->buffer(); int nshell = int2ev3->basis()->nshell(); for (i=0; ierep_2center(sh,sizes); ii = 0; for (k=0; k1.0e-15) cout << scprintf(" ((%d %d)|(%d %d)) = %15.11f\n", sh[0],k,sh[1],l, buffer[ii]); ii++; } } } } for (i=0; ierep_3center(sh,sizes); ii = 0; for (k=0; k1.0e-15) cout << scprintf( " ((%d %d)|(%d %d)(%d %d)) = %15.11f\n", sh[0],k,sh[1],l,sh[2],n, buffer[ii]); ii++; } } } } } } } void init_shell_perm(const Ref &int2ev3, double *integrals, double buff[maxint][maxint][maxint][maxint], int sh[4], int sizes[4]) { int i, j, k, l; int oldp = int2ev3->permute(); int2ev3->set_permute(0); int2ev3->erep(sh, sizes); int2ev3->set_permute(oldp); for (i=0; i &int2ev3, double *integrals, double buff[maxint][maxint][maxint][maxint], int sh[4], int sizes[4], int p0, int p1, int p2, int p3) { int ip[4], p[4]; int psizes[4]; int psh[4]; int index = 0; int i[4]; p[0] = p0; p[1] = p1; p[2] = p2; p[3] = p3; ip[p0] = 0; ip[p1] = 1; ip[p2] = 2; ip[p3] = 3; psh[0] = sh[p0]; psh[1] = sh[p1]; psh[2] = sh[p2]; psh[3] = sh[p3]; int oldp = int2ev3->permute(); int2ev3->set_permute(0); int2ev3->erep(psh, psizes); int2ev3->set_permute(oldp); for (i[0]=0; i[0] 1.0e-13) { cout << scprintf("perm %d %d %d %d failed:" "((%d %d)(%d %d)|(%d %d)(%d %d)):" " %18.14f != %18.14f " "((%d %d)(%d %d)|(%d %d)(%d %d))\n", p0, p1, p2, p3, sh[0],i[0], sh[1],i[1], sh[2],i[2], sh[3],i[3], buff[i[ip[0]]][i[ip[1]]][i[ip[2]]][i[ip[3]]], integrals[index], psh[0],i[p[0]], psh[1],i[p[1]], psh[2],i[p[2]], psh[3],i[p[3]]); } index++; } } } } } void do_shell_quartet_test(const Ref &int2ev3, int print, int printbounds, int bounds, int permute, const Ref& keyval, int i, int j, int k, int l) { int sh[4], sizes[4]; int ibuf; int ii, jj, kk, ll; sh[0] = i; sh[1] = j; sh[2] = k; sh[3] = l; double maxintegral, integralbound; int boundijkl; if (bounds) { integralbound = int2ev3->logbound_to_bound( (boundijkl = int2ev3->erep_4bound(i,j,k,l)) ); } double *buffer = int2ev3->buffer(); int2ev3->erep(sh,sizes); ibuf = 0; maxintegral = 0.0; for (ii=0; ii maxintegral) { maxintegral = absint; } if (bounds && absint > integralbound) { cout << scprintf("((%d %d)(%d %d)|(%d %d)(%d %d)) = %15.11f, " "bound = %15.11f\n", sh[0], ii, sh[1], jj, sh[2], kk, sh[3], ll, buffer[ibuf], integralbound); abort(); } if (print && (absint > 1.0e-9 || (bounds && integralbound > 1.0e-9))) { cout << scprintf(" ((%d %d)(%d %d)|(%d %d)(%d %d))" " = %15.11f", sh[0],ii, sh[1],jj, sh[2],kk, sh[3],ll, buffer[ibuf]); if (bounds) { cout << scprintf(" (%2d%% of bound)", (int)(100*(absint/integralbound))); } cout << scprintf("\n"); } ibuf++; } } } } if (permute) { double buff1[maxint][maxint][maxint][maxint]; sh[0] = i; sh[1] = j; sh[2] = k; sh[3] = l; init_shell_perm(int2ev3, buffer, buff1, sh, sizes); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 0, 1, 2, 3); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 1, 0, 2, 3); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 0, 1, 3, 2); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 1, 0, 3, 2); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 2, 3, 0, 1); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 2, 3, 1, 0); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 3, 2, 0, 1); check_shell_perm(int2ev3, buffer, buff1, sh, sizes, 3, 2, 1, 0); } if (bounds) { int boundij = int2ev3->erep_4bound(i,j,-1,-1); int boundkl = int2ev3->erep_4bound(-1,-1,k,l); int badbound = 0; if (boundij < boundijkl || boundkl < boundijkl) { badbound = 1; } if (badbound || printbounds) { cout << scprintf("max(%d,%d,%d,%d)=%7.4f, bnd=%7.4f, " "bnd(%d,%d,*,*)=%7.4f, bnd(*,*,%d,%d)=%7.4f\n", i, j, k, l, maxintegral, integralbound, i,j, int2ev3->logbound_to_bound(boundij), k,l, int2ev3->logbound_to_bound(boundkl)); } if (badbound) { cout << scprintf("ERROR: bad bound\n"); abort(); } } } void do_4_center_test(const Ref &int2ev3, int print, int printbounds, int bounds, int permute, const Ref& keyval) { int ii,jj,kk,ll, i,j,k,l, ibuf; int nshell = int2ev3->basis()->nshell(); int unique = keyval->booleanvalue("unique"); int timestats = keyval->booleanvalue("timestats"); Ref timer = new RegionTimer(); if (!timestats) { for (i=0; iget_cpu_time(); for (j=0; j= nshell) ish = nshell-1; if (jsh >= nshell) jsh = nshell-1; if (ksh >= nshell) ksh = nshell-1; if (lsh >= nshell) lsh = nshell-1; sh[0] = ish; sh[1] = jsh; sh[2] = ksh; sh[3] = lsh; int2ev3->erep(sh,sizes); } double t2 = timer->get_cpu_time(); times[i] = t2-t1; } double ave = 0.0; for (i=0; i& keyval, const Ref &int2ev3) { int i; cout << scprintf("4 center test:\n"); cout << scprintf(" on entry int2ev3 used %d bytes\n", int2ev3->used_storage()); int2ev3->set_permute(0); int2ev3->set_redundant(1); int storage = keyval->intvalue("storage") - int2ev3->used_storage(); if (storage < 0) storage = 0; if (keyval->booleanvalue("store_integrals")) storage = 0; int niter = keyval->intvalue("niter"); int print = keyval->booleanvalue("print"); int bounds = keyval->booleanvalue("bounds"); int permute = keyval->booleanvalue("permute"); int printbounds = keyval->booleanvalue("printbounds"); cout << scprintf(" storage = %d\n", storage); cout << scprintf(" niter = %d\n", niter); cout << scprintf(" print = %d\n", print); cout << scprintf(" bounds = %d\n", bounds); cout << scprintf(" permute = %d\n", permute); cout << scprintf("printbounds = %d\n", printbounds); if (bounds) int2ev3->init_bounds(); int2ev3->init_storage(storage); for (i=0; icount("quartet") == 4) { do_shell_quartet_test(int2ev3, print, printbounds, bounds, permute, keyval, keyval->intvalue("quartet", 0), keyval->intvalue("quartet", 1), keyval->intvalue("quartet", 2), keyval->intvalue("quartet", 3)); } int2ev3->done_storage(); int2ev3->done_bounds(); } void do_shell_quartet_der_test(const Ref &int2ev3, double* buffer, int print, int printbounds, int bounds, int permute, const Ref& keyval, int i, int j, int k, int l) { int ii,jj,kk,ll, ibuf, ider, xyz; der_centersv3_t dercenters; int sh[4], sizes[4]; sh[0] = i; sh[1] = j; sh[2] = k; sh[3] = l; double maxintegral = 0.0, integralbound; int boundijkl; if (bounds) { integralbound = int2ev3->logbound_to_bound( (boundijkl = int2ev3->erep_4bound_1der(i,j,k,l)) ); } int2ev3->erep_all1der(sh,sizes,&dercenters); ibuf = 0; for (ider=0; ider maxintegral) { maxintegral = absint; } if (bounds && absint > integralbound) { cout << scprintf("((%d %d)(%d %d)|(%d %d)(%d %d))" " = %15.11f, bound = %15.11f\n", sh[0], ii, sh[1], jj, sh[2], kk, sh[3], ll, buffer[ibuf], integralbound); abort(); } if (print && absint > 1.0e-15) { cout << scprintf(" ((%d %d)(%d %d)" "|(%d %d)(%d %d))(%d %d)" " = %15.11f\n", sh[0],ii, sh[1],jj, sh[2],kk, sh[3],ll, dercenters.num[ider], xyz, buffer[ibuf] ); } ibuf++; } } } } } } if (bounds) { int boundij = int2ev3->erep_4bound_1der(i,j,-1,-1); int boundkl = int2ev3->erep_4bound_1der(-1,-1,k,l); int badbound = 0; if (boundij < boundijkl || boundkl < boundijkl) { badbound = 1; } if (badbound || printbounds) { cout << scprintf("max(%d,%d,%d,%d)=%7.4f, bnd=%7.4f, " "bnd(%d,%d,*,*)=%8.4f, bnd(*,*,%d,%d)=%8.4f\n", i, j, k, l, maxintegral, integralbound, i,j, int2ev3->logbound_to_bound(boundij), k,l, int2ev3->logbound_to_bound(boundkl)); } if (badbound) { cout << scprintf("ERROR: bad bound\n"); abort(); } } } void do_test_4der_center(const Ref &int2ev3, double* buffer, int print, int printbounds, int bounds, int permute, const Ref& keyval) { int i,j,k,l; int nshell = int2ev3->basis()->nshell(); for (i=0; i& keyval, const Ref &int2ev3) { int i; int2ev3->set_permute(0); int2ev3->set_redundant(1); double *buffer = int2ev3->buffer(); int niter = keyval->intvalue("niter"); int print = keyval->booleanvalue("print"); int bounds = keyval->booleanvalue("bounds"); int printbounds = keyval->booleanvalue("printbounds"); int permute = keyval->booleanvalue("permute"); cout << scprintf("4 center derivative test:\n"); cout << scprintf(" niter = %d\n", niter); cout << scprintf(" print = %d\n", print); cout << scprintf(" bounds = %d\n", bounds); cout << scprintf("printbounds = %d\n", printbounds); cout << scprintf(" permute = %d\n", permute); if (bounds) int2ev3->init_bounds_1der(); for (i=0; icount("quartet") == 4) { do_shell_quartet_der_test(int2ev3, buffer, print, printbounds, bounds, permute, keyval, keyval->intvalue("quartet", 0), keyval->intvalue("quartet", 1), keyval->intvalue("quartet", 2), keyval->intvalue("quartet", 3)); } if (bounds) int2ev3->done_bounds_1der(); } ///////////////////////////////////////////////////////////////////////////// // Local Variables: // mode: c++ // c-file-style: "CLJ" // End: