cintstest.cc

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

					       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<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 = 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<nbf1*nbf2*nbf3*nbf4;ijkl++) {
	      double valuecints = buffercints[ijkl];
	      double valuev3 = bufferv3[ijkl];
	      if (fabs(valuecints-valuev3) > 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<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3)
{
  Ref<GaussianBasisSet> basis = tbcints->basis();
  const double *buffercints = tbcints->buffer();
  const double *bufferv3 = tbv3->buffer();

  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=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<nbf1; i++) {
	      for (int j=0; j<nbf2; j++) {
		for (int k=0; k<nbf3; k++) {
		  for (int l=0; l<nbf4; l++) {
		    sum_cints += buffercints[index];
		    sum_v3 += bufferv3[index];
		    /*		    cout << scprintf("index = %d TEIntegral(cints) = %20.15lf\n",
				     index,buffercints[index]);
		    cout << scprintf("index = %d TEIntegral(V3) = %20.15lf\n\n",
				     index,bufferv3[index]);
		    */
		    index++;
		  }
		}
	      }
	    }
	    
	    if (fabs(sum_cints-sum_v3) > 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<TwoBodyInt>& tbcints, Ref<TwoBodyInt>& tbv3)
{
  Ref<GaussianBasisSet> basis = tbcints->basis();
  const double *buffercints = tbcints->buffer();
  const double *bufferv3 = tbv3->buffer();

  for (int sh1=0; sh1<basis->nshell(); 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<nbf1; i++) {
	      int jmax = e12 ? i : nbf2-1;
	      for (int j=0; j<=jmax; j++) {
		int kmax = e13e24 ? i : nbf3-1;
		for (int k=0; k<=kmax; k++) {
		  int lmax = e34 ? ( (e13e24&&(i==k)) ? j : k) : ( (e13e24&&(i==k)) ? j : nbf4-1);
		  for (int l=0; l<=lmax; l++) {
		    sum_cints += buffercints[index];
		    sum_v3 += bufferv3[index];
		    /*		    cout << scprintf("index = %d TEIntegral(cints) = %20.15lf\n",
				     index,buffercints[index]);
		    cout << scprintf("index = %d TEIntegral(V3) = %20.15lf\n\n",
				     index,bufferv3[index]);
		    */
		    index++;
		  }
		}
	      }
	    }
	    
	    if (fabs(sum_cints-sum_v3) > 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<TwoBodyInt>& tbcints)
{
  Ref<GaussianBasisSet> 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; ush1<basis->nshell(); 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;uq<num;uq++) {
	      int sh1 = S1[uq];
	      int sh2 = S2[uq];
	      int sh3 = S3[uq];
	      int sh4 = S4[uq];

	      tbcints->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();
	    
	      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<nbf1; i++) {
		int jmax = e12 ? i : nbf2-1;
		for (int j=0; j<=jmax; j++) {
		  int kmax = e13e24 ? i : nbf3-1;
		  for (int k=0; k<=kmax; k++) {
		    int lmax = e34 ? ( (e13e24&&(i==k)) ? j : k) : ( (e13e24&&(i==k)) ? j : nbf4-1);
		    for (int l=0; l<=lmax; l++) {
		      
		      double integral = buffer[0][index];
		      if (fabs(integral) > 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<Integral>& cints)
{
  Ref<TwoBodyInt> tb1 = cints->electron_repulsion();
  Ref<TwoBodyInt> tb2 = cints->electron_repulsion();
  Ref<GaussianBasisSet> 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<nbf1*nbf2*nbf3*nbf4; index++)
    if (fabs(buffer1[index]-buffer2[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> &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<na*nb; a++) {
      buf[a] = buffer[a];
    }
  (int1ev3->*int_shell_1e)(j, i);
  for (a=0; a<na*nb; a++) {
      pbuf[a] = buffer[a];
    }
  for (a=0; a<na; a++) {
      for (int b=0; b<nb; b++) {
          if (fabs(buf[ii] - pbuf[a + na*b]) > 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>& keyval, const Ref<Int1eV3> &int1ev3,
                  void (Int1eV3::*int_shell_1e)(int,int),
                  int permute)
{
  int flags = 0;
  Ref<GaussianBasisSet> 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; i<nshell; i++) {
      int na = basis->shell(i).nfunction();
      for (int j=0; j<nshell; j++) {
          int nb = basis->shell(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>& keyval, const Ref<Int2eV3> &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; i<nshell; i++) {
      for (j=0; j<nshell; j++) {
          int sh[2], sizes[2];
          sh[0] = i;
          sh[1] = j;
          int2ev3->erep_2center(sh,sizes);
          ii = 0;
          for (k=0; k<sizes[0]; k++) {
              for (l=0; l<sizes[1]; l++) {
                  if (fabs(buffer[ii])>1.0e-15)
                      cout << scprintf(" ((%d %d)|(%d %d)) = %15.11f\n",
                                       sh[0],k,sh[1],l, buffer[ii]);
                  ii++;
                }
            }
        }
    }

  for (i=0; i<nshell; i++) {
      for (j=0; j<nshell; j++) {
          for (m=0; m<nshell; m++) {
              int sh[3], sizes[3];
              sh[0] = i;
              sh[1] = j;
              sh[2] = m;
              int2ev3->erep_3center(sh,sizes);
              ii = 0;
              for (k=0; k<sizes[0]; k++) {
                  for (l=0; l<sizes[1]; l++) {
                      for (n=0; n<sizes[2]; n++) {
                          if (fabs(buffer[ii])>1.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++;
                        }
                    }
                }
            }
        }
    }

}