tformv3.cc

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

  int i;
  double *tmp, *tmp2;

  /* Allocate more temporary space if needed. */
  source_space(nsource);

  tmp = source;
  tmp2 = integrals;
  for (i=0; i<nsource; i++) *tmp++ = *tmp2++;
}

void
Int1eV3::do_transform_1e(Integral *integ,
                         double *integrals,
                         GaussianShell *sh1, GaussianShell *sh2,
                         int chunk)
{
  int i, j;
  int ogc1, ogc2;
  int ogc1pure, ogc2pure;
  int am1, am2;
  int pure1 = sh1->has_pure();
  int pure2 = sh2->has_pure();
  int ncart1 = sh1->ncartesian();
  int ncart2 = sh2->ncartesian();
  int nfunc1 = sh1->nfunction();
  int nfunc2 = sh2->nfunction();
  int nfunci, nfuncj;

  if (!pure1 && !pure2) return;

  /* Loop through the generalized general contractions,
   * transforming the first index. */
  if (pure1) {
      copy_to_source(integrals, ncart1*ncart2*chunk);
      memset(integrals, 0, sizeof(double)*sh1->nfunction()*ncart2*chunk);

      ogc1 = 0;
      ogc1pure = 0;
      for (i=0; i<sh1->ncontraction(); i++) {
          am1 = sh1->am(i);
          nfunci = sh1->nfunction(i);
          ogc2 = 0;
          for (j=0; j<sh2->ncontraction(); j++) {
              am2 = sh2->am(j);
              nfuncj = sh2->nfunction(j);

              if (sh1->is_pure(i)) {
                  SphericalTransformIter
                      trans(integ->spherical_transform(sh1->am(i)));
                  do_sparse_transform11(source, integrals, chunk,
                                        trans,
                                        ogc1,
                                        ogc1pure,
                                        INT_NCART(am2), ncart2, ogc2);
                }
              else {
                  do_copy1(source, integrals, chunk,
                           nfunci, nfunc1, ogc1pure,
                           INT_NCART(am2), ncart2, ogc2);
                }
              ogc2 += INT_NCART(am2);
            }
          ogc1 += INT_NCART(am1);
          ogc1pure += INT_NPURE(am1);
        }
    }

  if (pure2) {
      copy_to_source(integrals, nfunc1*ncart2*chunk);
      memset(integrals, 0,
             sizeof(double)*sh1->nfunction()*sh2->nfunction()*chunk);

      ogc1 = 0;
      for (i=0; i<sh1->ncontraction(); i++) {
          am1 = sh1->am(i);
          nfunci = sh1->nfunction(i);
          ogc2 = 0;
          ogc2pure = 0;
          for (j=0; j<sh2->ncontraction(); j++) {
              am2 = sh2->am(j);
              nfuncj = sh2->nfunction(j);

              if (sh2->is_pure(j)) {
                  SphericalTransformIter
                      trans(integ->spherical_transform(sh2->am(j)));
                  do_sparse_transform12(source, integrals, chunk,
                                        trans,
                                        INT_NPURE(am1), ogc1,
                                        ncart2, ogc2,
                                        sh2->nfunction(), ogc2pure);
                }
              else {
                  do_copy1(source, integrals, chunk,
                           nfunci, nfunc1, ogc1,
                           nfuncj, nfunc2, ogc2pure);
                }
              ogc2 += INT_NCART(am2);
              ogc2pure += INT_NPURE(am2);
            }
          ogc1 += INT_NPURE(am1);
        }
    }
}

/* it is ok for integrals and target to overlap */
void
Int1eV3::transform_1e(Integral *integ,
                      double *integrals, double *target,
                      GaussianShell *sh1, GaussianShell *sh2, int chunk)
{
  int ntarget;

  do_transform_1e(integ, integrals, sh1, sh2, chunk);

  /* copy the integrals to the target, if necessary */
  ntarget = sh1->nfunction() * sh2->nfunction();
  if (integrals != target) {
      memmove(target, integrals, ntarget*sizeof(double)*chunk);
    }
}

/* it is not ok for integrals and target to overlap */
void
Int1eV3::accum_transform_1e(Integral *integ,
                            double *integrals, double *target,
                            GaussianShell *sh1, GaussianShell *sh2, int chunk)
{
  int i, ntarget;

  do_transform_1e(integ, integrals, sh1, sh2, chunk);

  /* accum the integrals to the target */
  ntarget = sh1->nfunction() * sh2->nfunction() * chunk;
  for (i=0; i<ntarget; i++) target[i] += integrals[i];
}

void
Int1eV3::transform_1e(Integral*integ,
                            double *integrals, double *target,
                            GaussianShell *sh1, GaussianShell *sh2)
{
  transform_1e(integ, integrals, target, sh1, sh2, 1);
}

void
Int1eV3::accum_transform_1e(Integral*integ,
                                  double *integrals, double *target,
                                  GaussianShell *sh1, GaussianShell *sh2)
{
  accum_transform_1e(integ, integrals, target, sh1, sh2, 1);
}

void
Int1eV3::transform_1e_xyz(Integral*integ,
                          double *integrals, double *target,
                          GaussianShell *sh1, GaussianShell *sh2)
{
  transform_1e(integ, integrals, target, sh1, sh2, 3);
}

void
Int1eV3::accum_transform_1e_xyz(Integral*integ,
                                double *integrals, double *target,
                                GaussianShell *sh1, GaussianShell *sh2)
{
  accum_transform_1e(integ, integrals, target, sh1, sh2, 3);
}

void
Int2eV3::do_gencon_sparse_transform_2e(Integral*integ,
                                       double *integrals, double *target,
                                       int index,
                                       GaussianShell *sh1, GaussianShell *sh2,
                                       GaussianShell *sh3, GaussianShell *sh4)
{
  int ncart[4];
  int nfunc[4];
  int nfunci, nfuncj, nfunck, nfuncl;
  int ncarti, ncartj, ncartk, ncartl;
  int i, j, k, l;
  int ogccart[4];
  int ogcfunc[4];
  int am1, am2, am3, am4;

  int ntarget1;
  int ntarget2;
  int ntarget3;
  int ntarget4;
              
  int nsource1;
  int nsource2;
  int nsource3;
  int nsource4;

  int *ni = &ncarti;
  int *nj = &ncartj;
  int *nk = &ncartk;
  int *nl = &ncartl;

  int *ogc1 = &ogccart[0];
  int *ogc2 = &ogccart[1];
  int *ogc3 = &ogccart[2];
  int *ogc4 = &ogccart[3];

  GaussianShell *shell;

  int *tgencon;

  ncart[0] = sh1->ncartesian();
  ncart[1] = sh2->ncartesian();
  ncart[2] = sh3->ncartesian();
  ncart[3] = sh4->ncartesian();
  nfunc[0] = sh1->nfunction();
  nfunc[1] = sh2->nfunction();
  nfunc[2] = sh3->nfunction();
  nfunc[3] = sh4->nfunction();

  ntarget1 = ncart[0];
  ntarget2 = ncart[1];
  ntarget3 = ncart[2];
  ntarget4 = ncart[3];
  
  nsource1 = ncart[0];
  nsource2 = ncart[1];
  nsource3 = ncart[2];
  nsource4 = ncart[3];

  if (index >= 0) {
      ntarget1 = nfunc[0];
      if (index >= 1) {
          ntarget2 = nfunc[1];
          nsource1 = nfunc[0];
          ni = &nfunci;
          ogc1 = &ogcfunc[0];
          if (index >= 2) {
              ntarget3 = nfunc[2];
              nsource2 = nfunc[1];
              nj = &nfuncj;
              ogc2 = &ogcfunc[1];
              if (index >= 3) {
                  ntarget4 = nfunc[3];
                  nsource3 = nfunc[2];
                  nk = &nfunck;
                  ogc3 = &ogcfunc[2];
                }
            }
        }
    }

  switch (index) {
  case 0:
      shell = sh1;
      tgencon = &i;
      break;
  case 1:
      shell = sh2;
      tgencon = &j;
      break;
  case 2:
      shell = sh3;
      tgencon = &k;
      break;
  case 3:
      shell = sh4;
      tgencon = &l;
      break;
  default:
      shell = 0;
      tgencon = 0;
      break;
    }

#if PRINT
    {
      double *tmp = integrals;
      ExEnv::outn() << scprintf("Before transform of index %d (%dx%dx%dx%d)\n",
             index, nsource1, nsource2, nsource3, nsource4);
      for (i=0; i<nsource1; i++) {
          for (j=0; j<nsource2; j++) {
              for (k=0; k<nsource3; k++) {
                  for (l=0; l<nsource4; l++) {
                      if (fabs(*tmp)>1.e-15) {
                          ExEnv::outn() << scprintf("(%d %d|%d %d) = %15.11lf\n",i,j,k,l,*tmp);
                        }
                      tmp++;
                    }
                }
            }
        }
    }
#endif

  copy_to_source(integrals, nsource1*nsource2*nsource3*nsource4);
  memset(target, 0, sizeof(double)*ntarget1*ntarget2*ntarget3*ntarget4);

  ogccart[0] = 0;
  ogcfunc[0] = 0;
  for (i=0; i<sh1->ncontraction(); i++) {
      am1 = sh1->am(i);
      nfunci = sh1->nfunction(i);
      ncarti = INT_NCART(am1);
      ogccart[1] = 0;
      ogcfunc[1] = 0;
      for (j=0; j<sh2->ncontraction(); j++) {
          am2 = sh2->am(j);
          nfuncj = sh2->nfunction(j);
          ncartj = INT_NCART(am2);
          ogccart[2] = 0;
          ogcfunc[2] = 0;
          for (k=0; k<sh3->ncontraction(); k++) {
              am3 = sh3->am(k);
              nfunck = sh3->nfunction(k);
              ncartk = INT_NCART(am3);
              ogccart[3] = 0;
              ogcfunc[3] = 0;
              for (l=0; l<sh4->ncontraction(); l++) {
                  am4 = sh4->am(l);
                  nfuncl = sh4->nfunction(l);
                  ncartl = INT_NCART(am4);

                  if (shell->is_pure(*tgencon)) {
                      SphericalTransformIter
                        trans(integ->spherical_transform(shell->am(*tgencon)));
                      do_sparse_transform2(source, target,
                                           index, trans,
                                           ncart[index], nfunc[index],
                                           ogccart[index], ogcfunc[index],
                                           *ni, nsource1, *ogc1,
                                           *nj, nsource2, *ogc2,
                                           *nk, nsource3, *ogc3,
                                           *nl, nsource4, *ogc4);
                    }
                  else {
                      do_copy2(source, integrals,
                               *ni, nsource1, *ogc1,
                               *nj, nsource2, *ogc2,
                               *nk, nsource3, *ogc3,
                               *nl, nsource4, *ogc4);
                    }
                  ogccart[3] += ncartl;
                  ogcfunc[3] += nfuncl;
                }
              ogccart[2] += ncartk;
              ogcfunc[2] += nfunck;
            }
          ogccart[1] += ncartj;
          ogcfunc[1] += nfuncj;
        }
      ogccart[0] += ncarti;
      ogcfunc[0] += nfunci;
    }

  
#if PRINT
    {
      double *tmp = integrals;
      ExEnv::outn() << scprintf("After transform of index %d (%dx%dx%dx%d)\n",
             index, ntarget1, ntarget2, ntarget3, ntarget4);
      for (i=0; i<ntarget1; i++) {
          for (j=0; j<ntarget2; j++) {
              for (k=0; k<ntarget3; k++) {
                  for (l=0; l<ntarget4; l++) {
                      if (fabs(*tmp)>1.e-15) {
                          ExEnv::outn()
                            << scprintf("(%d %d|%d %d) = %15.11lf\n",
                                        i,j,k,l,*tmp);
                        }
                      tmp++;
                    }
                }
            }
        }
    }
#endif
}