csgrad.cc

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

         << scprintf(" npass  rest  nbasis  nshell  nfuncmax") << endl;
    ExEnv::out0() << indent
         << scprintf("  %-4i   %-3i   %-5i    %-4i     %-3i",
                     npass,rest,nbasis,nshell,nfuncmax)
         << endl;
    ExEnv::out0() << indent
         << scprintf(" nocc   nvir   nfzc   nfzv") << endl;
    ExEnv::out0() << indent
         << scprintf("  %-4i   %-4i   %-4i   %-4i",
                     nocc,nvir,nfzc,nfzv)
         << endl;
    }

  int nijmax = 0;
  index = 0;
  for (i=0; i<ni; i++) {
      for (j=0; j<nocc; j++) {
          if (index++ % nproc == me) nijmax++;
        }
    }

  ////////////////////////////////////////////////
  // The scf vector is distributed between nodes;
  // put a copy of the scf vector on each node;
  ////////////////////////////////////////////////

  escf = reference_->energy();
  hf_energy_ = escf;

  RefDiagSCMatrix occ;
  RefSCMatrix Scf_Vec;
  RefDiagSCMatrix evalmat;
  eigen(evalmat, Scf_Vec, occ);

  if (debug_ > 1) {
    evalmat.print("eigenvalues");
    Scf_Vec.print("eigenvectors");
    }

  double *scf_vector_dat = new double[nbasis*noso];
  Scf_Vec.t()->convert(scf_vector_dat);

  evals = new double[noso];
  double** scf_vector = new double*[nbasis];
  for (i=0; i<nbasis; i++) {
    scf_vector[i] = &scf_vector_dat[i*noso];
    }
  for (i=0; i<noso; i++) {
      evals[i] = evalmat(i);
    }

  Scf_Vec = 0;
  evalmat = 0;

  //////////////////////////////////////////////////////////////
  // Allocate storage for various arrays needed for gradients
  // (Pkj and Pab are symmetric, so store only lower triangle)
  //////////////////////////////////////////////////////////////

  if (dograd) {
    Pkj            = (double*) malloc((nocc*(nocc+1)/2)*sizeof(double));
    Pab            = (double*) malloc((nvir*(nvir+1)/2)*sizeof(double));
    Wkj            = (double*) malloc(nocc*nocc*sizeof(double));
    Wab            = (double*) malloc(nvir*nvir*sizeof(double));
    Waj            = (double*) malloc(nvir*nocc*sizeof(double));
    if (do_d2_) {
      d2occ_mat =  new double[nocc_act*(nocc_act+1)/2];
      d2vir_mat =  new double[nvir_act*(nvir_act+1)/2];
      }

    gradient_dat = new double[natom*3];
    gradient = new double*[natom];
    for (i=0; i<natom; i++) {
      gradient[i] = &gradient_dat[i*3];
      }

    hf_gradient_dat = new double[natom*3];
    hf_gradient = new double*[natom];
    for (i=0; i<natom; i++) {
      hf_gradient[i] = &hf_gradient_dat[i*3];
      }

    ginter = new double*[natom];
    for (i=0; i<natom; i++) {
      ginter[i] = new double[3];
      for (xyz=0; xyz<3; xyz++) ginter[i][xyz] = 0;
      }

    hf_ginter = new double*[natom];
    for (i=0; i<natom; i++) {
      hf_ginter[i] = new double[3];
      for (xyz=0; xyz<3; xyz++) hf_ginter[i][xyz] = 0;
      }

    //////////////////////////////
    // Initialize various arrays
    //////////////////////////////

    bzerofast(Pkj,nocc*(nocc+1)/2);
    bzerofast(Wkj,nocc*nocc);
    bzerofast(Pab,nvir*(nvir+1)/2);
    bzerofast(Wab,nvir*nvir);
    bzerofast(Waj,nvir*nocc);
    if (do_d2_) {
      bzerofast(d2occ_mat,nocc_act*(nocc_act+1)/2);
      bzerofast(d2vir_mat,nvir_act*(nvir_act+1)/2);
      }

    if (me == 0) zero_gradients(gradient, natom, 3);
    if (me == 0) zero_gradients(hf_gradient, natom, 3);
    }

  if (dograd || do_d1_) {
    Laj = (double*) malloc(nvir*nocc*sizeof(double));
    bzerofast(Laj,nvir*nocc);
    }

  if (debug_ > 2 && me == 0) {
    for (j=0; j<noso; j++) {
      ExEnv::out0() << indent
           << scprintf("eigenvalue[%3d] = %15.10lf", j, evals[j]);
      if (j < nfzc) ExEnv::out0() << " (frozen docc)";
      else if (j < nocc_act + nfzc) ExEnv::out0() << " (active docc)";
      else if (j < nvir_act + nocc_act + nfzc) ExEnv::out0() << " (active uocc)";
      else ExEnv::out0() << " (frozen uocc)";
      ExEnv::out0() << endl;
      }
    }

  /////////////////////////////////////
  //  Begin MP2 loops
  /////////////////////////////////////

  // debug print
  if (debug_ && me == 0) {
    ExEnv::out0() << indent
         << scprintf("node %i, begin loop over i-batches",me) << endl;
    }
  // end of debug print

  // Initialize the integrals
  integral()->set_storage(mem_remaining);
  tbints_ = new Ref<TwoBodyInt>[thr_->nthread()];
  for (i=0; i<thr_->nthread(); i++) {
      tbints_[i] = integral()->electron_repulsion();
    }
  if (dograd || do_d1_) {
    tbintder_ = new Ref<TwoBodyDerivInt>[thr_->nthread()];
    for (i=0; i<thr_->nthread(); i++) {
      tbintder_[i] = integral()->electron_repulsion_deriv();
      }
    }

  int mem_integral_intermediates = integral()->storage_used();
  int mem_integral_storage = (mem_remaining - mem_integral_intermediates) / thr_->nthread();
  if (mem_integral_storage<0) mem_integral_storage = 0;
  for (i=0; i<thr_->nthread(); i++) {
      tbints_[i]->set_integral_storage(mem_integral_storage);
    }

  ExEnv::out0() << endl << indent
       << scprintf("Memory used for integral intermediates: %i Bytes",
                   mem_integral_intermediates)
       << endl;
  ExEnv::out0() << indent
       << scprintf("Memory used for integral storage:       %i Bytes",
                   mem_integral_storage)
       << endl;

  if (mem.null()) {
      ExEnv::errn() << "MBPT2: memory group not initialized" << endl;
      abort();
    }

  mem->set_localsize(nijmax*nbasis*nbasis*sizeof(double));
  ExEnv::out0() << indent
       << "Size of global distributed array:       "
       << mem->totalsize()
       << " Bytes"
       << endl;

  MemoryGrpBuf<double> membuf_remote(mem);

  int usep4 = !dograd;

  Ref<ThreadLock> lock = thr_->new_lock();
  CSGradErep12Qtr** e12thread = new CSGradErep12Qtr*[thr_->nthread()];
  for (i=0; i<thr_->nthread(); i++) {
    e12thread[i] = new CSGradErep12Qtr(i, thr_->nthread(), me, nproc,
                                       mem, msg_, lock, basis(), tbints_[i],
                                       nocc, scf_vector, tol, debug_,
                                       dynamic_, usep4);
    }

    CSGrad34Qbtr** qbt34thread;
    if (dograd || do_d1_) {
      qbt34thread = new CSGrad34Qbtr*[thr_->nthread()];
      for (i=0; i<thr_->nthread(); i++) {
        qbt34thread[i] = new CSGrad34Qbtr(i, thr_->nthread(), me, nproc,
                                          mem, msg_, lock, basis(), tbints_[i],
                                          tbintder_[i], nocc, nfzc, scf_vector,
                                          tol, debug_, dynamic_, dograd, natom);
        }
      }

  tim_enter("mp2 passes");
  for (pass=0; pass<npass; pass++) {

    if (me == 0) {
      ExEnv::out0() << indent << "Beginning pass " << pass+1 << endl;
      }

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

    // Compute number of of i,j pairs on each node for
    // two-el integrals and non-sep 2PDM elements 
    index = 0;
    nij = 0;
    for (i=0; i<ni; i++) {
      for (j=0; j<nocc; j++) {
        if (index++ % nproc == me) nij++;
        }
      }

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

    mem->sync(); // This must be here or gamma non-sep will be wrong when running
                 // on multiple processors with more than one pass

    r_offset = 0;

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

    // Allocate (and initialize) some arrays

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

    bzerofast(integral_iqjs, nij*nbasis*nbasis);

    integral_iqjs = 0;

    mem->sync();

    index = 0;

    if (me == 0) {
      ExEnv::out0() << indent
           << scprintf("Begin loop over shells (erep, 1.+2. q.t.)") << endl;
      }

    // Do the two eletron integrals and the first two quarter transformations
    tim_enter("erep+1.qt+2.qt");
    for (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("erep+1.qt+2.qt");

    if (me == 0) {
      ExEnv::out0() << indent << "End of loop over shells" << endl;
      }

    mem->sync();  // Make sure iqjs is complete on each node before continuing

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

#if PRINT2Q
    if (me == 0) {
      int index = 0;
      int ij_index = 0;
      for (int i = 0; i<ni; i++) {
	for (int j = 0; j<nocc; j++) {
	  if (index++ % nproc == me) {
	    if (j >= nfzc) {
	      double *integral_ij_offset = integral_iqjs + nbasis*nbasis*ij_index;
	      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,q,j,s,*integral_ijsq_ptr);
		  *integral_ijsq_ptr++;
		}
	      }
	    }
	    ij_index++;
	  }
	}
      }
    }
#endif

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

    if (me == 0) {
      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 or frz, and x any MO (act or frz)
    index = 0;
    ij_index = 0;
    for (i=0; i<ni; i++) {
      for (j=0; j<nocc; j++) {
        if (index++ % nproc == me) {

          for (s=0; s<nbasis; s++) {

            bzerofast(ixjs_tmp, nbasis);
            for (q=0; q<nbasis; q++) {
              integral_iqjs_ptr = &integral_iqjs[q + nbasis*(s + nbasis*ij_index)];
              ixjs_ptr = ixjs_tmp;
              c_qx = scf_vector[q];
              tmpval = *integral_iqjs_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 (x=0; x<noso; x++) {
                *ixjs_ptr++ += *c_qx++ * tmpval;
                }
              }   // exit q loop

            // Put ixjs into integral_iqjs, while overwriting what was there;
            // i.e., integral_iqjs will now contain three-quarter transformed
            // integrals ixjs
            integral_iqjs_ptr = &integral_iqjs[nbasis*(s + nbasis*ij_index)];
            ixjs_ptr = ixjs_tmp;
            for (x=0; x<noso; x++) {
#if PRINT_BIGGEST_INTS
              if (x>=nocc) {
                biggest_ints_3a.insert(*ixjs_ptr,i+i_offset,j,s,x-nocc);
                }
#endif
              *integral_iqjs_ptr++ = *ixjs_ptr++;
              }
            }   // exit s loop
          ij_index++;
          }     // endif
        }       // exit j loop
      }         // exit i loop
    // end of third quarter transformation
    tim_exit("3. q.t.");

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

    delete[] ixjs_tmp;

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

#if PRINT3Q
    if (me == 0) {
      int index = 0;
      int ij_index = 0;
      for (int i = 0; i<ni; i++) {
	for (int j = 0; j<nocc; j++) {
	  if (index++ % nproc == me) {
	    if (j >= nfzc) {
	      double *integral_ij_offset = integral_ixjs + nbasis*nbasis*ij_index;
	      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,x,j,s,*integral_ijsx_ptr);
		  *integral_ijsx_ptr++;
		}
	      }
	    }
	    ij_index++;
	  }
	}