comp1e.cc

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

 * cartesian exponents for centers 1 and 2.  The shell1 and shell2
 * globals are used. */
void
Int1eV3::comp_shell_block_nuclear(int gc1, int a, int gc2, int b,
                                  int gcsize2, int gcoff1, int gcoff2,
                                  double coef, double *buffer)
{
  int i,j,k,xyz;
  double Pi;
  double oozeta;
  double AmB,AmB2;
  double PmC2;
  double auxcoef;
  double tmp;
  int am = a + b;
  int size1 = INT_NCART_NN(a);
  int size2 = INT_NCART_NN(b);

#if DEBUG_NUC_SHELL
  double *shellints = new double[size1*size2];
  memset(shellints,0,sizeof(double)*size1*size2);
#endif

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

      /* Compute the intermediates. */
      zeta = gshell1->exponent(i) + gshell2->exponent(j);
      oozeta = 1.0/zeta;
      oo2zeta = 0.5*oozeta;
      AmB2 = 0.0;
      PmC2 = 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];
        PmC[xyz] = Pi - C[xyz];
        AmB = A[xyz] - B[xyz];
        AmB2 += AmB*AmB;
        PmC2 += PmC[xyz]*PmC[xyz];
        }

      /* The auxillary integral coeficients. */
      auxcoef =   2.0 * 3.141592653589793/(gshell1->exponent(i)
                                           +gshell2->exponent(j))
           * exp(- oozeta * gshell1->exponent(i)
                 * gshell2->exponent(j) * AmB2);

      /* The Fm(U) intermediates. */
      fjttable_ = fjt_->values(am,zeta*PmC2);

      /* Convert the Fm(U) intermediates into the auxillary
       * nuclear attraction integrals. */
      for (k=0; k<=am; k++) {
        fjttable_[k] *= auxcoef;
        }

      /* Compute the primitive nuclear attraction integral. */
      comp_prim_block_nuclear(a,b);

      tmp =  gshell1->coefficient_unnorm(gc1,i)
             * gshell2->coefficient_unnorm(gc2,j)
             * coef;

      if (exponent_weighted == 0) tmp *= gshell1->exponent(i);
      else if (exponent_weighted == 1) tmp *= gshell2->exponent(j);

#if DEBUG_NUC_SHELL
      double *tprimbuffer = inter(a,b,0);
      double *tmpshellints = shellints;
      for (int ip=0; ip<size1; ip++) {
        for (int jp=0; jp<size2; jp++) {
          *tmpshellints++ += tmp * *tprimbuffer++ / coef;
          }
        }
#endif
      double *primbuffer = inter(a,b,0);
      for (int ip=0; ip<size1; ip++) {
        for (int jp=0; jp<size2; jp++) {
          //ExEnv::outn() << scprintf("buffer[%d] += %18.15f",
          //                 (ip+gcoff1)*gcsize2+jp+gcoff2,
          //                 tmp * *primbuffer)
          //     << endl;
          buffer[(ip+gcoff1)*gcsize2+jp+gcoff2] += tmp * *primbuffer++;
          }
        }
      }
    }

#if DEBUG_NUC_SHELL
#  if DEBUG_NUC_SHELL > 1
  ExEnv::outn() << scprintf("GC = (%d %d), A = %d, B = %d", gc1, gc2, a, b)
       << endl;
#  endif
  int i1,j1,k1;
  int i2,j2,k2;
  int ip = 0;
  double *tmpshellints = shellints;
  FOR_CART(i1,j1,k1,a) {
    int jp = 0;
    FOR_CART(i2,j2,k2,b) {
      double fast = *tmpshellints++;
      double slow = comp_shell_nuclear(gc1, i1, j1, k1,
                                       gc2, i2, j2, k2);
      int bad = fabs(fast-slow)>1.0e-12;
      if (DEBUG_NUC_SHELL > 1 || bad) {
        ExEnv::outn() << scprintf("NUC SHELL: (%d %d %d) (%d %d %d): ",
                         i1,j1,k1, i2,j2,k2)
             << scprintf(" % 20.15f % 20.15f",
                         fast, slow);
        }
      if (bad) {
        ExEnv::outn() << " ****" << endl;
        }
      else if (DEBUG_NUC_SHELL > 1) {
        ExEnv::outn() << endl;
        }
      jp++;
      } END_FOR_CART;
    ip++;
    } END_FOR_CART;
  delete[] shellints;
#endif
  }

void
Int1eV3::comp_prim_block_nuclear(int a, int b)
{
  int im, ia, ib;
  int l = a + b;

  // fill in the ia+ib=0 integrals
  for (im=0; im<=l; im++) {
#if DEBUG_NUC_PRIM > 1
    ExEnv::outn() << "BUILD: M = " << im
         << " A = " << 0
         << " B = " << 0
         << endl;
#endif
    inter(0,0,im)[0] = fjttable_[im];
    }

  for (im=l-1; im>=0; im--) {
    int lm = l-im;
    // build the integrals for a = 0
    for (ib=1; ib<=lm && ib<=b; ib++) {
#if DEBUG_NUC_PRIM > 1
      ExEnv::outn() << "BUILD: M = " << im
           << " A = " << 0
           << " B = " << ib
           << endl;
#endif
      comp_prim_block_nuclear_build_b(ib,im);
      }
    for (ia=1; ia<=lm && ia<=a; ia++) {
      for (ib=0; ib<=lm-ia && ib<=b; ib++) {
#if DEBUG_NUC_PRIM > 1
        ExEnv::outn() << "BUILD: M = " << im
             << " A = " << ia
             << " B = " << ib
             << endl;
#endif
        comp_prim_block_nuclear_build_a(ia,ib,im);
        }
      }
    }

#if DEBUG_NUC_PRIM
  for (im=0; im<=l; im++) {
    int lm = l-im;
    for (ia=0; ia<=lm && ia<=a; ia++) {
      int na = INT_NCART_NN(a);
      for (ib=0; ib<=lm-ia && ib<=b; ib++) {
        int nb = INT_NCART_NN(b);
#if DEBUG_NUC_PRIM > 1
        ExEnv::outn() << "M = " << im
             << " A = " << ia
             << " B = " << ib
             << endl;
#endif
        double *buf = inter(ia,ib,im);
        int i1,j1,k1, i2,j2,k2;
        FOR_CART(i1,j1,k1,ia) {
          FOR_CART(i2,j2,k2,ib) {
            double fast = *buf++;
            double slow = comp_prim_nuclear(i1, j1, k1,
                                            i2, j2, k2, im);
            if (fast > 999.0) fast = 999.0;
            if (fast < -999.0) fast = -999.0;
            if (fabs(fast-slow)>1.0e-12) {
              ExEnv::outn() << scprintf("(%d %d %d) (%d %d %d) (%d): ",
                               i1,j1,k1, i2,j2,k2, im)
                   << scprintf(" % 20.15f % 20.15f",
                               fast, slow)
                   << endl;
              }
            } END_FOR_CART;
          } END_FOR_CART;
        }
      }
    }
#endif
}

void
Int1eV3::comp_prim_block_nuclear_build_a(int a, int b, int m)
{
  double *I000 = inter(a,b,m);
  double *I100 = inter(a-1,b,m);
  double *I101 = inter(a-1,b,m+1);
  double *I200 = (a>1?inter(a-2,b,m):0);
  double *I201 = (a>1?inter(a-2,b,m+1):0);
  double *I110 = (b?inter(a-1,b-1,m):0);
  double *I111 = (b?inter(a-1,b-1,m+1):0);
  int i1,j1,k1;
  int i2,j2,k2;
  int carta=0;
  int sizeb = INT_NCART_NN(b);
  int sizebm1 = INT_NCART_DEC(b,sizeb);
  FOR_CART(i1,j1,k1,a) {
    int cartb=0;
    FOR_CART(i2,j2,k2,b) {
      double result = 0.0;
      if (i1) {
        int am1 = INT_CARTINDEX(a-1,i1-1,j1);
        result  = PmA[0] * I100[am1*sizeb+cartb];
        result -= PmC[0] * I101[am1*sizeb+cartb];
        if (i1>1) {
          int am2 = INT_CARTINDEX(a-2,i1-2,j1);
          result += oo2zeta * (i1-1)
                    *(I200[am2*sizeb+cartb] - I201[am2*sizeb+cartb]);
          }
        if (i2) {
          int bm1 = INT_CARTINDEX(b-1,i2-1,j2);
          result += oo2zeta * i2
                    *(I110[am1*sizebm1+bm1] - I111[am1*sizebm1+bm1]);
          }
        }
      else if (j1) {
        int am1 = INT_CARTINDEX(a-1,i1,j1-1);
        result  = PmA[1] * I100[am1*sizeb+cartb];
        result -= PmC[1] * I101[am1*sizeb+cartb];
        if (j1>1) {
          int am2 = INT_CARTINDEX(a-2,i1,j1-2);
          result += oo2zeta * (j1-1)
                    *(I200[am2*sizeb+cartb] - I201[am2*sizeb+cartb]);
          }
        if (j2) {
          int bm1 = INT_CARTINDEX(b-1,i2,j2-1);
          result += oo2zeta * j2
                    *(I110[am1*sizebm1+bm1] - I111[am1*sizebm1+bm1]);
          }
        }
      else if (k1) {
        int am1 = INT_CARTINDEX(a-1,i1,j1);
        result  = PmA[2] * I100[am1*sizeb+cartb];
        result -= PmC[2] * I101[am1*sizeb+cartb];
        if (k1>1) {
          int am2 = INT_CARTINDEX(a-2,i1,j1);
          result += oo2zeta * (k1-1)
                    *(I200[am2*sizeb+cartb] - I201[am2*sizeb+cartb]);
          }
        if (k2) {
          int bm1 = INT_CARTINDEX(b-1,i2,j2);
          result += oo2zeta * k2
                    *(I110[am1*sizebm1+bm1] - I111[am1*sizebm1+bm1]);
          }
        }
      I000[carta*sizeb+cartb] = result;
      cartb++;
      } END_FOR_CART;
    carta++;
    } END_FOR_CART;
}

void
Int1eV3::comp_prim_block_nuclear_build_b(int b, int m)
{
  double *I000 = inter(0,b,m);
  double *I010 = inter(0,b-1,m);
  double *I011 = inter(0,b-1,m+1);
  double *I020 = (b>1?inter(0,b-2,m):0);
  double *I021 = (b>1?inter(0,b-2,m+1):0);
  int i2,j2,k2;
  int cartb=0;
  FOR_CART(i2,j2,k2,b) {
    double result = 0.0;

    if (i2) {
      int bm1 = INT_CARTINDEX(b-1,i2-1,j2);
      result  = PmB[0] * I010[bm1];
      result -= PmC[0] * I011[bm1];
      if (i2>1) {
        int bm2 = INT_CARTINDEX(b-2,i2-2,j2);
        result += oo2zeta * (i2-1) * (I020[bm2] - I021[bm2]);
        }
      }
    else if (j2) {
      int bm1 = INT_CARTINDEX(b-1,i2,j2-1);
      result  = PmB[1] * I010[bm1];
      result -= PmC[1] * I011[bm1];
      if (j2>1) {
        int bm2 = INT_CARTINDEX(b-2,i2,j2-2);
        result += oo2zeta * (j2-1) * (I020[bm2] - I021[bm2]);
        }
      }
    else if (k2) {
      int bm1 = INT_CARTINDEX(b-1,i2,j2);
      result  = PmB[2] * I010[bm1];
      result -= PmC[2] * I011[bm1];
      if (k2>1) {
        int bm2 = INT_CARTINDEX(b-2,i2,j2);
        result += oo2zeta * (k2-1) * (I020[bm2] - I021[bm2]);
        }
      }

    I000[cartb] = result;
    cartb++;
    } END_FOR_CART;
}

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

  if (i1) {
    result  = PmA[0] * comp_prim_nuclear(i1-1,j1,k1,i2,j2,k2,m);
    result -= PmC[0] * comp_prim_nuclear(i1-1,j1,k1,i2,j2,k2,m+1);
    if (i1>1) result += oo2zeta * (i1-1)
                       * (  comp_prim_nuclear(i1-2,j1,k1,i2,j2,k2,m)
                          - comp_prim_nuclear(i1-2,j1,k1,i2,j2,k2,m+1));
    if (i2) result += oo2zeta * i2
                     * (  comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m)
                        - comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m+1));
    }
  else if (j1) {
    result  = PmA[1] * comp_prim_nuclear(i1,j1-1,k1,i2,j2,k2,m);
    result -= PmC[1] * comp_prim_nuclear(i1,j1-1,k1,i2,j2,k2,m+1);
    if (j1>1) result += oo2zeta * (j1-1)
                       * (  comp_prim_nuclear(i1,j1-2,k1,i2,j2,k2,m)
                          - comp_prim_nuclear(i1,j1-2,k1,i2,j2,k2,m+1));
    if (j2) result += oo2zeta * j2
                     * (  comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m)
                        - comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m+1));
    }
  else if (k1) {
    result  = PmA[2] * comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2,m);
    result -= PmC[2] * comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2,m+1);
    if (k1>1) result += oo2zeta * (k1-1)
                       * (  comp_prim_nuclear(i1,j1,k1-2,i2,j2,k2,m)
                          - comp_prim_nuclear(i1,j1,k1-2,i2,j2,k2,m+1));
    if (k2) result += oo2zeta * k2
                     * (  comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m)
                        - comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m+1));
    }
  else if (i2) {
    result  = PmB[0] * comp_prim_nuclear(i1,j1,k1,i2-1,j2,k2,m);
    result -= PmC[0] * comp_prim_nuclear(i1,j1,k1,i2-1,j2,k2,m+1);
    if (i2>1) result += oo2zeta * (i2-1)
                       * (  comp_prim_nuclear(i1,j1,k1,i2-2,j2,k2,m)
                          - comp_prim_nuclear(i1,j1,k1,i2-2,j2,k2,m+1));
    if (i1) result += oo2zeta * i1
                     * (  comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m)
                        - comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m+1));
    }
  else if (j2) {
    result  = PmB[1] * comp_prim_nuclear(i1,j1,k1,i2,j2-1,k2,m);
    result -= PmC[1] * comp_prim_nuclear(i1,j1,k1,i2,j2-1,k2,m+1);
    if (j2>1) result += oo2zeta * (j2-1)
                       * (  comp_prim_nuclear(i1,j1,k1,i2,j2-2,k2,m)
                          - comp_prim_nuclear(i1,j1,k1,i2,j2-2,k2,m+1));
    if (j1) result += oo2zeta * j1
                     * (  comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m)
                        - comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m+1));
    }
  else if (k2) {
    result  = PmB[2] * comp_prim_nuclear(i1,j1,k1,i2,j2,k2-1,m);
    result -= PmC[2] * comp_prim_nuclear(i1,j1,k1,i2,j2,k2-1,m+1);
    if (k2>1) result += oo2zeta * (k2-1)
                       * (  comp_prim_nuclear(i1,j1,k1,i2,j2,k2-2,m)
                          - comp_prim_nuclear(i1,j1,k1,i2,j2,k2-2,m+1));
    if (k1) result += oo2zeta * k1
                     * (  comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m)
                        - comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m+1));
    }
  else result = fjttable_[m];

  return result;
  }

/* Compute the electric field integral for the shell.  The arguments are the
 * the electric field vector, the
 * cartesian exponents for centers 1 and 2.  The shell1 and shell2
 * globals are used. */
void
Int1eV3::comp_shell_efield(double *efield,
                           int gc1, int i1, int j1, int k1,
                           int gc2, int i2, int j2, int k2)
{
  int i,j,k,xyz;
  double result[3];
  double Pi;
  double oozeta;
  double AmB,AmB2;
  double PmC2;
  double auxcoef;
  int am;

  am = i1+j1+k1+i2+j2+k2;

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

      /* Compute the intermediates. */