build2e.cc

来自「大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CH」· CC 代码 · 共 1,562 行 · 第 1/4 页

//
// build2e.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[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 <assert.h>
#include <math.h>

#include <scconfig.h>
#include <util/misc/formio.h>
#include <chemistry/qc/intv3/macros.h>
#include <chemistry/qc/intv3/fjt.h>
#include <chemistry/qc/intv3/utils.h>
#include <chemistry/qc/intv3/int2e.h>

using namespace std;
using namespace sc;

#define CHECK_STACK_ALIGNMENT 0
#if CHECK_STACK_ALIGNMENT
static void
stack_alignment_error(void *ptr, const char *where)
{
  ExEnv::outn() << "UNALIGNED STACK: " << where << ": " << ptr << endl;
}
static inline void
stack_alignment_check(void *ptr, const char *where)
{
  if ((unsigned)ptr & 7) stack_alignment_error(ptr,where);
}
#else
#  define stack_alignment_check(ptr,where)
#endif

  /* MG is the maximum angular momentum for which we will use
   * the generated build routines. It is defined in oint3/build.h */
#define MINA(x) (((x)<MG)?(x):MG)

static inline void 
iswtch(int *i, int *j)
{
  int tmp;

  tmp = *i;
  *i = *j;
  *j = tmp;
}

static void
fail()
{
  ExEnv::errn() << scprintf("failing module:\n%s",__FILE__) << endl;
  abort();
  }

/* This initializes the build routines.  It is called from
 * int_initialize_erep.  This allocates storage for the
 * intermediate integrals. */
void
Int2eV3::int_init_buildgc(int order,
                          int am1, int am2, int am3, int am4,
                          int nc1, int nc2, int nc3, int nc4)
{
  int *jmax_for_con;
  int am12;
  int am34;
  int am;
  int i,j,k,l,m;
  int ci,cj,ck,cl;
  int e0f0_con_int_bufsize;
  double *e0f0_con_int_buf;
  int int_v_bufsize, int_v0_bufsize;
  double *int_v_buf, *int_v0_buf;

  used_storage_build_ = 0;

  /* Convert the am1-4 to their canonical ordering. */
  if (am2>am1) {
    iswtch(&am1,&am2);
    iswtch(&nc1,&nc2);
    }
  if (am4>am3) {
    iswtch(&am3,&am4);
    iswtch(&nc3,&nc4);
    }
  if ((am3 > am1)||((am3 == am1)&&(am4 > am2))) {
    iswtch(&am1,&am3);
    iswtch(&nc1,&nc3);
    iswtch(&am2,&am4);
    iswtch(&nc2,&nc4);
    }

  /* If the center permutation 1<->3 and 2<->4 is performed, then
   * we may need the am for center 2 to be as big as for center 4. */
  if (am4 > am2) am2 = am4;

  /* As far as this routine knows the biggest nc can end up anywhere. */
  if (nc2>nc1) nc1 = nc2;
  if (nc3>nc1) nc1 = nc3;
  if (nc4>nc1) nc1 = nc4;
  nc2 = nc3 = nc4 = nc1;

  jmax_for_con = (int *) malloc(sizeof(int)*nc1);
  // storage for jmax_for_con is not counted since it is freed below
  for (i=0; i<nc1; i++) {
    int tmp;
    jmax_for_con[i] = bs1_->max_am_for_contraction(i);
    if (  (bs2_ != bs1_)
        &&((tmp=bs2_->max_am_for_contraction(i))>jmax_for_con[i]))
      jmax_for_con[i] = tmp;
    if (  (bs3_ != bs1_) && (bs3_ != bs2_)
        &&((tmp=bs3_->max_am_for_contraction(i))>jmax_for_con[i]))
      jmax_for_con[i] = tmp;
    if (  (bs4_ != bs1_) && (bs4_ != bs2_) && (bs4_ != bs3_)
        &&((tmp=bs4_->max_am_for_contraction(i))>jmax_for_con[i]))
      jmax_for_con[i] = tmp;
    }

  /* If derivatives are needed, then am1 can be bigger. */
  if (order==1) am1++;
  /* To compute derivative integral bounds, am3 can be bigger also. */
  if (order==1 && int_derivative_bounds) am3++;

  am12 = am1 + am2;
  am34 = am3 + am4;
  am = am12 + am34;

  /* Allocate the intlist. */
  contract_length.set_dim(am12+1,am34+1,am34+1);
  build.int_v_list.set_dim(am12+1,am34+1,am+1);
  used_storage_build_ += contract_length.nbyte();
  used_storage_build_ += build.int_v_list.nbyte();
#if CHECK_INTEGRAL_ALGORITHM
  ExEnv::outn() << "contract_length: " << contract_length.nbyte() << endl;
  ExEnv::outn() << "int_v_list: " << build.int_v_list.nbyte() << endl;
#endif

  /* Set all slots to 0 */
  for (i=0; i<=am12; i++) {
    for (j=0; j<=am34; j++) {
      for (k=0; k<=am12+am34; k++) {
        build.int_v_list(i,j,k) = 0;
        }
      }
    }

  for (i=0; i<=am12; i++) {
      for (j=0; j<=am34; j++) {
          for (k=0; k<=am34; k++) {
              contract_length(i,j,k) = 0;
              for (l=j; l<=k; l++) {
                  contract_length(i,j,k) += INT_NCART(i)*INT_NCART(l);
                }
            }
        }
    }

  /* Compute the size of the buffer for the primitive integrals. */
  int_v_bufsize = 0;
  int_v0_bufsize = 0;
  for (i=0; i<=am12; i++) {
    for (j=0; j<=am34; j++) {
      int_v0_bufsize += INT_NCART(i)*INT_NCART(j);
      for (k=0; k<=am12+am34-i-j; k++) {
        int_v_bufsize += INT_NCART(i)*INT_NCART(j);
        }
      }
    }

  int_v0_buf = (double*) malloc(sizeof(double)*int_v_bufsize);
  used_storage_build_ += sizeof(double)*int_v_bufsize;
  if (!int_v0_buf) {
    ExEnv::errn() << scprintf("couldn't allocate all integral intermediates\n");
    fail();
    }
  add_store(int_v0_buf);
  int_v_buf = &int_v0_buf[int_v0_bufsize];

  /* Allocate storage for the needed slots. */
  for (i=0; i<=am12; i++) {
    for (j=0; j<=am34; j++) {
      build.int_v_list(i,j,0) = int_v0_buf;
      int_v0_buf += INT_NCART(i)*INT_NCART(j);
      for (k=1; k<=am12+am34-i-j; k++) {
        build.int_v_list(i,j,k) = int_v_buf;
        int_v_buf += INT_NCART(i)*INT_NCART(j);
        }
      }
    }


  /* Allocate storage for the contracted integrals (these are the output
   * of the build routines). */
  /* The ci, etc, indices refer to which set of contraction
   * coefficients we are using. */
  e0f0_con_int_bufsize = 0;
  e0f0_con_ints_array = new IntV3Arraydoublep2***[nc1];
  used_storage_build_ += sizeof(IntV3Arraydoublep2***)*nc1;
  for (ci=0; ci<nc1; ci++) {
    e0f0_con_ints_array[ci] = new IntV3Arraydoublep2**[nc2];
    used_storage_build_ += sizeof(IntV3Arraydoublep2**)*nc2;
    for (cj=0; cj<nc2; cj++) {
      e0f0_con_ints_array[ci][cj] = new IntV3Arraydoublep2*[nc3];
      used_storage_build_ += sizeof(IntV3Arraydoublep2*)*nc3;
      for (ck=0; ck<nc3; ck++) {
        e0f0_con_ints_array[ci][cj][ck] = new IntV3Arraydoublep2[nc4];
        used_storage_build_ += sizeof(IntV3Arraydoublep2)*nc4;
        for (cl=0; cl<nc4; cl++) {
  int am12_for_con;
  int am34_for_con;

  am12_for_con = jmax_for_con[ci] + jmax_for_con[cj] + order;
  if ((jmax_for_con[ck]!=am3)||(jmax_for_con[cl]!=am4)) {
    am34_for_con = jmax_for_con[ck] + jmax_for_con[cl] + order;
    }
  else {
    am34_for_con = jmax_for_con[ck] + jmax_for_con[cl];
    }

#if CHECK_INTEGRAL_ALGORITHM
  ExEnv::outn() << "am12_for_con: " << am12_for_con << endl;
  ExEnv::outn() << "am34_for_con: " << am34_for_con << endl;
#endif

  e0f0_con_ints_array[ci][cj][ck][cl].set_dim(am12+1,am34+1);
  used_storage_build_ += e0f0_con_ints_array[ci][cj][ck][cl].nbyte();
#if CHECK_INTEGRAL_ALGORITHM
  ExEnv::outn() << "e0f0_con_ints_array: "
       << e0f0_con_ints_array[ci][cj][ck][cl].nbyte()
       << endl;
#endif

  /* Count how much storage for the integrals is needed. */
  for (i=0; i<=am12_for_con; i++) {
    for (k=0; k<=am34_for_con; k++) {
      int s =  INT_NCART(i)
               * INT_NCART(k);
      e0f0_con_int_bufsize += s;
      }
    }
          }
        }
      }
    }
  e0f0_con_int_buf = (double*) malloc(sizeof(double)*e0f0_con_int_bufsize);
  used_storage_build_ += e0f0_con_int_bufsize * sizeof(double);
#if CHECK_INTEGRAL_ALGORITHM
  ExEnv::outn() << "e0f0_int_buf: " << e0f0_con_int_bufsize * sizeof(double) << endl;
#endif
  if (!e0f0_con_int_buf) {
    ExEnv::errn() << scprintf("couldn't allocate contracted integral storage\n");
    fail();
    }
  add_store(e0f0_con_int_buf);
  /* Allocate storage for the integrals which will be used by the shift
   * routine. */
  for (ci=0; ci<nc1; ci++) {
    for (cj=0; cj<nc2; cj++) {
      for (ck=0; ck<nc3; ck++) {
        for (cl=0; cl<nc4; cl++) {
  int am12_for_con;
  int am34_for_con;

  am12_for_con = jmax_for_con[ci] + jmax_for_con[cj] + order;
  if ((jmax_for_con[ck]!=am3)||(jmax_for_con[cl]!=am4)) {
    am34_for_con = jmax_for_con[ck] + jmax_for_con[cl] + order;
    }
  else {
    am34_for_con = jmax_for_con[ck] + jmax_for_con[cl];
    }

  for (i=0; i<=am12; i++) {
    for (k=0; k<=am34; k++) {
      e0f0_con_ints_array[ci][cj][ck][cl](i,k) = 0;
      }
    }
  for (i=0; i<=am12_for_con; i++) {
    for (k=0; k<=am34_for_con; k++) {
/* If there are Pople style s=p shells and the shells are ordered
 * first s and then p and there are no p or d shells on the molecule,
 * then this algorithm would will allocate a little more storage
 * than needed.  General contraction should be ordered high to
 * low angular momentum for this reason. */
      e0f0_con_ints_array[ci][cj][ck][cl](i,k)
        = e0f0_con_int_buf;
      e0f0_con_int_buf +=  INT_NCART(i)
                           * INT_NCART(k);
      }
    }
          }
        }
      }
    }

  /* Initialize the build_routine array. */
  for (i=0; i<4; i++) {
    for (j=0; j<4; j++) {
      for (k=0; k<4; k++) {
        for (l=0; l<4; l++) {
          for (m=0; m<2; m++) {
    build_routine[i][j][k][l][m] = &BuildIntV3::impossible_integral;
            }
          }
        }
      }
    }

#define ASSIGN_BUILD(ii,j,k,l) \
  build_routine[ii][j][k][l][0]= &BuildIntV3::i ## ii ## j ## k ## l ;\
  build_routine[ii][j][k][l][1]= &BuildIntV3::i ## ii ## j ## k ## l ## eAB;
#if (MG == 1) || (MG == 2) || (MG == 3) || (MG == 4)
  ASSIGN_BUILD(0,1,0,0)
  ASSIGN_BUILD(0,1,0,1)
  ASSIGN_BUILD(0,1,1,1)
  ASSIGN_BUILD(1,1,0,0)
  ASSIGN_BUILD(1,1,1,1)
#endif

#if (MG == 2) || (MG == 3) || (MG == 4)
  ASSIGN_BUILD(0,2,0,0)
  ASSIGN_BUILD(0,2,0,1)
  ASSIGN_BUILD(0,2,0,2)
  ASSIGN_BUILD(0,2,1,1)
  ASSIGN_BUILD(0,2,1,2)
  ASSIGN_BUILD(0,2,2,2)
  ASSIGN_BUILD(1,2,0,0)
  ASSIGN_BUILD(1,2,0,1)
  ASSIGN_BUILD(1,2,1,1)
  ASSIGN_BUILD(1,2,1,2)
  ASSIGN_BUILD(1,2,2,2)
  ASSIGN_BUILD(2,2,0,0)
  ASSIGN_BUILD(2,2,0,1)
  ASSIGN_BUILD(2,2,1,1)
  ASSIGN_BUILD(2,2,2,2)
#endif

#if (MG == 3) || (MG == 4)
  ASSIGN_BUILD(0,3,0,0)
  ASSIGN_BUILD(0,3,0,1)
  ASSIGN_BUILD(0,3,0,2)
  ASSIGN_BUILD(0,3,0,3)
  ASSIGN_BUILD(0,3,1,1)
  ASSIGN_BUILD(0,3,1,2)
  ASSIGN_BUILD(0,3,1,3)
  ASSIGN_BUILD(0,3,2,2)
  ASSIGN_BUILD(0,3,2,3)
  ASSIGN_BUILD(0,3,3,3)
  ASSIGN_BUILD(1,3,0,0)
  ASSIGN_BUILD(1,3,0,1)
  ASSIGN_BUILD(1,3,0,2)
  ASSIGN_BUILD(1,3,1,1)
  ASSIGN_BUILD(1,3,1,2)
  ASSIGN_BUILD(1,3,1,3)
  ASSIGN_BUILD(1,3,2,2)
  ASSIGN_BUILD(1,3,2,3)
  ASSIGN_BUILD(1,3,3,3)
  ASSIGN_BUILD(2,3,0,0)
  ASSIGN_BUILD(2,3,0,1)
  ASSIGN_BUILD(2,3,0,2)
  ASSIGN_BUILD(2,3,1,1)
  ASSIGN_BUILD(2,3,1,2)
  ASSIGN_BUILD(2,3,2,2)
  ASSIGN_BUILD(2,3,2,3)
  ASSIGN_BUILD(2,3,3,3)
  ASSIGN_BUILD(3,3,0,0)
  ASSIGN_BUILD(3,3,0,1)
  ASSIGN_BUILD(3,3,0,2)
  ASSIGN_BUILD(3,3,1,1)