build2e.cc

来自「大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CH」· CC 代码 · 共 1,562 行 · 第 1/4 页

  ASSIGN_BUILD(3,3,1,2)
  ASSIGN_BUILD(3,3,2,2)
  ASSIGN_BUILD(3,3,3,3)
#endif

  free(jmax_for_con);
  saved_am12 = am12;
  saved_am34 = am34;
  saved_ncon = nc1;

  used_storage_ += used_storage_build_;
#if CHECK_INTEGRAL_ALGORITHM
  ExEnv::outn() << "used_storage_build: " << used_storage_build_ << endl;
#endif
  }

void
Int2eV3::int_done_buildgc()
{
  int ci,cj,ck;

  used_storage_ -= used_storage_build_;
  used_storage_build_ = 0;

  free_store();

  for (ci=0; ci<saved_ncon; ci++) {
    for (cj=0; cj<saved_ncon; cj++) {
      for (ck=0; ck<saved_ncon; ck++) {
        delete[] e0f0_con_ints_array[ci][cj][ck];
        }
      delete[] e0f0_con_ints_array[ci][cj];
      }
    delete[] e0f0_con_ints_array[ci];
    }
  delete[] e0f0_con_ints_array;

  }

/* add_store maintains a list of free storage allocated by int_init_buildgc */
void
Int2eV3::add_store(void *p)
{
  if (!store) {
    store = (store_list_t*) malloc(sizeof(store_list_t));
    assert(store);
    store->p = 0;
    n_store_last = 0;
    }
  if (n_store_last >= STORAGE_CHUNK) {
    store_list_t* tmp = (store_list_t*) malloc(sizeof(store_list_t));
    assert(tmp);
    tmp->p = store;
    store = tmp;
    n_store_last = 0;
    }
  store->data[n_store_last++] = p;
  }

/* free_store frees the memory that add_store keeps track of */
void
Int2eV3::free_store()
{
  _free_store(store,n_store_last);
  store = 0;
  }

void
Int2eV3::_free_store(store_list_t* s, int n)
{
  int i;
  if (!s) return;
  for (i=0; i<n; i++) {
    free(s->data[i]);
    }
  _free_store(s->p,STORAGE_CHUNK);
  free(s);
  }


void
Int2eV3::int_buildgcam(int minam1, int minam2, int minam3, int minam4,
                       int maxam1, int maxam2, int maxam3, int maxam4,
                       int dam1, int dam2, int dam3, int dam4,
                       int sh1, int sh2, int sh3, int sh4,
                       int eAB)
{
  int k,m,n;
  int ci,cj,ck,cl;
  int maxam12,maxam34;
  int nc1,nc2,nc3,nc4;

  if (maxam1<0 || maxam2<0 || maxam3<0 || maxam4<0) return;
  if (minam1<0) minam1=0;
  if (minam2<0) minam2=0;
  if (minam3<0) minam3=0;
  if (minam4<0) minam4=0;

  maxam12 = maxam1 + maxam2;
  maxam34 = maxam3 + maxam4;

  /* Compute the offset shell numbers. */
  osh1 = sh1 + bs1_shell_offset_;
  if (!int_unit2) osh2 = sh2 + bs2_shell_offset_;
  osh3 = sh3 + bs3_shell_offset_;
  if (!int_unit4) osh4 = sh4 + bs4_shell_offset_;

  nc1 = bs1_->shell(sh1).ncontraction();
  if (int_unit2) nc2 = 1;
  else nc2 = bs2_->shell(sh2).ncontraction();
  nc3 = bs3_->shell(sh3).ncontraction();
  if (int_unit4) nc4 = 1;
  else nc4 = bs4_->shell(sh4).ncontraction();

  /* Zero the target contracted integrals that the build routine
   * will accumulate into. */
  for (m=minam1; m<=maxam12; m++) {
    for (n=minam3; n<=maxam34; n++) {
  int nm_cart = INT_NCART(m)*INT_NCART(n);
  for (ci=0; ci<nc1; ci++) {
    if (m < int_shell1->am(ci)+dam1) continue;
    for (cj=0; cj<nc2; cj++) {
      if (int_shell1->am(ci)+dam1+int_shell2->am(cj)+dam2 < m)
        continue;
      for (ck=0; ck<nc3; ck++) {
        if (n < int_shell3->am(ck)+dam3) continue;
        for (cl=0; cl<nc4; cl++) {
          if (int_shell3->am(ck)+dam3 +int_shell4->am(cl)+dam4 < n)
            continue;
          double *tmp = e0f0_con_ints_array[ci][cj][ck][cl](m,n);
          for (int ii=0; ii<nm_cart; ii++) tmp[ii] = 0.0;
          }
        }
      }
    }
      }
    }

  gen_shell_intermediates(sh1,sh2,sh3,sh4);

  /* If enough of the functions come from generalized contractions
   * to make it workwhile, then don't do redundant primitives
   * at the additional cost of slower normalization computations.
   */
  if (nc1 + nc2 + nc3 + nc4 > 4)
    build_using_gcs(nc1,nc2,nc3,nc4,
                    minam1,minam3,maxam12,maxam34,dam1,dam2,dam3,dam4,eAB);
  else
    build_not_using_gcs(nc1,nc2,nc3,nc4,
                    minam1,minam3,maxam12,maxam34,dam1,dam2,dam3,dam4,eAB);
  }

void
Int2eV3::build_not_using_gcs(int nc1, int nc2, int nc3, int nc4,
                             int minam1, int minam3, int maxam12, int maxam34,
                             int dam1, int dam2, int dam3, int dam4, int eAB)
{
  int i,j,k,l,m;
  int ci,cj,ck,cl;
  double *bufferprim;

#if 0
  ExEnv::outn() << scprintf("not_gcs: %d%d%d%d\n",
                   int_expweight1,
                   int_expweight2,
                   int_expweight3,
                   int_expweight4
      );
#endif

          /* Sum thru all possible contractions. */
  for (ci=0; ci<nc1; ci++) {
    int mlower = int_shell1->am(ci) + dam1;
    if (mlower < 0) continue;
    IntV3Arraydoublep2 ***e0f0_i = e0f0_con_ints_array[ci];
    for (cj=0; cj<nc2; cj++) {
      int mupper = mlower + int_shell2->am(cj) + dam2;
      if (mupper < mlower) continue;
      if (mlower < minam1) mlower = minam1;
      if (mupper > maxam12) mupper = maxam12;
      IntV3Arraydoublep2 **e0f0_ij = e0f0_i[cj];
      for (ck=0; ck<nc3; ck++) {
        int nlower = int_shell3->am(ck) + dam3;
        if (nlower < 0) continue;
        IntV3Arraydoublep2 *e0f0_ijk = e0f0_ij[ck];
        for (cl=0; cl<nc4; cl++) {
          int nupper = nlower + int_shell4->am(cl) + dam4;
          if (nupper < nlower) continue;
          if (nlower < minam3) nlower = minam3;
          if (nupper > maxam34) nupper = maxam34;

  /* Loop over the primitives. */
  for (i=0; i<int_shell1->nprimitive(); i++) {
    double coef0;
    coef0 = int_shell1->coefficient_unnorm(ci,i);
    if (int_expweight1) coef0 = coef0
                                    * int_shell1->exponent(i);
    /* This factor of two comes from the derivative integral formula. */
    if (int_expweight1) coef0 *= 2.0;
    if (int_expweight2) coef0 *= 2.0;
    if (int_expweight3) coef0 *= 2.0;
    if (int_expweight4) coef0 *= 2.0;
    if (int_store1) opr1 = int_shell_to_prim[osh1] + i;
    for (j=0; j<int_shell2->nprimitive(); j++) {
      double coef1;
      coef1 = int_shell2->coefficient_unnorm(cj,j);
      if (int_expweight2) coef1 *=  coef0
                                      * int_shell2->exponent(j);
      else                     coef1 *= coef0;
      if (int_store1) opr2 = int_shell_to_prim[osh2] + j;
      for (k=0; k<int_shell3->nprimitive(); k++) {
        double coef2;
        coef2 = int_shell3->coefficient_unnorm(ck,k);
        if (int_expweight3) coef2 *=  coef1
                                        * int_shell3->exponent(k);
        else                     coef2 *= coef1;
        if (int_store1) opr3 = int_shell_to_prim[osh3] + k;
        for (l=0; l<int_shell4->nprimitive(); l++) {
          double coef3;
          coef3 = int_shell4->coefficient_unnorm(cl,l);
          if (int_expweight4) coef3 *=  coef2
                                          * int_shell4->exponent(l);
          else                     coef3 *= coef2;
          if (int_store1) opr4 = int_shell_to_prim[osh4] + l;

          /* Produce the remaining intermediates. */
          gen_prim_intermediates_with_norm(i,j,k,l, maxam12+maxam34,coef3);

          /* Generate the target integrals. */
          if ((maxam12 == 0) && (maxam34 == 0)) {
            /* Do nothing: gen_prim_intermediates has set everything up. */
            }
          else if ((minam1<=MG)&&(minam3<=MG)&&(maxam12<=MG)&&(maxam34<=MG)) {
            if (build_routine[minam1]
                             [maxam12]
                             [minam3]
                             [maxam34][eAB]==&BuildIntV3::impossible_integral){
              ExEnv::errn() << scprintf("trying to build with int2v%d%d%d%d (exact)\n",
                      minam1,maxam12,minam3,maxam34);
              }
            if (!(build.*build_routine[minam1]
                                      [maxam12]
                                      [minam3]
                                      [maxam34][eAB])()) {
              ExEnv::outn() << "build2e.cc: did not succeed in building all integrals"
                   << endl;
              abort();
              }
            }
          else {
            blockbuildprim(minam1,maxam12,minam3,maxam34);
            }

          /* Contract the primitive target integrals. */
          /* Throw out all unneeded contractions. */
          if (i||j||k||l) {
            for (m=mlower; m<=mupper; m++) {
              int o;
              int sizec = contract_length(m,nlower,nupper);
              double *restrictxx con_ints = e0f0_ijk[cl](m,nlower);
              bufferprim = build.int_v_list(m,nlower,0);

              for (o=sizec; o!=0; o--) {
                *con_ints++ += *bufferprim++;
                }

              }
            }
          else {
            // for the first primitive write to con_ints rather
            // than accumulate into it
            for (m=mlower; m<=mupper; m++) {
              int o;
              int sizec = contract_length(m,nlower,nupper);
              double *restrictxx con_ints = e0f0_ijk[cl](m,nlower);
              bufferprim = build.int_v_list(m,nlower,0);

              for (o=sizec; o!=0; o--) {
                *con_ints++ = *bufferprim++;
                }

              }
            }

          }
        }
      }
    }

          }
        }
      }
    }

  }

void
Int2eV3::build_using_gcs(int nc1, int nc2, int nc3, int nc4,
                         int minam1, int minam3, int maxam12, int maxam34,
                         int dam1, int dam2, int dam3, int dam4, int eAB)
{
  int i,j,k,l,m;
  int ci,cj,ck,cl;
  int maxam1234=maxam12+maxam34;
  double coef0,coef1,coef2,coef3;
  double ishl1expi=1.0, ishl2expj=1.0, ishl3expk=1.0;
  double *bufferprim;
  double c0scale;

  /* Loop over the primitives. */
  for (i=0; i<int_shell1->nprimitive(); i++) {
    if (int_store1) opr1 = int_shell_to_prim[osh1] + i;
    if (int_expweight1) ishl1expi=2.0*int_shell1->exponent(i);

    for (j=0; j<int_shell2->nprimitive(); j++) {
      if (int_store1) opr2 = int_shell_to_prim[osh2] + j;
      ishl2expj = (int_expweight2) ? 
                        2.0*int_shell2->exponent(j)*ishl1expi : ishl1expi;

      for (k=0; k<int_shell3->nprimitive(); k++) {
        if (int_store1) opr3 = int_shell_to_prim[osh3] + k;
        ishl3expk = (int_expweight3) ? 
                        2.0*int_shell3->exponent(k)*ishl2expj : ishl2expj;

        for (l=0; l<int_shell4->nprimitive(); l++) {
          if (int_store1) opr4 = int_shell_to_prim[osh4] + l;
          c0scale = (int_expweight4) ? 
                        2.0*int_shell4->exponent(l)*ishl3expk : ishl3expk;

          /* Produce the remaining intermediates. */
          gen_prim_intermediates(i,j,k,l, maxam1234);

          /* Generate the target integrals. */
          if (!maxam1234) {
            /* Do nothing: gen_prim_intermediates has set everything up. */
            }
          else if ((minam1<=MG)&&(minam3<=MG)&&(maxam12<=MG)&&(maxam34<=MG)) {
            intfunc brptr=build_routine[minam1][maxam12][minam3][maxam34][eAB];
            if (brptr == &BuildIntV3::impossible_integral) {
              ExEnv::errn() << scprintf("trying to build with int2v%d%d%d%d (exact)\n",
                      minam1,maxam12,minam3,maxam34);
              }
            if (!(build.*brptr)()) {
              ExEnv::outn() << "build2e.cc: did not succeed in building all integrals"
                   << endl;
              abort();
              }
            }
          else {
            blockbuildprim(minam1,maxam12,minam3,maxam34);
            }

          /* Sum thru all possible contractions.
           * Throw out all unneeded contractions. */

  for (ci=0; ci<nc1; ci++) {
    int mlower = int_shell1->am(ci) + dam1;
    if (mlower < 0) continue;
    coef0 = int_shell1->coefficient_unnorm(ci,i)*c0scale;
    IntV3Arraydoublep2 ***e0f0_i = e0f0_con_ints_array[ci];
    for (cj=0; cj<nc2; cj++) {
      int mupper = mlower + int_shell2->am(cj) + dam2;
      if (mupper < mlower) continue;
      if (mlower < minam1) mlower = minam1;
      if (mupper > maxam12) mupper = maxam12;
      coef1 = int_shell2->coefficient_unnorm(cj,j)*coef0;
      IntV3Arraydoublep2 **e0f0_ij = e0f0_i[cj];
      for (ck=0; ck<nc3; ck++) {
        int nlower = int_shell3->am(ck) + dam3;
        if (nlower < 0) continue;
        coef2 = int_shell3->coefficient_unnorm(ck,k)*coef1;
        IntV3Arraydoublep2 *e0f0_ijk = e0f0_ij[ck];
        for (cl=0; cl<nc4; cl++) {
          int nupper = nlower + int_shell4->am(cl) + dam4;
          if (nupper < nlower) continue;
          if (nlower < minam3) nlower = minam3;
          if (nupper > maxam34) nupper = maxam34;
          coef3 = int_shell4->coefficient_unnorm(cl,l)*coef2;

          /* Contract the primitive target integrals. */
          if (i||j||k||l) {
            for (m=mlower; m<=mupper; m++) {
              int o;
              int sizec = contract_length(m,nlower,nupper);
              double *restrictxx con_ints = e0f0_ijk[cl](m,nlower);
              bufferprim = build.int_v_list(m,nlower,0);
              /* Sum the integrals into the contracted integrals. */
#ifdef SUNMOS
              for (o=0; o < sizec; o++) {
                con_ints[o] += coef3 * bufferprim[o];
                }