gaussshval.cc

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

                  if (j.c()) cart_g_values[i_cart] += j.c() * norm_precon
                    * xs[j.a()] * ys[j.b()] * zs[j.c()-1];
                  i_cart++;
                }
              SphericalTransformIter *ti = sivec[l[i]];
              int n = ti->n();
              memset(&g_values[i_grad], 0, sizeof(double)*n*3);
              for (ti->start(); ti->ready(); ti->next()) {
                  double coef = ti->coef();
                  int pi = ti->pureindex();
                  int ci = ti->cartindex();
                  for (int xyz=0; xyz<3; xyz++) {
                      g_values[i_grad + pi*3 + xyz]
                          += coef * cart_g_values[ci*3 + xyz];
                    }
                }
              i_grad += 3*n;
            }
        }
    }

  // compute the hessian of the shell values
  if (h_values) {
      int i_hess=0;                // Basis function counter
      for (i=0; i<ncon; i++) {
          // handle s functions with a special case to speed things up
          if (l[i] == 0) {
              double norm_precon_g = precon_g[i];
              double norm_precon_h = precon_h[i];
              // xx
              h_values[i_hess] = norm_precon_h*xs[2] - norm_precon_g;
              i_hess++;
              // yx
              h_values[i_hess] = norm_precon_h*xs[1]*ys[1];
              i_hess++;
              // yy
              h_values[i_hess] = norm_precon_h*ys[2] - norm_precon_g;
              i_hess++;
              // zx
              h_values[i_hess] = norm_precon_h*zs[1]*xs[1];
              i_hess++;
              // zy
              h_values[i_hess] = norm_precon_h*zs[1]*ys[1];
              i_hess++;
              // zz
              h_values[i_hess] = norm_precon_h*zs[2] - norm_precon_g;
              i_hess++;
            }
          else {
              double *cart_h;
              double tmp_cart_h[6*((MAX_AM+1)*(MAX_AM+2))/2];
              if (!puream[i]) {
                  cart_h = &h_values[i_hess];
                }
              else {
                  cart_h = tmp_cart_h;
                }
              CartesianIter *jp = civec[l[i]];
              CartesianIter& j = *jp;
              int i_cart = 0;
              for (j.start(); j; j.next()) {
                  double pre = precon[i];
                  double pre_g = - precon_g[i];
                  double pre_h = precon_h[i];
                  int a = j.a();
                  int b = j.b();
                  int c = j.c();
                  // xx
                  cart_h[i_cart] = pre_h * xs[a+2]*ys[b]*zs[c]
                                 + pre_g * (a+1) * xs[a]*ys[b]*zs[c];
                  if (a>0) {
                      cart_h[i_cart] += pre_g * a*xs[a]*ys[b]*zs[c];
                      if (a>1) cart_h[i_cart] += pre * a*(a-1)
                                               * xs[a-2]*ys[b]*zs[c];
                    }
                  i_cart++;

                  // yx
                  cart_h[i_cart] = pre_h * xs[a+1]*ys[b+1]*zs[c];
                  if (a>0)
                      cart_h[i_cart] += pre_g * a * xs[a-1]*ys[b+1]*zs[c];
                  if (b>0)
                      cart_h[i_cart] += pre_g * b * xs[a+1]*ys[b-1]*zs[c];
                  if (a>0 && b>0)
                      cart_h[i_cart] += pre * a*b * xs[a-1]*ys[b-1]*zs[c];
                  i_cart++;

                  // yy
                  cart_h[i_cart] = pre_h * xs[a]*ys[b+2]*zs[c]
                                 + pre_g * (b+1) * xs[a]*ys[b]*zs[c];
                  if (b>0) {
                      cart_h[i_cart] += pre_g * b*xs[a]*ys[b]*zs[c];
                      if (b>1) cart_h[i_cart] += pre * b*(b-1)
                                               * xs[a]*ys[b-2]*zs[c];
                    }
                  i_cart++;

                  // zx
                  cart_h[i_cart] = pre_h * xs[a+1]*ys[b]*zs[c+1];
                  if (a>0)
                      cart_h[i_cart] += pre_g * a * xs[a-1]*ys[b]*zs[c+1];
                  if (c>0)
                      cart_h[i_cart] += pre_g * c * xs[a+1]*ys[b]*zs[c-1];
                  if (a>0 && c>0)
                      cart_h[i_cart] += pre * a*c * xs[a-1]*ys[b]*zs[c-1];
                  i_cart++;

                  // zy
                  cart_h[i_cart] = pre_h * xs[a]*ys[b+1]*zs[c+1];
                  if (c>0)
                      cart_h[i_cart] += pre_g * c * xs[a]*ys[b+1]*zs[c-1];
                  if (b>0)
                      cart_h[i_cart] += pre_g * b * xs[a]*ys[b-1]*zs[c+1];
                  if (c>0 && b>0)
                      cart_h[i_cart] += pre * c*b * xs[a]*ys[b-1]*zs[c-1];
                  i_cart++;

                  // zz
                  cart_h[i_cart] = pre_h * xs[a]*ys[b]*zs[c+2]
                                 + pre_g * (c+1) * xs[a]*ys[b]*zs[c];
                  if (c>0) {
                      cart_h[i_cart] += pre_g * c*xs[a]*ys[b]*zs[c];
                      if (c>1) cart_h[i_cart] += pre * c*(c-1)
                                               * xs[a]*ys[b]*zs[c-2];
                    }
                  i_cart++;
                }
              if (puream[i]) {
                  SphericalTransformIter *ti = sivec[l[i]];
                  int n = ti->n();
                  memset(&h_values[i_hess], 0, sizeof(double)*n*6);
                  for (ti->start(); ti->ready(); ti->next()) {
                      double coef = ti->coef();
                      int pi = ti->pureindex();
                      int ci = ti->cartindex();
                      for (int xyz2=0; xyz2<6; xyz2++) {
                          h_values[i_hess + pi*6 + xyz2]
                              += coef * cart_h[ci*6 + xyz2];
                        }
                    }
                  i_hess += 6*n;
                }
              else {
                  i_hess += 3*(l[i]+1)*(l[i]+2);
                }
            }
        }
    }

  return i_basis;
}

int
GaussianShell::test_monobound(double &r, double &bound) const
{
  // compute the maximum angular momentum component of the shell
  // add one since derivatives will be needed
  int maxam = max_am() + 1;

  // check limitations
  if (nprim > MAX_NPRIM || ncon > MAX_NCON || maxam >= MAX_AM) {
      ExEnv::err0() << indent
           << "GaussianShell::gaussshval: limit exceeded:\n"
           << indent
           << scprintf(
               "ncon = %d (%d max) nprim = %d (%d max) maxam = %d (%d max)\n",
               ncon,MAX_NCON,nprim,MAX_NPRIM,maxam,MAX_AM-1);
      abort();
    }

  // loop variables
  int i,j;

  // precompute powers of r
  double rs[MAX_AM+1];
  rs[0] = 1.0;
  if (maxam>0) {
      rs[1] = r;
    }
  for (i=2; i<=maxam; i++) {
      rs[i] = rs[i-1]*r;
    }

  // precompute r*r
  double r2 = r*r;

  // precompute exponentials
  double exps[MAX_NPRIM];
  for (i=0; i<nprim; i++) {
      exps[i]=::exp(-r2*exp[i]);
    }

  // precompute contractions over exponentials
  double precon[MAX_NCON];
  for (i=0; i<ncon; i++) {
      precon[i] = 0.0;
      for (j=0; j<nprim; j++) {
          // using fabs since we want a monotonically decreasing bound
          precon[i] += fabs(coef[i][j]) * exps[j];
        }
    }

  // precompute contractions over exponentials with exponent weighting
  double precon_w[MAX_NCON];
  for (i=0; i<ncon; i++) {
      precon_w[i] = 0.0;
      for (j=0; j<nprim; j++) {
          precon_w[i] += exp[j] * fabs(coef[i][j]) * exps[j];
        }
    }

  double max_bound = 0.0;
  bound = 0.0;
  for (i=0; i<ncon; i++) {
      // using r^l since r^l >= x^a y^b z^c
      double component_bound = rs[l[i]]*precon[i];
      if (l[i] > 0) {
          double d1 = -2.0*rs[l[i]+1]*precon_w[i];
          double d2 = l[i]*rs[l[i]-1]*precon[i];
          if (d1+d2 > 0) {
              // This bound is no good since the contraction is increasing
              // at this position.  Move r out and return to let the driver
              // call again.
              double rold = r;
              r = sqrt(l[i]*precon[i]/(2.0*precon_w[i]));
              if (r<rold+0.01) r = rold+0.01;
              //ExEnv::outn() << "rejected at " << rold << " trying again at "
              //     << r << endl;
              return 1;
            }
        }
      if (component_bound > max_bound) {
          max_bound = component_bound;
      }
    }

  bound = max_bound;
  return 0;
}

double
GaussianShell::monobound(double r) const
{
  // doesn't work at r <= zero
  if (r<=0.001) r = 0.001;
  double b;
  while (test_monobound(r, b));
  return b;
}

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

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