build2e.cc

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

#else
              for (o=sizec; o; o--) {
                *con_ints++ += coef3 * *bufferprim++;
                }
#endif
              }

            }
          else {
            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);
              /* Write the integrals to the contracted integrals. */
#ifdef SUNMOS
              for (o=0; o < sizec; o++) {
                con_ints[o] = coef3 * bufferprim[o];
                }
#else
              for (o=sizec; o; o--) {
                *con_ints++ = coef3 * *bufferprim++;
                }
#endif
              }
            }
          }
        }
      }
    }


          }
        }
      }
    }
  }

/* This routine constructs intermediates needed for each quartet of
 * primitives.  It is given the total angular momentum as the argument
 * and requires that the global primitive offsets and other global
 * constants be initialized. */
void
Int2eV3::gen_prim_intermediates(int pr1, int pr2, int pr3, int pr4, int am)
{
  int i;
  double T;
  double pmq,pmq2;
  double AmB,AmB2;
  /* This is 2^(1/2) * pi^(5/4) */
  const double sqrt2pi54 = 5.9149671727956129;
  double conv_to_s;

  if (int_store2 && !int_unit2 && !int_unit4) {
    double *tmp;
    build.int_v_zeta12 = int_prim_zeta(opr1,opr2);
    build.int_v_zeta34 = int_prim_zeta(opr3,opr4);
    build.int_v_oo2zeta12 = int_prim_oo2zeta(opr1,opr2);
    build.int_v_oo2zeta34 = int_prim_oo2zeta(opr3,opr4);
    tmp = int_prim_p(opr1,opr2);
    build.int_v_p120 = *tmp++;
    build.int_v_p121 = *tmp++;
    build.int_v_p122 = *tmp;
    tmp = int_prim_p(opr3,opr4);
    build.int_v_p340 = *tmp++;
    build.int_v_p341 = *tmp++;
    build.int_v_p342 = *tmp;
    build.int_v_k12 = int_prim_k(opr1,opr2);
    build.int_v_k34 = int_prim_k(opr3,opr4);
    }
  else {
    build.int_v_zeta12 = int_shell1->exponent(pr1) + int_shell2->exponent(pr2);
    build.int_v_zeta34 = int_shell3->exponent(pr3) + int_shell4->exponent(pr4);
    build.int_v_oo2zeta12 = 1.0/build.int_v_zeta12;
    build.int_v_oo2zeta34 = 1.0/build.int_v_zeta34;
    build.int_v_p120 = build.int_v_oo2zeta12
                     * ( int_shell1->exponent(pr1) * build.int_v_r10
                         + int_shell2->exponent(pr2) * build.int_v_r20 );
    build.int_v_p121 = build.int_v_oo2zeta12
                     * ( int_shell1->exponent(pr1) * build.int_v_r11
                         + int_shell2->exponent(pr2) * build.int_v_r21 );
    build.int_v_p122 = build.int_v_oo2zeta12
                     * ( int_shell1->exponent(pr1) * build.int_v_r12
                         + int_shell2->exponent(pr2) * build.int_v_r22 );
    build.int_v_p340 = build.int_v_oo2zeta34
                     * ( int_shell3->exponent(pr3) * build.int_v_r30
                         + int_shell4->exponent(pr4) * build.int_v_r40 );
    build.int_v_p341 = build.int_v_oo2zeta34
                     * ( int_shell3->exponent(pr3) * build.int_v_r31
                         + int_shell4->exponent(pr4) * build.int_v_r41 );
    build.int_v_p342 = build.int_v_oo2zeta34
                     * ( int_shell3->exponent(pr3) * build.int_v_r32
                         + int_shell4->exponent(pr4) * build.int_v_r42 );

    /* Compute AmB^2 for shell 1 and 2. */
    AmB = build.int_v_r20 - build.int_v_r10;
    AmB2 = AmB*AmB;
    AmB = build.int_v_r21 - build.int_v_r11;
    AmB2 += AmB*AmB;
    AmB = build.int_v_r22 - build.int_v_r12;
    AmB2 += AmB*AmB;

    build.int_v_k12 =    sqrt2pi54
                 * build.int_v_oo2zeta12
                 * exp( -   int_shell1->exponent(pr1)*int_shell2->exponent(pr2)
                          * build.int_v_oo2zeta12
                          * AmB2 );

    /* Compute AmB^2 for shells 3 and 4. */
    AmB = build.int_v_r40 - build.int_v_r30;
    AmB2 = AmB*AmB;
    AmB = build.int_v_r41 - build.int_v_r31;
    AmB2 += AmB*AmB;
    AmB = build.int_v_r42 - build.int_v_r32;
    AmB2 += AmB*AmB;

    build.int_v_k34 =    sqrt2pi54
                 * build.int_v_oo2zeta34
                 * exp( -   int_shell3->exponent(pr3)*int_shell4->exponent(pr4)
                          * build.int_v_oo2zeta34
                          * AmB2 );

    build.int_v_oo2zeta12 *= 0.5;
    build.int_v_oo2zeta34 *= 0.5;
    }

  build.int_v_ooze = 1.0/(build.int_v_zeta12 + build.int_v_zeta34);

  build.int_v_W0 = build.int_v_ooze*(  build.int_v_zeta12 * build.int_v_p120
                         + build.int_v_zeta34 * build.int_v_p340 );
  build.int_v_W1 = build.int_v_ooze*(  build.int_v_zeta12 * build.int_v_p121
                         + build.int_v_zeta34 * build.int_v_p341 );
  build.int_v_W2 = build.int_v_ooze*(  build.int_v_zeta12 * build.int_v_p122
                         + build.int_v_zeta34 * build.int_v_p342 );

  pmq = build.int_v_p120 - build.int_v_p340;
  pmq2 = pmq*pmq;
  pmq = build.int_v_p121 - build.int_v_p341;
  pmq2 += pmq*pmq;
  pmq = build.int_v_p122 - build.int_v_p342;
  pmq2 += pmq*pmq;

  T =   build.int_v_zeta12
      * build.int_v_zeta34
      * build.int_v_ooze * pmq2;

  double *fjttable = fjt_->values(am,T);

  /* Convert the fjttable produced by int_fjt into the S integrals */
  conv_to_s = sqrt(build.int_v_ooze) * build.int_v_k12 * build.int_v_k34;
  for (i=0; i<=am; i++) {
    build.int_v_list(0,0,i)[0] =   fjttable[i] * conv_to_s;
    }

  }

/* This is like gen_prim_intermediates, except the normalization is
 * put into the ssss integrals. */
void
Int2eV3::gen_prim_intermediates_with_norm(int pr1, int pr2, int pr3, int pr4,
                                          int am, double norm)
{
  int i;
  double T;
  double pmq,pmq2;
  double AmB,AmB2;
  /* This is 2^(1/2) * pi^(5/4) */
  const double sqrt2pi54 = 5.9149671727956129;
  double conv_to_s;

  if (int_store2) {
    build.int_v_zeta12 = int_prim_zeta(opr1,opr2);
    build.int_v_zeta34 = int_prim_zeta(opr3,opr4);
    build.int_v_oo2zeta12 = int_prim_oo2zeta(opr1,opr2);
    build.int_v_oo2zeta34 = int_prim_oo2zeta(opr3,opr4);
    build.int_v_p120 = int_prim_p(opr1,opr2,0);
    build.int_v_p121 = int_prim_p(opr1,opr2,1);
    build.int_v_p122 = int_prim_p(opr1,opr2,2);
    build.int_v_p340 = int_prim_p(opr3,opr4,0);
    build.int_v_p341 = int_prim_p(opr3,opr4,1);
    build.int_v_p342 = int_prim_p(opr3,opr4,2);
    build.int_v_k12 = int_prim_k(opr1,opr2);
    build.int_v_k34 = int_prim_k(opr3,opr4);
    }
  else {
    build.int_v_zeta12 = int_shell1->exponent(pr1) + int_shell2->exponent(pr2);
    build.int_v_zeta34 = int_shell3->exponent(pr3) + int_shell4->exponent(pr4);
    build.int_v_oo2zeta12 = 1.0/build.int_v_zeta12;
    build.int_v_oo2zeta34 = 1.0/build.int_v_zeta34;
    build.int_v_p120 = build.int_v_oo2zeta12
                     * ( int_shell1->exponent(pr1) * build.int_v_r10
                         + int_shell2->exponent(pr2) * build.int_v_r20 );
    build.int_v_p121 = build.int_v_oo2zeta12
                     * ( int_shell1->exponent(pr1) * build.int_v_r11
                         + int_shell2->exponent(pr2) * build.int_v_r21 );
    build.int_v_p122 = build.int_v_oo2zeta12
                     * ( int_shell1->exponent(pr1) * build.int_v_r12
                         + int_shell2->exponent(pr2) * build.int_v_r22 );
    build.int_v_p340 = build.int_v_oo2zeta34
                     * ( int_shell3->exponent(pr3) * build.int_v_r30
                         + int_shell4->exponent(pr4) * build.int_v_r40 );
    build.int_v_p341 = build.int_v_oo2zeta34
                     * ( int_shell3->exponent(pr3) * build.int_v_r31
                         + int_shell4->exponent(pr4) * build.int_v_r41 );
    build.int_v_p342 = build.int_v_oo2zeta34
                     * ( int_shell3->exponent(pr3) * build.int_v_r32
                         + int_shell4->exponent(pr4) * build.int_v_r42 );

    /* Compute AmB^2 for shell 1 and 2. */
    AmB = build.int_v_r20 - build.int_v_r10;
    AmB2 = AmB*AmB;
    AmB = build.int_v_r21 - build.int_v_r11;
    AmB2 += AmB*AmB;
    AmB = build.int_v_r22 - build.int_v_r12;
    AmB2 += AmB*AmB;

    build.int_v_k12 =    sqrt2pi54
                 * build.int_v_oo2zeta12
                 * exp( -   int_shell1->exponent(pr1)*int_shell2->exponent(pr2)
                          * build.int_v_oo2zeta12
                          * AmB2 );

    /* Compute AmB^2 for shells 3 and 4. */
    AmB = build.int_v_r40 - build.int_v_r30;
    AmB2 = AmB*AmB;
    AmB = build.int_v_r41 - build.int_v_r31;
    AmB2 += AmB*AmB;
    AmB = build.int_v_r42 - build.int_v_r32;
    AmB2 += AmB*AmB;

    build.int_v_k34 =    sqrt2pi54
                 * build.int_v_oo2zeta34
                 * exp( -   int_shell3->exponent(pr3)*int_shell4->exponent(pr4)
                          * build.int_v_oo2zeta34
                          * AmB2 );

    build.int_v_oo2zeta12 *= 0.5;
    build.int_v_oo2zeta34 *= 0.5;
    }

  build.int_v_ooze = 1.0/(build.int_v_zeta12 + build.int_v_zeta34);

  build.int_v_W0 = build.int_v_ooze*(  build.int_v_zeta12 * build.int_v_p120
                         + build.int_v_zeta34 * build.int_v_p340 );
  build.int_v_W1 = build.int_v_ooze*(  build.int_v_zeta12 * build.int_v_p121
                         + build.int_v_zeta34 * build.int_v_p341 );
  build.int_v_W2 = build.int_v_ooze*(  build.int_v_zeta12 * build.int_v_p122
                         + build.int_v_zeta34 * build.int_v_p342 );

  pmq = build.int_v_p120 - build.int_v_p340;
  pmq2 = pmq*pmq;
  pmq = build.int_v_p121 - build.int_v_p341;
  pmq2 += pmq*pmq;
  pmq = build.int_v_p122 - build.int_v_p342;
  pmq2 += pmq*pmq;

  T =   build.int_v_zeta12
      * build.int_v_zeta34
      * build.int_v_ooze * pmq2;

  double *fjttable = fjt_->values(am,T);

  /* Convert the fjttable produced by int_fjt into the S integrals */
  conv_to_s = sqrt(build.int_v_ooze)
            * build.int_v_k12 * build.int_v_k34 * norm;
  for (i=0; i<=am; i++) {
    build.int_v_list(0,0,i)[0] =   fjttable[i] * conv_to_s;
    }

  }


/* This routine computes the shell intermediates. */
void
Int2eV3::gen_shell_intermediates(int sh1, int sh2, int sh3, int sh4)
{
  if (int_store1 && !int_unit2 && !int_unit4) {
    build.int_v_r10 = int_shell_r(osh1,0);
    build.int_v_r11 = int_shell_r(osh1,1);
    build.int_v_r12 = int_shell_r(osh1,2);
    build.int_v_r20 = int_shell_r(osh2,0);
    build.int_v_r21 = int_shell_r(osh2,1);
    build.int_v_r22 = int_shell_r(osh2,2);
    build.int_v_r30 = int_shell_r(osh3,0);
    build.int_v_r31 = int_shell_r(osh3,1);
    build.int_v_r32 = int_shell_r(osh3,2);
    build.int_v_r40 = int_shell_r(osh4,0);
    build.int_v_r41 = int_shell_r(osh4,1);
    build.int_v_r42 = int_shell_r(osh4,2);
    }
  else {
    build.int_v_r10 = bs1_->r(bs1_->shell_to_center(sh1),0);
    build.int_v_r11 = bs1_->r(bs1_->shell_to_center(sh1),1);
    build.int_v_r12 = bs1_->r(bs1_->shell_to_center(sh1),2);
    if (int_unit2) {
        build.int_v_r20 = build.int_v_r10;
        build.int_v_r21 = build.int_v_r11;
        build.int_v_r22 = build.int_v_r12;
      }
    else {
        build.int_v_r20 = bs2_->r(bs2_->shell_to_center(sh2),0);
        build.int_v_r21 = bs2_->r(bs2_->shell_to_center(sh2),1);
        build.int_v_r22 = bs2_->r(bs2_->shell_to_center(sh2),2);
      }
    build.int_v_r30 = bs3_->r(bs3_->shell_to_center(sh3),0);
    build.int_v_r31 = bs3_->r(bs3_->shell_to_center(sh3),1);
    build.int_v_r32 = bs3_->r(bs3_->shell_to_center(sh3),2);
    if (int_unit4) {
        build.int_v_r40 = build.int_v_r30;
        build.int_v_r41 = build.int_v_r31;
        build.int_v_r42 = build.int_v_r32;
      }
    else {
        build.int_v_r40 = bs4_->r(bs4_->shell_to_center(sh4),0);
        build.int_v_r41 = bs4_->r(bs4_->shell_to_center(sh4),1);
        build.int_v_r42 = bs4_->r(bs4_->shell_to_center(sh4),2);
      }
    }
  }

void
Int2eV3::blockbuildprim(int minam1,int maxam12,int minam3,int maxam34)
{
  int m, b;
  int l=maxam12+maxam34;

  // the (0,0,m) integrals have already been initialized

  // compute (0,b,m) integrals
  for (m=l-1; m>=0; m--) {
    int bmax = l-m;
    if (bmax>maxam34) bmax=maxam34;
    blockbuildprim_3(1,bmax,m);
    }

  // compute (a,b,m) integrals
  for (m=maxam12-1; m>=0; m--) {
    for (b=0; b<=maxam34; b++) {
// This is how the code was for a long while,
// but at some point it started giving the wrong
// answers and seems wrong from inspection.  Valgrind
// flags that uninitialized I10i integrals are being
// used, which results from amin > 1.  I have switched
// to the correctly behaving amin = 1.
//       int amin = minam1-m;
//       if (amin<1) amin=1;
//       int amax = maxam12-m;
//       blockbuildprim_1(amin,amax,b,m);
      int amax = maxam12-m;
      blockbuildprim_1(1,amax,b,m);
      }
    }
}

void
Int2eV3::blockbuildprim_1(int amin,int amax,int am34,int m)
{
  double *I00;
  double *I10; /* = [a0|c0](m) */
  double *I11; /* = [a0|c0](m+1) */
  double *I20; /* = [a-1 0|c0](m) */
  double *I21; /* = [a-1 0|c0](m+1) */
  double *I31; /* = [a0|c-1 0](m+1) */
  int cartindex12;
  int cartindex34;
  int cartindex1234;
  int size34=0,size34m1;
  int i12, j12, k12;
  int i34, j34, k34;

  double ***vlist1;
  double **vlist10;
  double **vlist11;
  double ***vlist2;
  double **vlist20;

  vlist1 = build.int_v_list(amin-1);
  vlist10 = vlist1[am34];

  if (am34) {
    vlist11 = vlist1[am34-1];
    }

  if (amin>1) {
    vlist2 = build.int_v_list(amin-2);
    vlist20 = vlist2[am34];
    }

  /* The size of the am34 group of primitives. */
  size34 = INT_NCART_NN(am34);
  /* The size of the group of primitives with ang. mom. = am34 - 1 */