comp2e.cc

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

      ogc3m += INT_NCART(am3-1);
      }
    index += size1234;
    }
  else if (dercenter == 3) {
    int ogc4,ogc4m,gc4,i4,k4,f123,size123;
    size123 = size1*size2*size3;

#ifdef DER_TIMING
  tim_change("- 3");
#endif
    /* The center 3 d/dx integrals */
    ogc4 = 0;
    ogc4m = 0;
    for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
      int am4 = shell4->am(gc4);
      // only integrals with x^n n>0 on center 3 contribute
      // so skip over n==0
      int c4 = am4+1;
      for (i4=1; i4<=am4; i4++) {
        double factor = -i4;
        for (k4=0; k4<=am4-i4; k4++) {
          int c4xm1 = c4-am4-1;//=INT_CARTINDEX(am4-1,i4-1,j4)
          int buffer_index = ogc4+c4;
          int int_buffer_index = ogc4m+c4xm1;
          for (f123=0; f123<size123; f123++) {
            buffer[buffer_index] += factor * int_buffer[int_buffer_index];
            buffer_index += size4;
            int_buffer_index += sizem4;
            }
          c4++;
          }
        }
      ogc4 += c4;
      ogc4m += INT_NCART(am4-1);
      }
    index += size1234;

     /* The center 3 d/dy integrals */
    ogc4 = 0;
    ogc4m = 0;
    for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
      int am4 = shell4->am(gc4);
      // only integrals with y^n n>0 on center 3 contribute
      // so skip over n==0
      int c4 = 0;
      for (i4=0; i4<=am4; i4++) {
        for (k4=0; k4<=am4-i4-1; k4++) {
          double factor = -(am4-k4-i4);
          int c4ym1 = c4-i4;//=INT_CARTINDEX(am4-1,i4,j4-1)
          int buffer_index = size1234 + ogc4+c4;
          int int_buffer_index = ogc4m+c4ym1;
          for (f123=0; f123<size123; f123++) {
            buffer[buffer_index] += factor * int_buffer[int_buffer_index];
            buffer_index += size4;
            int_buffer_index += sizem4;
            }
          c4++;
          }
        // account for the y^n n==0 case by an extra increment
        c4++;
        }
      ogc4 += c4;
      ogc4m += INT_NCART(am4-1);
      }
    index += size1234;

    /* The center 3 d/dz integrals */
    ogc4 = 0;
    ogc4m = 0;
    for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
      int am4 = shell4->am(gc4);
      // only integrals with z^n n>0 on center 3 contribute
      // so skip over n==0
      int c4 = 0;
      for (i4=0; i4<=am4; i4++) {
        // account for the z^n n==0 case by an extra increment
        c4++;
        for (k4=1; k4<=am4-i4; k4++) {
          double factor = -k4;
          int c4zm1 = c4-i4-1;//=INT_CARTINDEX(am4-1,i4,j4-1)
          int buffer_index = size1234+size1234+ogc4+c4;
          int int_buffer_index = ogc4m+c4zm1;
          for (f123=0; f123<size123; f123++) {
            buffer[buffer_index] += factor * int_buffer[int_buffer_index];
            buffer_index += size4;
            int_buffer_index += sizem4;
            }
          c4++;
          }
        }
      ogc4 += c4;
      ogc4m += INT_NCART(am4-1);
      }
    index += size1234;
    }

#ifdef DER_TIMING
  tim_change("+ erep");
#endif

  /* 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. */
  if (dercenter==0) int_expweight1 = 1;
  else if (dercenter==1) int_expweight2 = 1;
  else if (dercenter==2) int_expweight3 = 1;
  else if (dercenter==3) int_expweight4 = 1;
  old_perm = permute();
  set_permute(0);
  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);
  if (dercenter==0) int_expweight1 = 0;
  else if (dercenter==1) int_expweight2 = 0;
  else if (dercenter==2) int_expweight3 = 0;
  else if (dercenter==3) int_expweight4 = 0;

  /* Place the contributions into the user integral buffer. */
  index = 0;
  if (dercenter==0) {
    int ogc1,ogc1p,gc1,i1,k1,f234,size234;
    size234=size2*size3*size4;

#ifdef DER_TIMING
  tim_change("+ 0");
#endif
    /* The center 0 d/dx integrals */
    ogc1 = 0;
    ogc1p = 0;
    for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
      int am1 = shell1->am(gc1);
      int c1 = 0;
      for (i1=0; i1<=am1; i1++) {
        for (k1=0; k1<=am1-i1; k1++) {
          int c1xp1 = c1+am1+2;//=INT_CARTINDEX(am1+1,i1+1,j1)
          double *tmp_buffer=&buffer[index];
          double *tmp_int_buffer=&int_buffer[(ogc1p+c1xp1)*size234];
          for (f234=0; f234<size234; f234++) {
            tmp_buffer[f234] += tmp_int_buffer[f234];
            }
          index+=size234;
          c1++;
          }
        }
      ogc1 += c1;
      ogc1p += INT_NCART(am1+1);
      }

    /* The center 0 d/dy integrals */
    ogc1 = 0;
    ogc1p = 0;
    for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
      int am1 = shell1->am(gc1);
      int c1 = 0;
      for (i1=0; i1<=am1; i1++) {
        for (k1=0; k1<=am1-i1; k1++) {
          int c1yp1 = c1+i1;//=INT_CARTINDEX(am1+1,i1,j1+1)
          double *tmp_buffer=&buffer[index];
          double *tmp_int_buffer=&int_buffer[(ogc1p+c1yp1)*size234];
          for (f234=0; f234<size234; f234++) {
            tmp_buffer[f234] += tmp_int_buffer[f234];
            }
          index+=size234;
          c1++;
          }
        }
      ogc1 += c1;
      ogc1p += INT_NCART(am1+1);
      }

    /* The center 0 d/dz integrals */
    ogc1 = 0;
    ogc1p = 0;
    for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
      int am1 = shell1->am(gc1);
      int c1 = 0;
      for (i1=0; i1<=am1; i1++) {
        for (k1=0; k1<=am1-i1; k1++) {
          int c1zp1 = c1+i1+1;//=INT_CARTINDEX(am1+1,i1,j1)
          double *tmp_buffer=&buffer[index];
          double *tmp_int_buffer=&int_buffer[(ogc1p+c1zp1)*size234];
          for (f234=0; f234<size234; f234++) {
            tmp_buffer[f234] += tmp_int_buffer[f234];
            }
          index+=size234;
          c1++;
          }
        }
      ogc1 += c1;
      ogc1p += INT_NCART(am1+1);
      }
    }
  else if (dercenter == 1) {
    int ogc2,ogc2p,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;
    ogc2p = 0;
    for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
      int am2 = shell2->am(gc2);
      int c2=0;
      for (i2=0; i2<=am2; i2++) {
        for (k2=0; k2<=am2-i2; k2++) {
          int c2xp1 = c2+am2+2;//=INT_CARTINDEX(am2+1,i2+1,j2)
          int buffer_index = (ogc2+c2)*size34;
          int int_buffer_index = (ogc2p+c2xp1)*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] += tmp_int_buffer[f34];
              }
            buffer_index += size234;
            int_buffer_index += sizep234;
            }
          c2++;
          }
        }
      ogc2 += c2;
      ogc2p += INT_NCART(am2+1);
      }
    index += size1234;

     /* The center 1 d/dy integrals */
    ogc2 = 0;
    ogc2p = 0;
    for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
      int am2 = shell2->am(gc2);
      int c2 = 0;
      for (i2=0; i2<=am2; i2++) {
        for (k2=0; k2<=am2-i2; k2++) {
          int c2yp1 = c2+i2;//=INT_CARTINDEX(am2+1,i2,j2+1)
          int buffer_index = size1234 + (ogc2+c2)*size34;
          int int_buffer_index = (ogc2p+c2yp1)*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] += tmp_int_buffer[f34];
              }
            buffer_index += size234;
            int_buffer_index += sizep234;
            }
          c2++;
          }
        }
      ogc2 += c2;
      ogc2p += INT_NCART(am2+1);
      }
    index += size1234;

    /* The center 1 d/dz integrals */
    ogc2 = 0;
    ogc2p = 0;
    for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
      int am2 = shell2->am(gc2);
      int c2 = 0;
      for (i2=0; i2<=am2; i2++) {
        for (k2=0; k2<=am2-i2; k2++) {
          int c2zp1 = c2+i2+1;//=INT_CARTINDEX(am2+1,i2,j2+1)
          int buffer_index = size1234+size1234+(ogc2+c2)*size34;
          int int_buffer_index = (ogc2p+c2zp1)*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] += tmp_int_buffer[f34];
              }
            buffer_index += size234;
            int_buffer_index += sizep234;
            }
          c2++;
          }
        }
      ogc2 += c2;
      ogc2p += INT_NCART(am2+1);
      }
    index += size1234;
    }
  else if (dercenter == 2) {
    int ogc3,ogc3p,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;
    ogc3p = 0;
    for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
      int am3 = shell3->am(gc3);
      int c3 = 0;
      for (i3=0; i3<=am3; i3++) {
        for (k3=0; k3<=am3-i3; k3++) {
          int c3xp1 = c3+am3+2;//=INT_CARTINDEX(am3+1,i3+1,j3)
          int buffer_index = (ogc3+c3)*size4;
          int int_buffer_index = (ogc3p+c3xp1)*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] += tmp_int_buffer[f4];
              }
            buffer_index += size34;
            int_buffer_index += sizep34;
            }
          c3++;
          }
        }
      ogc3 += c3;
      ogc3p += INT_NCART(am3+1);
      }
    index += size1234;

     /* The center 2 d/dy integrals */
    ogc3 = 0;
    ogc3p = 0;
    for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
      int am3 = shell3->am(gc3);
      int c3 = 0;
      for (i3=0; i3<=am3; i3++) {
        for (k3=0; k3<=am3-i3; k3++) {
          int c3yp1 = c3+i3;//=INT_CARTINDEX(am3+1,i3,j3+1)
          int buffer_index = size1234 + (ogc3+c3)*size4;
          int int_buffer_index = (ogc3p+c3yp1)*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] += tmp_int_buffer[f4];
              }
            buffer_index += size34;
            int_buffer_index += sizep34;
            }
          c3++;
          }
        }
      ogc3 += c3;
      ogc3p += INT_NCART(am3+1);
      }
    index += size1234;

     /* The center 2 d/dz integrals */
    ogc3 = 0;
    ogc3p = 0;
    for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
      int am3 = shell3->am(gc3);
      int c3 = 0;
      for (i3=0; i3<=am3; i3++) {
        for (k3=0; k3<=am3-i3; k3++) {
          int c3zp1 = c3+i3+1;//=INT_CARTINDEX(am3+1,i3,j3)
          int buffer_index = size1234+size1234+(ogc3+c3)*size4;
          int int_buffer_index = (ogc3p+c3zp1)*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] += tmp_int_buffer[f4];
              }
            buffer_index += size34;
            int_buffer_index += sizep34;
            }
          c3++;
          }
        }
      ogc3 += c3;
      ogc3p += INT_NCART(am3+1);
      }
    index += size1234;
    }
  else if (dercenter == 3) {
    int ogc4,ogc4p,gc4,i4,k4,f123,size123;
    size123 = size1*size2*size3;

#ifdef DER_TIMING
  tim_change("+ 3");
#endif
    /* The center 3 d/dx 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 c4xp1 = c4+am4+2;//=INT_CARTINDEX(am4+1,i4+1,j4)
          int buffer_index = ogc4+c4;
          int int_buffer_index = ogc4p+c4xp1;
          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/dy 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++) {