compute_sbs_a.cc

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

	    if (i != j && k != l) {
	      Xaa_ij[kl_aa] += Xaa_ijkl;
	      Xab_ji[lk_ab] += Xab_jilk;
	    }
            Tab_ij[kl_ab] += Tab_ijkl;
	    if (i != j)
	      Tab_ji[kl_ab] += Tab_jikl;
	    if (k != l)
	      Tab_ij[lk_ab] += Tab_ijlk;
	    if (i != j && k != l) {
	      Taa_ij[kl_aa] += Taa_ijkl;
	      Tab_ji[lk_ab] += Tab_jilk;
	    }
            tim_exit("MO ints contraction");

#if PRINT_R12_INTERMED
	    if (i != j && k != l)
	      printf("Vaa[%d][%d] = %lf\n",ij_aa,kl_aa,Vaa_ij[kl_aa]);
            printf("Vab[%d][%d] = %lf\n",ij_ab,kl_ab,Vab_ij[kl_ab]);
	    if (i != j)
	      printf("Vab[%d][%d] = %lf\n",ji_ab,kl_ab,Vab_ji[kl_ab]);
	    if (k != l)
	      printf("Vab[%d][%d] = %lf\n",ij_ab,lk_ab,Vab_ij[lk_ab]);
	    if (i != j && k != l)
	      printf("Vab[%d][%d] = %lf\n",ji_ab,lk_ab,Vab_ji[lk_ab]);
	    if (i != j && k != l)
	      printf("Xaa[%d][%d] = %lf\n",ij_aa,kl_aa,Xaa_ij[kl_aa]);
            printf("Xab[%d][%d] = %lf\n",ij_ab,kl_ab,Xab_ij[kl_ab]);
	    if (i != j)
	      printf("Xab[%d][%d] = %lf\n",ji_ab,kl_ab,Xab_ji[kl_ab]);
	    if (k != l)
	      printf("Xab[%d][%d] = %lf\n",ij_ab,lk_ab,Xab_ij[lk_ab]);
	    if (i != j && k != l)
	      printf("Xab[%d][%d] = %lf\n",ji_ab,lk_ab,Xab_ji[lk_ab]);
	    if (i != j && k != l)
	      printf("Taa[%d][%d] = %lf\n",ij_aa,kl_aa,Taa_ij[kl_aa]);
            printf("Tab[%d][%d] = %lf\n",ij_ab,kl_ab,Tab_ij[kl_ab]);
	    if (i != j)
	      printf("Tab[%d][%d] = %lf\n",ji_ab,kl_ab,Tab_ji[kl_ab]);
	    if (k != l)
	      printf("Tab[%d][%d] = %lf\n",ij_ab,lk_ab,Tab_ij[lk_ab]);
	    if (i != j && k != l)
	      printf("Tab[%d][%d] = %lf\n",ji_ab,lk_ab,Tab_ji[lk_ab]);
#endif
            r12intsacc->release_pair_block(i,j,R12IntsAcc::r12);
          }
        r12intsacc->release_pair_block(k,l,R12IntsAcc::eri);
        r12intsacc->release_pair_block(k,l,R12IntsAcc::r12);
        r12intsacc->release_pair_block(k,l,R12IntsAcc::r12t1);
	r12intsacc->release_pair_block(l,k,R12IntsAcc::r12t1);
      }
  }
  else {
    // Tasks that don't do any work here still need to create these timers
    tim_enter("MO ints retrieve");
    tim_exit("MO ints retrieve");
    tim_enter("MO ints contraction");
    tim_exit("MO ints contraction");
  }
  delete[] proc_with_ints;
  tim_exit("mp2-r12a intermeds");
  r12intsacc->deactivate();
  if (debug_)
    ExEnv::out0() << indent << "Computed intermediates V, X, and T" << endl;
  
  // If running in distributed environment use Message group to collect intermediates and pair energies on node 0
  if (nproc > 1) {
    // collect contributions from the nodes that computed the intermediates and broadcast to all nodes
    msg->sum(Vaa_ijkl,naa*naa,0,-1);
    msg->sum(Vab_ijkl,nab*nab,0,-1);
    msg->sum(Xaa_ijkl,naa*naa,0,-1);
    msg->sum(Xab_ijkl,nab*nab,0,-1);
    msg->sum(Taa_ijkl,naa*naa,0,-1);
    msg->sum(Tab_ijkl,nab*nab,0,-1);

    int naa = emp2_aa.dim().n();
    int nab = emp2_ab.dim().n();
    double* epair_aa = new double[naa];
    double* epair_ab = new double[nab];
    bzerofast(epair_aa,naa);
    bzerofast(epair_ab,nab);
    emp2_aa.convert(epair_aa);
    emp2_ab.convert(epair_ab);
    msg->sum(epair_aa,naa,0,-1);
    msg->sum(epair_ab,nab,0,-1);
    msg->sync();
    emp2_aa.assign(epair_aa);
    emp2_ab.assign(epair_ab);
    delete[] epair_aa;
    delete[] epair_ab;
  }
  
  if (debug_)
    ExEnv::out0() << indent << "Gathered intermediates V, X, and T and MP2 pair energies" << endl;

  // Initialize global intermediates
  Vaa->unit();
  Vab->unit();
  Baa->unit();
  Bab->unit();
  Xaa.assign(0.0);
  Xab.assign(0.0);
  
  // Add intermediates contribution to their global values
  for(int ij=0;ij<naa;ij++)
    for(int kl=0;kl<=ij;kl++) {
      int ijkl = ij*naa+kl;
      int klij = kl*naa+ij;
      double velem = Vaa->get_element(ij,kl) + Vaa_ijkl[ijkl];
      Vaa->set_element(ij,kl,velem);
      if (ij != kl) {
        velem = Vaa->get_element(kl,ij) + Vaa_ijkl[klij];
        Vaa->set_element(kl,ij,velem);
      }
      double xelem = Xaa->get_element(ij,kl) + Xaa_ijkl[ijkl];
      Xaa->set_element(ij,kl,xelem);
      if (ij != kl) {
        xelem = Xaa->get_element(kl,ij) + Xaa_ijkl[klij];
        Xaa->set_element(kl,ij,xelem);
      }
      double belem = Baa->get_element(ij,kl) + 0.5*(Taa_ijkl[ijkl] + Taa_ijkl[klij]);
      Baa->set_element(ij,kl,belem);
      Baa->set_element(kl,ij,belem);
    }

  for(int ij=0;ij<nab;ij++)
    for(int kl=0;kl<=ij;kl++) {
      int ijkl = ij*nab+kl;
      int klij = kl*nab+ij;
      double velem = Vab->get_element(ij,kl) + Vab_ijkl[ijkl];
      Vab->set_element(ij,kl,velem);
      if (ij != kl) {
        velem = Vab->get_element(kl,ij) + Vab_ijkl[klij];
        Vab->set_element(kl,ij,velem);
      }
      double xelem = Xab->get_element(ij,kl) + Xab_ijkl[ijkl];
      Xab->set_element(ij,kl,xelem);
      if (ij != kl) {
        xelem = Xab->get_element(kl,ij) + Xab_ijkl[klij];
        Xab->set_element(kl,ij,xelem);
      }
      double belem = Bab->get_element(ij,kl) + 0.5*(Tab_ijkl[ijkl] + Tab_ijkl[klij]);
      Bab->set_element(ij,kl,belem);
      Bab->set_element(kl,ij,belem);
    }
  msg->sum(aoint_computed);

#if PRINT_BIGGEST_INTS
  biggest_ints_1.combine(msg);
  biggest_ints_2.combine(msg);
  biggest_ints_2s.combine(msg);
  biggest_ints_3a.combine(msg);
  biggest_ints_3.combine(msg);
#endif

#if PRINT_BIGGEST_INTS
  ExEnv::out0() << "biggest 1/4 transformed ints" << endl;
  for (int i=0; i<biggest_ints_1.ncontrib(); i++) {
    ExEnv::out0() << scprintf("%3d %3d %3d %3d %16.12f",
                              biggest_ints_1.indices(i)[0],
                              biggest_ints_1.indices(i)[1],
                              biggest_ints_1.indices(i)[2],
                              biggest_ints_1.indices(i)[3],
                              biggest_ints_1.val(i)
                              )
    << endl;
  }
  ExEnv::out0() << "biggest 2/4 transformed ints" << endl;
  for (int i=0; i<biggest_ints_2.ncontrib(); i++) {
    ExEnv::out0() << scprintf("%3d %3d %3d %3d %16.12f",
                              biggest_ints_2.indices(i)[0],
                              biggest_ints_2.indices(i)[1],
                              biggest_ints_2.indices(i)[2],
                              biggest_ints_2.indices(i)[3],
                              biggest_ints_2.val(i)
                              )
    << endl;
  }
  ExEnv::out0() << "restricted 2/4 transformed ints" << endl;
  for (int i=0; i<biggest_ints_2s.ncontrib(); i++) {
    ExEnv::out0() << scprintf("%3d %3d %3d %3d %16.12f",
                              biggest_ints_2s.indices(i)[0],
                              biggest_ints_2s.indices(i)[1],
                              biggest_ints_2s.indices(i)[2],
                              biggest_ints_2s.indices(i)[3],
                              biggest_ints_2s.val(i)
                              )
    << endl;
  }
  ExEnv::out0() << "biggest 3/4 transformed ints (in 3.)" << endl;
  for (int i=0; i<biggest_ints_3a.ncontrib(); i++) {
    ExEnv::out0() << scprintf("%3d %3d %3d %3d %16.12f",
                              biggest_ints_3a.indices(i)[0],
                              biggest_ints_3a.indices(i)[1],
                              biggest_ints_3a.indices(i)[2],
                              biggest_ints_3a.indices(i)[3],
                              biggest_ints_3a.val(i)
                              )
    << endl;
  }
  ExEnv::out0() << "biggest 3/4 transformed ints (in 4.)" << endl;
  for (int i=0; i<biggest_ints_3.ncontrib(); i++) {
  ExEnv::out0() << scprintf("%3d %3d %3d %3d %16.12f",
                            biggest_ints_3.indices(i)[0],
                            biggest_ints_3.indices(i)[1],
                            biggest_ints_3.indices(i)[2],
                            biggest_ints_3.indices(i)[3],
                            biggest_ints_3.val(i)
                            )
    << endl;
  }
#endif

  // Print out various energies etc.

  if (debug_) {
    ExEnv::out0() << indent << "Number of shell quartets for which AO integrals\n"
    << indent << "would have been computed without bounds checking: "
    << npass*nshell*nshell*(nshell+1)*(nshell+1)/2
    << endl;

    ExEnv::out0() << indent << "Number of shell quartets for which AO integrals\n"
    << indent << "were computed: " << aoint_computed
    << endl;
  }
  
  /*--------
    Cleanup
   --------*/
  delete[] Vaa_ijkl;
  delete[] Taa_ijkl;
  delete[] Xaa_ijkl;
  delete[] Vab_ijkl;
  delete[] Tab_ijkl;
  delete[] Xab_ijkl;
  
  for (int i=0; i<thr->nthread(); i++) {
    delete e12thread[i];
  }
  delete[] e12thread;
  
  
  delete[] tbints_; tbints_ = 0;
  delete[] scf_vector;
  delete[] scf_vector_dat;
  delete[] evals;
  tim_exit("r12a-sbs-mem");

  evaluated_ = true;

  if (me == 0 && mole->if_to_checkpoint()) {
    StateOutBin stateout(mole->checkpoint_file());
    SavableState::save_state(mole,stateout);
    ExEnv::out0() << indent << "Checkpointed the wave function" << endl;
  }
  
  return;
}



///////////////////////////////////////////////////////
// Compute the batchsize for the transformation
//
// Only arrays allocated before exiting the loop over
// i-batches are included here  - only these arrays
// affect the batch size.
///////////////////////////////////////////////////////
int
R12IntEval_sbs_A::compute_transform_batchsize_(size_t mem_alloc, size_t mem_static, int nocc_act, const int num_te_types)
{
  // Check is have enough for even static objects
  size_t mem_dyn = 0;
  if (mem_alloc <= mem_static)
    return 0;
  else
    mem_dyn = mem_alloc - mem_static;

  // Determine if calculation is possible at all (i.e., if ni=1 possible)
  int ni = 1;
  distsize_t maxdyn = compute_transform_dynamic_memory_(ni, nocc_act, num_te_types);
  if (maxdyn > mem_dyn) {
    return 0;
  }

  ni = 2;
  while (ni<=nocc_act) {
    maxdyn = compute_transform_dynamic_memory_(ni, nocc_act, num_te_types);
    if (maxdyn >= mem_dyn) {
      ni--;
      break;
    }
    ni++;
  }
  if (ni > nocc_act) ni = nocc_act;

  return ni;
}

//////////////////////////////////////////////////////
// Compute required (dynamic) memory
// for a given batch size of the transformation
//
// Only arrays allocated before exiting the loop over
// i-batches are included here  - only these arrays
// affect the batch size.
//////////////////////////////////////////////////////
distsize_t
R12IntEval_sbs_A::compute_transform_dynamic_memory_(int ni, int nocc_act, const int num_te_types)
{
  int nproc = r12info()->msg()->n();

  ///////////////////////////////////////
  // the largest memory requirement will
  // occur just before
  // the end of the i-batch loop (mem)
  ///////////////////////////////////////

  // compute nij as nij on node 0, since nij on node 0 is >= nij on other nodes
  int index = 0;
  int nij = 0;
  for (int i=0; i<ni; i++) {
    for (int j=0; j<nocc_act; j++) {
      if (index++ % nproc == 0) nij++;
    }
  }

  int nbasis = r12info()->basis()->nbasis();
  int nfuncmax = r12info()->basis()->max_nfunction_in_shell();
  int nthread = r12info()->thr()->nthread();

  distsize_t memsize = sizeof(double)*(num_te_types*((distsize_t)nthread * ni * nbasis * nfuncmax * nfuncmax // iqrs
						     + (distsize_t)nij * 2 * nbasis * nfuncmax  // iqjs and iqjr buffers
						     + (distsize_t)nij * nbasis * nbasis // iqjs_contrib - buffer of half and higher
						     // transformed integrals
						     )
				       );
  return memsize;
}



////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "CLJ-CONDENSED"
// End: