shift2e.cc

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

      }
      // swap the buffers.
      double *tmp = buf12_1;
      buf12_1 = buf12_2;
      buf12_2 = tmp;
    }

  /* Construct the target integrals. */
  return buf12_2;
  }

/* Shift angular momentum from center 1 to center 2.
 * I0100 are the target integrals.
 * am1-4 is the angular momentum on each of the centers in the target set.
 */
void
Int2eV3::shiftam_12(double *I0100, double *I1000, double *I0000,
                    int am1, int am2, int am3, int am4)
{
  int i;
  int i1,k1;
  int size2, size2m134, size34;

#if CHECK_INTEGRAL_ALGORITHM > 1
  ExEnv::outn() << "(" << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
       << " <- "
       << "(" << am1+1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
       << "(" << am1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
       << endl;
#endif

  size2m134 = INT_NCART(am2-1)*INT_NCART(am3)*INT_NCART(am4);
  size34    = INT_NCART(am3)*INT_NCART(am4);
  size2     = INT_NCART(am2);

  int size_zcontrib = am2*size34;
  int size_xcontrib = (size2-(am2+1))*size34;

  double AmB0 = AmB[0];
  double AmB1 = AmB[1];
  double AmB2 = AmB[2];

  /* Loop over the target integrals. */
  double *restrictxx I0100i=I0100;
  int cartindex1 = 0;
  for (i1=0; i1<=am1; i1++) {
    for (k1=0; k1<=am1-i1; k1++) {
      //int j1 = am1 - i1 - k1;
      int ci1x1 = (cartindex1 + am1 + 2) * size2m134;
      int ci1y1 = (cartindex1 + i1) * size2m134;
      int ci1z1 = (cartindex1 + i1 + 1)  * size2m134;
      //note:
      //ci1x1 = INT_CARTINDEX(am1+1,i1+1,j1) * size2m134;
      //ci1y1 = INT_CARTINDEX(am1+1,i1,j1+1) * size2m134;
      //ci1z1 = INT_CARTINDEX(am1+1,i1,j1) * size2m134;
      int ci1 = cartindex1 * size2m134;
      // i2 == 0, k2 == 0, j2 == am2 (>0)
      double *I1000i=&I1000[ci1y1];
      double *I0000i=&I0000[ci1];
      for (i=0; i<size34; i++) {
        I0100i[i] = I1000i[i] + I0000i[i] * AmB1;
        }
      I0100i=&I0100i[size34];
      // i2 == 0, k2 > 0
      I1000i=&I1000[ci1z1];
      I0000i=&I0000[ci1];
      for (i=0; i<size_zcontrib; i++) {
        I0100i[i] = I1000i[i] + I0000i[i] * AmB2;
        }
      I0100i=&I0100i[size_zcontrib];
      // i2 >= 1
      I1000i=&I1000[ci1x1];
      I0000i=&I0000[ci1];
      for (i=0; i<size_xcontrib; i++) {
        I0100i[i] = I1000i[i] + I0000i[i] * AmB0;
        }
      I0100i=&I0100i[size_xcontrib];

      cartindex1++;
      }
    }
  }


/* Shift angular momentum from center 1 to center 2 when centers
 * one and two are the same.
 * I0100 are the target integrals.
 * am1-4 is the angular momentum on each of the centers in the target set.
 */
void
Int2eV3::shiftam_12eAB(double *I0100, double *I1000, double *I0000,
                       int am1, int am2, int am3, int am4)
{
  int i;
  int i1,k1;
  int size2, size2m134, size34;

#if CHECK_INTEGRAL_ALGORITHM > 1
  ExEnv::outn() << "(" << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
       << " <- "
       << "(" << am1+1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
       << "(" << am1 << "," << am2-1 << "," << am3 << "," << am4 << ")"
       << endl;
#endif

  size2m134 = INT_NCART(am2-1)*INT_NCART(am3)*INT_NCART(am4);
  size34    = INT_NCART(am3)*INT_NCART(am4);
  size2     = INT_NCART(am2);

  int size_zcontrib = am2*size34;
  int size_xcontrib = (size2-(am2+1))*size34;

  /* Loop over the target integrals. */
  double *restrictxx I0100i=I0100;
  int cartindex1 = 0;
  for (i1=0; i1<=am1; i1++) {
    for (k1=0; k1<=am1-i1; k1++) {
      //int j1 = am1 - i1 - k1;
      int ci1x1 = (cartindex1 + am1 + 2) * size2m134;
      int ci1y1 = (cartindex1 + i1) * size2m134;
      int ci1z1 = (cartindex1 + i1 + 1)  * size2m134;
      //note:
      //ci1x1 = INT_CARTINDEX(am1+1,i1+1,j1) * size2m134;
      //ci1y1 = INT_CARTINDEX(am1+1,i1,j1+1) * size2m134;
      //ci1z1 = INT_CARTINDEX(am1+1,i1,j1) * size2m134;
      // i2 == 0, k2 == 0, j2 == am2 (>0)
      double *I1000i=&I1000[ci1y1];
      for (i=0; i<size34; i++) {
        I0100i[i] = I1000i[i];
        }
      I0100i=&I0100i[size34];
      // i2 == 0, k2 > 0
      I1000i=&I1000[ci1z1];
      for (i=0; i<size_zcontrib; i++) {
        I0100i[i] = I1000i[i];
        }
      I0100i=&I0100i[size_zcontrib];
      // i2 >= 1
      I1000i=&I1000[ci1x1];
      for (i=0; i<size_xcontrib; i++) {
        I0100i[i] = I1000i[i];
        }
      I0100i=&I0100i[size_xcontrib];

      cartindex1++;
      }
    }
  }

void
Int2eV3::shiftam_34(double *restrictxx I0001, double *I0010, double *I0000,
                    int am1, int am2, int am3, int am4)
{
  int i1,k1,cartindex1;
  int i2,k2,cartindex2;
  int i3,k3,cartindex3;
  int i4,k4,cartindex4;
  int cartindex1234;
  int size23p14m1,size3p14m1,size4m1,size234m1,size34m1;

#if CHECK_INTEGRAL_ALGORITHM > 1
  ExEnv::outn() << "(" << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
       << " <- "
       << "(" << am1 << "," << am2 << "," << am3+1 << "," << am4-1 << ")"
       << "(" << am1 << "," << am2 << "," << am3 << "," << am4-1 << ")"
       << endl;
#endif

  size23p14m1 = INT_NCART(am2)*INT_NCART(am3+1)*INT_NCART(am4-1);
  size3p14m1 = INT_NCART(am3+1)*INT_NCART(am4-1);
  size4m1 = INT_NCART(am4-1);

  size234m1 = INT_NCART(am2)*INT_NCART(am3)*INT_NCART(am4-1);
  size34m1 = INT_NCART(am3)*INT_NCART(am4-1);

  double CmD0 = CmD[0];
  double CmD1 = CmD[1];
  double CmD2 = CmD[2];

  /* Loop over the target integrals. */
  cartindex1234 = 0;
  cartindex1 = 0;
  for (i1=0; i1<=am1; i1++) {
    for (k1=0; k1<=am1-i1; k1++) {
      //int j1 = am1 - i1 - k1;
      int ci1_I0010 = cartindex1 * size23p14m1;
      int ci1_I0000 = cartindex1 * size234m1;
      cartindex2 = 0;
      for (i2=0; i2<=am2; i2++) {
        for (k2=0; k2<=am2-i2; k2++) {
          //int j2 = am2 - i2 - k2;
          int ci2_I0010 = ci1_I0010 + cartindex2 * size3p14m1;
          int ci2_I0000 = ci1_I0000 + cartindex2 * size34m1;
          cartindex3 = 0;
          for (i3=0; i3<=am3; i3++) {
            for (k3=0; k3<=am3-i3; k3++) {
              //int j3 = am3 - i3 - k3;
              //note: cartindex3 + am3 + 2 = INT_CARTINDEX(am3+1,i3+1,j3)
              int ci3_I0010 = ci2_I0010
                              + (cartindex3 + am3 + 2)*size4m1;
              int ci3_I0000 = ci2_I0000 + cartindex3*size4m1;
              //cartindex4 = 0;
              // this routine called only when am4 > 0
              ///// CASE 1: i4 = 0 k4 = 0 j4 = am4; shift on y
              //note: j4 = am4;
              //note: cartindex4 - i4 = INT_CARTINDEX(am4-1,i4,j4-1)
              //note: cartindex3 - i3 = INT_CARTINDEX(am3+1,i3,j3+1)
              int ci3 = cartindex3 + i3;
              I0001[cartindex1234]
                =   I0010[ci2_I0010 + ci3 * size4m1]
                + I0000[ci3_I0000]
                * CmD1;
              cartindex1234++;
              //cartindex4++;
              ///// CASE 2: i4 = 0 k4 > 0; shift on z
              ci3++;
              for (int ci4=0; ci4<am4; ci4++) {
                //note: j4 = am4 - i4 - k4;
                //note: cartindex4 - i4 - 1 = INT_CARTINDEX(am4-1,i4,j4)
                //note: ci4 = cartindex4 - i4 - 1;
                //note: cartindex3 - i3 - 1 = INT_CARTINDEX(am3+1,i3,j3)
                I0001[cartindex1234]
                  =   I0010[ci2_I0010 + ci3 * size4m1 + ci4 ]
                  + I0000[ci3_I0000 + ci4 ]
                  * CmD2;
                cartindex1234++;
                //cartindex4++;
                }
              ///// CASE 3: i4 > 0; shift on x
              int ncart_remain = INT_NCART(am4) - (am4+1);
              for (int ci4=0; ci4<ncart_remain; ci4++) {
                  //note: j4 = am4 - i4 - k4;
                  //note: cartindex4 - am4 - 1 = INT_CARTINDEX(am4-1,i4-1,j4)
                  //note: ci4 = cartindex4 - am4 - 1;
                  I0001[cartindex1234]
                    =   I0010[ci3_I0010 + ci4]
                    + I0000[ci3_I0000 + ci4]
                    * CmD0;
                  cartindex1234++;
                  //cartindex4++;
                  }
              cartindex3++;
              }
            }
          cartindex2++;
          }
        }
      cartindex1++;
      }
    }
  }


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

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