csgrad.cc

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

//
// csgrad.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 <stdlib.h>
#include <math.h>
#include <limits.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 <math/scmat/matrix.h>
#include <math/scmat/blocked.h>
#include <chemistry/molecule/molecule.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/mbpt/bzerofast.h>
#include <chemistry/qc/mbpt/mbpt.h>
#include <chemistry/qc/mbpt/util.h>
#include <chemistry/qc/mbpt/csgrade12.h>
#include <chemistry/qc/mbpt/csgrad34qb.h>
#include <chemistry/qc/mbpt/csgrads2pdm.h>

using namespace std;
using namespace sc;

#define SINGLE_THREAD_E12   0
#define SINGLE_THREAD_QBT34 0
#define SINGLE_THREAD_S2PDM 0

#define PRINT2Q 0
#define PRINT3Q 0
#define PRINT4Q 0
#if PRINT_BIGGEST_INTS
BiggestContribs biggest_ints_1(4,40);
#endif

#define WRITE_DOUBLES 0

static void sum_gradients(const Ref<MessageGrp>& msg, double **f, int n1, int n2);
static void zero_gradients(double **f, int n1, int n2);
static void accum_gradients(double **g, double **f, int n1, int n2);

#define PRINT1Q 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

void
MBPT2::compute_cs_grad()
{

  // New version of MP2 gradient program which uses the full
  // permutational symmetry of the two-electron integral derivatives

  Ref<SCMatrixKit> kit = basis()->matrixkit();

  int do_d2_ = 1;  // if true, compute d2 diagnostic

  int nij;        // number of i,j pairs on a node (for e.g., mo_int)
  double *mo_int; // MO integrals of type (ov|ov)
                  // (and these integrals divided by
                  // orbital energy denominators)
  double *integral_iqjs; // half-transformed integrals

  int nocc_act, nvir_act;
  int i, j, k;
  int ii, bb;
  int x, y;
  int a, b, c;
  int nshell;
  int offset;
  int ik_offset;
  int i_offset; 
  int npass, pass;
  int tmpint;
  int np, nq, nr, ns; 
  int P, Q, R, S;
  int p, q, r, s;
  int bf1, bf2, bf3, bf4;
  int index;
  int me;
  int nproc;
  int rest;
  int p_offset, q_offset, r_offset, s_offset;

  int aoint_computed = 0; 
  int aointder_computed = 0; 
  int xyz;
  int natom = molecule()->natom();     // the number of atoms
  int int_index;
  size_t mem_static;    // static memory in bytes
  int ij_proc;          // the processor which has ij pair
  int ij_index;         // of the ij pairs on a proc, this ij pair is number ij_index
                        // (i.e., ij_index < nij)
  int ik_proc;          // the processor which has ik pair
  int ik_index;
  int jloop, kloop;

  int ni;

  double *evals;              // scf eigenvalues
  double *iajb_ptr, *ibja_ptr, *iakb_ptr, *ibka_ptr;
  double *iajc_ptr, *ibjc_ptr, *icjb_ptr, *icja_ptr;
  double *ijkb_ptr, *ibkj_ptr;
  double pqrs;
  double *c_sa, c_rj;
  double *c_pi, *c_qi, *c_sj;
  double *c_qx, *c_qa, *c_sb, *c_pa, *c_pq, *c_sy;
  double delta_ijab, delta_ijbc, delta_ijac;
  double ecorr_mp2 = 0.0;
  double escf;
  double emp2=0.0;
  int tol;                    // log2 of the erep tolerance
                              // (erep < 2^tol => discard)
  double *Wkj=0,*Wab=0,*Waj=0;// occ-occ, vir-vir and vir-occ parts of 
                              // second order correction to MP2
                              // energy weighted density matrix
  double *Pkj=0,*Pab=0;       // occ-occ and vir-vir parts of second order
                              // correction to MP2 density matrix
  double *d2occ_mat, *d2vir_mat; // matrices for computation of D2 diagnostic
  double *Laj=0;              // MP2 Lagrangian
  double *Lpi;                // contrib to MP2 Lagrangian partially in AO basis
  double *pkj_ptr=0, *pab_ptr;
  double *d2occ_mat_ptr;
  double *d2vir_mat_ptr;
  double *wkj_ptr, *wjk_ptr, *wab_ptr, *wba_ptr, *waj_ptr=0;
  double *laj_ptr, *lpi_ptr, *lqi_ptr;
  double *gamma_iajs, *gamma_iajs_tmp; 
                              // partially back-transformed non-sep 2PDM's
  double *gamma_iqjs_tmp;
  double *gamma_iajs_ptr;
  double *gamma_iqjs_ptr;
  double *gammabuf;           // buffer used for sending elements of gamma_iqjs
  double *mo_intbuf;          // buffer used for sending mo integrals
  double tmpval, tmpval1;
  double *P2AO, *W2AO;
  double *p2ao_ptr, *w2ao_ptr;
  double *PHF, *WHF;
  double *phf_ptr, *whf_ptr;
  double *PMP2, *WMP2;
  double *pmp2_ptr, *wmp2_ptr;

  double *ixjs_tmp;      // three-quarter transformed two-el integrals
  double *integral_ixjs;  // all three-quarter transformed two-el integrals
  double *integral_iajy; // mo integrals (y = any MO)
  double *integral_ikja; // mo integrals
  double *integral_iqjs_ptr;
  double *iajy_ptr;
  double *ixjs_ptr;
  double *ikja_ptr;
  double *iajs_ptr, *ikjs_ptr;

  double **gradient=0, *gradient_dat=0;  // The MP2 gradient
  double **hf_gradient=0, *hf_gradient_dat=0;  // The HF gradient
  double **ginter=0;    // Intermediates for the MP2 gradient
  double **hf_ginter=0;    // Intermediates for the HF gradient
  double d2o, d2v, d2_diag;

  BiggestContribs biggest_coefs(5,10);
  CharacterTable ct = molecule()->point_group()->char_table();

#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

  int dograd = gradient_needed();

  tim_enter("mp2-mem");

  nfuncmax = basis()->max_nfunction_in_shell();

  nshell = basis()->nshell();

  me = msg_->me();

  if (me == 0) {
    ExEnv::out0() << endl << indent
         << "Entered memgrp based MP2 routine" << endl;
    }
  
  nproc = msg_->n();
  if (me == 0)
    ExEnv::out0() << indent << scprintf("nproc = %i", nproc) << endl;

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

  nocc = 0;
  for (i=0; i<oso_dimension()->n(); i++) {
    if (reference_->occupation(i) == 2.0) nocc++;
    }

  nocc_act = nocc - nfzc;
  nvir  = noso - nocc;
  nvir_act = nvir - nfzv;

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

  if (nocc_act <= 0) {
    if (me == 0) {
      ExEnv::err0() << "There are no active occupied orbitals; program exiting" << endl;
      }
    abort();
    }

  if (nvir_act <= 0) {
    if (me == 0) {
      ExEnv::err0() << "There are no active virtual orbitals; program exiting" << endl;
      }
    abort();
    }
    
  if (restart_orbital_memgrp_) {
    if (!dograd && !do_d1_ && !do_d2_) {
      ExEnv::out0() << indent
           << scprintf("Restarting at orbital %d with partial energy %18.14f",
                       restart_orbital_memgrp_, restart_ecorr_)
           << endl;
      ecorr_mp2 = restart_ecorr_;
      }
    else {
      ExEnv::out0() << indent
           << "Restart requested but not possible with gradients, D1, or D2"
           << endl;
      restart_ecorr_ = 0.0;
      restart_orbital_memgrp_ = 0;
      }
    }
  else {
      restart_ecorr_ = 0.0;
    }

  ////////////////////////////////////////////////////////
  // Compute batch size ni for mp2 loops;
  //
  // The following arrays are kept throughout (all of type double):
  //   scf_vector, gradient, ginter, Pkj, Pab, Wkj, Wab, Waj, Laj
  // and memory allocated for these arrays  and integral evaluators
  // is called mem_static
  //
  ////////////////////////////////////////////////////////
  if (me == 0) {
    mem_static = nbasis*noso; // scf vector
    mem_static += 2*nbasis*nfuncmax; // iqjs & iqjr
    if (dograd) {
      mem_static += 9*natom; // gradient & ginter & hf_ginter
      mem_static += (nocc*(nocc+1))/2; // Pkj
      mem_static += (nvir*(nvir+1))/2; // Pab
      mem_static += nocc*nocc; // Wkj
      mem_static += nvir*nvir; // Wab
      mem_static += 2*nocc*nvir; // Waj & Laj
      if (do_d2_) {
        mem_static += (nocc_act*(nocc_act+1))/2; // d2occ_mat
        mem_static += (nvir_act*(nvir_act+1))/2; // d2vir_mat
        }
      }
    else if (do_d1_) {
      mem_static += nocc*nvir; // partial Laj
      }
    mem_static *= sizeof(double);
    int nthreads = thr_->nthread();
    mem_static += nthreads * integral()->storage_required_eri(basis()); // integral evaluators
    ni = compute_cs_batchsize(mem_static, nocc_act-restart_orbital_memgrp_); 
    }

  if (max_norb_ > 0 && ni > max_norb_) {
      ExEnv::out0() << indent
           << "\"max_norb\" set: could have done "
           << ni << " orbitals per pass otherwise."
           << endl;
      ni = max_norb_;
    }

  // Send value of ni and mem_static to other nodes
  msg_->bcast(ni);
  double dmem_static = mem_static;
  msg_->bcast(dmem_static);
  mem_static = size_t(dmem_static);

  // Compute the storage to be used by the integral routines (required plus optional)
  size_t dyn_mem = distsize_to_size(compute_cs_dynamic_memory(ni,nocc_act));
  int mem_remaining;
  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_eri(basis());

  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_cs_dynamic_memory(nocc_act,nocc_act)+mem_static
       << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Minimum memory required:        "
       << compute_cs_dynamic_memory(1,nocc_act)+mem_static
       << " Bytes"
       << endl;
  ExEnv::out0() << indent
       << "Batch size:                     " << ni
       << endl;

  if (ni == 0) {
    ExEnv::err0() << "Batch size is 0: more memory or processors are needed"
         << endl;
    abort();
    }

  if (dynamic_) {
    ExEnv::out0() << indent << "Using dynamic load balancing." << endl;
    }

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

  if (me == 0) {
    ExEnv::out0() << indent