csgrad.cc

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

      }
    }
#endif

    integral_iajy = new double[noso];
    // in iajy: i act; a,j act or frz; y act or frz occ or virt.
    integral_ikja = new double[nvir_act];
    // in ikja: i,j act; k act or frz; a act.

    if (me == 0) {
      ExEnv::out0() << indent << "Begin fourth q.t." << endl;
      }

    // Begin fourth quarter transformation
    // generating MO integrals (ov|ov), (ov|oo) and (oo|ov)
    tim_enter("4. q.t.");
    index = 0;
    ij_index = 0;
    for (i=0; i<ni; i++) {
      for (j=0; j<nocc; j++) {
        if (index++ % nproc == me) {

          for (a=0; a<nvir; a++) {
            bzerofast(integral_iajy, noso);
            iajs_ptr = &integral_ixjs[a+nocc + nbasis*nbasis*ij_index];
            for (s=0; s<nbasis; s++) {
              c_sy = scf_vector[s];
              iajy_ptr = integral_iajy;
              tmpval = *iajs_ptr;
#if PRINT_BIGGEST_INTS
              biggest_ints_3.insert(tmpval,i+i_offset,j,s,a);
              if ((i+i_offset==105 && j == 2 && s == 170 && a == 3)
                  ||(i+i_offset==102 && j == 2 && s == 170 && a == 2)) {
                ExEnv::outn() << scprintf("3/4: %3d %3d %3d %3d: %16.10f",
                                 i+i_offset, j, s, x-nocc)
                     << endl;
                }
#endif
              for (y=0; y<noso; y++) {
                *iajy_ptr++ += *c_sy++ * tmpval;
                } // exit y loop
              iajs_ptr += nbasis;
              }   // exit s loop
            // Put integral_iajy into ixjs for one i,a,j while
            // overwriting elements of ixjs
            iajs_ptr = &integral_ixjs[a+nocc + nbasis*nbasis*ij_index];
            iajy_ptr = integral_iajy;
            for (y=0; y<noso; y++) {
              *iajs_ptr = *iajy_ptr++;
              iajs_ptr += nbasis;
              } // exit y loop
            }   // exit a loop

          // this is only needed for gradients
          if ((dograd || do_d1_) && j >= nfzc) {
            for (k=0; k<nocc; k++) {
              bzerofast(integral_ikja, nvir_act);
              ikjs_ptr = &integral_ixjs[k + nbasis*nbasis*ij_index];
              for (s=0; s<nbasis; s++) {
                c_sa = &scf_vector[s][nocc];
                ikja_ptr = integral_ikja;
                tmpval = *ikjs_ptr;
                for (a=0; a<nvir_act; a++) {
                  *ikja_ptr++ += *c_sa++ * tmpval;
                  } // exit a loop 
                ikjs_ptr += nbasis;
                }   // exit s loop 
              // Put integral_ikja into ixjs for one i,k,j while
              // overwriting elements of ixjs
              ikjs_ptr = &integral_ixjs[k + nbasis*(nocc + nbasis*ij_index)];
              ikja_ptr = integral_ikja;
              for (a=0; a<nvir_act; a++) {
                *ikjs_ptr = *ikja_ptr++;
                ikjs_ptr += nbasis;
                } // exit a loop 
              }   // exit k loop 
            }     //endif

          ij_index++;
          }   // endif
        }     // exit j loop
      }       // exit i loop
    // end of fourth quarter transformation
    tim_exit("4. q.t.");

    if (me == 0) {
      ExEnv::out0() << indent << "End of fourth q.t." << endl;
      }

    // The array integral_ixjs has now been overwritten by MO integrals
    // iajy and ikja, so rename the array mo_int
    mo_int = integral_ixjs;

    delete[] integral_iajy;
    delete[] integral_ikja;

    // Divide the (ia|jb) MO integrals by the term 
    // evals[i]+evals[j]-evals[a]-evals[b]
    // and keep these integrals in mo_int;
    // do this only for active i, j, a, and b
    tim_enter("divide (ia|jb)'s");

    index = 0;
    ij_index = 0;
    for (i=0; i<ni; i++) {
      ii = i+i_offset;
      for (j=0; j<nocc; j++) {
        if (index++ % nproc == me) {
          if (j>=nfzc) {
            for (b=0; b<nvir_act; b++) {
              iajb_ptr = &mo_int[nocc + nbasis*(b+nocc + nbasis*ij_index)];
              bb = b+nocc;
              for (a=0; a<nvir_act; a++) {
#if PRINT4Q
	      printf("4Q: (%d %d|%d %d) = %12.8f\n",
		     i,a+nocc,j,b+nocc,*iajb_ptr);
#endif
		// Zero out nonsymmetric integral, else divide by denominators
	        if (usep4 && ( symorb_irrep_[ii] ^
		      symorb_irrep_[j] ^
		      symorb_irrep_[a+nocc] ^
		      symorb_irrep_[bb]) ) {
		  *iajb_ptr++ = 0.0;
		}
		else
		  *iajb_ptr++ /= evals[ii]+evals[j]-evals[a+nocc]-evals[bb];
                } // exit a loop
              }   // exit b loop
            }     // endif
          ij_index++;
          }       // endif
        }         // exit j loop
      }           // exit i loop
    tim_exit("divide (ia|jb)'s");

    // We now have the fully transformed integrals (ia|jb)
    // divided by the proper orbital energy denominators
    // for one batch of i, all j<nocc, and all a<nvir and b<nvir,
    // where i, j, a, and b are all active;
    // compute contribution to the MP2 correlation energy
    // from these integrals 

#if WRITE_DOUBLES
    if (nproc > 1 || npass > 1) {
      ExEnv::outn() << "csgrad.cc: WRITE_DOUBLES set but case not allowed" << endl;
      abort();
      }
    ExEnv::outn() << "csgrad.cc: WRITING DOUBLES: CHECK ORDER" << endl;
    char *doutname = SCFormIO::fileext_to_filename(".mp2");
    FILE *dout = fopen(doutname,"w");
    delete[] doutname;
    fwrite(&nocc_act, sizeof(int), 1, dout);
    fwrite(&nvir_act, sizeof(int), 1, dout);
    for (j=nfzc; j<nocc; j++) {
      for (b=0; b<nvir_act; b++) {
        for (i=0; i<ni; i++) {
          ij_index = nocc*i + j;
          iajb_ptr = &mo_int[nocc + nbasis*(b+nocc + nbasis*ij_index)];
          fwrite(iajb_ptr, sizeof(double), nvir_act, dout);
          for (a=0; a<nvir_act; a++) {
            if (fabs(iajb_ptr[a])>1.0e-8) {
              ExEnv::outn() << scprintf(" Djbia(%2d %2d %2d %2d) = %12.8f",
                               j+1-nfzc,b+1,i+1,a+1,iajb_ptr[a])
                   << endl;
              }
            }
          }
        }
      }
    fclose(dout);
#endif

    tim_enter("compute ecorr");

    index = 0;
    ij_index = 0;
    for (i=0; i<ni; i++) {
      double ecorr_i = 0.0;
      for (j=0; j<nocc; j++) {
        double ecorr_ij = 0.0;
        if (index++ % nproc == me) {

          if (j>=nfzc) {
            for (b=0; b<nvir_act; b++) {
              iajb_ptr = &mo_int[nocc + nbasis*(b+nocc + nbasis*ij_index)];
              ibja_ptr = &mo_int[b+nocc + nbasis*(nocc + nbasis*ij_index)];
              for (a=0; a<nvir_act; a++) {
                delta_ijab = evals[i_offset+i]+evals[j]-evals[nocc+a]-evals[nocc+b];
                // only include determinants with unique coefficients
                if (a>=b && i_offset+i>=j) {
                  if (a>b && i_offset+i>j) {
                    // aaaa or bbbb
                    biggest_coefs.insert(*iajb_ptr - *ibja_ptr,
                                         i_offset+i,j,a,b,1111);
                    // aabb or bbaa or abba or baab
                    biggest_coefs.insert(*ibja_ptr,i_offset+i,j,b,a,1212);
                    } // endif
                  // aabb or bbaa or abba or baab
                  biggest_coefs.insert(*iajb_ptr,i_offset+i,j,a,b,1212);
                  } // endif

                tmpval = *iajb_ptr*(2**iajb_ptr - *ibja_ptr)*delta_ijab;
                ecorr_mp2 += tmpval;
                if (debug_) ecorr_ij += tmpval;
                iajb_ptr++;
                ibja_ptr += nbasis;;
                } // exit a loop
              }   // exit b loop
            }     // endif
          ij_index++;
          }       // endif
        if (debug_) {
          msg_->sum(ecorr_ij);
          ecorr_i += ecorr_ij;
          ExEnv::out0() << indent
               << scprintf("correlation energy for pair %3d %3d = %16.12f",
                           i+i_offset, j, ecorr_ij)
               << endl;
          }
        }         // exit j loop
      if (debug_) {
        ExEnv::out0() << indent
             << scprintf("correlation energy for orbital %3d = %16.12f",
                         i+i_offset, ecorr_i)
             << endl;
        }
      }           // exit i loop
    tim_exit("compute ecorr");

    // debug print
    if (debug_ && me == 0) {
      ExEnv::out0() << indent << "End of ecorr" << endl;
      }
    // end of debug print

    if (npass > 1 && pass < npass - 1) {
      double passe = ecorr_mp2;
      msg_->sum(passe);
      ExEnv::out0() << indent
           << "Partial correlation energy for pass " << pass << ":" << endl;
      ExEnv::out0() << indent
           << scprintf("  restart_ecorr          = %18.14f", passe)
           << endl;
      ExEnv::out0() << indent
           << scprintf("  restart_orbital_memgrp = %d", ((pass+1) * ni))
           << endl;
      }

    integral_iqjs = 0;
    mem->sync(); // Make sure MO integrals are complete on all nodes before continuing

    // don't go beyond this point if only the energy is needed
    if (!dograd && !do_d1_) continue;

    mo_int = (double*) mem->localdata();

    if (!dograd) goto compute_L;

    // Update the matrices Pkj and Wkj with
    // contributions from (occ vir|occ vir) integrals
    index = 0;
    ij_index = 0;
    tim_enter("Pkj and Wkj");
    for (i=0; i<ni; i++) {
      for (j=0; j<nocc; j++) {
        if (index++ % nproc == me) {

          if (j>=nfzc) {
            for (kloop=me; kloop<me+nocc; kloop++) {
              // stagger k's to minimize contention
              k = kloop%nocc;
              if (k<=j) pkj_ptr = &Pkj[j*(j+1)/2 + k];
              if (do_d2_) {
                if (k<=j && k>=nfzc) {
                  d2occ_mat_ptr = &d2occ_mat[(j-nfzc)*(j-nfzc+1)/2 + k-nfzc];
                  }
                }
              wjk_ptr = &Wkj[j*nocc + k];
              // Send for iakb, if necessary
              ik_index = (i*nocc + k)/nproc;
              ik_proc = (i*nocc + k)%nproc;
              ik_offset = nocc + nocc*nbasis + nbasis*nbasis*ik_index;
              mo_intbuf = (double*) membuf_remote.readonly_on_node(ik_offset,
                                                                   nbasis*nvir-nocc,
                                                                   ik_proc);
              for (a=0; a<nvir_act; a++) {
                ibja_ptr = &mo_int[nocc + nbasis*(a+nocc + nbasis*ij_index)];
                iajb_ptr = &mo_int[a+nocc + nbasis*(nocc + nbasis*ij_index)];
                iakb_ptr = &mo_intbuf[a];
                for (b=0; b<nvir_act; b++) {
                  tmpval1 = *iakb_ptr * *iajb_ptr;
                  tmpval = 2**iakb_ptr * *ibja_ptr++ - 4*tmpval1;
                  // tmpval = 2**iakb_ptr * (*ibja_ptr++ - 2 * *iajb_ptr);
                  iakb_ptr += nbasis;
                  iajb_ptr += nbasis;
                  if (k<nfzc && k<=j) *pkj_ptr += tmpval/(evals[k]-evals[j]);
                  else if (k<=j) {
                    *pkj_ptr += tmpval;
                    if (do_d2_) *d2occ_mat_ptr += tmpval1;
                    }
                  if (k>=nfzc) {
                    delta_ijab = evals[i_offset+i]+evals[j]-evals[nocc+a]-evals[nocc+b];
                    *wjk_ptr += tmpval*delta_ijab;
                    } 
                  } // exit b loop
                }   // exit a loop
              mo_intbuf = 0;
              membuf_remote.release();
              }     // end kloop loop
            }       // endif

          ij_index++;
          }         // endif
        }           // exit j loop
      }             // exit i loop
    tim_exit("Pkj and Wkj");

    // debug print
    if (debug_ && me == 0) {
      ExEnv::out0() << indent << "End of Pkj and Wkj" << endl;
      }
    // end of debug print

    // Update the matrices Pab and Wab with
    // contributions from (occ vir|occ vir) integrals
    tim_enter("Pab and Wab");
    index = 0;
    ij_index = 0;
    for (i=0; i<ni; i++) {
      for (j=0; j<nocc; j++) {
        if (index++ % nproc == me) {
          if (j>=nfzc) {

            offset = nocc + nocc*nbasis + nbasis*nbasis*ij_index;
            for (a=0; a<nvir_act; a++) {
              pab_ptr = &Pab[a*(a+1)/2];
              if (do_d2_) d2vir_mat_ptr = &d2vir_mat[a*(a+1)/2];
              for (b=0; b<=a; b++) {  // active-active part of Pab and Wab
                wab_ptr = &Wab[a*nvir + b];
                wba_ptr = &Wab[b*nvir + a];
                ibjc_ptr = &mo_int[offset + b];
                icjb_ptr = &mo_int[offset + b*nbasis];
                iajc_ptr = &mo_int[offset + a];
                icja_ptr = &mo_int[offset + a*nbasis];
                for (c=0; c<nvir_act; c++) {
                  tmpval = *iajc_ptr**ibjc_ptr;
                  if (do_d2_) *d2vir_mat_ptr += tmpval;
                  *pab_ptr += 4*tmpval - 2**iajc_ptr * *icjb_ptr;
                  // *pab_ptr += 2**iajc_ptr * (2 * *ibjc_ptr - *icjb_ptr);

                  if (a == b) {
                    delta_ijac = evals[i_offset+i]+evals[j]-evals[nocc+a]-evals[nocc+c];
                    *wab_ptr += 2**ibjc_ptr*(*icja_ptr - 2 * *iajc_ptr)*delta_ijac;
                    }
                  else {
                    delta_ijbc = evals[i_offset+i]+evals[j]-evals[nocc+b]-evals[nocc+c];
                    delta_ijac = evals[i_offset+i]+evals[j]-evals[nocc+a]-evals[nocc+c];
                    *wab_ptr += 2**ibjc_ptr * (*icja_ptr - 2 * *iajc_ptr)*delta_ijac;
                    *wba_ptr += 2**iajc_ptr * (*icjb_ptr - 2 * *ibjc_ptr)*delta_ijbc;
                    }
                  iajc_ptr += nbasis;
                  ibjc_ptr += nbasis;
                  icja_ptr++;
                  icjb_ptr++;
                  } // exit c loop
                pab_ptr++;
                d2vir_mat_ptr++;
                }   // exit b loop
              }     // exit a loop

            for (a=0; a<nfzv; a++) {        // active-frozen part of Pab
              pab_ptr = &Pab[(a+nvir_act)*(a+nvir_act+1)/2];
              for (b=0; b<nvir_act; b++) {
                tmpval = evals[nocc+b] - evals[nocc+nvir_act+a];
                ibjc_ptr = &mo_int[offset+b];
                iajc_ptr = &mo_int[offset+a+nvir_act];
                icja_ptr = &mo_int[offset+(a+nvir_act)*nbasis];
                for (c=0; c<nvir_act; c++) {
                  *pab_ptr += 2**ibjc_ptr*(2**iajc_ptr - *icja_ptr++)/tmpval;
                  ibjc_ptr += nbasis;
                  iajc_ptr += nbasis;
                  }  // exit c loop
                pab_ptr++;
                }    // exit b loop
              }      // exit a loop

            }        // endif
          ij_index++;
          }     // endif
        }       // exit j loop
      }         // exit i loop
    tim_exit("Pab and Wab");

    // debug print