compute_abs_a.cc

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

//
// compute_abs_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 <limits.h>

#include <scconfig.h>
#include <util/misc/formio.h>
#include <util/misc/timer.h>
#include <util/group/memory.h>
#include <util/group/message.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/basis/orthog.h>
#include <chemistry/qc/mbpt/util.h>
#include <chemistry/qc/mbpt/bzerofast.h>
#include <chemistry/qc/mbptr12/trans123_r12a_abs.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_abs_a.h>

using namespace std;
using namespace sc;

#define SINGLE_THREAD_E123   0
#define PRINT3Q 0
#define PRINT4Q 0
#define PRINT4Q_MP2 0
#define PRINT_NUM_TE_TYPES 4
#define PRINT_R12_INTERMED 0
#define LINDEP_TOL 1.e-6

#define USE_GLOBAL_ORTHOG 1

#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_abs_A::compute(RefSCMatrix& Vaa, RefSCMatrix& Xaa, RefSCMatrix& Baa,
			  RefSCMatrix& Vab, RefSCMatrix& Xab, RefSCMatrix& Bab)
{
  if (evaluated_)
    return;
  
  int debug_ = r12info()->debug_level();

  MolecularEnergy* mole = r12info()->mole();
  Ref<Integral> integral = r12info()->integral();
  Ref<GaussianBasisSet> bs = r12info()->basis();
  Ref<GaussianBasisSet> bs_aux = r12info()->basis_aux();
  bool two_basis_form = (bs != bs_aux);
  if (!two_basis_form)
    throw std::runtime_error("R12IntEval_abs_A::compute called when basis sets are identical");
  integral->set_basis(bs,bs,bs,bs_aux);

  Ref<SCMatrixKit> matrixkit_aux = bs_aux->matrixkit();
  Ref<SCMatrixKit> so_matrixkit_aux = bs_aux->so_matrixkit();
  Ref<MessageGrp> msg = r12info()->msg();
  Ref<MemoryGrp> mem = r12info()->mem();
  Ref<ThreadGrp> thr = r12info()->thr();
  const int num_te_types = 4;
  enum te_types {eri=0, r12=1, r12t1=2, r12t2=3};

  // 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-abs-mem");

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

  int nbasis = bs->nbasis();
  int nbasis_aux = bs_aux->nbasis();
  int nfuncmax = bs->max_nfunction_in_shell();
  int nfuncmax_aux = bs_aux->max_nfunction_in_shell();
  int nshell = bs->nshell();
  int nshell_aux = bs_aux->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;

  ExEnv::out0() << endl << indent
	       << "Entered ABS 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 + nbasis_aux*nbasis_aux; // scf vector + aux_basis orthogonalizer
    mem_static *= sizeof(double);
    int nthreads = thr->nthread();
    mem_static += nthreads * integral->storage_required_grt(bs,bs,bs,bs_aux); // 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,bs,bs,bs_aux);
  
  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_abs_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  nbasis(ABS) nshell(ABS) nfuncmax(ABS)") << endl;
  ExEnv::out0() << indent
		<< scprintf("  %-4i   %-3i   %-5i    %-4i     %-3i   %-5i        %-4i         %-3i",
			    npass,rest,nbasis,nshell,nfuncmax,nbasis_aux,nshell_aux,nfuncmax_aux)
		<< 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;


  //////////////////////////////////////////////////
  // For the auxiliary basis form an orthogonalizer
  // and distribute to nodes
  // Use canonical orthogonalization
  // ( does Wavefunction call that symmetric? )
  ////////////////////////////////////////////////

  RefSCMatrix orthog_aux;
  RefSCDimension osodim_aux;
  {
    //
    // Compute overlap matrix
    //
    
    // Make an Integral and initialize with bs_aux
    Ref<Integral> integral_aux = integral->clone();
    integral_aux->set_basis(bs_aux);
    Ref<PetiteList> pl_aux = integral_aux->petite_list();
    Ref<OneBodyInt> ov_aux_engine = integral_aux->overlap();

    // form skeleton s matrix
    RefSymmSCMatrix s(bs_aux->basisdim(), matrixkit_aux);
    Ref<SCElementOp> ov =
      new OneBodyIntOp(new SymmOneBodyIntIter(ov_aux_engine, pl_aux));
    
    s.assign(0.0);
    s.element_op(ov);
    ov=0;
    if (debug_ > 1) s.print("AO skeleton overlap (auxiliary basis)");
    
    // then symmetrize it
    RefSCDimension sodim_aux = pl_aux->SO_basisdim();
    RefSymmSCMatrix overlap_aux(sodim_aux, so_matrixkit_aux);
    pl_aux->symmetrize(s,overlap_aux);
    
    // and clean up a bit
    ov_aux_engine = 0;
    s = 0;
    integral_aux = 0;

    //
    // Compute orthogonalizer for the auxiliary basis
    //
#if USE_GLOBAL_ORTHOG
    OverlapOrthog orthog(OverlapOrthog::Canonical,
                         overlap_aux,
                         so_matrixkit_aux,
                         LINDEP_TOL,
                         0);

    RefSCMatrix orthog_aux_so = orthog.basis_to_orthog_basis();
    orthog_aux_so = orthog_aux_so.t();
    osodim_aux = orthog_aux_so.coldim();
#else    
    RefSCMatrix overlap_aux_eigvec;
    RefDiagSCMatrix overlap_isqrt_eigval;
    RefDiagSCMatrix overlap_sqrt_eigval;
    
    // diagonalize a copy of overlap_
    RefSymmSCMatrix M = overlap_aux.copy();
    RefSCMatrix U(M.dim(), M.dim(), M.kit());
    RefDiagSCMatrix m(M.dim(), M.kit());
    M.diagonalize(m,U);
    M = 0;
    Ref<SCElementMaxAbs> maxabsop = new SCElementMaxAbs;
    m.element_op(maxabsop.pointer());
    double maxabs = maxabsop->result();
    double s_tol = LINDEP_TOL * maxabs;
    
    double minabs = maxabs;
    BlockedDiagSCMatrix *bm = dynamic_cast<BlockedDiagSCMatrix*>(m.pointer());
    if (bm == 0) {
      ExEnv::out0() << "R12A_intermed_spinorb_abs: orthog_aux: expected blocked overlap" << endl;
    }
    
    double *pm_sqrt = new double[bm->dim()->n()];
    double *pm_isqrt = new double[bm->dim()->n()];
    int *pm_index = new int[bm->dim()->n()];
    int *nfunc = new int[bm->nblocks()];
    int nfunctot = 0;
    int nlindep = 0;
    for (int i=0; i<bm->nblocks(); i++) {
      nfunc[i] = 0;
      if (bm->block(i).null()) continue;
      int n = bm->block(i)->dim()->n();
      double *pm = new double[n];
      bm->block(i)->convert(pm);
      for (int j=0; j<n; j++) {
	if (pm[j] > s_tol) {
	  if (pm[j] < minabs) { minabs = pm[j]; }
	  pm_sqrt[nfunctot] = sqrt(pm[j]);
	  pm_isqrt[nfunctot] = 1.0/pm_sqrt[nfunctot];
	  pm_index[nfunctot] = j;
	  nfunc[i]++;
	  nfunctot++;
	}
	else {
	  nlindep++;
	}
      }
      delete[] pm;
    }
    
    if (debug_)
      ExEnv::out0() << indent << "Removed " << nlindep << " linearly dependent functions from the auxiliary basis" << endl;

    // make sure all nodes end up with exactly the same data
    msg->bcast(nfunctot);
    msg->bcast(nfunc, bm->nblocks());
    msg->bcast(pm_sqrt,nfunctot);
    msg->bcast(pm_isqrt,nfunctot);
    msg->bcast(pm_index,nfunctot);
    osodim_aux = new SCDimension(nfunctot, bm->nblocks(),
				 nfunc, "ortho aux SO (canonical)");
    for (int i=0; i<bm->nblocks(); i++) {
      osodim_aux->blocks()->set_subdim(i, new SCDimension(nfunc[i]));
    }

    overlap_aux_eigvec = so_matrixkit_aux->matrix(sodim_aux, osodim_aux);
    BlockedSCMatrix *bev
      = dynamic_cast<BlockedSCMatrix*>(overlap_aux_eigvec.pointer());
    BlockedSCMatrix *bU
      = dynamic_cast<BlockedSCMatrix*>(U.pointer());

    int ifunc = 0;
    for (int i=0; i<bev->nblocks(); i++) {