hsosv1.cc

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

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

typedef int dmt_matrix;

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

#include <util/misc/formio.h>
#include <util/misc/timer.h>
#include <util/class/class.h>
#include <util/state/state.h>
#include <util/group/message.h>
#include <math/scmat/matrix.h>
#include <chemistry/molecule/molecule.h>
#include <chemistry/qc/scf/scf.h>
#include <chemistry/qc/mbpt/mbpt.h>
#include <chemistry/qc/mbpt/bzerofast.h>
#include <chemistry/qc/mbpt/hsosv1e1.h>

using namespace std;
using namespace sc;

static distsize_t
compute_v1_memory(int ni,
                  int nfuncmax, int nbasis, int noso,
                  int a_number, int nshell,
                  int ndocc, int nsocc, int nvir,
                  int nfzc, int nfzv,
                  int nproc)
{
  distsize_t mem = 0;
  int nocc = ndocc + nsocc;
  int dim_ij = nocc*ni - (ni*(ni-1))/2;
  mem += nproc*sizeof(int);
  mem += (noso+nsocc-nfzc-nfzv)*sizeof(double);
  mem += nfuncmax*nfuncmax*nbasis*ni*sizeof(double);
  mem += nfuncmax*nfuncmax*nbasis*ni*sizeof(double);
  mem += (distsize_t)nbasis*a_number*dim_ij*sizeof(double);
  mem += nvir*a_number*sizeof(double);
  mem += nvir*nvir*sizeof(double);
  if (nsocc) {
    mem += nsocc*sizeof(double);
    mem += ndocc*nsocc*(nvir-nsocc)*sizeof(double);
    mem += ndocc*nsocc*(nvir-nsocc)*sizeof(double);
    }
  mem += sizeof(double*)*(nbasis);
  mem += sizeof(double)*((nocc+nvir)*nbasis);
  return mem;
}

void
MBPT2::compute_hsos_v1()
{
  int i, j;
  int s1, s2;
  int a, b;
  int isocc, asocc;   /* indices running over singly occupied orbitals */
  int nfuncmax = basis()->max_nfunction_in_shell();
  int nvir;
  int nocc=0;
  int ndocc=0,nsocc=0;
  int i_offset; 
  int npass, pass;
  int ni;             /* batch size */
  int nr, ns; 
  int R, S;
  int q, r, s;
  int bf3,bf4;
  int docc_index, socc_index, vir_index;
  int me;
  int nproc;
  int rest;
  int a_rest;
  int a_number;              /* number of a-values processed by each node */
  int a_offset;
  int *a_vector;           /* each node's # of iajb integrals for one i,j */
  int compute_index;
  int tmp_index;
  int dim_ij;
  int nshell;
  double *evals_open;    /* reordered scf eigenvalues                      */
  double *trans_int1;    /* partially transformed integrals                */
  double *trans_int2;    /* partially transformed integrals                */
  double *trans_int3;    /* partially transformed integrals                */
  double *trans_int4_node;/* each node's subset of fully transf. integrals */
  double *trans_int4;    /* fully transformed integrals                    */
  double *mo_int_do_so_vir=0;/*mo integral (is|sa); i:d.o.,s:s.o.,a:vir     */
  double *mo_int_tmp=0;  /* scratch array used in global summations        */
  double *socc_sum=0;    /* sum of 2-el integrals involving only s.o.'s    */
  double *iqrs;
  double *iars_ptr, *iajs_ptr, *iajr_ptr;
  double iajr;
  double iars;
  double *iajb;
  double *c_qa;
  double *c_rb, *c_rj, *c_sj;
  double delta_ijab;
  double delta;
  double contrib1, contrib2;
  double ecorr_opt2=0,ecorr_opt1=0;
  double ecorr_zapt2;
  double ecorr_opt2_contrib=0, ecorr_zapt2_contrib=0;
  double escf;
  double eopt2,eopt1,ezapt2;
  double tol;     /* log2 of the erep tolerance (erep < 2^tol => discard) */
  int ithread;

  me = msg_->me();
  
  ExEnv::out0() << indent << "Just entered OPT2 program (opt2_v1)" << endl;

  tol = (int) (-10.0/log10(2.0));  /* discard ereps smaller than 10^-10 */

  nproc = msg_->n();
  ExEnv::out0() << indent << "nproc = " << nproc << endl;

  ndocc = nsocc = 0;
  const double epsilon = 1.0e-4;
  for (i=0; i<oso_dimension()->n(); i++) {
    if      (reference_->occupation(i) >= 2.0 - epsilon) ndocc++;
    else if (reference_->occupation(i) >= 1.0 - epsilon) nsocc++;
    }

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

  if (ndocc == 0 && nsocc == 0) {
    ExEnv::err0() << "There are no occupied orbitals; program exiting" << endl;
    abort();
    }

  if (nfzc > ndocc) {
    ExEnv::err0()
         << "The number of frozen core orbitals exceeds the number" << endl
         << "of doubly occupied orbitals; program exiting" << endl;
    abort();
    }

  if (nfzv > noso - ndocc - nsocc) {
    ExEnv::err0()
         << "The number of frozen virtual orbitals exceeds the number" << endl
         << "of unoccupied orbitals; program exiting" << endl;
    abort();
    }

  ndocc = ndocc - nfzc;
  /* nvir = # of unocc. orb. + # of s.o. orb. - # of frozen virt. orb. */
  nvir  = noso - ndocc - nfzc - nfzv; 
  /* nocc = # of d.o. orb. + # of s.o. orb - # of frozen d.o. orb. */
  nocc  = ndocc + nsocc;


  /* compute number of a-values (a_number) processed by each node */

  a_number = nvir/nproc; 
  a_rest = nvir%nproc;
  if (me < a_rest) a_number++;

  if (me == 0 && a_number < nsocc) { 
    ExEnv::err0() << "not enough memory allocated" << endl;
    /* must have all socc's on node 0 for computation of socc_sum*/
    abort();
    }

  if (me < a_rest) a_offset = me*a_number; /* a_offset for each node */
  else a_offset = a_rest*(a_number + 1) + (me - a_rest)*a_number;

  /* fill in elements of a_vector for gcollect */

  a_vector = (int*) malloc(nproc*sizeof(int));
  if (!a_vector) {
    ExEnv::errn() << "could not allocate storage for a_vector" << endl;
    abort();
    }
  for (i=0; i<nproc; i++) {
    a_vector[i] = nvir*(nvir/nproc)*sizeof(double);
    }
  for (i=0; i<a_rest; i++) {
    a_vector[i] += nvir*sizeof(double); /* first a_rest nodes hold an extra a */
    }

  // Cannot restart when singly occupied orbitals are present
  if (nsocc) {
    restart_orbital_v1_ = 0;
    }
  else if (restart_orbital_v1_) {
    ExEnv::out0() << indent
         << scprintf("Restarting at orbital %d with partial energy %18.14f",
                     restart_orbital_v1_, restart_ecorr_)
         << endl;
    }

  /* compute batch size ni for opt2 loops                                 *
   * need to store the following arrays: trans_int1-4, trans_int4_node,   *
   * scf_vector, evals_open, socc_sum, mo_int_do_so_vir, mo_int_tmp and   *
   * a_vector;                                                            *
   * since a_number is not the same on all nodes, use node 0's a_number   *
   * (which is >= all other a_numbers) and broadcast ni afterwords        */

  nshell = basis()->nshell();
  size_t memused = 0;
  ni = 0;
  for (i=1; i<=nocc-restart_orbital_v1_; i++) {
    distsize_t tmpmem = compute_v1_memory(i,
                                          nfuncmax, nbasis, noso,
                                          a_number, nshell,
                                          ndocc, nsocc, nvir,
                                          nfzc, nfzv, nproc);
    if (tmpmem > mem_alloc) break;
    ni = i;
    memused = distsize_to_size(tmpmem);
    }

  size_t mem_remaining = mem_alloc - memused;

  /* set ni equal to the smallest batch size for any node */
  msg_->min(ni);
  msg_->bcast(ni);

  ExEnv::out0() << indent
       << "Memory available per node:      " << mem_alloc << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Total memory used per node:     " << memused << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Memory required for one pass:   "
       << compute_v1_memory(nocc-restart_orbital_v1_,
                            nfuncmax, nbasis, noso, a_number, nshell,
                            ndocc, nsocc, nvir, nfzc, nfzv, nproc)
       << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Minimum memory required:        "
       << compute_v1_memory(1,
                            nfuncmax, nbasis, noso, a_number, nshell,
                            ndocc, nsocc, nvir, nfzc, nfzv, nproc)
       << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Batch size:                     " << ni
       << endl;

  if (ni < nsocc) {
    ExEnv::out0() << indent << "Not enough memory allocated to handle"
         << " SOCC orbs in first pass" << endl;
    abort();
    }

  if (ni < 1) {
    ExEnv::out0() << indent << "Not enough memory allocated" << endl;
    abort();
    }

  rest = (nocc-restart_orbital_v1_)%ni;
  npass = (nocc - restart_orbital_v1_ - rest)/ni + 1;
  if (rest == 0) npass--;

  if (me == 0) {
    ExEnv::out0() << indent << " npass  rest  nbasis  nshell  nfuncmax"
         << "  ndocc  nsocc  nvir  nfzc  nfzv" << endl;
    ExEnv::out0() << indent << scprintf("   %-4i   %-3i   %-5i   %-4i     %-3i"
                     "     %-3i     %-3i   %-3i    %-3i   %-3i",
                     npass,rest,nbasis,nshell,nfuncmax,ndocc,nsocc,nvir,nfzc,nfzv)
         << endl;
    }

  /* the scf vector might be distributed between the nodes, but for OPT2 *