cintstest.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

//
// cintstest.cc
//
// Copyright (C) 2001 Edward Valeev
//
// Author: Edward Valeev <[e-mail]>
// 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 <stdlib.h>
#include <string.h>

#include <util/misc/formio.h>
#include <util/misc/regtime.h>
#include <util/keyval/keyval.h>
#include <util/group/message.h>
#include <util/group/pregtime.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/intv3/int1e.h>
#include <chemistry/qc/intv3/int2e.h>
#include <chemistry/qc/intv3/intv3.h>
#include <chemistry/qc/intv3/cartitv3.h>
#include <chemistry/qc/cints/cints.h>
#include <chemistry/qc/cints/int2e.h>
#include <chemistry/qc/cints/cartit.h>


using namespace std;
using namespace sc;

#define CINTS

void compare_1e_cints_vs_v3(Ref<OneBodyInt>& obcints, Ref<OneBodyInt>& obv3);
void compare_2e_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
void compare_2e_puream_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
void compare_2e_bufsum_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
void compare_2e_unique_bufsum_cints_vs_v3(Ref<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3);
void print_grt_ints(Ref<TwoBodyInt>& tbcints);
void compare_2e_permute(Ref<Integral>& cints);
void test_int_shell_1e(const Ref<KeyVal>&, const Ref<Int1eV3> &int1ev3,
                       void (Int1eV3::*int_shell_1e)(int,int),
                       int permute);
void test_3_center(const Ref<KeyVal>&, const Ref<Int2eV3> &);
void test_4_center(const Ref<KeyVal>& keyval, const Ref<Int2eV3> &int2ev3);
void test_4der_center(const Ref<KeyVal>&, const Ref<Int2eV3> &int2ev3);

#define maxint 9

void
testint(const Ref<OneBodyInt>& in)
{
  if (in.null()) {
      cout << "null integral generator" << endl;
      abort();
    }
}

void
testint(const Ref<OneBodyDerivInt>& in)
{
  if (in.null()) {
      cout << "null integral generator" << endl;
      abort();
    }
}

void
testint(const Ref<TwoBodyInt>& in)
{
  if (in.null()) {
      cout << "null integral generator" << endl;
      abort();
    }
}

void
testint(const Ref<TwoBodyDerivInt>& in)
{
  if (in.null()) {
      cout << "null integral generator" << endl;
      abort();
    }
}

/*void
do_bounds_stats(const Ref<KeyVal>& keyval,
                const Ref<Int2eV3> &int2ev3)
{
  int i,j;
  int nshell = int2ev3->basis()->nshell();
  int eps = -10;
  int *nonzero = new int[nshell];
  for (i=0; i<nshell; i++) {
      if (i==0) nonzero[i] = 0;
      else nonzero[i] = nonzero[i-1];
      for (j=0; j<=i; j++) {
          if (int2ev3->erep_4bound(i,j,-1,-1) > eps) {
              nonzero[i]++;
            }
        }
    }
  for (i=0; i<nshell; i++) {
      int natom = 1 + int2ev3->basis()->shell_to_center(i);
      int npq = (i*(i+1))/2;
      cout<<scprintf("nsh=%2d nat=%2d npq=%4d npq>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<MessageGrp> msg = MessageGrp::initial_messagegrp(argc,argv);
  if (msg.null()) msg = new ProcMessageGrp();
  MessageGrp::set_default_messagegrp(msg);

  Ref<RegionTimer> 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<KeyVal> pkv(new ParsedKeyVal(infile));
  Ref<KeyVal> tkeyval(new PrefixKeyVal(":test", pkv));

  Ref<GaussianBasisSet> basis = require_dynamic_cast<GaussianBasisSet*>(
    tkeyval->describedclassvalue("basis").pointer(),"main\n");
  Ref<Molecule> 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<Integral> intgrlv3 = new IntegralV3(basis,basis,basis,basis);
  Ref<Int1eV3> int1ev3 = new Int1eV3(intgrlv3.pointer(),basis,basis,1);
  Ref<Int2eV3> 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> integral = new IntegralV3(basis);
#ifdef CINTS
  Ref<Integral> integralcints = new IntegralCints(basis);
#endif

  Ref<OneBodyInt> overlapv3 = integral->overlap();
  Ref<OneBodyInt> kineticv3 = integral->kinetic();
  Ref<OneBodyInt> nuclearv3 = integral->nuclear();
  Ref<OneBodyInt> hcorev3 = integral->hcore();

#ifdef CINTS
  Ref<OneBodyInt> overlapcints = integralcints->overlap();
  testint(overlapcints);
  Ref<OneBodyInt> kineticcints = integralcints->kinetic();
  testint(kineticcints);
  Ref<OneBodyInt> nuclearcints = integralcints->nuclear();
  testint(nuclearcints);
  Ref<OneBodyInt> hcorecints = integralcints->hcore();
  testint(hcorecints);
#endif

  Ref<TwoBodyInt> 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<TwoBodyInt> 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<TwoBodyInt> 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<OneBodyInt>& obcints, Ref<OneBodyInt>& obv3)
{
  Ref<GaussianBasisSet> basis = obcints->basis();
  for (int sh1=4; sh1<basis->nshell(); sh1++)
    for (int sh2=0; sh2<basis->nshell(); 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; gc1<basis->shell(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; gc2<basis->shell(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<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3)
{
  const double *buffercints = tbcints->buffer();
  const double *bufferv3 = tbv3->buffer();

  Ref<GaussianBasisSet> basis = tbcints->basis();
  for (int sh1=0; sh1<basis->nshell(); sh1++)
    for (int sh2=0; sh2<basis->nshell(); sh2++)
      for (int sh3=0; sh3<basis->nshell(); sh3++)
	for (int sh4=0; sh4<basis->nshell(); 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;gc1<basis->shell(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;gc2<basis->shell(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;gc3<basis->shell(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;gc4<basis->shell(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",