comp1e.cc

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

{
  int i;
  int gc1,gc2;
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int index1,index2;
  double tmp[3];
  double *ctmp = cartesianbuffer;

  c1 = bs1_->shell_to_center(ish);
  c2 = bs2_->shell_to_center(jsh);
  for (int xyz=0; xyz<3; xyz++) {
      A[xyz] = bs1_->r(c1,xyz);
      B[xyz] = bs2_->r(c2,xyz);
    }
  gshell1 = &bs1_->shell(ish);
  gshell2 = &bs2_->shell(jsh);
  FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
      if ((bs1_==bs2_)&&(c1==c2)) {
        if (    three_center
             && !((bs1_==third_centers)&&(c1==third_centernum))
             && ((bs1_==dercs)&&(c1==centernum))) {
          for (i=0; i<3; i++) {
            /* Derivative wrt first shell. */
            exponent_weighted = 0;
            tmp[i] = 2.0 *
               (this->*shell_function)(gc1,i1+IN(i,0),j1+IN(i,1),k1+IN(i,2),gc2,i2,j2,k2);
            exponent_weighted = -1;
            if (SELECT(i1,j1,k1,i)) {
              tmp[i] -= SELECT(i1,j1,k1,i) *
                (this->*shell_function)(gc1,i1-IN(i,0),j1-IN(i,1),k1-IN(i,2),gc2,i2,j2,k2);
              }
            /* Derviative wrt second shell. */
            exponent_weighted = 1;
            tmp[i] += 2.0 *
               (this->*shell_function)(gc1,i1,j1,k1,gc2,i2+IN(i,0),j2+IN(i,1),k2+IN(i,2));
            exponent_weighted = -1;
            if (SELECT(i2,j2,k2,i)) {
              tmp[i] -= SELECT(i2,j2,k2,i) *
                (this->*shell_function)(gc1,i1,j1,k1,gc2,i2-IN(i,0),j2-IN(i,1),k2-IN(i,2));
              }
            }
	  }
        else {
          /* If there are two centers and they are the same, then we
           * use translational invariance to get a net contrib of 0.0 */
          for (i=0; i<3; i++) tmp[i] = 0.0;
          }
        }
      else if ((bs1_==dercs)&&(c1==centernum)) {
        for (i=0; i<3; i++) {
          exponent_weighted = 0;
          tmp[i] = 2.0 *
             (this->*shell_function)(gc1,i1+IN(i,0),j1+IN(i,1),k1+IN(i,2),gc2,i2,j2,k2);
          exponent_weighted = -1;
          if (SELECT(i1,j1,k1,i)) {
            tmp[i] -= SELECT(i1,j1,k1,i) *
              (this->*shell_function)(gc1,i1-IN(i,0),j1-IN(i,1),k1-IN(i,2),gc2,i2,j2,k2);
            }
          }
        }
      else if ((bs2_==dercs)&&(c2==centernum)) {
        for (i=0; i<3; i++) {
          exponent_weighted = 1;
          tmp[i] = 2.0 *
             (this->*shell_function)(gc1,i1,j1,k1,gc2,i2+IN(i,0),j2+IN(i,1),k2+IN(i,2));
          exponent_weighted = -1;
          if (SELECT(i2,j2,k2,i)) {
            tmp[i] -= SELECT(i2,j2,k2,i) *
              (this->*shell_function)(gc1,i1,j1,k1,gc2,i2-IN(i,0),j2-IN(i,1),k2-IN(i,2));
            }
          }
        }
      else {
        for (i=0; i<3; i++) tmp[i] = 0.0;
        }

      if (scale_shell_result) {
        for (i=0; i<3; i++) tmp[i] *= result_scale_factor;
        }

      for (i=0; i<3; i++) ctmp[i] = tmp[i];

      /* Increment the pointer to the xyz for the next atom. */
      ctmp += 3;
      END_FOR_GCCART_GS(index2)
    END_FOR_GCCART_GS(index1)

  accum_transform_1e_xyz(integral_,
                         cartesianbuffer, buff, gshell1, gshell2);
  }

/* This computes the basis function part of 
 * generic derivative integrals between functions
 * in two shells.  The result is accumulated in the buffer which is ordered
 * atom 0 x, y, z, atom 1, ... .
 * The function used to compute the nonderivative integrals is shell_function.
 */
void
Int1eV3::accum_shell_block_1der(double *buff, int ish, int jsh,
                                Ref<GaussianBasisSet> dercs, int centernum,
                                void (Int1eV3::*shell_block_function)
                                  (int gc1, int a, int gc2, int b,
                                   int gcsize2, int gcoff1, int gcoff2,
                                   double coef, double *buffer))
{
  int i;
  int gc1,gc2;
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int index1,index2;

  c1 = bs1_->shell_to_center(ish);
  c2 = bs2_->shell_to_center(jsh);

  int docenter1=0, docenter2=0;
  int equiv12 = (bs1_==bs2_)&&(c1==c2);
  int der1 = (bs1_==dercs)&&(c1==centernum);
  int der2 = (bs2_==dercs)&&(c2==centernum);
  if (!equiv12) {
    docenter1 = der1;
    docenter2 = der2;
    }
  else if (three_center) {
    int equiv123 = (bs1_==third_centers)&&(c1==third_centernum);
    if (!equiv123) {
      docenter1 = der1;
      docenter2 = der2;
      }
    }

  gshell1 = &bs1_->shell(ish);
  gshell2 = &bs2_->shell(jsh);
  int gcsize1 = gshell1->ncartesian();
  int gcsize2 = gshell2->ncartesian();
  memset(cartesianbuffer,0,sizeof(double)*gcsize1*gcsize2*3);

  if (!docenter1 && !docenter2) return;

  double coef;
  if (scale_shell_result) {
    coef = result_scale_factor;
    }
  else coef = 1.0;

  for (int xyz=0; xyz<3; xyz++) {
      A[xyz] = bs1_->r(c1,xyz);
      B[xyz] = bs2_->r(c2,xyz);
    }
  int gcoff1 = 0;
  for (gc1=0; gc1<gshell1->ncontraction(); gc1++) {
    int a = gshell1->am(gc1);
    int sizea = INT_NCART_NN(a);
    int sizeap1 = INT_NCART_NN(a+1);
    int sizeam1 = INT_NCART(a-1);
    int gcoff2 = 0;
    for (gc2=0; gc2<gshell2->ncontraction(); gc2++) {
      int b = gshell2->am(gc2);
      int sizeb = INT_NCART_NN(b);
      int sizebp1 = INT_NCART_NN(b+1);
      int sizebm1 = INT_NCART(b-1);
      /* Derivative wrt first shell. */
      if (docenter1) {
        exponent_weighted = 0;
        memset(cartesianbuffer_scratch,0,sizeof(double)*sizeap1*sizeb);
        (this->*shell_block_function)(gc1, a+1, gc2, b,
                                      sizeb, 0, 0,
                                      coef, cartesianbuffer_scratch);
        index1=0;
        FOR_CART(i1,j1,k1,a) {
          index2=0;
          FOR_CART(i2,j2,k2,b) {
            double *ctmp = &cartesianbuffer[((index1+gcoff1)*gcsize2
                                             +index2+gcoff2)*3];
            for (i=0; i<3; i++) {
              int ind = INT_CARTINDEX(a+1,i1+IN(i,0),j1+IN(i,1));
              *ctmp++ += 2.0*cartesianbuffer_scratch[ind*sizeb+index2];
              }
            index2++;
            } END_FOR_CART;
          index1++;
          } END_FOR_CART;

        if (a) {
          exponent_weighted = -1;
          memset(cartesianbuffer_scratch,0,sizeof(double)*sizeam1*sizeb);
          (this->*shell_block_function)(gc1, a-1, gc2, b,
                                        sizeb, 0, 0,
                                        coef, cartesianbuffer_scratch);
          index1=0;
          FOR_CART(i1,j1,k1,a) {
            index2=0;
            FOR_CART(i2,j2,k2,b) {
              double *ctmp = &cartesianbuffer[((index1+gcoff1)*gcsize2
                                               +index2+gcoff2)*3];
              for (i=0; i<3; i++) {
                int sel = SELECT(i1,j1,k1,i);
                if (sel) {
                  int ind = INT_CARTINDEX(a-1,i1-IN(i,0),j1-IN(i,1));
                  ctmp[i] -= sel * cartesianbuffer_scratch[ind*sizeb+index2];
                  }
                }
              index2++;
              } END_FOR_CART;
            index1++;
            } END_FOR_CART;
          }
        }
      if (docenter2) {
        /* Derviative wrt second shell. */
        exponent_weighted = 1;
        memset(cartesianbuffer_scratch,0,sizeof(double)*sizea*sizebp1);
        (this->*shell_block_function)(gc1, a, gc2, b+1,
                                      sizebp1, 0, 0,
                                      coef, cartesianbuffer_scratch);
        index1=0;
        FOR_CART(i1,j1,k1,a) {
          index2=0;
          FOR_CART(i2,j2,k2,b) {
            double *ctmp = &cartesianbuffer[((index1+gcoff1)*gcsize2
                                             +index2+gcoff2)*3];
            for (i=0; i<3; i++) {
              int ind = INT_CARTINDEX(b+1,i2+IN(i,0),j2+IN(i,1));
              *ctmp++ += 2.0*cartesianbuffer_scratch[index1*sizebp1+ind];
              }
            index2++;
            } END_FOR_CART;
          index1++;
          } END_FOR_CART;

        if (b) {
          exponent_weighted = -1;
          memset(cartesianbuffer_scratch,0,sizeof(double)*sizea*sizebm1);
          (this->*shell_block_function)(gc1, a, gc2, b-1,
                                        sizebm1, 0, 0,
                                        coef, cartesianbuffer_scratch);
          index1=0;
          FOR_CART(i1,j1,k1,a) {
            index2=0;
            FOR_CART(i2,j2,k2,b) {
              double *ctmp = &cartesianbuffer[((index1+gcoff1)*gcsize2
                                               +index2+gcoff2)*3];
              for (i=0; i<3; i++) {
                int sel = SELECT(i2,j2,k2,i);
                if (sel) {
                  int ind = INT_CARTINDEX(b-1,i2-IN(i,0),j2-IN(i,1));
                  ctmp[i] -= sel * cartesianbuffer_scratch[index1*sizebm1+ind];
                  }
                }
              index2++;
              } END_FOR_CART;
            index1++;
            } END_FOR_CART;
          }
        }
      gcoff2 += INT_NCART_NN(b);
      }
    gcoff1 += INT_NCART_NN(a);
    }

  accum_transform_1e_xyz(integral_,
                         cartesianbuffer, buff, gshell1, gshell2);

#if DEBUG_NUC_SHELL_DER
  double *fastbuff = cartesianbuffer;
  cartesianbuffer = new double[gcsize1*gcsize2*3];
  double *junkbuff = new double[gcsize1*gcsize2*3];
  memset(junkbuff,0,sizeof(double)*gcsize1*gcsize2*3);
  accum_shell_1der(junkbuff,ish,jsh,dercs,centernum,
                   &Int1eV3::comp_shell_nuclear);
  delete[] junkbuff;
  int index = 0;
  for (i=0; i<gcsize1; i++) {
      for (int j=0; j<gcsize2; j++, index++) {
          ExEnv::outn() << scprintf(" (%d %d): % 18.15f % 18.15f",
                           i,j,cartesianbuffer[index],fastbuff[index]);
          if (fabs(cartesianbuffer[index]-fastbuff[index])>1.0e-12) {
              ExEnv::outn() << " **";
            }
          ExEnv::outn() << endl;
        }
    }
  delete[] cartesianbuffer;
  cartesianbuffer = fastbuff;
#endif
  }

/* Compute the kinetic energy for the shell.  The arguments are the
 * cartesian exponents for centers 1 and 2.  The shell1 and shell2
 * globals are used. */
double
Int1eV3::comp_shell_kinetic(int gc1, int i1, int j1, int k1,
                          int gc2, int i2, int j2, int k2)
{
  int i,j,xyz;
  double result;
  double Pi;
  double oozeta;
  double AmB,AmB2;
  double tmp;

  /* Loop over the primitives in the shells. */
  result = 0.0;
  for (i=0; i<gshell1->nprimitive(); i++) {
    for (j=0; j<gshell2->nprimitive(); j++) {

      /* Compute the intermediates. */
      oo2zeta_a = 0.5/gshell1->exponent(i);
      oo2zeta_b = 0.5/gshell2->exponent(j);
      oozeta = 1.0/(gshell1->exponent(i) + gshell2->exponent(j));
      oo2zeta = 0.5*oozeta;
      xi = oozeta * gshell1->exponent(i) * gshell2->exponent(j);
      AmB2 = 0.0;
      for (xyz=0; xyz<3; xyz++) {
        Pi = oozeta*(gshell1->exponent(i) * A[xyz]
                     + gshell2->exponent(j) * B[xyz]);
        PmA[xyz] = Pi - A[xyz];
        PmB[xyz] = Pi - B[xyz];
        AmB = A[xyz] - B[xyz];
        AmB2 += AmB*AmB;
        }
      /* The s integral kinetic energy. */
      ss =   pow(3.141592653589793/(gshell1->exponent(i)
                                    +gshell2->exponent(j)),1.5)
           * exp(- xi * AmB2);
      sTs =  ss
           * xi
           * (3.0 - 2.0 * xi * AmB2);
      tmp     =  gshell1->coefficient_unnorm(gc1,i)
               * gshell2->coefficient_unnorm(gc2,j)
               * comp_prim_kinetic(i1,j1,k1,i2,j2,k2);
      if (exponent_weighted == 0) tmp *= gshell1->exponent(i);
      else if (exponent_weighted == 1) tmp *= gshell2->exponent(j);
      result += tmp;
      }
    }

  return result;
  }

double
Int1eV3::comp_prim_kinetic(int i1, int j1, int k1,
                         int i2, int j2, int k2)
{
  double tmp;
  double result;

  if (i1) {
    result = PmA[0] * comp_prim_kinetic(i1-1,j1,k1,i2,j2,k2);
    if (i1>1) result += oo2zeta*(i1-1)*comp_prim_kinetic(i1-2,j1,k1,i2,j2,k2);
    if (i2) result += oo2zeta*i2*comp_prim_kinetic(i1-1,j1,k1,i2-1,j2,k2);
    tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
    if (i1>1) tmp -= oo2zeta_a*(i1-1)*comp_prim_overlap(i1-2,j1,k1,i2,j2,k2);
    result += 2.0 * xi * tmp;
    return result;
    }
  if (j1) {
    result = PmA[1] * comp_prim_kinetic(i1,j1-1,k1,i2,j2,k2);
    if (j1>1) result += oo2zeta*(j1-1)*comp_prim_kinetic(i1,j1-2,k1,i2,j2,k2);
    if (j2) result += oo2zeta*j2*comp_prim_kinetic(i1,j1-1,k1,i2,j2-1,k2);
    tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
    if (j1>1) tmp -= oo2zeta_a*(j1-1)*comp_prim_overlap(i1,j1-2,k1,i2,j2,k2);
    result += 2.0 * xi * tmp;
    return result;
    }
  if (k1) {
    result = PmA[2] * comp_prim_kinetic(i1,j1,k1-1,i2,j2,k2);
    if (k1>1) result += oo2zeta*(k1-1)*comp_prim_kinetic(i1,j1,k1-2,i2,j2,k2);
    if (k2) result += oo2zeta*k2*comp_prim_kinetic(i1,j1,k1-1,i2,j2,k2-1);
    tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
    if (k1>1) tmp -= oo2zeta_a*(k1-1)*comp_prim_overlap(i1,j1,k1-2,i2,j2,k2);
    result += 2.0 * xi * tmp;
    return result;
    }
  if (i2) {
    result = PmB[0] * comp_prim_kinetic(i1,j1,k1,i2-1,j2,k2);
    if (i2>1) result += oo2zeta*(i2-1)*comp_prim_kinetic(i1,j1,k1,i2-2,j2,k2);
    if (i1) result += oo2zeta*i1*comp_prim_kinetic(i1-1,j1,k1,i2-1,j2,k2);
    tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
    if (i2>1) tmp -= oo2zeta_b*(i2-1)*comp_prim_overlap(i1,j1,k1,i2-2,j2,k2);
    result += 2.0 * xi * tmp;
    return result;
    }
  if (j2) {
    result = PmB[1] * comp_prim_kinetic(i1,j1,k1,i2,j2-1,k2);
    if (j2>1) result += oo2zeta*(j2-1)*comp_prim_kinetic(i1,j1,k1,i2,j2-2,k2);
    if (j1) result += oo2zeta*i1*comp_prim_kinetic(i1,j1-1,k1,i2,j2-1,k2);
    tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
    if (j2>1) tmp -= oo2zeta_b*(j2-1)*comp_prim_overlap(i1,j1,k1,i2,j2-2,k2);
    result += 2.0 * xi * tmp;
    return result;
    }
  if (k2) {
    result = PmB[2] * comp_prim_kinetic(i1,j1,k1,i2,j2,k2-1);
    if (k2>1) result += oo2zeta*(k2-1)*comp_prim_kinetic(i1,j1,k1,i2,j2,k2-2);
    if (k1) result += oo2zeta*i1*comp_prim_kinetic(i1,j1,k1-1,i2,j2,k2-1);
    tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
    if (k2>1) tmp -= oo2zeta_b*(k2-1)*comp_prim_overlap(i1,j1,k1,i2,j2,k2-2);
    result += 2.0 * xi * tmp;
    return result;
    }

  return sTs;
  }

/* --------------------------------------------------------------- */
/* ------------- Routines for the nuclear attraction ------------- */
/* --------------------------------------------------------------- */

/* This computes the nuclear attraction derivative integrals between functions
 * in two shells.  The result is accumulated in the buffer which is ordered
 * atom 0 x, y, z, atom 1, ... .
 */
void
Int1eV3::int_accum_shell_nuclear_1der(int ish, int jsh,
                                      Ref<GaussianBasisSet> dercs, int centernum)
{
  int_accum_shell_nuclear_hf_1der(ish,jsh,dercs,centernum);
  int_accum_shell_nuclear_nonhf_1der(ish,jsh,dercs,centernum);
  }

/* A correction to the Hellman-Feynman part is computed which
 * is not included in the original HF routine.  This is only needed