compute_sbs_a.cc

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

  ExEnv::out0() << indent
		<< scprintf("Memory used for integral storage:       %i Bytes",
			    integral->storage_used())
		<< endl;


  if (mem.null())
    throw std::runtime_error("R12IntEval_sbs_A: memory group not initialized");
  mem->set_localsize(num_te_types*nijmax*nbasis*nbasis*sizeof(double));
  ExEnv::out0() << indent
	       << "Size of global distributed array:       "
	       << mem->totalsize()
	       << " Bytes"
	       << endl;


  Ref<ThreadLock> lock = thr->new_lock();
  R12A_GRT_12Qtr** e12thread = new R12A_GRT_12Qtr*[thr->nthread()];
  for (int i=0; i<thr->nthread(); i++) {
    e12thread[i] = new R12A_GRT_12Qtr(i, thr->nthread(), me, nproc,
				      mem, msg, lock, bs, bs, tbints_[i],
				      nocc, nocc_act, scf_vector, tol, debug_,
				      r12info()->dynamic());
  }

  
  ///////////////////////////////////////////////////////////
  // Figure out which integrals accumulator should be used
  ///////////////////////////////////////////////////////////

  Ref<R12IntsAcc> r12intsacc;
  R12IntEvalInfo::StoreMethod ints_method = r12info()->ints_method();
  char *r12ints_file = r12info()->ints_file();
  bool restart = (restart_orbital_ > 0);

  switch (ints_method) {

  case R12IntEvalInfo::mem_only:
    if (restart)
      throw std::runtime_error("R12IntEval_sbs_A::compute -- cannot use MemoryGrp-based accumulator when restarting");
    ExEnv::out0() << indent << "Will hold transformed integrals in memory" << endl;
    r12intsacc = new R12IntsAcc_MemoryGrp(mem,num_te_types,nbasis,nbasis,nocc,nfzc);
    break;

  case R12IntEvalInfo::mem_posix:
    if (npass == 1) {
      ExEnv::out0() << indent << "Will hold transformed integrals in memory" << endl;
      r12intsacc = new R12IntsAcc_MemoryGrp(mem,num_te_types,nbasis,nbasis,nocc,nfzc);
      break;
    }
    // else use the next case
      
  case R12IntEvalInfo::posix:
    ExEnv::out0() << indent << "Will use POSIX I/O on node 0 to handle transformed integrals" << endl;
    r12intsacc = new R12IntsAcc_Node0File(mem,r12ints_file,num_te_types,nbasis,nbasis,nocc,nfzc,restart);
    break;

#if HAVE_MPIIO
  case R12IntEvalInfo::mem_mpi:
    if (npass == 1) {
      ExEnv::out0() << indent << "Will hold transformed integrals in memory" << endl;
      r12intsacc = new R12IntsAcc_MemoryGrp(mem,num_te_types,nbasis,nbasis,nocc,nfzc);
      break;
    }
    // else use the next case

  case R12IntEvalInfo::mpi:
    ExEnv::out0() << indent << "Will use MPI-IO (individual I/O) to handle transformed integrals" << endl;
    r12intsacc = new R12IntsAcc_MPIIOFile_Ind(mem,r12ints_file,num_te_types,nbasis,nbasis,nocc,nfzc,restart);
    break;
#endif
  
  default:
    throw std::runtime_error("R12IntEval_sbs_A::compute -- invalid integrals store method");
  }
  free(r12ints_file);


  /////////////////////////////////////////////////
  // zero out data arrays prior to the FIRST pass
  /////////////////////////////////////////////////

  if (!restart) {
    emp2_aa.assign(0.0);
    emp2_ab.assign(0.0);
  }
  
  /*-----------------------------------

    Start the integrals transformation

   -----------------------------------*/
  tim_enter("mp2-r12/a passes");
  if (me == 0 && mole->if_to_checkpoint() && r12intsacc->can_restart()) {
    StateOutBin stateout(mole->checkpoint_file());
    SavableState::save_state(mole,stateout);
    ExEnv::out0() << indent << "Checkpointed the wave function" << endl;
  }

  for (int pass=0; pass<npass; pass++) {

    ExEnv::out0() << indent << "Beginning pass " << pass+1 << endl;

    int i_offset = restart_orbital_ + pass*ni + nfzc;
    if ((pass == npass - 1) && (rest != 0)) ni = rest;

    // Compute number of of i,j pairs on each node during current pass for
    // two-el integrals
    int nij = 0;
    {
      int index = 0;
      for (int i=0; i<ni; i++) {
	for (int j=0; j<nocc_act; j++) {
	  if (index++ % nproc == me) nij++;
	}
      }
    }

    // debug print
    if (debug_)
      ExEnv::outn() << indent << "node " << me << ", nij = " << nij << endl;
    // end of debug print

    // Allocate and initialize some arrays
    // (done here to avoid having these arrays
    // overlap with arrays allocated later)

    // Allocate (and initialize) some arrays

    double* integral_ijsq = (double*) mem->localdata();
    bzerofast(integral_ijsq, (num_te_types*nij*nbasis*nbasis));
    integral_ijsq = 0;
    mem->sync();
    ExEnv::out0() << indent
		  << scprintf("Begin loop over shells (grt, 1.+2. q.t.)") << endl;

    // Do the two electron integrals and the first two quarter transformations
    tim_enter("grt+1.qt+2.qt");
    for (int i=0; i<thr->nthread(); i++) {
      e12thread[i]->set_i_offset(i_offset);
      e12thread[i]->set_ni(ni);
      thr->add_thread(i,e12thread[i]);
#     if SINGLE_THREAD_E12
      e12thread[i]->run();
#     endif
    }
#   if !SINGLE_THREAD_E12
    thr->start_threads();
    thr->wait_threads();
#   endif
    tim_exit("grt+1.qt+2.qt");
    ExEnv::out0() << indent << "End of loop over shells" << endl;

    mem->sync();  // Make sure ijsq is complete on each node before continuing
    integral_ijsq = (double*) mem->localdata();

#if PRINT2Q
    if (me == 0) {
      int index = 0;
      int ij_index = 0;
      int j_offset = nfzc;
      for (int i = 0; i<ni; i++) {
	for (int j = 0; j<nocc_act; j++) {
	  if (index++ % nproc == me) {
	    double *integral_ij_offset = integral_ijsq + num_te_types*nbasis*nbasis*ij_index;
	    for(int te_type=0; te_type<PRINT_NUM_TE_TYPES; te_type++,integral_ij_offset+=nbasis*nbasis) {
	      for (int s = 0; s<nbasis; s++) {
		double *integral_ijsq_ptr = integral_ij_offset + s*nbasis;
		for (int q = 0; q<nbasis; q++) {
		  printf("2Q: (%d %d|%d %d) = %12.8f\n",
		       i+i_offset,q,j+j_offset,s,*integral_ijsq_ptr);
		integral_ijsq_ptr++;
		}
	      }
	    }
	    ij_index++;
	  }
	}
      }
    }
#endif

    // Allocate and initialize some arrays
    double* ijsx_tmp = new double[nbasis];

    ExEnv::out0() << indent << "Begin third q.t." << endl;
    tim_enter("3. q.t.");
    // Begin third quarter transformation;
    // generate (ix|js) for i act, j act, and x any MO
    int index = 0;
    int ij_index = 0;
    for (int i=0; i<ni; i++) {
      for (int j=0; j<nocc_act; j++) {
        if (index++ % nproc == me) {
	  double* integral_ij_offset = integral_ijsq + num_te_types*nbasis*nbasis*ij_index;
	  for(int te_type=0; te_type<num_te_types; te_type++,integral_ij_offset+=nbasis*nbasis) {
	    for (int s=0; s<nbasis; s++) {
	      double* integral_ijsq_ptr = integral_ij_offset + s*nbasis;
	      bzerofast(ijsx_tmp, nbasis);
	      for (int q=0; q<nbasis; q++,integral_ijsq_ptr++) {
		double* ijsx_ptr = ijsx_tmp;
		double* c_qx = scf_vector[q];
		double tmpval = *integral_ijsq_ptr;
#if PRINT_BIGGEST_INTS
		biggest_ints_2.insert(tmpval,i+i_offset,j,s,q);
		if ((i+i_offset==104 && j == 1)
		    ||(i+i_offset==104 && j == 2)) {
		  biggest_ints_2s.insert(tmpval,i+i_offset,j,s,q);
		}
#endif
		for (int x=0; x<noso; x++) {
		  *ijsx_ptr++ += *c_qx++ * tmpval;
		}
              }   // exit q loop

	      // Put ixjs into integral_ijsq, while overwriting what was there;
	      // i.e., integral_ijsq will now contain three-quarter transformed
	      // integrals ijsx
	      integral_ijsq_ptr = integral_ij_offset + s*nbasis;
	      double* ijsx_ptr = ijsx_tmp;
	      for (int x=0; x<noso; x++) {
#if PRINT_BIGGEST_INTS
		if (x>=nocc) {
		  biggest_ints_3a.insert(*ijsx_ptr,i+i_offset,j,s,x-nocc);
                }
#endif
		*integral_ijsq_ptr++ = *ijsx_ptr++;
              }
            }   // exit s loop
	  }
	  ij_index++;
	}     // endif
      }       // exit j loop
    }         // exit i loop
    // end of third quarter transformation
    tim_exit("3. q.t.");
    ExEnv::out0() << indent << "End of third q.t." << endl;


    delete[] ijsx_tmp;

    // The array of half-transformed integrals integral_ijsq has now
    // been overwritten by three-quarter transformed integrals ixjs;
    // rename the array integral_ixjs, where x = any MO
    double* integral_ijsx = integral_ijsq;

#if PRINT3Q
    if (me == 0) {
      int index = 0;
      int ij_index = 0;
      int j_offset = nfzc;
      for (int i = 0; i<ni; i++) {
	for (int j = 0; j<nocc_act; j++) {
	  if (index++ % nproc == me) {
	    double *integral_ij_offset = integral_ijsx + num_te_types*nbasis*nbasis*ij_index;
	    for(int te_type=0; te_type<PRINT_NUM_TE_TYPES; te_type++,integral_ij_offset+=nbasis*nbasis) {
	      for (int s = 0; s<nbasis; s++) {
		double *integral_ijsx_ptr = integral_ij_offset + s*nbasis;
		for (int x = 0; x<noso; x++) {
		  printf("3Q: (%d %d|%d %d) = %12.8f\n",
		       i+i_offset,x,j+j_offset,s,*integral_ijsx_ptr);
		integral_ijsx_ptr++;
		}
	      }
	    }
	    ij_index++;
	  }
	}
      }
    }
#endif

    double* ijyx_tmp = new double[noso];
    // in ijyx: i act; x any MO, j act; y any MO.

    // Begin fourth quarter transformation
    ExEnv::out0() << indent << "Begin fourth q.t." << endl;
    tim_enter("4. q.t.");
    index = 0;
    ij_index = 0;
    for (int i=0; i<ni; i++) {
      for (int j=0; j<nocc_act; j++) {
        if (index++ % nproc == me) {
	  double *integral_ij_offset = integral_ijsx + num_te_types*nbasis*nbasis*ij_index;
	  for(int te_type=0; te_type<num_te_types; te_type++,integral_ij_offset+=nbasis*nbasis) {
	    for (int x=0; x<noso; x++) {
	      bzerofast(ijyx_tmp, noso);
	      double* ijsx_ptr = integral_ij_offset + x;
	      for (int s=0; s<nbasis; s++) {
		double* c_sy = scf_vector[s];
		double* ijyx_ptr = ijyx_tmp;
		double tmpval = *ijsx_ptr;
#if PRINT_BIGGEST_INTS
		biggest_ints_3.insert(tmpval,i+i_offset,j,s,x);
		if ((i+i_offset==105 && j == 2 && s == 170 && x == 3)
		    ||(i+i_offset==102 && j == 2 && s == 170 && x == 2)) {
		  ExEnv::outn() << scprintf("3/4: %3d %3d %3d %3d: %16.10f",
					   i+i_offset, j, s, x-nocc)
			       << endl;
                }
#endif
		for (int y=0; y<noso; y++) {
		  *ijyx_ptr++ += *c_sy++ * tmpval;
		} // exit y loop
		ijsx_ptr += nbasis;
	      }   // exit s loop
	      // Put ijyx_tmp into ijsx for one i,x,j while
	      // overwriting elements of ijsx
	      ijsx_ptr = integral_ij_offset + x;
	      double* ijyx_ptr = ijyx_tmp;
	      for (int y=0; y<noso; y++) {
		*ijsx_ptr = *ijyx_ptr++;
		ijsx_ptr += nbasis;
	      } // exit y loop
	    }   // exit x loop
	  }
          ij_index++;
	}   // endif
      }     // exit j loop
    }       // exit i loop
    // end of fourth quarter transformation
    tim_exit("4. q.t.");
    ExEnv::out0() << indent << "End of fourth q.t." << endl;
    
    // The array integral_ijsx has now been overwritten by MO integrals ijyx
    // rename the array mo_int
    double* mo_int = integral_ijsx;
    delete[] ijyx_tmp;

    // Zero out nonsymmetric integrals
    {
    int index = 0;
    int ij_index = 0;
    int j_offset = nfzc;
    for (int i = 0; i<ni; i++) {
      for (int j = 0; j<nocc_act; j++) {
	if (index++ % nproc == me) {
	  double *integral_ij_offset = mo_int + num_te_types*nbasis*nbasis*ij_index;
	  for(int te_type=0; te_type<PRINT_NUM_TE_TYPES; te_type++,integral_ij_offset+=nbasis*nbasis) {
	    for (int y = 0; y < noso; y++) {
	      double *integral_ijyx_ptr = integral_ij_offset + y*nbasis;
	      for (int x = 0; x<noso; x++) {
		if (( mo_irrep[i+i_offset] ^
		      mo_irrep[j+j_offset] ^
		      mo_irrep[x] ^
		      mo_irrep[y]) ) {
		  *integral_ijyx_ptr = 0.0;
		}
		integral_ijyx_ptr++;