mp2r12_energy.cc

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

//
// mp2r12_energy.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.
//

#ifdef __GNUG__
#pragma implementation
#endif

#include <stdexcept>
#include <util/misc/formio.h>
#include <util/misc/timer.h>
#include <util/ref/ref.h>
#include <math/scmat/local.h>
#include <chemistry/qc/mbpt/bzerofast.h>
#include <chemistry/qc/mbptr12/mp2r12_energy.h>
#include <chemistry/qc/mbptr12/vxb_eval_info.h>

using namespace std;
using namespace sc;

inline int max(int a,int b) { return (a > b) ? a : b;}

/*-------------
  MP2R12Energy
 -------------*/
static ClassDesc MP2R12Energy_cd(
  typeid(MP2R12Energy),"MP2R12Energy",1,"virtual public SavableState",
  0, 0, create<MP2R12Energy>);

MP2R12Energy::MP2R12Energy(Ref<R12IntEval>& r12eval, LinearR12::StandardApproximation stdapp, int debug)
{
  r12eval_ = r12eval;
  stdapprox_ = stdapp;
  if (debug >= 0)
    debug_ = debug;
  else
    debug_ = 0;
  evaluated_ = false;
}

MP2R12Energy::MP2R12Energy(StateIn& si) : SavableState(si)
{
  r12eval_ << SavableState::restore_state(si);

  RefSCDimension dim_aa = r12eval_->dim_aa();
  RefSCDimension dim_ab = r12eval_->dim_ab();
  Ref<SCMatrixKit> kit = r12eval_->r12info()->matrixkit();
  er12_aa_ = kit->vector(dim_aa);
  er12_ab_ = kit->vector(dim_ab);
  emp2r12_aa_ = kit->vector(dim_aa);
  emp2r12_ab_ = kit->vector(dim_ab);
 
  er12_aa_.restore(si);
  er12_ab_.restore(si);
  emp2r12_aa_.restore(si);
  emp2r12_ab_.restore(si);
  
  int stdapprox;
  si.get(stdapprox);
  stdapprox_ = (LinearR12::StandardApproximation) stdapprox;
  si.get(debug_);
  int evaluated;
  si.get(evaluated);
  evaluated_ = (bool) evaluated;
}

MP2R12Energy::~MP2R12Energy()
{
  r12eval_ = 0;
}

void MP2R12Energy::save_data_state(StateOut& so)
{
  SavableState::save_state(r12eval_.pointer(),so);

  er12_aa_.save(so);
  er12_ab_.save(so);
  emp2r12_aa_.save(so);
  emp2r12_ab_.save(so);
  
  so.put((int)stdapprox_);
  so.put(debug_);
  so.put((int)evaluated_);
}

void MP2R12Energy::obsolete()
{
  evaluated_ = false;
}

Ref<R12IntEval> MP2R12Energy::r12eval() const { return r12eval_; };
LinearR12::StandardApproximation MP2R12Energy::stdapp() const { return stdapprox_; };
void MP2R12Energy::set_debug(int debug) { debug_ = debug; };
int MP2R12Energy::get_debug() const { return debug_; };

double MP2R12Energy::energy()
{
  double value = emp2tot_aa_() + emp2tot_ab_() + er12tot_aa_() + er12tot_ab_();
  return value;
}

double MP2R12Energy::emp2tot_aa_() const
{
  RefSCVector emp2_aa = r12eval_->emp2_aa();
  int nij = emp2_aa.dim().n();
  double value = 0;
  for(int ij=0; ij<nij; ij++)
    value += emp2_aa.get_element(ij);

  return value;
}

double MP2R12Energy::emp2tot_ab_() const
{
  RefSCVector emp2_ab = r12eval_->emp2_ab();
  int nij = emp2_ab.dim().n();
  double value = 0;
  for(int ij=0; ij<nij; ij++)
    value += emp2_ab.get_element(ij);

  return value;
}

double MP2R12Energy::er12tot_aa_()
{
  compute();
  int nij = er12_aa_.dim().n();
  double value = 0;
  for(int ij=0; ij<nij; ij++)
    value += er12_aa_.get_element(ij);

  return value;
}

double MP2R12Energy::er12tot_ab_()
{
  compute();
  int nij = er12_ab_.dim().n();
  double value = 0;
  for(int ij=0; ij<nij; ij++)
    value += er12_ab_.get_element(ij);

  return value;
}

void MP2R12Energy::compute()
{
  if (evaluated_)
    return;
  
  Ref<R12IntEvalInfo> r12info = r12eval_->r12info();
  Ref<MessageGrp> msg = r12info->msg();
  int me = msg->me();
  int ntasks = msg->n();

  //
  // Evaluate pair energies:
  // distribute workload among nodes by pair index
  //

  // Need eigenvalues
  int nocc = r12info->nocc();
  int nfzc = r12info->nfzc();
  int nocc_act = r12info->nocc_act();
  RefDiagSCMatrix evalmat = r12eval_->evals();
  double* evals = new double[nocc_act];
  for(int i=nfzc; i<nocc; i++)
    evals[i-nfzc] = evalmat(i);
  evalmat = 0;

  // Get the intermediates
  RefSCMatrix Vaa = r12eval_->V_aa();
  RefSCMatrix Xaa = r12eval_->X_aa();
  RefSCMatrix Baa = r12eval_->B_aa();
  RefSCMatrix Vab = r12eval_->V_ab();
  RefSCMatrix Xab = r12eval_->X_ab();
  RefSCMatrix Bab = r12eval_->B_ab();
  RefSCVector emp2_aa = r12eval_->emp2_aa();
  RefSCVector emp2_ab = r12eval_->emp2_ab();

  // Prepare total and R12 pairs
  Ref<SCMatrixKit> localkit = Vaa.kit();
  RefSCDimension dim_aa = r12eval_->dim_aa();
  RefSCDimension dim_ab = r12eval_->dim_ab();
  int naa = dim_aa.n();
  int nab = dim_ab.n();
  emp2r12_aa_ = localkit->vector(dim_aa);
  emp2r12_ab_ = localkit->vector(dim_ab);
  er12_aa_ = localkit->vector(dim_aa);
  er12_ab_ = localkit->vector(dim_ab);
  double* er12_aa_vec = new double[naa];
  double* er12_ab_vec = new double[nab];
  bzerofast(er12_aa_vec,naa);
  bzerofast(er12_ab_vec,nab);

  //
  // Alpha-alpha pairs
  //
  
  // Allocate the B matrix:
  // 1) in MP2-R12/A the B matrix is the same for all pairs
  // 2) int MP2-R12/A' the B matrix is pair-specific
  RefSCMatrix Baa_inv = Baa.clone();
  if (stdapprox_ == LinearR12::StdApprox_A) {
    Baa_inv->assign(Baa);
    Baa_inv->gen_invert_this();
    if (debug_ > 1)
      Baa_inv.print("Inverse alpha-alpha MP2-R12/A B matrix");
  }
  
  int ij=0;
  for(int i=0; i<nocc_act; i++)
    for(int j=0; j<i; j++, ij++) {

      if (ij%ntasks != me)
        continue;

      RefSCVector Vaa_ij = Vaa.get_column(ij);

      // In MP2-R12/A' matrices B are pair-specific:
      // Form B(ij)kl,ow = Bkl,ow + 1/2(ek + el + eo + ew - 2ei - 2ej)Xkl,ow
      if (stdapprox_ == LinearR12::StdApprox_Ap) {
        Baa_inv.assign(Baa);
        int kl=0;
        for(int k=0; k<nocc_act; k++)
          for(int l=0; l<k; l++, kl++) {
            int ow=0;
            for(int o=0; o<nocc_act; o++)
              for(int w=0; w<o; w++, ow++) {
                double fx = 0.5 * (evals[k] + evals[l] + evals[o] + evals[w] - 2.0*evals[i] - 2.0*evals[j]) *
                Xaa.get_element(kl,ow);
                Baa_inv.accumulate_element(kl,ow,fx);
              }
          }
        if (debug_ > 1)
          Baa_inv.print("Alpha-alpha MP2-R12/A' B matrix");

        Baa_inv->gen_invert_this();

        if (debug_ > 1)
          Baa_inv.print("Inverse alpha-alpha MP2-R12/A' B matrix");
      }

      double eaa_ij = -2.0*Vaa_ij.dot(Baa_inv * Vaa_ij);
      er12_aa_vec[ij] = eaa_ij;
    }
  Baa_inv = 0;
  msg->sum(er12_aa_vec,naa,0,-1);
  er12_aa_->assign(er12_aa_vec);
  emp2r12_aa_->assign(emp2_aa);
  emp2r12_aa_->accumulate(er12_aa_);