mbpt.cc

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

//
// mbpt.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Ida Nielsen <[e-mail]>
// Maintainer: LPS
//
// 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 __GNUC__
#pragma implementation
#endif

#include <util/misc/formio.h>
#include <util/misc/exenv.h>
#include <util/state/stateio.h>
#include <math/scmat/blocked.h>
#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/mbpt/mbpt.h>

using namespace std;
using namespace sc;

/////////////////////////////////////////////////////////////////
// Function dquicksort performs a quick sort (smaller -> larger) 
// of the double data in item by the integer indices in index;
// data in item remain unchanged
/////////////////////////////////////////////////////////////////
static void
dqs(double *item,int *index,int left,int right)
{
  register int i,j;
  double x;
  int y;
 
  i=left; j=right;
  x=item[index[(left+right)/2]];
 
  do {
    while(item[index[i]]<x && i<right) i++;
    while(x<item[index[j]] && j>left) j--;
 
    if (i<=j) {
      if (item[index[i]] != item[index[j]]) {
        y=index[i];
        index[i]=index[j];
        index[j]=y;
        }
      i++; j--;
      }
    } while(i<=j);
       
  if (left<j) dqs(item,index,left,j);
  if (i<right) dqs(item,index,i,right);
}
static void
dquicksort(double *item,int *index,int n)
{
  int i;
  if (n<=0) return;
  for (i=0; i<n; i++) {
    index[i] = i;
    }
  dqs(item,index,0,n-1);
  }

///////////////////////////////////////////////////////////////////////////
// MBPT2

static ClassDesc MBPT2_cd(
  typeid(MBPT2),"MBPT2",8,"public Wavefunction",
  0, create<MBPT2>, create<MBPT2>);

MBPT2::MBPT2(StateIn& s):
  SavableState(s),
  Wavefunction(s)
{
  reference_ << SavableState::restore_state(s);
  s.get(nfzc);
  s.get(nfzv);
  if (s.version(::class_desc<MBPT2>()) >= 8) {
      double dmem_alloc;
      s.get(dmem_alloc);
      mem_alloc = size_t(dmem_alloc);
    }
  else {
      unsigned int imem_alloc;
      s.get(imem_alloc);
      mem_alloc = imem_alloc;
    }
  s.getstring(method_);
  s.getstring(algorithm_);

  if (s.version(::class_desc<MBPT2>()) <= 6) {
      int debug_old;
      s.get(debug_old);
    }

  if (s.version(::class_desc<MBPT2>()) >= 2) {
      s.get(do_d1_);
    }
  else {
      do_d1_ = 0;
    }

  if (s.version(::class_desc<MBPT2>()) >= 3) {
      s.get(dynamic_);
    }
  else {
      dynamic_ = 0;
    }

  if (s.version(::class_desc<MBPT2>()) >= 4) {
      s.get(cphf_epsilon_);
    }
  else {
      cphf_epsilon_ = 1.0e-8;
    }

  if (s.version(::class_desc<MBPT2>()) >= 5) {
      s.get(max_norb_);
    }
  else {
      max_norb_ = 0;
    }

  if (s.version(::class_desc<MBPT2>()) >= 6) {
      s.get(do_d2_);
    }
  else {
      do_d2_ = 1;
    }

  hf_energy_ = 0.0;

  symorb_irrep_ = 0;
  symorb_num_ = 0;

  eliminate_in_gmat_ = 1;

  mem = MemoryGrp::get_default_memorygrp();
  msg_ = MessageGrp::get_default_messagegrp();
  thr_ = ThreadGrp::get_default_threadgrp();

  restart_ecorr_ = 0.0;
  restart_orbital_v1_ = 0;
  restart_orbital_memgrp_ = 0;
}

MBPT2::MBPT2(const Ref<KeyVal>& keyval):
  Wavefunction(keyval)
{
  reference_ << keyval->describedclassvalue("reference");
  if (reference_.null()) {
      ExEnv::err0() << "MBPT2::MBPT2: no reference wavefunction" << endl;
      abort();
    }
  copy_orthog_info(reference_);
  nfzc = keyval->intvalue("nfzc");
  char *nfzc_charval = keyval->pcharvalue("nfzc");
  if (nfzc_charval && !strcmp(nfzc_charval, "auto")) {
      if (molecule()->max_z() > 30) {
          ExEnv::err0()
               << "MBPT2: cannot use \"nfzc = auto\" for Z > 30" << endl;
          abort();
        }
      nfzc = molecule()->n_core_electrons()/2;
      ExEnv::out0() << indent
           << "MBPT2: auto-freezing " << nfzc << " core orbitals" << endl;
    }
  delete[] nfzc_charval;
  nfzv = keyval->intvalue("nfzv");
  mem_alloc = keyval->sizevalue("memory");
  if (keyval->error() != KeyVal::OK) {
      // by default, take half of the memory
      mem_alloc = ExEnv::memory()/2;
      if (mem_alloc == 0) {
          mem_alloc = DEFAULT_SC_MEMORY;
        }
    }
  mem << keyval->describedclassvalue("memorygrp");
  if (mem.null()) {
      mem = MemoryGrp::get_default_memorygrp();
    }
  msg_ = MessageGrp::get_default_messagegrp();
  thr_ = ThreadGrp::get_default_threadgrp();

  method_ = keyval->pcharvalue("method");

  algorithm_ = keyval->pcharvalue("algorithm");

  do_d1_ = keyval->booleanvalue("compute_d1");
  do_d2_ = keyval->booleanvalue("compute_d2",KeyValValueboolean(1));

  restart_ecorr_ = keyval->doublevalue("restart_ecorr");
  restart_orbital_v1_ = keyval->intvalue("restart_orbital_v1");
  restart_orbital_memgrp_ = keyval->intvalue("restart_orbital_memgrp");

  KeyValValueint default_dynamic(0);
  dynamic_ = keyval->booleanvalue("dynamic", default_dynamic);

  cphf_epsilon_ = keyval->doublevalue("cphf_epsilon",KeyValValuedouble(1.e-8));

  max_norb_ = keyval->intvalue("max_norb",KeyValValueint(-1));

  hf_energy_ = 0.0;

  symorb_irrep_ = 0;
  symorb_num_ = 0;

  eliminate_in_gmat_ = 1;
}

MBPT2::~MBPT2()
{
  delete[] method_;
  delete[] algorithm_;
  delete[] symorb_irrep_;
  delete[] symorb_num_;
}

void
MBPT2::save_data_state(StateOut& s)
{
  Wavefunction::save_data_state(s);
  SavableState::save_state(reference_.pointer(),s);
  s.put(nfzc);
  s.put(nfzv);
  double dmem_alloc = mem_alloc;
  s.put(dmem_alloc);
  s.putstring(method_);
  s.putstring(algorithm_);
  s.put(do_d1_);
  s.put(dynamic_);
  s.put(cphf_epsilon_);
  s.put(max_norb_);
  s.put(do_d2_);
}

void
MBPT2::print(ostream&o) const
{
  o << indent << "MBPT2:" << endl;
  o << incindent;
  Wavefunction::print(o);
  o << indent << "Reference Wavefunction:" << endl;
  o << incindent; reference_->print(o); o << decindent << endl;
  o << decindent;
}

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

int
MBPT2::spin_polarized()
{
  return reference_->spin_polarized();
}

RefSymmSCMatrix
MBPT2::density()
{
  return 0;
}

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

void
MBPT2::compute()
{
  init_variables();

  reference_->set_desired_value_accuracy(desired_value_accuracy()
                                         / ref_to_mp2_acc);
  if (gradient_needed()) {
      if (nsocc) {
          ExEnv::errn() << "MBPT2: cannot compute open shell gradients" << endl;
          abort();
        }
      compute_cs_grad();
    }
  else {
      if (nsocc && algorithm_ && !strcmp(algorithm_,"memgrp")) {
          ExEnv::errn() << "MBPT2: memgrp algorithm cannot compute open shell energy"
               << endl;
          abort();
        }
      if (!nsocc && (!algorithm_ || !strcmp(algorithm_,"memgrp"))) {
          compute_cs_grad();
        }
      else if ((!algorithm_ && msg_->n() <= 32)
               || (algorithm_ && !strcmp(algorithm_,"v1"))) {
          compute_hsos_v1();
        }
      else if (!algorithm_ || !strcmp(algorithm_,"v2")) {
          compute_hsos_v2();
        }
      else if (!strcmp(algorithm_,"v2lb")) {
          compute_hsos_v2_lb();
        }
      else {
          ExEnv::errn() << "MBPT2: unknown algorithm: " << algorithm_ << endl;
          abort();
        }
    }
}

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

void
MBPT2::obsolete()