build2e.cc

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

  size34m1 = INT_NCART_DEC(am34,size34);

  // Some local intermediates
  double half_ooze = 0.5 * build.int_v_ooze;
  double zeta34_ooze = build.int_v_zeta34 * build.int_v_ooze;
  double W0_m_p120 = build.int_v_W0 - build.int_v_p120;
  double p120_m_r10 = build.int_v_p120 - build.int_v_r10;
  double oo2zeta12 = build.int_v_oo2zeta12;
  double p121_m_r11 = build.int_v_p121 - build.int_v_r11;
  double W1_m_p121 = build.int_v_W1 - build.int_v_p121;
  double p122_m_r12 = build.int_v_p122 - build.int_v_r12;
  double W2_m_p122 = build.int_v_W2 - build.int_v_p122;

  stack_alignment_check(&half_ooze, "buildprim_1: half_ooze");

  for (int am12=amin; am12<=amax; am12++) {
    /* Construct the needed intermediate integrals. */
    double ***vlist0 = build.int_v_list(am12);
    double **vlist00 = vlist0[am34];
    I00 = vlist00[m];
    I10 = vlist10[m];
    I11 = vlist10[m+1];
    //I00 = build.int_v_list(am12,am34,m);
    //I10 = build.int_v_list(am12-1,am34,m);
    //I11 = build.int_v_list(am12-1,am34,m+1);
    if (am34) {
      I31 = vlist11[m+1];
      //I31 = build.int_v_list(am12 - 1, am34 - 1, m + 1);
      vlist11 = vlist0[am34-1];
      }
    if (am12>1) {
      I20 = vlist20[m];
      I21 = vlist20[m+1];
      //I20 = build.int_v_list(am12 - 2, am34, m);
      //I21 = build.int_v_list(am12 - 2, am34, m + 1);
      }
    vlist20 = vlist10;
    vlist10 = vlist00;

    /* Construct the new integrals. */
    cartindex12 = 0;
    cartindex1234 = 0;
    // the i12==0, k12==0, j12=am12 case (build on y)
    i12 = 0;
    j12 = am12;
    k12 = 0;

    int i12y1 = 0;  //= INT_CARTINDEX(am12-1,i12,j12-1);
    int i12y1s34 = i12y1*size34;
    int i12y1s34m1 = i12y1*size34m1;
    double *I10i = &I10[i12y1s34];
    double *I11i = &I11[i12y1s34];
    double *restrictxx I00i = &I00[cartindex1234];
    if (j12==1) {
      for (cartindex34=0; cartindex34<size34; cartindex34++) {
        I00i[cartindex34]
          = I10i[cartindex34] * p121_m_r11
          + I11i[cartindex34] * W1_m_p121;
        }
      }
    else { // j12 > 1
      int i12y2s34 = 0; // = INT_CARTINDEX(am12-2,i12,j12-2)*size34;
      double *I20i = &I20[i12y2s34];
      double *I21i = &I21[i12y2s34];
      for (cartindex34=0; cartindex34<size34; cartindex34++) {
        I00i[cartindex34]
          = I10i[cartindex34] * p121_m_r11
          + I11i[cartindex34] * W1_m_p121
          + (j12 - 1) * oo2zeta12 * (I20i[cartindex34]
                                     - I21i[cartindex34] * zeta34_ooze);
        }
      }
    if (am34) {
      double *I31i = &I31[i12y1s34m1];
      cartindex34 = 0;
      for (i34=0; i34<=am34; i34++) {
        //note: k34 < am34-i34 instead of <= am34-i34, so j34 > 0
        int i34y1 = cartindex34-i34;//=INT_CARTINDEX(am34-1,i34,j34-1)
        j34 = am34 - i34;
        double j34_half_ooze = j34 * half_ooze;
        for (k34=0; k34<am34-i34; k34++) {

          I00i[cartindex34] +=  j34_half_ooze * I31i[i34y1];

          j34_half_ooze -= half_ooze;
          i34y1++;
          /* cartindex34 == INT_CARTINDEX(am34,i34,j34) */
          cartindex34++;
          }
        // increment cartindex34 here since the last k was skipped
        cartindex34++;
        }
      }
    cartindex12++;
    cartindex1234+=size34;

    // the i12==0, j12==am12-1, k12==1 case (build on z)
    i12 = 0;
    j12 = am12 - 1;
    k12 = 1;
    int i12z1 = 0;//= INT_CARTINDEX(am12-1,i12,j12);
    int i12z1s34 = i12z1*size34;
    int i12z1s34m1 = i12z1*size34m1;
    I10i = &I10[i12z1s34];
    I11i = &I11[i12z1s34];
    I00i = &I00[cartindex1234];
    for (cartindex34=0; cartindex34<size34; cartindex34++) {
      I00i[cartindex34]
        = I10i[cartindex34] * p122_m_r12
        + I11i[cartindex34] * W2_m_p122;
      }
    if (am34) {
      double *I31i = &I31[i12z1s34m1];
      cartindex34 = 0;
      for (i34=0; i34<=am34; i34++) {
        // skip k34 == 0
        cartindex34++;
        int i34z1 = cartindex34-i34-1;//=INT_CARTINDEX(am34-1,i34,j34)
        double k34_half_ooze = half_ooze;
        for (k34=1; k34<=am34-i34; k34++) {
          I00i[cartindex34] += k34_half_ooze * I31i[i34z1];
          k34_half_ooze += half_ooze;
          i34z1++;
          cartindex34++;
          }
        }
      }
    cartindex12++;
    cartindex1234+=size34;
    // the i12==0, j12==am12-k12, k12>1 case (build on z)
    double k12m1_oo2zeta12 = oo2zeta12;
    for (k12=2; k12<=am12-i12; k12++) {
      j12 = am12 - k12;
      i12z1 = cartindex12-i12-1;//=INT_CARTINDEX(am12-1,i12,j12);
      i12z1s34 = i12z1*size34;
      i12z1s34m1 = i12z1*size34m1;
      int i12z2s34 = (cartindex12-i12-i12-2)*size34;
      //=INT_CARTINDEX(am12-2,i12,j12)*size34;
      I10i = &I10[i12z1s34];
      I11i = &I11[i12z1s34];
      double *I20i = &I20[i12z2s34];
      double *I21i = &I21[i12z2s34];
      I00i = &I00[cartindex1234];
      for (cartindex34=0; cartindex34<size34; cartindex34++) {
        I00i[cartindex34]
          = I10i[cartindex34] * p122_m_r12
          + I11i[cartindex34] * W2_m_p122
          + k12m1_oo2zeta12 * (I20i[cartindex34]
                               - I21i[cartindex34] * zeta34_ooze);
        }
      if (am34) {
        double *I31i = &I31[i12z1s34m1];
        cartindex34 = 0;
        for (i34=0; i34<=am34; i34++) {
          // skip k34 == 0
          cartindex34++;
          int i34z1 = cartindex34-i34-1;//=INT_CARTINDEX(am34-1,i34,j34)
          double k34_half_ooze = half_ooze;
          for (k34=1; k34<=am34-i34; k34++) {
            I00i[cartindex34]
              +=  k34_half_ooze * I31i[i34z1];
            k34_half_ooze += half_ooze;
            i34z1++;
            cartindex34++;
            }
          }
        }
      cartindex12++;
      cartindex1234+=size34;
      k12m1_oo2zeta12 += oo2zeta12;
      }

    // the i12==1 case (build on x)
    i12 = 1;
    int i12x1 = cartindex12-am12-1;//=INT_CARTINDEX(am12-1,i12-1,am12-i12)
    int i12x1s34 = i12x1*size34;
    int i12x1s34m1 = i12x1*size34m1;
    I00i = &I00[cartindex1234];
    I10i = &I10[i12x1s34];
    I11i = &I11[i12x1s34];
    //for (k12=0; k12<=am12-i12; k12++)
    int k12_cartindex34;
    int nk12_size34 = am12*size34;
    for (k12_cartindex34=0; k12_cartindex34<nk12_size34; k12_cartindex34++) {
      *I00i++ = *I10i++ * p120_m_r10 + *I11i++ * W0_m_p120;
      }
    I00i = &I00[cartindex1234];
    if (am34) {
      double *I31i = &I31[i12x1s34m1];
      for (k12=0; k12<am12; k12++) {
        // skip over i34==0
        double *I00is=&I00i[am34+1];
        double i34_half_ooze = half_ooze;
        for (i34=1; i34<=am34; i34++) {
          for (k34=i34; k34<=am34; k34++) { // index_k34 = true_k34 + i34
            *I00is++ +=  i34_half_ooze * *I31i++;
            }
          i34_half_ooze += half_ooze;
          }
        I00i += size34;
        }
      }
    cartindex12 += am12;
    cartindex1234 += am12*size34;
    // the i12>1 case (build on x)
    if (am12<2) continue;
    double i12m1_oo2zeta12 = oo2zeta12;
    i12x1 = cartindex12-am12-1;
    i12x1s34 = i12x1*size34;
    i12x1s34m1 = i12x1*size34m1;
    int i12x2s34 = (cartindex12-am12-am12-1)*size34;
    I10i = &I10[i12x1s34];
    I11i = &I11[i12x1s34];
    double *I20i = &I20[i12x2s34];
    double *I21i = &I21[i12x2s34];
    I00i = &I00[cartindex1234];
    for (i12=2; i12<=am12; i12++) {
      int sizek12_size34 = (am12-i12+1)*size34;
      int k12_c34;
      for (k12_c34=0; k12_c34<sizek12_size34; k12_c34++) {
        *I00i++
          = *I10i++ * p120_m_r10
          + *I11i++ * W0_m_p120
          + i12m1_oo2zeta12 * (*I20i++
                               - *I21i++
                               * zeta34_ooze);
        }
      i12m1_oo2zeta12 += oo2zeta12;
      }
    if (am34) {
      double *I31i = &I31[i12x1s34m1];
      I00i = &I00[cartindex1234];
      for (i12=2; i12<=am12; i12++) {
        for (k12=0; k12<=am12-i12; k12++) {
          // skip over i34==0
          double *I00is=&I00i[am34+1];
          double i34_half_ooze = half_ooze;
          for (i34=1; i34<=am34; i34++) {
            for (k34=0; k34<=am34-i34; k34++) {
              *I00is++ += i34_half_ooze * *I31i++;
              }
            i34_half_ooze += half_ooze;
            }
          I00i += size34;
          }
        }
      }
    }
}

void
Int2eV3::blockbuildprim_3(int bmin,int bmax,int m)
{
  double *I00;
  double *I10; /* = [a0|c0](m) */
  double *I11; /* = [a0|c0](m+1) */
  double *I20; /* = [a0|c-1 0](m) */
  double *I21; /* = [a0|c-1 0](m+1) */
  int ci34m1,ci34m2;
  int size34,size34m1,size34m2;
  int i34, k34;

  // These temporaries point to subblocks within the integrals arrays.
  double *I10o,*I11o,*I20o,*I21o;

  double ***vlist0;
  double **vlist01;
  double **vlist02;

  vlist0 = build.int_v_list(0);
  vlist01 = vlist0[bmin-1];
  if (bmin>1) {
    vlist02 = vlist0[bmin-2];
    }

  for (int am34=bmin; am34<=bmax; am34++) {

    /* Construct the needed intermediate integrals. */
    double **vlist00 = vlist0[am34];
    I00 = vlist00[m];
    I10 = vlist01[m];
    I11 = vlist01[m+1];
    //I00 = build.int_v_list(0, am34, m);
    //I10 = build.int_v_list(0, am34 - 1, m);
    //I11 = build.int_v_list(0, am34 - 1, m + 1);
    if (am34>1) {
      I20 = vlist02[m];
      I21 = vlist02[m+1];
      //I20 = build.int_v_list(0, am34 - 2, m);
      //I21 = build.int_v_list(0, am34 - 2, m + 1);
      }
    vlist02 = vlist01;
    vlist01 = vlist00;

    /* The size of the group of primitives with ang. mom. = am34 - 1 */
    size34 = INT_NCART_NN(am34);
    size34m1 = INT_NCART_DEC(am34,size34);
    size34m2 = INT_NCART(am34-2);

    // Useful constants
    double p340_m_r30 = build.int_v_p340 - build.int_v_r30;
    double W0_m_p340 = build.int_v_W0 - build.int_v_p340;
    double p341_m_r31 = build.int_v_p341 - build.int_v_r31;
    double W1_m_p341 = build.int_v_W1 - build.int_v_p341;
    double p342_m_r32 = build.int_v_p342 - build.int_v_r32;
    double W2_m_p342 = build.int_v_W2 - build.int_v_p342;
    double oo2zeta34 = build.int_v_oo2zeta34;
    double zeta12_ooze = build.int_v_zeta12 * build.int_v_ooze;

    stack_alignment_check(&p340_m_r30, "buildprim_3: p340_m_r30");

    /* Construct the new integrals. */
    double *restrictxx I00o = I00; // points the current target integral
    I10o = I10;
    I11o = I11;
    //int cartindex34 = 0;
    // i34 == 0, k34 == 0, j34 = am34
    /* ------------------ Build from the y position. */
    /* I10 I11 and I21 */
    *I00o = *I10o * p341_m_r31 + *I11o * W1_m_p341;
    if (am34>1) {
      I20o = I20;
      I21o = I21;
      *I00o += (am34 - 1) * oo2zeta34 * (*I20o
                                         - *I21o * zeta12_ooze);
      }
    //cartindex34++;
    // i34 == 0, k34 >= 1
    // loop over a portion of the l=am34-1 integrals
    I00o = &I00o[1];
    for (ci34m1=0; ci34m1<am34; ci34m1++) {
      /* ------------------ Build from the z position. */
      //note: ci34m1 = cartindex34 - i34 - 1;//=INT_CARTINDEX(am34-1,i34,j34)
      /* I10 and I11 */
      I00o[ci34m1] = I10o[ci34m1] * p342_m_r32 + I11o[ci34m1] * W2_m_p342;
      }
    // skip over i34 == 0, k34 == 1
    //cartindex34++;
    // i34 == 0, k34 > 1
    I00o = &I00o[1];
    // loop over a portion of the l=am34-2 integrals
    double k34m1_oo2zeta34 = oo2zeta34;
    for (ci34m2=0; ci34m2<am34-1; ci34m2++) {
      //note: k34 = 2+ci34m2
      /* ------------------ Build from the z position. */
      /* I20 and I21 */
      I00o[ci34m2]
        +=  k34m1_oo2zeta34 * (I20o[ci34m2] - I21o[ci34m2] * zeta12_ooze);
      k34m1_oo2zeta34 += oo2zeta34;
      }
    //cartindex34+=am34-1;
    // i34 >= 1
    I00o = &I00o[am34-1];
    //note: ci34m1 = INT_CARTINDEX(am34-1,i34-1,j34)
    for (ci34m1=0; ci34m1<size34m1; ci34m1++) {
      /* I10 and I11 contrib */
      /* ------------------ Build from the x position. */
      I00o[ci34m1] = I10o[ci34m1] * p340_m_r30 + I11o[ci34m1] * W0_m_p340;
      }
    // skip past i34 == 1
    //cartindex34 += am34;
    // i34 > 1
    I00o = &I00o[am34];
    //note: ci34m2=INT_CARTINDEX(am34-2,i34-2,j34)
    ci34m2=0;
    double i34m1_oo2zeta34 = oo2zeta34;
    for (i34=2; i34<=am34; i34++) {
      for (k34=0; k34<=am34-i34; k34++) {
        /* I20 and I21 contrib */
        /* ------------------ Build from the x position. */
        I00o[ci34m2]
          +=  i34m1_oo2zeta34 * (I20o[ci34m2] - I21o[ci34m2] * zeta12_ooze);
        ci34m2++;
        }
      i34m1_oo2zeta34 += oo2zeta34;
      }
    //cartindex34 += size34m2;

    I00o = &I00o[size34m2];
    }
}

/////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "CLJ-CONDENSED"
// End: