comp2e.cc

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

      if (p12) {
        tmp=stride1; stride1=stride2; stride2=tmp;
        }
      if (p34) {
        tmp=stride3; stride3=stride4; stride4=tmp;
        }
      if (p13p24) {
        tmp=stride1; stride1=stride3; stride3=tmp;
        tmp=stride2; stride2=stride4; stride4=tmp;
        }
      int newindex1 = newindexoffset;
      for (int ci1=0; ci1<tsize1; ci1++) {
        int newindex12 = newindex1;
        for (int ci2=0; ci2<tsize2; ci2++) {
          int newindex123 = newindex12;
          for (int ci3=0; ci3<tsize3; ci3++) {
            double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
            int newindex1234 = newindex123;
            for (int ci4=0; ci4<tsize4; ci4++) {
              int_buffer[newindex1234] = tmp_shiftbuffer[ci4];
              newindex1234 += stride4;
              }
            redundant_index+=tsize4;
            newindex123 += stride3;
            }
          newindex12 += stride2;
          }
        newindex1 += stride1;
        }
      }
    else if (nc3 == 1 && nc4 == 1) {
      // this is the p12 only case w/o gen contractions on 3 & 4
      // this special case collapses the 3rd and 4th indices together
      for (int ci1=0; ci1<tsize1; ci1++) {
        for (int ci2=0; ci2<tsize2; ci2++) {
          int pci1=ci1;
          int pci2=ci2;
          if (p12) {
            int tmp=pci1; pci1=pci2; pci2=tmp;
            }
          int newindex123 = newindexoffset + pci1*psize234 + pci2*psize34;
          double *tmp_int_buffer = &int_buffer[newindex123];
          double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
          for (int ci34=0; ci34<psize34; ci34++)
            tmp_int_buffer[ci34] = tmp_shiftbuffer[ci34];
          redundant_index += psize34;
          }
        }
      }
    else {
      // this is the p12 only case w/gen. contr. on 3 & 4
      for (int ci1=0; ci1<tsize1; ci1++) {
        for (int ci2=0; ci2<tsize2; ci2++) {
          int pci1=ci1;
          int pci2=ci2;
          if (p12) {
            int tmp=pci1; pci1=pci2; pci2=tmp;
            }
          int newindex123 = newindexoffset + pci1*psize234 + pci2*psize34;
          for (int ci3=0; ci3<tsize3; ci3++) {
            double *tmp_int_buffer = &int_buffer[newindex123];
            double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
            for (int ci4=0; ci4<tsize4; ci4++) {
              tmp_int_buffer[ci4] = tmp_shiftbuffer[ci4];
              }
            redundant_index += tsize4;
            newindex123 += psize4;
            }
          }
        }
      }
    }
  else if (nc3 == 1 && nc4 == 1) {
    // this special case collapses the 3rd and 4th indices together
    int size34 =  size3 * size4;
    int size234 = size2 * size34;
    double* redund_ints = shiftbuffer;
    redundant_index = 0;
    int newindex1 = ogc1*size234 + ogc2*size34 + ogc3*size4 + ogc4;
    for (int ci1=0; ci1<tsize1; ci1++) {
      int newindex12 = newindex1;
      for (int ci2=0; ci2<tsize2; ci2++) {
        double *tmp_int_buffer = &int_buffer[newindex12];
        double *tmp_redund_ints = &redund_ints[redundant_index];
        for (int ci34=0; ci34<size34; ci34++)
          tmp_int_buffer[ci34] = tmp_redund_ints[ci34];
        redundant_index += size34;
        newindex12 += size34;
        }
      newindex1 += size234;
      }
    }
  else {
    int size34 =  size3 * size4;
    int size234 = size2 * size34;
    double* redund_ints = shiftbuffer;
    redundant_index = 0;
    int newindex1 = ogc1*size234 + ogc2*size34 + ogc3*size4 + ogc4;
    for (int ci1=0; ci1<tsize1; ci1++) {
      int newindex12 = newindex1;
      for (int ci2=0; ci2<tsize2; ci2++) {
        int newindex123 = newindex12;
        for (int ci3=0; ci3<tsize3; ci3++) {
          double *tmp_int_buffer = &int_buffer[newindex123];
          double *tmp_redund_ints = &redund_ints[redundant_index];
          for (int ci4=0; ci4<tsize4; ci4++) {
            tmp_int_buffer[ci4] = tmp_redund_ints[ci4];
            }
          redundant_index += tsize4;
          newindex123 += size4;
          }
        newindex12 += size34;
        }
      newindex1 += size234;
      }
    }

    /* End loop over generalized contractions. */
          ogc4 += tsize4;
          }
        ogc3 += tsize3;
        }
      ogc2 += tsize2;
      }
    ogc1 += tsize1;
    }

  if (   !int_unit2
      && !int_unit4
      && int_integral_storage) {
#ifdef EREP_TIMING
      tim_change("maybe store");
#endif
      int_store_integral(sh1,sh2,sh3,sh4,p12,p34,p13p24,size);
    }

  /* We branch here if an integral was precomputed and the int_buffer
   * has been already filled. */
  post_computation:

#ifdef EREP_TIMING
  tim_change("post");
#endif

  /* Unpermute all of the permuted quantities. */
  if ((!permute_)&&(p12||p34||p13p24)) {
    if (p13p24) {
      iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
      iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
      iswtch(&int_unit2,&int_unit4);
      iswtch(&am1,&am3);
      iswtch(&am2,&am4);
      iswtch(&am12,&am34);
      sswtch(&int_shell1,&int_shell3);
      swtch(pbs1,pbs3);
      sswtch(&int_shell2,&int_shell4);
      swtch(pbs2,pbs4);
      iswtch(&int_expweight1,&int_expweight3);
      iswtch(&int_expweight2,&int_expweight4);
      }
    if (p34) {
      iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
      iswtch(&am3,&am4);
      sswtch(&int_shell3,&int_shell4);
      swtch(pbs3,pbs4);
      iswtch(&int_expweight3,&int_expweight4);
      }
    if (p12) {
      iswtch(&sh1,&sh2);iswtch(psh1,psh2);iswtch(&osh1,&osh2);
      iswtch(&am1,&am2);
      sswtch(&int_shell1,&int_shell2);
      swtch(pbs1,pbs2);
      iswtch(&int_expweight1,&int_expweight2);
      }
    }

  /* Transform to pure am (if requested in the centers structure). */
  if (!(flags&INT_NOPURE)) {
      transform_2e(integral_, int_buffer, int_buffer,
                   &bs1_->shell(sh1),
                   int_unit2?int_unit_shell:&bs2_->shell(sh2),
                   &bs3_->shell(sh3),
                   int_unit4?int_unit_shell:&bs4_->shell(sh4));
    }

  /* Remove the redundant integrals, unless redundant_ is set. */
  if (!redundant_) {
    int redundant_offset = 0;
    int nonredundant_offset = 0;
    if ((osh1 == osh4)&&(osh2 == osh3)&&(osh1 != osh2)) {
      ExEnv::errn() << scprintf("nonredundant integrals cannot be generated\n");
      fail();
      }
    e12 = (int_unit2?0:(osh1 == osh2));
    e13e24 = ((osh1 == osh3)
              && ((int_unit2 && int_unit4)
                  || ((int_unit2||int_unit4)?0:(osh2 == osh4))));
    e34 = (int_unit4?0:(osh3 == osh4));
    if (e12||e34||e13e24) {
      nonredundant_erep(int_buffer,e12,e34,e13e24,
                        int_shell1->nfunction(),
                        int_shell2->nfunction(),
                        int_shell3->nfunction(),
                        int_shell4->nfunction(),
                        &redundant_offset,
                        &nonredundant_offset);
      }
    }
    
#ifdef EREP_TIMING
  tim_exit("post");
  tim_exit(section);
#endif
  }

/* This computes the two electron derivatives for all unique
 * centers in the passed shell quartet.  One center in
 * the set of unique centers is not included.  This can
 * be computed as minus the sum of the other derivatives.
 * The list of centers for which integrals were computed can
 * be determined from the contents of der_centers.
 * The results are put into the global integral buffer in the
 * format:
 * +------------------+
 * | dercenter1 +---+ |
 * |            | x | |
 * |            +---+ |
 * |            | y | |
 * |            +---+ |
 * |            | z | |
 * |            +---+ |
 * +------------------+
 * | dercenter2 +---+ |
 * |            | x | |
 * |            +---+ |
 * |            | y | |
 * |            +---+ |
 * |            | z | |
 * |            +---+ |
 * +------------------+
 * | dercenter3 +---+ |
 * |            | x | |
 * |            +---+ |
 * |            | y | |
 * |            +---+ |
 * |            | z | |
 * |            +---+ |
 * +------------------+
 */

void
Int2eV3::erep_all1der(int &psh1, int &psh2, int &psh3, int &psh4,
                      der_centersv3_t *der_centers)
{
  double *current_buffer;
  int nints;
  double *user_int_buffer;
  int omit;
  GaussianBasisSet *cs[4];
  int sh[4];
  int n_unique;
  int i,j;
  GaussianShell *shell1,*shell2,*shell3,*shell4;
  GaussianBasisSet *ucs[4]; /* The centers struct for the unique centers. */
  int unum[4];        /* The number of times that this unique center occurs. */
  int uam[4];         /* The total angular momentum on each unique center. */
  int am[4];
  int osh[4];
  int cen[4];
  int ucen[4];
  int ncart;
  double *current_pure_buffer;

  cs[0] = bs1_.pointer();
  cs[1] = bs2_.pointer();
  cs[2] = bs3_.pointer();
  cs[3] = bs4_.pointer();

  sh[0] = psh1;
  sh[1] = psh2;
  sh[2] = psh3;
  sh[3] = psh4;

  /* Set up pointers to the current shells. */
  shell1 = &bs1_->shell(psh1);
  shell2 = &bs2_->shell(psh2);
  shell3 = &bs3_->shell(psh3);
  shell4 = &bs4_->shell(psh4);

  /* Number of cartesian and pure integrals. */
  ncart = shell1->ncartesian()*shell2->ncartesian()
         *shell3->ncartesian()*shell4->ncartesian();
  nints = shell1->nfunction()*shell2->nfunction()
         *shell3->nfunction()*shell4->nfunction();

  am[0] = shell1->max_am();
  am[1] = shell2->max_am();
  am[2] = shell3->max_am();
  am[3] = shell4->max_am();

  /* Compute the offset shell numbers. */
  osh[0] = psh1 + bs1_shell_offset_;
  osh[1] = psh2 + bs2_shell_offset_;
  osh[2] = psh3 + bs3_shell_offset_;
  osh[3] = psh4 + bs4_shell_offset_;

  for (i=0; i<4; i++) cen[i] = cs[i]->shell_to_center(sh[i]);

  /* This macro returns true if two shell centers are the same. */
#define SC(cs1,shc1,cs2,shc2) (((cs1)==(cs2))&&((shc1)==(shc2)))

  /* Build the list of unique centers structures and shells. */
  n_unique = 0;
  for (i=0; i<4; i++) {
    int unique = 1;
    for (j=0; j<n_unique; j++) {
      if (SC(ucs[j],ucen[j],cs[i],cen[i])) {
        unique = 0;
        uam[j] += am[i];
        unum[j]++;
        break;
        }
      }
    if (unique) {
      ucs[n_unique] = cs[i];
      ucen[n_unique] = cen[i];
      uam[n_unique] = am[i];
      unum[n_unique] = 1;
      n_unique++;
      }
    }

  /* Choose which unique center is to be omitted from the calculation. */
  if (n_unique == 1) {
    omit = 0;
    }
  else if (n_unique == 2) {
    if (unum[0]==3) omit = 0;
    else if (unum[1]==3) omit = 1;
    else if (uam[1]>uam[0]) omit = 1;
    else omit = 0;
    }
  else if (n_unique == 3) {
    if (unum[0]==2) omit = 0;
    else if (unum[1]==2) omit = 1;
    else omit = 2;
    }
  else {
    int max = 0;
    omit = 0;
    for (i=0; i<4; i++) {
      if (uam[i]>max) {
        max = uam[i];
        omit = i;
        }
      }
    }

  /* Save the location of the int_buffer. */
  user_int_buffer = int_buffer;
  int_buffer = int_derint_buffer;

  /* Zero out the result integrals. */
  for (i=0; i<3*(n_unique-1)*ncart; i++) user_int_buffer[i] = 0.0;

  /* Loop thru the unique centers, computing the integrals and
   * skip the derivative on the unique center specified by omit. */
  der_centers->n = 0;
  current_buffer = user_int_buffer;
  for (i=0; i<n_unique; i++) {
    if (i==omit) continue;

    der_centers->cs[der_centers->n] = ucs[i];
    der_centers->num[der_centers->n] = ucen[i];
    der_centers->n++;

    for (j=0; j<4; j++) {
      if (SC(ucs[i],ucen[i],cs[j],cen[j])) {
        int old_perm = permute();
        set_permute(0);
        compute_erep_1der(0,current_buffer,
                          &psh1,&psh2,&psh3,&psh4,j);
        set_permute(old_perm);
        }
      }

    current_buffer = &current_buffer[3*ncart];
    }

  /* Put the information about the omitted center into der_centers. */
  der_centers->ocs = ucs[omit];
  der_centers->onum = ucen[omit];

  /* Transform to pure am. */
  current_buffer = user_int_buffer;
  current_pure_buffer = user_int_buffer;
  for (i=0; i<3*der_centers->n; i++) {
      transform_2e(integral_, current_buffer, current_pure_buffer,
                   shell1, shell2, shell3, shell4);
      current_buffer = &current_buffer[ncart];
      current_pure_buffer = &current_pure_buffer[nints];
    }

  /* Eliminate redundant integrals, unless flags specifies otherwise. */
  current_buffer = user_int_buffer;
  if (!redundant_) {
    int redundant_offset = 0;
    int nonredundant_offset = 0;
    int e12,e13e24,e34;
    int i;

    if ((osh[0] == osh[3])&&(osh[1] == osh[2])&&(osh[0] != osh[1])) {
      ExEnv::errn() << scprintf("nonredundant integrals cannot be generated (1der)\n");
      fail();
      }

    /* Shell equivalence information. */
    e12 = (osh[0] == osh[1]);
    e13e24 = ((osh[0] == osh[2]) && (osh[1] == osh[3]));
    e34 = (osh[2] == osh[3]);
    if (e12||e13e24||e34) {
      /* Repack x, y, and z integrals. */
      for (i=0; i<3*der_centers->n; i++) {