comp2e.cc

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

        for (k4=0; k4<=am4-i4; k4++) {
          int c4yp1 = c4+i4;//=INT_CARTINDEX(am4+1,i4,j4+1)
          int buffer_index = size1234 + ogc4+c4;
          int int_buffer_index = ogc4p+c4yp1;
          for (f123=0; f123<size123; f123++) {
            buffer[buffer_index] += int_buffer[int_buffer_index];
            buffer_index += size4;
            int_buffer_index += sizep4;
            }
          c4++;
          }
        }
      ogc4 += c4;
      ogc4p += INT_NCART(am4+1);
      }
    index += size1234;

    /* The center 3 d/dz integrals */
    ogc4 = 0;
    ogc4p = 0;
    for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
      int am4 = shell4->am(gc4);
      int c4 = 0;
      for (i4=0; i4<=am4; i4++) {
        for (k4=0; k4<=am4-i4; k4++) {
          int c4zp1 = c4+i4+1;//=INT_CARTINDEX(am4+1,i4,j4)
          int buffer_index = size1234+size1234+ogc4+c4;
          int int_buffer_index = ogc4p+c4zp1;
          for (f123=0; f123<size123; f123++) {
            buffer[buffer_index] += int_buffer[int_buffer_index];
            buffer_index += size4;
            int_buffer_index += sizep4;
            }
          c4++;
          }
        }
      ogc4 += c4;
      ogc4p += INT_NCART(am4+1);
      }
    index += size1234;
    }
#ifdef DER_TIMING
  tim_exit(0);
  tim_exit(0);
#endif
  }

void
Int2eV3::nonredundant_erep(double *buffer, int e12, int e34, int e13e24,
                           int n1, int n2, int n3, int n4,
                           int *red_off, int *nonred_off)
{
  int nonredundant_index;
  int i,j,k,l;

  double *redundant_ptr = &buffer[*red_off];
  double *nonredundant_ptr = &buffer[*nonred_off];

  nonredundant_index = 0;
  int n34 = n3*n4;
  for (i=0; i<n1; i++) {
    int jmax = INT_MAX2(e12,i,n2);
    for (j=0; j<=jmax; j++) {
      int kmax = INT_MAX3(e13e24,i,n3);
      for (k=0; k<=kmax; k++) {
        int lmax = INT_MAX4(e13e24,e34,i,j,k,n4);
        for (l=0; l<=lmax; l++) {
          nonredundant_ptr[l] = redundant_ptr[l];
          }
        redundant_ptr += n4;
        nonredundant_index += lmax+1;
        nonredundant_ptr += lmax+1;
        }
      redundant_ptr += (n3-(kmax+1))*n4;
      }
    redundant_ptr += (n2-(jmax+1))*n34;
    }
  *red_off += n1*n2*n34;
  *nonred_off += nonredundant_index;
  }

/* This is an alternate interface to int_erep_all1der.  It takes
 * as arguments the flags, an integer vector of shell numbers
 * and an integer vector which will be filled in with size
 * information, if it is non-NULL, and the dercenters pointer. */
void
Int2eV3::erep_all1der(int *shells, int  *sizes,
                      der_centersv3_t *dercenters)
{
  erep_all1der(shells[0],shells[1],shells[2],shells[3],
               dercenters);
  if (sizes) {
    sizes[0] = bs1_->shell(shells[0]).nfunction();
    sizes[1] = bs2_->shell(shells[1]).nfunction();
    sizes[2] = bs3_->shell(shells[2]).nfunction();
    sizes[3] = bs4_->shell(shells[3]).nfunction();
    }
  }

void
Int2eV3::int_erep_bound1der(int flags, int bsh1, int bsh2, int *size)
{
  double *current_buffer;
  int nints;
  double *user_int_buffer;
  int i;
  GaussianShell *shell1,*shell2,*shell3,*shell4;
  int osh[4];
  int sh1 = bsh1;
  int sh2 = bsh2;
  int sh3 = bsh1;
  int sh4 = bsh2;
  int *psh1 = &sh1;
  int *psh2 = &sh2;
  int *psh3 = &sh3;
  int *psh4 = &sh4;
  int ncart;
  double *current_pure_buffer;

  /* 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();

  /* 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_;

  /* 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*ncart; i++) user_int_buffer[i] = 0.0;

  /* Set the size so it is available to the caller. */
  *size = nints * 3;

  current_buffer = user_int_buffer;
  int old_perm = permute();
  compute_erep_bound1der(flags|INT_NOPURE,current_buffer,
                          psh1,psh2,psh3,psh4);
  set_permute(old_perm);

  /* Transform to pure am. */
  current_buffer = user_int_buffer;
  current_pure_buffer = user_int_buffer;
  for (i=0; i<3; i++) {
      transform_2e(integral_, current_buffer, current_pure_buffer,
                   &bs1_->shell(sh1),
                   &bs2_->shell(sh2),
                   &bs3_->shell(sh3),
                   &bs4_->shell(sh4));
      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]);
    /* Repack x, y, and z integrals. */
    for (i=0; i<3; i++) {
      nonredundant_erep(current_buffer,e12,e34,e13e24,
                             shell1->nfunction(),
                             shell2->nfunction(),
                             shell3->nfunction(),
                             shell4->nfunction(),
                             &redundant_offset,
                             &nonredundant_offset);
      }
    }

  /* Return the integral buffers to their normal state. */
  int_derint_buffer = int_buffer;
  int_buffer = user_int_buffer;
  }

/* This routine computes a quantity needed to compute the
 * derivative integral bounds.
 * It fills int_buffer with (sh1+i sh2|sh1+i sh2).
 */
void
Int2eV3::compute_erep_bound1der(int flags, double *buffer,
                                int *psh1, int *psh2, int *psh3, int *psh4)
{
  int oc1,oc2,oc3,oc4;
  int ii;
  int c1,c2,c3,c4;
  int i[4],j[4],k[4],am[4];
  int index;
  int sizem234,sizem34,sizem2,sizem3,sizem4;
  int sizep234,sizep34,sizep2,sizep3,sizep4;
  int sizepm234,sizepm34,sizepm2,sizepm3,sizepm4;
  GaussianShell *shell1,*shell2,*shell3,*shell4;

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

#define DCT1(n) ((((n)==0)||((n)==2))?1:0)
#define DCT2(n) ((((n)==0)||((n)==2))?((((n)==0))?1:-1):0)
  /* Offsets for the intermediates with angular momentum decremented. */
  for (ii=sizem2=0; ii<shell2->ncontraction(); ii++) 
    sizem2 += INT_NCART(shell2->am(ii)-DCT1(1));
  for (ii=sizem3=0; ii<shell3->ncontraction(); ii++) 
    sizem3 += INT_NCART(shell3->am(ii)-DCT1(2));
  for (ii=sizem4=0; ii<shell4->ncontraction(); ii++) 
    sizem4 += INT_NCART(shell4->am(ii)-DCT1(3));
  sizem34 = sizem4 * sizem3;
  sizem234 = sizem34 * sizem2;

  /* Offsets for the intermediates with angular momentum incremented. */
  for (ii=sizep2=0; ii<shell2->ncontraction(); ii++) 
    sizep2 += INT_NCART(shell2->am(ii)+DCT1(1));
  for (ii=sizep3=0; ii<shell3->ncontraction(); ii++) 
    sizep3 += INT_NCART(shell3->am(ii)+DCT1(2));
  for (ii=sizep4=0; ii<shell4->ncontraction(); ii++) 
    sizep4 += INT_NCART(shell4->am(ii)+DCT1(3));
  sizep34 = sizep4 * sizep3;
  sizep234 = sizep34 * sizep2;

  /* Offsets for the intermediates with angular momentum incremented and
   * decremented. */
  for (ii=sizepm2=0; ii<shell2->ncontraction(); ii++) 
    sizepm2 += INT_NCART(shell2->am(ii)+DCT2(1));
  for (ii=sizepm3=0; ii<shell3->ncontraction(); ii++) 
    sizepm3 += INT_NCART(shell3->am(ii)+DCT2(2));
  for (ii=sizepm4=0; ii<shell4->ncontraction(); ii++) 
    sizepm4 += INT_NCART(shell4->am(ii)+DCT2(3));
  sizepm34 = sizepm4 * sizepm3;
  sizepm234 = sizepm34 * sizepm2;

  /* END_DERLOOP must be redefined here because it previously depended
   * on the DCTEST macro */
#undef END_DERLOOP
#define END_DERLOOP(index,indexp1,sign) \
       END_FOR_CART\
     oc##indexp1 += INT_NCART(am[index] sign DCT1(index));\
     }

#undef END_ALLDERLOOPS
#define END_ALLDERLOOPS(sign)\
            END_DERLOOP(3,4,sign)\
          END_DERLOOP(2,3,sign)\
        END_DERLOOP(1,2,sign)\
      END_DERLOOP(0,1,sign)

  int old_perm = permute();
  set_permute(0);
  int old_red = redundant();
  set_redundant(1);
  compute_erep(flags,psh1,psh2,psh3,psh4,
                   -DCT1(0),
                   -DCT1(1),
                   -DCT1(2),
                   -DCT1(3));
  set_permute(old_perm);
  set_redundant(old_red);

   /* Place the contributions into the user integral buffer. */
   index = 0;
   /* The d/dx integrals */
  ALLDERLOOPS
    if (i[0]>0 && i[2]>0) {
      buffer[index] += i[0] * i[2] * int_buffer[
        (oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0]-DCT1(0),j[0])) * sizem234
       +(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1]-DCT1(1),j[1])) * sizem34
       +(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2]-DCT1(2),j[2])) * sizem4
       +(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3]-DCT1(3),j[3]))
       ];
      }
    index++;
    END_ALLDERLOOPS(-)

   /* The d/dy integrals */
  ALLDERLOOPS
    if (j[0]>0 && j[2]>0) {
    buffer[index] += j[0] * j[2] * int_buffer[
         (oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0],j[0]-DCT1(0))) * sizem234
        +(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1],j[1]-DCT1(1))) * sizem34
        +(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2],j[2]-DCT1(2))) * sizem4
        +(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3],j[3]-DCT1(3)))
        ];
      }
    index++;
  END_ALLDERLOOPS(-)

   /* The d/dz integrals */
  ALLDERLOOPS
    if (k[0]>0 && k[2]>0) {
    buffer[index] += k[0] * k[2] * int_buffer[
         (oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0],j[0])) * sizem234
        +(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1],j[1])) * sizem34
        +(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2],j[2])) * sizem4
        +(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3],j[3]))
        ];
      }
    index++;
  END_ALLDERLOOPS(-)

  /* Compute the next contribution to the integrals. */
  /* Tell the build routine that we need an exponent weighted contraction
   * with the exponents taken from the dercenter and adjust the
   * angular momentum of dercenter to the needed value. */
  int_expweight1 = 1;
  int_expweight3 = 1;
  old_perm = permute();
  set_permute(0);
  old_red = redundant();
  set_redundant(1);
  compute_erep(flags,psh1,psh2,psh3,psh4,
                     DCT1(0),
                     DCT1(1),
                     DCT1(2),
                     DCT1(3));
  set_permute(old_perm);
  set_redundant(old_red);
  int_expweight1 = 0;
  int_expweight3 = 0;

  /* Place the contributions into the user integral buffer. */
  index = 0;
  /* The d/dx integrals */
  ALLDERLOOPS
          buffer[index] += int_buffer[
             (oc1+INT_CARTINDEX(am[0]+DCT1(0),i[0]+DCT1(0),j[0]))*sizep234
            +(oc2+INT_CARTINDEX(am[1]+DCT1(1),i[1]+DCT1(1),j[1]))*sizep34
            +(oc3+INT_CARTINDEX(am[2]+DCT1(2),i[2]+DCT1(2),j[2]))*sizep4
            +(oc4+INT_CARTINDEX(am[3]+DCT1(3),i[3]+DCT1(3),j[3]))
            ];
    index++;
    END_ALLDERLOOPS(+)

  /* The d/dy integrals */
  ALLDERLOOPS
          buffer[index] += int_buffer[
             (oc1+INT_CARTINDEX(am[0]+DCT1(0),i[0],j[0]+DCT1(0)))*sizep234
            +(oc2+INT_CARTINDEX(am[1]+DCT1(1),i[1],j[1]+DCT1(1)))*sizep34
            +(oc3+INT_CARTINDEX(am[2]+DCT1(2),i[2],j[2]+DCT1(2)))*sizep4
            +(oc4+INT_CARTINDEX(am[3]+DCT1(3),i[3],j[3]+DCT1(3)))
            ];
          index++;
    END_ALLDERLOOPS(+)

  /* The d/dz integrals */
  ALLDERLOOPS
          buffer[index] += int_buffer[
               (oc1 + INT_CARTINDEX(am[0]+DCT1(0),i[0],j[0])) * sizep234
              +(oc2 + INT_CARTINDEX(am[1]+DCT1(1),i[1],j[1])) * sizep34
              +(oc3 + INT_CARTINDEX(am[2]+DCT1(2),i[2],j[2])) * sizep4
              +(oc4 + INT_CARTINDEX(am[3]+DCT1(3),i[3],j[3]))
              ];
          index++;
    END_ALLDERLOOPS(+)

  /* END_DERLOOP must be redefined here because it previously depended
   * on the DCT1 macro */
#undef END_DERLOOP
#define END_DERLOOP(index,indexp1,sign) \
       END_FOR_CART\
     oc##indexp1 += INT_NCART(am[index] sign DCT2(index));\
     }

#undef END_ALLDERLOOPS
#define END_ALLDERLOOPS(sign)\
            END_DERLOOP(3,4,sign)\
          END_DERLOOP(2,3,sign)\
        END_DERLOOP(1,2,sign)\
      END_DERLOOP(0,1,sign)

  /* Compute the next contribution to the integrals. */
  /* Tell the build routine that we need an exponent weighted contraction
   * with the exponents taken from the dercenter and adjust the
   * angular momentum of dercenter to the needed value. */
  int_expweight1 = 1;
  old_perm = permute();
  set_permute(0);
  old_red = redundant();
  set_redundant(1);
  compute_erep(flags,psh1,psh2,psh3,psh4,
                     DCT2(0),
                     DCT2(1),
                     DCT2(2),
                     DCT2(3));
  set_permute(old_perm);
  set_redundant(old_red);
  int_expweight1 = 0;

  /* Place the contributions into the user integral buffer. */
  index = 0;
  /* The d/dx integrals */
  ALLDERLOOPS
    if (i[2] > 0)  {
          buffer[