comp1e.cc

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

 * if the real Hellman-Feynman forces are desired, because the sum
 * of the hf_1der and nonhf_1der forces are still correct.
 */
void
Int1eV3::int_accum_shell_nuclear_hfc_1der(int ish, int jsh,
                                          Ref<GaussianBasisSet> dercs,
                                          int centernum)
{
  /* If both ish and jsh are not on the der center,
   * then there's no correction. */
  if (!(  (bs1_==dercs)
        &&(bs2_==dercs)
        &&(bs1_->shell_to_center(ish)==centernum)
        &&(bs2_->shell_to_center(jsh)==centernum))) {
    return;
    }

  /* Compute the nuc attr part of the nuclear derivative for three equal
   * centers. */
  scale_shell_result = 1;
  result_scale_factor = -bs1_->molecule()->charge(centernum);
  for (int xyz=0; xyz<3; xyz++) {
      C[xyz] = bs1_->r(centernum,xyz);
    }
  accum_shell_efield(buff,ish,jsh);
  scale_shell_result = 0;

  }

/* 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, ... .  Only the Hellman-Feynman part is computed.
 */
void
Int1eV3::int_accum_shell_nuclear_hf_1der(int ish, int jsh,
                                         Ref<GaussianBasisSet> dercs,
                                         int centernum)
{

  /* If both ish and jsh are on the der center, then the contrib is zero. */
  if (  (bs1_==dercs)
      &&(bs2_==dercs)
      &&(bs1_->shell_to_center(ish)==centernum)
      &&(bs2_->shell_to_center(jsh)==centernum)) {
    return;
    }

  /* Compute the nuc attr part of the nuclear derivative. */
  if (bs1_ == dercs) {
    scale_shell_result = 1;
    result_scale_factor= -bs1_->molecule()->charge(centernum);
    for (int xyz=0; xyz<3; xyz++) {
        C[xyz] = bs1_->r(centernum,xyz);
      }
    //accum_shell_efield(buff,ish,jsh);
    accum_shell_block_efield(buff,ish,jsh);
    scale_shell_result = 0;
    }
  else if (bs2_ == dercs) {
    scale_shell_result = 1;
    result_scale_factor= -bs2_->molecule()->charge(centernum);
    for (int xyz=0; xyz<3; xyz++) {
        C[xyz] = bs2_->r(centernum,xyz);
      }
    //accum_shell_efield(buff,ish,jsh);
    accum_shell_block_efield(buff,ish,jsh);
    scale_shell_result = 0;
    }

  }

/* 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, ... .  Only the non Hellman-Feynman part is computed.
 */
void
Int1eV3::int_accum_shell_nuclear_nonhf_1der(int ish, int jsh,
                                            Ref<GaussianBasisSet> dercs,
                                            int centernum)
{
  int i;

  /* Get the basis function part of the nuclear derivative. */
  three_center = 1;
  third_centers = bs1_;
  for (i=0; i<bs1_->ncenter(); i++) {
    third_centernum = i;
    for (int xyz=0; xyz<3; xyz++) {
        C[xyz] = bs1_->r(i,xyz);
      }
    scale_shell_result = 1;
    result_scale_factor = -bs1_->molecule()->charge(i);
    //accum_shell_1der(buff,ish,jsh,dercs,centernum,
    //                 &Int1eV3::comp_shell_nuclear);
    accum_shell_block_1der(buff,ish,jsh,dercs,centernum,
                           &Int1eV3::comp_shell_block_nuclear);
    scale_shell_result = 0;
    }
  if (bs2_!=bs1_) {
    third_centers = bs2_;
    for (i=0; i<bs2_->ncenter(); i++) {
      third_centernum = i;
      for (int xyz=0; xyz<3; xyz++) {
          C[xyz] = bs2_->r(i,xyz);
        }
      scale_shell_result = 1;
      result_scale_factor = -bs2_->molecule()->charge(i);
      //accum_shell_1der(buff,ish,jsh,dercs,centernum,
      //                 &Int1eV3::comp_shell_nuclear);
      accum_shell_block_1der(buff,ish,jsh,dercs,centernum,
                             &Int1eV3::comp_shell_block_nuclear);
      scale_shell_result = 0;
      }
    }
  three_center = 0;

  }

/* This computes the efield integrals between functions in two shells.
 * The result is accumulated in the buffer in the form bf1 x y z, bf2
 * x y z, etc.
 */
void
Int1eV3::int_accum_shell_efield(int ish, int jsh,
                                double *position)
{
  scale_shell_result = 0;
  for (int xyz=0; xyz<3; xyz++) {
      C[xyz] = position[xyz];
    }
  accum_shell_efield(buff,ish,jsh);
}

/* This computes the efield integrals between functions in two shells.
 * The result is accumulated in the buffer in the form bf1 x y z, bf2
 * x y z, etc.  The globals scale_shell_result, result_scale_factor,
 * and C must be set before this is called.
 */
void
Int1eV3::accum_shell_efield(double *buff, int ish, int jsh)
{
  int i;
  int c1,i1,j1,k1,c2,i2,j2,k2;
  double efield[3];
  int gc1,gc2;
  int index1,index2;
  double *tmp = cartesianbuffer;

  if (!(init_order >= 1)) {
    ExEnv::errn() << scprintf("accum_shell_efield: 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);

  FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
      comp_shell_efield(efield,gc1,i1,j1,k1,gc2,i2,j2,k2);
      if (scale_shell_result) {
        for (i=0; i<3; i++) efield[i] *= result_scale_factor;
        }
      for (i=0; i<3; i++) tmp[i] = efield[i];
      tmp += 3;
      END_FOR_GCCART_GS(index2)
    END_FOR_GCCART_GS(index1)

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

/* This computes the efield integrals between functions in two shells.
 * The result is accumulated in the buffer in the form bf1 x y z, bf2
 * x y z, etc.  The globals scale_shell_result, result_scale_factor,
 * and C must be set before this is called.
 */
void
Int1eV3::accum_shell_block_efield(double *buff, int ish, int jsh)
{
  int c1,c2;
  int gc1,gc2;

  if (!(init_order >= 1)) {
    ExEnv::errn() << scprintf("accum_shell_efield: 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);

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

  int gcsize1 = gshell1->ncartesian();
  int gcsize2 = gshell2->ncartesian();
  memset(cartesianbuffer,0,sizeof(double)*gcsize1*gcsize2*3);
  int gcoff1 = 0;
  for (gc1=0; gc1<gshell1->ncontraction(); gc1++) {
    int a = gshell1->am(gc1);
    int sizea = INT_NCART_NN(a);
    int gcoff2 = 0;
    for (gc2=0; gc2<gshell2->ncontraction(); gc2++) {
      int b = gshell2->am(gc2);
      int sizeb = INT_NCART_NN(b);
      comp_shell_block_efield(gc1,a,gc2,b,
                              gcsize2, gcoff1, gcoff2,
                              coef, cartesianbuffer);
      gcoff2 += sizeb;
      }
    gcoff1 += sizea;
    }

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

/* This computes the efield integrals between functions in two shells.
 * The result is placed in the buffer in the form bf1 x y z, bf2
 * x y z, etc.
 */
void
Int1eV3::efield(int ish, int jsh, double *position)
{
  scale_shell_result = 0;
  int xyz;

  for (xyz=0; xyz<3; xyz++) {
      C[xyz] = position[xyz];
    }

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

  if (!(init_order >= 1)) {
    ExEnv::errn() << scprintf("Int1eV3::efield one electron routines are not ready\n");
    exit(1);
    }

  c1 = bs1_->shell_to_center(ish);
  c2 = bs2_->shell_to_center(jsh);
  for (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)
      comp_shell_efield(efield,gc1,i1,j1,k1,gc2,i2,j2,k2);
      if (scale_shell_result) {
        for (i=0; i<3; i++) efield[i] *= result_scale_factor;
        }
      for (i=0; i<3; i++) tmp[i] = efield[i];
      tmp += 3;
      END_FOR_GCCART_GS(index2)
    END_FOR_GCCART_GS(index1)

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

/* This slowly computes the nuc rep energy integrals between functions in
 * two shells.  The result is placed in the buffer.  */
void
Int1eV3::nuclear_slow(int ish, int jsh)
{
  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("int_shell_nuclear: 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 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 nuc rep energy integrals between functions in two
 * shells.  The result is placed in the buffer.  */
void
Int1eV3::nuclear(int ish, int jsh)
{
  int i;
  int c1,c2;
  int gc1,gc2;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_nuclear: 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);

  int ni = gshell1->ncartesian();
  int nj = gshell2->ncartesian();
  memset(cartesianbuffer,0,sizeof(double)*ni*nj);

  int offi = 0;
  for (gc1=0; gc1<gshell1->ncontraction(); gc1++) {
    int a = gshell1->am(gc1);
    int offj = 0;
    for (gc2=0; gc2<gshell2->ncontraction(); gc2++) {
      int b = gshell2->am(gc2);
      /* Loop thru the centers on bs1_. */
      for (i=0; i<bs1_->ncenter(); i++) {
        double charge = bs1_->molecule()->charge(i);
        for (int xyz=0; xyz<3; xyz++) C[xyz] = bs1_->r(i,xyz);
        comp_shell_block_nuclear(gc1, a, gc2, b,
                                 nj, offi, offj,
                                 -charge, cartesianbuffer);
        }
      /* Loop thru the centers on bs2_ if necessary. */
      if (bs2_ != bs1_) {
        for (i=0; i<bs2_->ncenter(); i++) {
          double charge = bs2_->molecule()->charge(i);
          for (int xyz=0; xyz<3; xyz++) C[xyz] = bs2_->r(i,xyz);
          comp_shell_block_nuclear(gc1, a, gc2, b,
                                   nj, offi, offj,
                                   -charge, cartesianbuffer);
          }
        }
      offj += INT_NCART_NN(b);
      }
    offi += INT_NCART_NN(a);
    }

#if DEBUG_NUC_SHELL
  double *fastbuf = new double[ni*nj];
  memcpy(fastbuf,cartesianbuffer,sizeof(double)*ni*nj);
  nuclear_slow(ish,jsh);

  index = 0;
  int cart1, cart2;
  FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1) {
    FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2) {
      double fast = fastbuf[index];
      double slow = cartesianbuffer[index];
      if (fabs(fast-slow)>1.0e-12) {
        ExEnv::outn() << scprintf("NUC SHELL FINAL: %d (%d %d %d) %d (%d %d %d): ",
                         gc1, i1,j1,k1, gc2, i2,j2,k2)
             << scprintf(" % 20.15f % 20.15f",
                         fast, slow)
             << endl;
        }
      index++;
      } END_FOR_GCCART_GS(cart2);
    } END_FOR_GCCART_GS(cart1);

  memcpy(cartesianbuffer,fastbuf,sizeof(double)*ni*nj);
  delete[] fastbuf;
#endif

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

/* This computes the integrals between functions in two shells for
 * a point charge interaction operator.