compute_sbs_a.cc

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

//
// compute_sbs_a.cc
//
// Copyright (C) 2003 Edward Valeev
//
// Author: Edward Valeev <[e-mail]>
// Maintainer: EV
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#include <stdexcept>
#include <stdlib.h>
#include <math.h>
#include <limits.h>

#include <scconfig.h>
#include <util/misc/formio.h>
#include <util/misc/timer.h>
#include <util/class/class.h>
#include <util/state/state.h>
#include <util/state/state_text.h>
#include <util/state/state_bin.h>
#include <math/scmat/matrix.h>
#include <chemistry/molecule/molecule.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/basis/obint.h>
#include <chemistry/qc/basis/symmint.h>
#include <chemistry/qc/mbpt/util.h>
#include <chemistry/qc/mbpt/bzerofast.h>
#include <chemistry/qc/mbptr12/mbptr12.h>
#include <chemistry/qc/mbptr12/trans12_grt.h>
#include <chemistry/qc/mbptr12/r12ia.h>
#include <chemistry/qc/mbptr12/r12ia_memgrp.h>
#include <chemistry/qc/mbptr12/r12ia_node0file.h>
#ifdef HAVE_MPIIO
  #include <chemistry/qc/mbptr12/r12ia_mpiiofile.h>
#endif
#include <chemistry/qc/mbptr12/vxb_eval_info.h>
#include <chemistry/qc/mbptr12/vxb_eval_sbs_a.h>

using namespace std;
using namespace sc;

#define SINGLE_THREAD_E12   0

#define PRINT2Q 0
#define PRINT3Q 0
#define PRINT4Q 0
#define PRINT4Q_MP2 0
#define PRINT_NUM_TE_TYPES 3
#define PRINT_R12_INTERMED 0

#if PRINT_BIGGEST_INTS
BiggestContribs biggest_ints_1(4,40);
#endif

#define WRITE_DOUBLES 0

#if PRINT_CONTRIB
static void
sw(int&i,int&j)
{
  int tmp = i;
  i = j;
  j = tmp;
}

static void
print_contrib(double tmpval, int num, int onum,
              int P,int Q,int R,int S, int p,int q,int r,int s)
{

  printf("noncanon: z(%d)(%d %d %d %d)(%d %d %d %d) contrib = % 6.4f\n",
         num, P, Q, R, S, p, q, r, s, tmpval);
  printf("noncanon: z(%d)(%d %d %d %d)(%d %d %d %d) contrib = % 6.4f\n",
         onum, P, Q, R, S, p, q, r, s, -tmpval);

  if (p < q) {
      sw(p,q); sw(P,Q);
    }
  if (r < s) {
      sw(r,s); sw(R,S);
    }
  if (p < r || (p == r && q < s)) {
      sw(P,R); sw(p,r);
      sw(Q,S); sw(q,s);
    }

  printf("z(%d)(%d %d %d %d)(%d %d %d %d) contrib = % 6.4f\n",
         num, P, Q, R, S, p, q, r, s, tmpval);
  printf("z(%d)(%d %d %d %d)(%d %d %d %d) contrib = % 6.4f\n",
         onum, P, Q, R, S, p, q, r, s, -tmpval);
}
#endif

/*-------------------------------------
  Based on MBPT2::compute_mp2_energy()
 -------------------------------------*/
void
R12IntEval_sbs_A::compute(RefSCMatrix& Vaa, RefSCMatrix& Xaa, RefSCMatrix& Baa,
			  RefSCMatrix& Vab, RefSCMatrix& Xab, RefSCMatrix& Bab,
			  RefSCVector& emp2_aa, RefSCVector& emp2_ab)
{
  if (evaluated_)
    return;
  
  int debug_ = r12info()->debug_level();

  MolecularEnergy* mole = r12info()->mole();
  Ref<Integral> integral = r12info()->integral();
  Ref<GaussianBasisSet> bs = r12info()->basis();
  bool two_basis_form = (bs != r12info()->basis_aux());
  Ref<MessageGrp> msg = r12info()->msg();
  Ref<MemoryGrp> mem = r12info()->mem();
  Ref<ThreadGrp> thr = r12info()->thr();
  const int num_te_types = 3;
  enum te_types {eri=0, r12=1, r12t1=2};

  // log2 of the erep tolerance
  // (erep < 2^tol => discard)
  const int tol = (int) (-10.0/log10(2.0));  // discard ints smaller than 10^-20

  int aoint_computed = 0; 

  BiggestContribs biggest_coefs(5,10);

#if PRINT_BIGGEST_INTS
  BiggestContribs biggest_ints_2(4,40);
  BiggestContribs biggest_ints_2s(4,40);
  BiggestContribs biggest_ints_3a(4,40);
  BiggestContribs biggest_ints_3(4,40);
#endif

  tim_enter("r12a-sbs-mem");

  int me = msg->me();
  int nproc = msg->n();
  const size_t mem_alloc = r12info()->memory();

  int nbasis = bs->nbasis();
  int nfuncmax = bs->max_nfunction_in_shell();
  int nshell = bs->nshell();
  int nocc = r12info()->nocc();
  int nocc_act = r12info()->nocc_act();
  int nfzc = r12info()->nfzc();
  int nfzv = r12info()->nfzv();
  int noso = r12info()->noso();
  int nvir  = noso - nocc;

  double pfac_xy_1, pfac_xy_2;
  if (two_basis_form) {
    pfac_xy_1 = 0.5;
    pfac_xy_2 = -0.5;
  }
  else {
    pfac_xy_1 = 0.5;
    pfac_xy_2 = 0.5;
  }

  ExEnv::out0() << endl << indent
	       << "Entered SBS A intermediates evaluator" << endl;
  ExEnv::out0() << indent << scprintf("nproc = %i", nproc) << endl;


  // Do a few preliminary tests to make sure the desired calculation
  // can be done (and appears to be meaningful!)

  if (nocc_act <= 0)
    throw std::runtime_error("There are no active occupied orbitals; program exiting");

  if (nvir-nfzv <= 0)
    throw std::runtime_error("There are no active virtual orbitals; program exiting");
  
  if (restart_orbital_) {
    ExEnv::out0() << indent
		  << scprintf("Restarting at orbital %d",
			      restart_orbital_)
		  << endl;
  }

  ////////////////////////////////////////////////////////
  // Compute batch size ni for mp2 loops;
  //
  // The following arrays are kept throughout (all of type double):
  //   scf_vector
  // The following objects are kept throughout:
  //   integrals evaluators
  // memory allocated for these arrays and objects is
  // called mem_static
  //
  ////////////////////////////////////////////////////////
  size_t mem_static = 0;
  int ni = 0;
  if (me == 0) {
    mem_static = nbasis*noso; // scf vector
    mem_static *= sizeof(double);
    int nthreads = thr->nthread();
    mem_static += nthreads * integral->storage_required_grt(bs); // integral evaluators
    ni = compute_transform_batchsize_(mem_alloc, mem_static, nocc_act-restart_orbital_, num_te_types); 
  }

  int max_norb = nocc_act - restart_orbital_;
  if (ni > max_norb)
    ni = max_norb;

  // Send value of ni and mem_static to other nodes
  msg->bcast(ni);
  double mem_static_double = (double) mem_static;
  msg->bcast(mem_static_double);
  mem_static = (size_t) mem_static_double;
  
  // Compute the storage remaining for the integral routines
  size_t dyn_mem = distsize_to_size(compute_transform_dynamic_memory_(ni,nocc_act,num_te_types));
  size_t mem_remaining = 0;
  if (mem_alloc <= (dyn_mem + mem_static)) mem_remaining += 0;
  else mem_remaining += mem_alloc - dyn_mem - mem_static;
  mem_remaining += thr->nthread() * integral->storage_required_grt(bs);
  
  ExEnv::out0() << indent
		<< "Memory available per node:      " << mem_alloc << " Bytes"
		<< endl;
  ExEnv::out0() << indent
		<< "Static memory used per node:    " << mem_static << " Bytes"
		<< endl;
  ExEnv::out0() << indent
		<< "Total memory used per node:     " << dyn_mem+mem_static << " Bytes"
		<< endl;
  ExEnv::out0() << indent
		<< "Memory required for one pass:   "
		<< compute_transform_dynamic_memory_(nocc_act,nocc_act,num_te_types)+mem_static
		<< " Bytes"
		<< endl;
  ExEnv::out0() << indent
		<< "Minimum memory required:        "
		<< compute_transform_dynamic_memory_(1,nocc_act,num_te_types)+mem_static
		<< " Bytes"
		<< endl;
  ExEnv::out0() << indent
		<< "Batch size:                     " << ni
		<< endl;
  
  if (ni == 0)
    throw std::runtime_error("R12IntEval_sbs_A: batch size is 0: more memory or processors are needed");
  
  if (r12info()->dynamic()) {
    ExEnv::out0() << indent << "Using dynamic load balancing." << endl;
  }

  int npass = 0;
  int rest = 0;
  if (ni == nocc_act-restart_orbital_) {
    npass = 1;
    rest = 0;
  }
  else {
    rest = (nocc_act-restart_orbital_)%ni;
    npass = (nocc_act-restart_orbital_ - rest)/ni + 1;
    if (rest == 0) npass--;
  }

  ExEnv::out0() << indent
		<< 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;
  {
    int index = 0;
    for (int i=0; i<ni; i++) {
      for (int j=0; j<nocc_act; j++) {
	if (index++ % nproc == me) nijmax++;
      }
    }
  }
  

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

  RefSCMatrix Scf_Vec = r12info()->scf_vec();
  RefDiagSCMatrix evalmat = r12info()->evals();
  int *mo_irrep = r12info()->orbsym();
  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);

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





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

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

  // Initialize the integrals
  integral->set_storage(mem_remaining);
  Ref<TwoBodyInt>* tbints_ = new Ref<TwoBodyInt>[thr->nthread()];
  for (int i=0; i<thr->nthread(); i++) {
    tbints_[i] = integral->grt();
  }