comp1e.cc

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

 * The result is placed in the buffer.
 */
void
Int1eV3::int_accum_shell_point_charge(int ish, int jsh,
                                      int ncharge, const double* charge,
                                      const double*const* position)
{
  int i;
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int index;
  int gc1,gc2;
  int cart1,cart2;
  double tmp;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_pointcharge: one electron routines are not init'ed\n");
    exit(1);
    }

  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);
  index = 0;

  FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
      /* Loop thru the point charges. */
      tmp = 0.0;
      for (i=0; i<ncharge; i++) {
          for (int xyz=0; xyz<3; xyz++) {
              C[xyz] = position[i][xyz];
            }
        tmp -=  comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
                       * charge[i];
        }
      cartesianbuffer[index] = tmp;
      index++;
      END_FOR_GCCART_GS(cart2)
    END_FOR_GCCART_GS(cart1)

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

/* This computes the integrals between functions in two shells for
 * a point charge interaction operator.
 * The result is placed in the buffer.
 */
void
Int1eV3::point_charge(int ish, int jsh,
                      int ncharge,
                      const double* charge, const double*const* position)
{
  int i;
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int index;
  int gc1,gc2;
  int cart1,cart2;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("Int1eV3::point_charge: one electron routines are not init'ed\n");
    exit(1);
    }

  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);
  index = 0;

  FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
      cartesianbuffer[index] = 0.0;
      /* Loop thru the point charges. */
      for (i=0; i<ncharge; i++) {
        for (int xyz=0; xyz<3; xyz++) {
          C[xyz] = position[i][xyz];
          }
        cartesianbuffer[index] -= comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
                                * charge[i];
        }
      index++;
      END_FOR_GCCART_GS(cart2)
    END_FOR_GCCART_GS(cart1)

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


/* This computes the 1e Hamiltonian integrals between functions in two shells.
 * The result is placed in the buffer.
 */
void
Int1eV3::hcore(int ish, int jsh)
{
  int i;
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int index;
  int cart1,cart2;
  int gc1,gc2;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("hcore: one electron routines are not init'ed\n");
    exit(1);
    }

  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);

  index = 0;
  FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
      cartesianbuffer[index] = comp_shell_kinetic(gc1,i1,j1,k1,gc2,i2,j2,k2);
      /* Loop thru the centers on bs1_. */
      for (i=0; i<bs1_->ncenter(); i++) {
        for (int xyz=0; xyz<3; xyz++) {
          C[xyz] = bs1_->r(i,xyz);
          }
        cartesianbuffer[index] -= comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
                                * bs1_->molecule()->charge(i);
        }
      /* Loop thru the centers on bs2_ if necessary. */
      if (bs2_ != bs1_) {
        for (i=0; i<bs2_->ncenter(); i++) {
          for (int xyz=0; xyz<3; xyz++) {
            C[xyz] = bs2_->r(i,xyz);
            }
          cartesianbuffer[index]-=comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
                                * bs2_->molecule()->charge(i);
          }
        }
      index++;
      END_FOR_GCCART_GS(cart2)
    END_FOR_GCCART_GS(cart1)

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

/* This computes the 1e Hamiltonian deriv ints between functions in two shells.
 * The result is placed in the buffer.
 */
void
Int1eV3::hcore_1der(int ish, int jsh,
                    int idercs, int centernum)
{
  int i;
  int c1,c2;
  int ni,nj;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_hcore: one electron routines are not init'ed\n");
    exit(1);
    }

  Ref<GaussianBasisSet> dercs;
  if (idercs == 0) dercs = bs1_;
  else dercs = bs2_;

  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);

  ni = gshell1->nfunction();
  nj = gshell2->nfunction();

  for (i=0; i<ni*nj*3; i++) {
    buff[i] = 0.0;
    }

  int_accum_shell_nuclear_1der(ish,jsh,dercs,centernum);
  int_accum_shell_kinetic_1der(ish,jsh,dercs,centernum);
  }

/* This computes the kinetic deriv ints between functions in two shells.
 * The result is placed in the buffer.
 */
void
Int1eV3::kinetic_1der(int ish, int jsh,
                      int idercs, int centernum)
{
  int i;
  int c1,c2;
  int ni,nj;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_kinetic: one electron routines are not init'ed\n");
    exit(1);
    }

  Ref<GaussianBasisSet> dercs;
  if (idercs == 0) dercs = bs1_;
  else dercs = bs2_;

  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);

  ni = gshell1->nfunction();
  nj = gshell2->nfunction();

  for (i=0; i<ni*nj*3; i++) {
    buff[i] = 0.0;
    }

  int_accum_shell_kinetic_1der(ish,jsh,dercs,centernum);
  }

/* This computes the nuclear deriv ints between functions in two shells.
 * The result is placed in the buffer.
 */
void
Int1eV3::nuclear_1der(int ish, int jsh, int idercs, int centernum)
{
  int i;
  int c1,c2;
  int ni,nj;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_nuclear: one electron routines are not init'ed\n");
    exit(1);
    }

  Ref<GaussianBasisSet> dercs;
  if (idercs == 0) dercs = bs1_;
  else dercs = bs2_;

  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);

  ni = gshell1->nfunction();
  nj = gshell2->nfunction();

  for (i=0; i<ni*nj*3; i++) {
    buff[i] = 0.0;
    }

  int_accum_shell_nuclear_1der(ish,jsh,dercs,centernum);
  }

/* This computes the nuclear deriv ints between functions in two shells.
 * Only the Hellman-Feynman part is computed.
 * The result is placed in the buffer.
 */
void
Int1eV3::int_shell_nuclear_hf_1der(int ish, int jsh,
                                   Ref<GaussianBasisSet> dercs, int centernum)
{
  int i;
  int c1,c2;
  int ni,nj;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_nuclear_hf_1der: one electron routines are not init'ed\n");
    exit(1);
    }

  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);

  ni = gshell1->nfunction();
  nj = gshell2->nfunction();

  for (i=0; i<ni*nj*3; i++) {
    buff[i] = 0.0;
    }

  int_accum_shell_nuclear_hf_1der(ish,jsh,dercs,centernum);
  }

/* This computes the nuclear deriv ints between functions in two shells.
 * Only the non Hellman-Feynman part is computed.
 * The result is placed in the buffer.
 */
void
Int1eV3::int_shell_nuclear_nonhf_1der(int ish, int jsh,
                                      Ref<GaussianBasisSet> dercs, int centernum)
{
  int i;
  int c1,c2;
  int ni,nj;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_nuclear_nonhf_1der: one electron routines are not init'ed\n");
    exit(1);
    }

  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);

  ni = gshell1->nfunction();
  nj = gshell2->nfunction();

#if 0
  ExEnv::outn() << scprintf("int_shell_nuclear_nonhf_1der: zeroing %d doubles in buff\n",ni*nj*3);
#endif
  for (i=0; i<ni*nj*3; i++) {
    buff[i] = 0.0;
    }

  int_accum_shell_nuclear_nonhf_1der(ish,jsh,dercs,centernum);
  }

/* Compute the nuclear attraction 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_nuclear(int gc1, int i1, int j1, int k1,
                          int gc2, int i2, int j2, int k2)
{
  int i,j,k,xyz;
  double result;
  double Pi;
  double oozeta;
  double AmB,AmB2;
  double PmC2;
  double auxcoef;
  int am;
  double tmp;

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

  /* 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. */
      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 nuclear attraction integral. */
      tmp     =  gshell1->coefficient_unnorm(gc1,i)
               * gshell2->coefficient_unnorm(gc2,j)
               * comp_prim_nuclear(i1,j1,k1,i2,j2,k2,0);

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

      result += tmp;
      }
    }

  /* printf("comp_shell_nuclear(%d,%d,%d,%d,%d,%d): result = % 12.8lf\n",
   *         i1,j1,k1,i2,j2,k2,result);
   */
  return result;
  }

/* Compute the nuclear attraction for the shell.  The arguments are the