comp2e.cc

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

        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 will call compute_erep
 * to compute the derivatives in terms of order == 0 integrals.
 * flags are the flags to the int_comp_erep routine
 * psh1-4 are pointers to the shell numbers
 * dercenter is 0, 1, 2, or 3 -- the center within the integral
 *           which we are taking the derivative with respect to.
 * The results are accumulated in buffer, which cannot be the same
 * as the current int_buffer.
 */
void
Int2eV3::compute_erep_1der(int flags, double *buffer,
                           int *psh1, int *psh2, int *psh3, int *psh4,
                           int dercenter)
{
  int ii;
  int index;
  int size1,size2,size3,size4,size1234;
  int sizem234,sizem34,sizem2,sizem3,sizem4;
  int sizep234,sizep34,sizep2,sizep3,sizep4;
  GaussianShell *shell1,*shell2,*shell3,*shell4;

#ifdef DER_TIMING
  tim_enter("erep_1der");
#endif

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

  if ((dercenter<0) || (dercenter > 3)) {
    ExEnv::errn() << scprintf("illegal derivative center -- must be 0, 1, 2, or 3\n");
    fail();
    }

  /* Offsets for the intermediates with original angular momentum. */
  for (ii=size1=0; ii<shell1->ncontraction(); ii++)
    size1 += INT_NCART(shell1->am(ii));
  for (ii=size2=0; ii<shell2->ncontraction(); ii++)
    size2 += INT_NCART(shell2->am(ii));
  for (ii=size3=0; ii<shell3->ncontraction(); ii++)
    size3 += INT_NCART(shell3->am(ii));
  for (ii=size4=0; ii<shell4->ncontraction(); ii++)
    size4 += INT_NCART(shell4->am(ii));

  size1234 = size1*size2*size3*size4;

#define DCTEST(n) ((dercenter==n)?1:0)
  /* Offsets for the intermediates with angular momentum decremented. */
  for (ii=sizem2=0; ii<shell2->ncontraction(); ii++) 
    sizem2 += INT_NCART(shell2->am(ii)-DCTEST(1));
  for (ii=sizem3=0; ii<shell3->ncontraction(); ii++) 
    sizem3 += INT_NCART(shell3->am(ii)-DCTEST(2));
  for (ii=sizem4=0; ii<shell4->ncontraction(); ii++) 
    sizem4 += INT_NCART(shell4->am(ii)-DCTEST(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)+DCTEST(1));
  for (ii=sizep3=0; ii<shell3->ncontraction(); ii++) 
    sizep3 += INT_NCART(shell3->am(ii)+DCTEST(2));
  for (ii=sizep4=0; ii<shell4->ncontraction(); ii++) 
    sizep4 += INT_NCART(shell4->am(ii)+DCTEST(3));
  sizep34 = sizep4 * sizep3;
  sizep234 = sizep34 * sizep2;

#ifdef DER_TIMING
  tim_enter("- erep");
#endif

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

  /* Trouble if cpp is nonANSI. */
#define DERLOOP(index,indexp1) \
   oc##indexp1 = 0;\
   for ( c##indexp1 =0; c##indexp1 <shell##indexp1->ncontraction(); c##indexp1 ++) {\
     am[index] = shell##indexp1->am(c##indexp1);\
     FOR_CART(i[index],j[index],k[index],am[index])

#define END_DERLOOP(index,indexp1,sign) \
       END_FOR_CART\
     oc##indexp1 += INT_NCART(am[index] sign DCTEST(index));\
     }

#define ALLDERLOOPS\
    DERLOOP(0,1)\
      DERLOOP(1,2)\
        DERLOOP(2,3)\
          DERLOOP(3,4)

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

  /* Place the contributions into the user integral buffer. */
  index = 0;

  if (dercenter==0) {
    int ogc1,ogc1m,gc1,i1,k1,f234,size234;
    size234=size2*size3*size4;

#ifdef DER_TIMING
  tim_change("- 0");
#endif
    /* The center 0 d/dx integrals */
    ogc1 = 0;
    ogc1m = 0;
    for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
      int am1 = shell1->am(gc1);
      // only integrals with x^n n>0 on center 0 contribute
      // so skip over n==0
      index += (am1+1)*size234;
      int c1 = am1+1;
      for (i1=1; i1<=am1; i1++) {
        double factor = -i1;
        for (k1=0; k1<=am1-i1; k1++) {
          int c1xm1 = c1-am1-1;//=INT_CARTINDEX(am1-1,i1-1,j1)
          double *tmp_buffer=&buffer[index];
          double *tmp_int_buffer=&int_buffer[(ogc1m+c1xm1)*size234];
          for (f234=0; f234<size234; f234++) {
            tmp_buffer[f234] += factor * tmp_int_buffer[f234];
            }
          index+=size234;
          c1++;
          }
        }
      ogc1 += c1;
      ogc1m += INT_NCART(am1-1);
      }

    /* The center 0 d/dy integrals */
    ogc1 = 0;
    ogc1m = 0;
    for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
      int am1 = shell1->am(gc1);
      // only integrals with y^n n>0 on center 0 contribute
      // so skip over n==0 (by making i1+k1<am1)
      int c1 = 0;
      for (i1=0; i1<=am1; i1++) {
        for (k1=0; k1<=am1-i1-1; k1++) {
          double factor = -(am1-i1-k1);
          int c1ym1 = c1-i1;//=INT_CARTINDEX(am1-1,i1,j1-1)
          double *tmp_buffer=&buffer[index];
          double *tmp_int_buffer=&int_buffer[(ogc1m+c1ym1)*size234];
          for (f234=0; f234<size234; f234++) {
            tmp_buffer[f234] += factor * tmp_int_buffer[f234];
            }
          index+=size234;
          c1++;
          }
        // account for the y^n n==0 case by an extra increment
        c1++;
        index+=size234;
        }
      ogc1 += c1;
      ogc1m += INT_NCART(am1-1);
      }

    /* The center 0 d/dz integrals */
    ogc1 = 0;
    ogc1m = 0;
    for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
      int am1 = shell1->am(gc1);
      int c1 = 0;
      for (i1=0; i1<=am1; i1++) {
        // only integrals with z^n n>0 on center 0 contribute
        // so skip over n==0
        c1++;
        index+=size234;
        for (k1=1; k1<=am1-i1; k1++) {
          double factor = -k1;
          int c1zm1 = c1-i1-1;//=INT_CARTINDEX(am1-1,i1,j1)
          double *tmp_buffer=&buffer[index];
          double *tmp_int_buffer=&int_buffer[(ogc1m+c1zm1)*size234];
          for (f234=0; f234<size234; f234++) {
            tmp_buffer[f234] += factor * tmp_int_buffer[f234];
            }
          index+=size234;
          c1++;
          }
        }
      ogc1 += c1;
      ogc1m += INT_NCART(am1-1);
      }
    }
  else if (dercenter == 1) {
    int ogc2,ogc2m,gc2,i2,k2,f1,f34,size34,size234;
    size34 = size3*size4;
    size234 = size2*size3*size4;

#ifdef DER_TIMING
  tim_change("- 1");
#endif
    /* The center 1 d/dx integrals */
    ogc2 = 0;
    ogc2m = 0;
    for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
      int am2 = shell2->am(gc2);
      // only integrals with x^n n>0 on center 1 contribute
      // so skip over n==0
      int c2 = am2+1;
      for (i2=1; i2<=am2; i2++) {
        double factor = -i2;
        for (k2=0; k2<=am2-i2; k2++) {
          int c2xm1 = c2-am2-1;//=INT_CARTINDEX(am2-1,i2-1,j2)
          int buffer_index = (ogc2+c2)*size34;
          int int_buffer_index = (ogc2m+c2xm1)*size34;
          for (f1=0; f1<size1; f1++) {
            double *tmp_buffer=&buffer[buffer_index];
            double *tmp_int_buffer=&int_buffer[int_buffer_index];
            for (f34=0; f34<size34; f34++) {
                tmp_buffer[f34] += factor * tmp_int_buffer[f34];
              }
            buffer_index += size234;
            int_buffer_index += sizem234;
            }
          c2++;
          }
        }
      ogc2 += c2;
      ogc2m += INT_NCART(am2-1);
      }
    index += size1234;

     /* The center 1 d/dy integrals */
    ogc2 = 0;
    ogc2m = 0;
    for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
      int am2 = shell2->am(gc2);
      // only integrals with y^n n>0 on center 1 contribute
      // so skip over n==0
      int c2 = 0;
      for (i2=0; i2<=am2; i2++) {
        for (k2=0; k2<=am2-i2-1; k2++) {
          double factor = -(am2-k2-i2);
          int c2ym1 = c2-i2;//=INT_CARTINDEX(am2-1,i2,j2-1)
          int buffer_index = size1234 + (ogc2+c2)*size34;
          int int_buffer_index = (ogc2m+c2ym1)*size34;
          for (f1=0; f1<size1; f1++) {
            double *tmp_buffer=&buffer[buffer_index];
            double *tmp_int_buffer=&int_buffer[int_buffer_index];
            for (f34=0; f34<size34; f34++) {
                tmp_buffer[f34] += factor * tmp_int_buffer[f34];
              }
            buffer_index += size234;
            int_buffer_index += sizem234;
            }
          c2++;
          }
        // account for the y^n n==0 case by an extra increment
        c2++;
        }
      ogc2 += c2;
      ogc2m += INT_NCART(am2-1);
      }
    index += size1234;

     /* The center 1 d/dz integrals */
    ogc2 = 0;
    ogc2m = 0;
    for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
      int am2 = shell2->am(gc2);
      // only integrals with z^n n>0 on center 1 contribute
      // so skip over n==0
      int c2 = 0;
      for (i2=0; i2<=am2; i2++) {
        // account for the z^n n==0 case by an extra increment
        c2++;
        for (k2=1; k2<=am2-i2; k2++) {
          double factor = -k2;
          int c2zm1 = c2-i2-1;//=INT_CARTINDEX(am2-1,i2,j2-1)
          int buffer_index = size1234+size1234+(ogc2+c2)*size34;
          int int_buffer_index = (ogc2m+c2zm1)*size34;
          for (f1=0; f1<size1; f1++) {
            double *tmp_buffer=&buffer[buffer_index];
            double *tmp_int_buffer=&int_buffer[int_buffer_index];
            for (f34=0; f34<size34; f34++) {
                tmp_buffer[f34] += factor * tmp_int_buffer[f34];
              }
            buffer_index += size234;
            int_buffer_index += sizem234;
            }
          c2++;
          }
        }
      ogc2 += c2;
      ogc2m += INT_NCART(am2-1);
      }
    index += size1234;
    }
  else if (dercenter == 2) {
    int ogc3,ogc3m,gc3,i3,k3,f12,f4,size12,size34;
    size12 = size1*size2;
    size34 = size3*size4;

#ifdef DER_TIMING
  tim_change("- 2");
#endif
    /* The center 2 d/dx integrals */
    ogc3 = 0;
    ogc3m = 0;
    for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
      int am3 = shell3->am(gc3);
      // only integrals with x^n n>0 on center 2 contribute
      // so skip over n==0
      int c3 = am3+1;
      for (i3=1; i3<=am3; i3++) {
        double factor = -i3;
        for (k3=0; k3<=am3-i3; k3++) {
          int c3xm1 = c3-am3-1;//=INT_CARTINDEX(am3-1,i3-1,j3)
          int buffer_index = (ogc3+c3)*size4;
          int int_buffer_index = (ogc3m+c3xm1)*size4;
          for (f12=0; f12<size12; f12++) {
            double *tmp_buffer=&buffer[buffer_index];
            double *tmp_int_buffer=&int_buffer[int_buffer_index];
            for (f4=0; f4<size4; f4++) {
                tmp_buffer[f4] += factor * tmp_int_buffer[f4];
              }
            buffer_index += size34;
            int_buffer_index += sizem34;
            }
          c3++;
          }
        }
      ogc3 += c3;
      ogc3m += INT_NCART(am3-1);
      }
    index += size1234;

     /* The center 2 d/dy integrals */
    ogc3 = 0;
    ogc3m = 0;
    for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
      int am3 = shell3->am(gc3);
      // only integrals with y^n n>0 on center 2 contribute
      // so skip over n==0
      int c3 = 0;
      for (i3=0; i3<=am3; i3++) {
        for (k3=0; k3<=am3-i3-1; k3++) {
          double factor = -(am3-k3-i3);
          int c3ym1 = c3-i3;//=INT_CARTINDEX(am3-1,i3,j3-1)
          int buffer_index = size1234 + (ogc3+c3)*size4;
          int int_buffer_index = (ogc3m+c3ym1)*size4;
          for (f12=0; f12<size12; f12++) {
            double *tmp_buffer=&buffer[buffer_index];
            double *tmp_int_buffer=&int_buffer[int_buffer_index];
            for (f4=0; f4<size4; f4++) {
                tmp_buffer[f4] += factor * tmp_int_buffer[f4];
              }
            buffer_index += size34;
            int_buffer_index += sizem34;
            }
          c3++;
          }
        // account for the y^n n==0 case by an extra increment
        c3++;
        }
      ogc3 += c3;
      ogc3m += INT_NCART(am3-1);
      }
    index += size1234;

     /* The center 2 d/dz integrals */
    ogc3 = 0;
    ogc3m = 0;
    for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
      int am3 = shell3->am(gc3);
      // only integrals with z^n n>0 on center 2 contribute
      // so skip over n==0
      int c3 = 0;
      for (i3=0; i3<=am3; i3++) {
        // account for the z^n n==0 case by an extra increment
        c3++;
        for (k3=1; k3<=am3-i3; k3++) {
          double factor = -k3;
          int c3zm1 = c3-i3-1;//=INT_CARTINDEX(am3-1,i3,j3)
          int buffer_index = size1234+size1234+(ogc3+c3)*size4;
          int int_buffer_index = (ogc3m+c3zm1)*size4;
          for (f12=0; f12<size12; f12++) {
            double *tmp_buffer=&buffer[buffer_index];
            double *tmp_int_buffer=&int_buffer[int_buffer_index];
            for (f4=0; f4<size4; f4++) {
                tmp_buffer[f4] += factor * tmp_int_buffer[f4];
              }
            buffer_index += size34;
            int_buffer_index += sizem34;
            }
          c3++;
          }
        }
      ogc3 += c3;