hsosv2.cc

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

  if (debug_) {
    evals.print("eigenvalues");
    Scf_Vec.print("eigenvectors");
    }

  double *scf_vectort_dat = new double[nbasis*noso];
  Scf_Vec->convert(scf_vectort_dat);

  double** scf_vectort = new double*[nocc + nvir];

  int idoc = 0, ivir = 0, isoc = 0;
  for (i=nfzc; i<noso-nfzv; i++) {
    if (occ(i) >= 2.0 - epsilon) {
      evals_open[idoc+nsocc] = evals(i);
      scf_vectort[idoc+nsocc] = &scf_vectort_dat[i*nbasis];
      idoc++;
      }
    else if (occ(i) >= 1.0 - epsilon) {
      evals_open[isoc] = evals(i);
      scf_vectort[isoc] = &scf_vectort_dat[i*nbasis];
      evals_open[isoc+nocc] = evals(i);
      scf_vectort[isoc+nocc] = &scf_vectort_dat[i*nbasis];
      isoc++;
      }
    else {
      if (ivir < nvir) {
        evals_open[ivir+nocc+nsocc] = evals(i);
        scf_vectort[ivir+nocc+nsocc] = &scf_vectort_dat[i*nbasis];
        }
      ivir++;
      }
    }
  // need the transpose of the vector
  double **scf_vector = new double*[nbasis];
  double *scf_vector_dat = new double[(nocc+nvir)*nbasis];
  for (i=0; i<nbasis; i++) {
    scf_vector[i] = &scf_vector_dat[(nocc+nvir)*i];
    for (j=0; j<nocc+nvir; j++) {
      scf_vector[i][j] = scf_vectort[j][i];
      }
    }
  delete[] scf_vectort;
  delete[] scf_vectort_dat;

  /* allocate storage for various arrays */

  dim_ij = nocc*ni - (ni*(ni - 1))/2;

  trans_int1 = (double*) malloc(nfuncmax*nfuncmax*nbasis*ni*sizeof(double));
  trans_int2 = (double*) malloc(nvir*ni*sizeof(double));
  trans_int3 = (double*) malloc(nbf[me]*nvir*dim_ij*sizeof(double));
  trans_int4 = (double*) malloc(nvir*nvir*sizeof(double));
  trans_int4_tmp = (double*) malloc(nvir*nvir*sizeof(double));
  if (nsocc) socc_sum = (double*) malloc(nsocc*sizeof(double));
  if (nsocc) socc_sum_tmp = (double*) malloc(nsocc*sizeof(double));
  if (nsocc) mo_int_do_so_vir = 
                     (double*) malloc(ndocc*nsocc*(nvir-nsocc)*sizeof(double));
  if (nsocc) mo_int_tmp = 
                     (double*) malloc(ndocc*nsocc*(nvir-nsocc)*sizeof(double));

  if (nsocc) bzerofast(mo_int_do_so_vir,ndocc*nsocc*(nvir-nsocc));


  // create the integrals object
  integral()->set_storage(mem_remaining);
  tbint_ = integral()->electron_repulsion();
  intbuf = tbint_->buffer();

/**************************************************************************
 *   begin opt2 loops                                                     *
 **************************************************************************/


  for (pass=0; pass<npass; pass++) {
    i_offset = pass*ni;  
    if ((pass == npass - 1) && (rest != 0)) ni = rest;

    r_offset = 0;
    bzerofast(trans_int3,nbf[me]*nvir*dim_ij);

    tim_enter("RS loop");

    for (imyshell=0; imyshell<nRshell; imyshell++) {

      R = myshells[imyshell];
      nr = basis()->shell(R).nfunction();

      for (S = 0; S < nshell; S++) {
        ns = basis()->shell(S).nfunction();
        tim_enter("bzerofast trans_int1");
        bzerofast(trans_int1,nfuncmax*nfuncmax*nbasis*ni);
        tim_exit("bzerofast trans_int1");

        tim_enter("PQ loop");

        for (P = 0; P < nshell; P++) {
          np = basis()->shell(P).nfunction();

          for (Q = 0; Q <= P; Q++) {
            if (tbint_->log2_shell_bound(P,Q,R,S) < tol) {
              continue;                          /* skip ereps less than tol */
              }

            aoint_computed++;

            nq = basis()->shell(Q).nfunction();

            tim_enter("erep");
            tbint_->compute_shell(P,Q,R,S);
            tim_exit("erep");

            tim_enter("1. quart. tr.");

            index = 0;

            for (bf1 = 0; bf1 < np; bf1++) {
              p = basis()->shell_to_function(P) + bf1;
 
              for (bf2 = 0; bf2 < nq; bf2++) {
                q = basis()->shell_to_function(Q) + bf2;
                if (q > p) {
                  /* if q > p: want to skip the loops over bf3-4  */
                  /* and larger bf2 values, so increment bf1 by 1 */
                  /* ("break") and adjust the value of index      */
                  index = (bf1 + 1) * nq * nr * ns;
                  break;
                  }

                for (bf3 = 0; bf3 < nr; bf3++) {

                  for (bf4 = 0; bf4 < ns; bf4++,index++) {

                    if (fabs(intbuf[index]) > 1.0e-15) {
                      pqrs = intbuf[index];

                      iqrs = &trans_int1[((bf4*nr + bf3)*nbasis + q)*ni];
                      iprs = &trans_int1[((bf4*nr + bf3)*nbasis + p)*ni];
                    
                      if (p == q) pqrs *= 0.5;
                      col_index = i_offset;
                      c_pi = &scf_vector[p][col_index];
                      c_qi = &scf_vector[q][col_index];

                      for (i=ni; i; i--) {
                        *iqrs++ += pqrs * *c_pi++;
                        *iprs++ += pqrs * *c_qi++;
                        }
                      }
                    }   /* exit bf4 loop */
                  }     /* exit bf3 loop */
                }       /* exit bf2 loop */
              }         /* exit bf1 loop */
            tim_exit("1. quart. tr.");
            }           /* exit Q loop */
          }             /* exit P loop */
        tim_exit("PQ loop");

        /* begin second and third quarter transformations */

        for (bf3 = 0; bf3 < nr; bf3++) {
          r = r_offset + bf3;

          for (bf4 = 0; bf4 < ns; bf4++) {
            s = basis()->shell_to_function(S) + bf4;

            tim_enter("bzerofast trans_int2");
            bzerofast(trans_int2,nvir*ni);
            tim_exit("bzerofast trans_int2");

            tim_enter("2. quart. tr.");

            for (q = 0; q < nbasis; q++) {
              iars_ptr = trans_int2;
              iqrs = &trans_int1[((bf4*nr + bf3)*nbasis + q)*ni];
              c_qa = &scf_vector[q][nocc];

              for (a = 0; a < nvir; a++) {

                for (i=ni; i; i--) {
                  *iars_ptr++ += *c_qa * *iqrs++;
                  }

                iqrs -= ni;
                c_qa++;
                }
              }             /* exit q loop */
            tim_exit("2. quart. tr.");

            /* begin third quarter transformation */

            tim_enter("3. quart. tr.");

            for (i=0; i<ni; i++) {
              tmp_index = i*(i+1)/2 + i*i_offset;

              for (a=0; a<nvir; a++) {
                iars = trans_int2[a*ni + i];
                c_sj = scf_vector[s];
                iajr_ptr = &trans_int3[tmp_index + dim_ij*(a + nvir*r)];

                for (j=0; j<=i+i_offset; j++) {
                  *iajr_ptr++ += *c_sj++ * iars;
                  }
                }
              }   /* exit i loop */
            tim_exit("3. quart. tr.");

            } /* exit bf4 loop */
          }   /* exit bf3 loop */

        }         /* exit S loop */
      r_offset += nr;
      }           /* exit R loop */
    tim_exit("RS loop");


    /* begin fourth quarter transformation;                                *
     * first tansform integrals with only s.o. indices;                    *
     * these integrals are needed to compute the denominators              *
     * in the various terms contributing to the correlation energy         *
     * and must all be computed in the first pass;                         *
     * the integrals are summed into the array socc_sum:                   *
     * socc_sum[isocc] = sum over asocc of (isocc asocc|asocc isocc)       *
     * (isocc, asocc = s.o. and the sum over asocc runs over all s.o.'s)   *
     * the individual integrals are not saved here, only the sums are kept */


    if (pass == 0) {
      tim_enter("4. quart. tr.");
      if (nsocc) bzerofast(socc_sum,nsocc);
      for (isocc=0; isocc<nsocc; isocc++) {

        for (index=0; index<nbf[me]; index++) {
          i = 0;
          sum = basis()->shell(myshells[i]).nfunction();
          while (sum <= index) {
            i++;
            sum += basis()->shell(myshells[i]).nfunction();
            }
          sum -= basis()->shell(myshells[i]).nfunction();
          r = basis()->shell_to_function(myshells[i]) + index - sum;

          for (asocc=0; asocc<nsocc; asocc++) {
            socc_sum[isocc] += scf_vector[r][nocc+asocc]*
                           trans_int3[isocc*(isocc+1)/2 + isocc*i_offset
                                     + isocc + dim_ij*(asocc + nvir*index)];
            }
          }
        }       /* exit isocc loop */

      tim_exit("4. quart. tr.");

      /* sum socc_sum contributions from each node (only if nsocc > 0 *
       * since gop1 will fail if nsocc = 0)                           */
      if (nsocc > 0) {
        tim_enter("global sum socc_sum");
        msg_->sum(socc_sum,nsocc,socc_sum_tmp);
        tim_exit("global sum socc_sum");
        }

      } 

    /* now we have all the sums of integrals involving s.o.'s (socc_sum);   *
     * begin fourth quarter transformation for all integrals (including     *
     * integrals with only s.o. indices); use restriction j <= (i_offset+i) *
     * to save flops                                                        */

    compute_index = 0;

    for (i=0; i<ni; i++) {

      for (j=0; j <= (i_offset+i); j++) {

        tim_enter("4. quart. tr.");

        bzerofast(trans_int4,nvir*nvir);

        for (index=0; index<nbf[me]; index++) {
          k = 0;
          sum = basis()->shell(myshells[k]).nfunction();
          while (sum <= index) {
            k++;
            sum += basis()->shell(myshells[k]).nfunction();
            }
          sum -= basis()->shell(myshells[k]).nfunction();
          r = basis()->shell_to_function(myshells[k]) + index - sum;

          for (a=0; a<nvir; a++) {
            iajb = &trans_int4[a*nvir];
            iajr = trans_int3[i*(i+1)/2 + i*i_offset + j + dim_ij*(a+nvir*index)];
            c_rb = &scf_vector[r][nocc];

            for (b=0; b<nvir; b++) {
              *iajb++ += *c_rb++ * iajr;
              }
            }
          }

        tim_exit("4. quart. tr.");

        tim_enter("global sum trans_int4");
        msg_->sum(trans_int4,nvir*nvir,trans_int4_tmp);
        tim_exit("global sum trans_int4");

        /* we now have the fully transformed integrals (ia|jb)      *
         * for one i, one j (j <= i_offset+i), and all a and b;     *
         * compute contribution to the OPT1 and OPT2 correlation    *
         * energies; use restriction b <= a to save flops           */

        tim_enter("compute ecorr");

        for (a=0; a<nvir; a++) {
          for (b=0; b<=a; b++) {
            compute_index++;
            if (compute_index%nproc != me) continue;

            docc_index = ((i_offset+i) >= nsocc && (i_offset+i) < nocc) 
                        + (j >= nsocc && j < nocc);
            socc_index = ((i_offset+i)<nsocc)+(j<nsocc)+(a<nsocc)+(b<nsocc);
            vir_index = (a >= nsocc) + (b >= nsocc);

            if (socc_index >= 3) continue; /* skip to next b value */
 
            delta_ijab = evals_open[i_offset+i] + evals_open[j]