hsosv2.cc

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

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

#include <iostream>
#include <math.h>

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

using namespace std;
using namespace sc;

static void iqs(int *item,int *index,int left,int right);
static void iquicksort(int *item,int *index,int n);

void
MBPT2::compute_hsos_v2()
{
  int i, j, k;
  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 nshell;
  int nocc=0,ndocc=0,nsocc=0;
  int i_offset; 
  int npass, pass;
  int ni;
  int np, nq, nr, ns; 
  int P, Q, R, S;
  int p, q, r, s;
  int bf1, bf2, bf3, bf4;
  int index;
  int compute_index;
  int col_index;
  int tmp_index;
  int dim_ij;
  int docc_index, socc_index, vir_index;
  int me;
  int nproc;
  int rest;
  int r_offset;
  int sum;
  int min;
  int iproc;
  int nRshell;
  int imyshell;
  int *myshells;       /* the R indices processed by node me */
  int *shellsize;      /* size of each shell                 */
  int *sorted_shells;  /* sorted shell indices: large shells->small shells */
  int *nbf;            /* number of basis functions processed by each node */
  int *proc;           /* element k: processor which will process shell k  */
  int aoint_computed = 0; 
  double *evals_open;    /* reordered scf eigenvalues                      */
  const double *intbuf;  /* 2-electron AO integral buffer                  */
  double *trans_int1;    /* partially transformed integrals                */
  double *trans_int2;    /* partially transformed integrals                */
  double *trans_int3;    /* partially transformed integrals                */
  double *trans_int4;    /* fully transformed integrals                    */
  double *trans_int4_tmp; /* scratch array                                 */ 
  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 *socc_sum_tmp=0;/* scratch array                                  */ 
  double *iqrs, *iprs;
  double *iars_ptr;
  double iars;
  double iajr;
  double *iajr_ptr;
  double *iajb;
  double pqrs;
  double *c_qa;
  double *c_rb, *c_pi, *c_qi, *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) */

  me = msg_->me();

  ExEnv::out0() << indent << "Just entered OPT2 program (opt2_v2)" << endl;

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

  nproc = msg_->n();

  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;
  nshell = basis()->nshell();

  /* allocate storage for some arrays used for keeping track of which R   *
   * indices are processed by each node                                   */
  shellsize = (int*) malloc(nshell*sizeof(int));
  sorted_shells = (int*) malloc(nshell*sizeof(int));
  nbf = (int*) malloc(nproc*sizeof(int));
  proc = (int*) malloc(nshell*sizeof(int));


  /******************************************************
  * Begin distributing R shells between nodes so each   *
  * node gets ca. the same number of r basis functions  *
  ******************************************************/

  /* compute size of each shell */
  for (i=0; i<nshell; i++) {
    shellsize[i] = basis()->shell(i).nfunction();
    }

  /* do an index sort (large -> small) of shellsize to form sorted_shells */
  iquicksort(shellsize,sorted_shells,nshell);

  /* initialize nbf */
  for (i=0; i<nproc; i++) nbf[i] = 0;

  for (i=0; i<nshell; i++) {
    min = nbf[0];
    iproc = 0;
    for (j=1; j<nproc; j++) {
      if (nbf[j] < min) {
        iproc = j;
        min = nbf[j];
        }
      }
    proc[sorted_shells[i]] = iproc;
    nbf[iproc] += shellsize[sorted_shells[i]];
    }
  if (me == 0) {
    ExEnv::out0() << indent << "Distribution of basis functions between nodes:" << endl;
    for (i=0; i<nproc; i++) {
      if (i%12 == 0) ExEnv::out0() << indent;
      ExEnv::out0() << scprintf(" %4i",nbf[i]);
      if ((i+1)%12 == 0) ExEnv::out0() << endl;
      }
    ExEnv::out0() << endl;
    }

  /* determine which shells are to be processed by node me */
  nRshell = 0;
  for (i=0; i<nshell; i++) {
    if (proc[i] == me) nRshell++;
    }
  myshells = (int*) malloc(nRshell*sizeof(int));
  imyshell = 0;
  for (i=0; i<nshell; i++) {
    if (proc[i] == me) {
      myshells[imyshell] = i;
      imyshell++;
      }
    }

  /************************************************
  * End of distribution of R shells between nodes *
  ************************************************/


  /* compute batch size ni for opt2 loops;                                *
   * need to store the following arrays of type double : trans_int1-4,    *
   * trans_int4_tmp, scf_vector, evals_open, socc_sum, socc_sum_tmp,      *
   * mo_int_do_so_vir, mo_int_tmp,                                        *
   * and the following arrays of type int: myshells, shellsize,           *
   * sorted_shells, nbf, and proc                                         */
    
  size_t memused = 0;
  ni = 0;
  for (i=1; i<=nocc; i++) {
    distsize_t tmpmem = compute_v2_memory(i,
                                          nfuncmax, nbf[me], nshell,
                                          ndocc, nsocc, nvir, 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);

  int nbfmax = nbf[me];
  msg_->max(nbfmax);

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


  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_v2_memory(nocc,
                            nfuncmax, nbfmax, nshell,
                            ndocc, nsocc, nvir, nproc)
       << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Minimum memory required:        "
       << compute_v2_memory(1,
                            nfuncmax, nbfmax, nshell,
                            ndocc, nsocc, nvir, nproc)
       << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Batch size:                     " << ni
       << endl;

  if (ni < nsocc) {
    ExEnv::err0() << "Not enough memory allocated" << endl;
    abort();
    }

  if (ni < 1) {     /* this applies only to a closed shell case */
    ExEnv::err0() << "Not enough memory allocated" << endl;
    abort();
    }

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

  ExEnv::out0() << indent
       << "npass = " << npass
       << " rest = " << rest
       << endl;

  /* the scf vector might be distributed between the nodes, but for OPT2 *
   * each node needs its own copy of the vector;                         *
   * therefore, put a copy of the scf vector on each node;               * 
   * while doing this, duplicate columns corresponding to singly         *
   * occupied orbitals and order columns as [socc docc socc unocc]       */
  /* also rearrange scf eigenvalues as [socc docc socc unocc]            *
   * want socc first to get the socc's in the first batch                *
   * (need socc's to compute energy denominators - see                   *
   * socc_sum comment below)                                             */

  evals_open = (double*) malloc((noso+nsocc-nfzc-nfzv)*sizeof(double));

  RefDiagSCMatrix occ;
  RefDiagSCMatrix evals;
  RefSCMatrix Scf_Vec;
  eigen(evals, Scf_Vec, occ);