shift2e.cc

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

//
// shift2e.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 <util/misc/formio.h>
#include <chemistry/qc/intv3/macros.h>
#include <chemistry/qc/intv3/int2e.h>

using namespace std;
using namespace sc;

//#undef CHECK_INTEGRAL_ALGORITHM
//#define CHECK_INTEGRAL_ALGORITHM 1

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

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

/* This initializes the shift routines.  It is called by int_initialize_erep.
 * It is passed the maximum am to be found on each center.
 */
void
Int2eV3::int_init_shiftgc(int order, int am1, int am2, int am3, int am4)
{
  /* The intermediate integral arrays are allocated by the
   * build initialization routine. */

  used_storage_shift_ = 0;

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

  /* 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;

  /* If derivatives are needed am1 will need to be larger. */
  if (order==1) am1++;
  /* For derivative integral bounds am3 will need to be larger. */
  if (order==1 && int_derivative_bounds) am3++;

  // Set up the new intermediate arrays.
  int e, c, d;
  int ndata34_e = 0;
  for (e=am1; e<=am1+am2; e++) {
    int size_e = INT_NCART(e);
    ndata34_e += size_e;
    }
  int ndata34_f = 0;
  for (d=1; d<=am4; d++) {
    int size_d = INT_NCART(d);
    int size_dm1 = INT_NCART(d-1);
    int off_cp1_dm1 = INT_NCART(am3) * size_dm1;
    int off_c_d = 0;
    for (c=am3; c<=am3+am4-d; c++) {
      int size_c = INT_NCART(c);
      int size_cp1 = INT_NCART(c+1);
      off_c_d += size_c * size_d;
      off_cp1_dm1 += size_cp1 * size_dm1;
      }
    if (off_c_d > ndata34_f) ndata34_f = off_c_d;
    if (off_cp1_dm1 > ndata34_f) ndata34_f = off_cp1_dm1;
    }
  int ndata34 = ndata34_e * ndata34_f;

  int ndata12 = 0;
  int a, b;
  int size_c_d = INT_NCART(am3)*INT_NCART(am4);
  for (b=1; b<=am2; b++) {
    int size_b = INT_NCART(b);
    int size_bm1 = INT_NCART(b-1);
    int off_a_b = 0;
    int off_ap1_bm1 = INT_NCART(am1) * size_bm1 * size_c_d;
    for (a=am1; a<=am1+am2-b; a++) {
      int size_a = INT_NCART(a);
      int size_ap1 = INT_NCART(a+1);
      off_a_b += size_a * size_b * size_c_d;
      off_ap1_bm1 += size_ap1 * size_bm1 * size_c_d;
      }
    if (off_a_b > ndata12) ndata12 = off_a_b;
    if (off_ap1_bm1 > ndata12) ndata12 = off_ap1_bm1;
    }
  int ndatamax = (ndata12>ndata34?ndata12:ndata34);
  buf34 = new double[ndata34];
  buf12 = new double[ndata12];
  bufshared = new double[ndatamax];

  used_storage_shift_ += sizeof(double)*(ndata34+ndata12+ndatamax);

  used_storage_ += used_storage_shift_;
  }

void
Int2eV3::int_done_shiftgc()
{
  used_storage_ -= used_storage_shift_;
  delete[] buf12;
  delete[] buf34;
  delete[] bufshared;
  }

/* This is the principle entry point for the am shifting routines.
 * tam1-4 is the target angular momentum on centers 1-4
 * sh1-4 are the shell numbers on centers 1-4
 */
double *
Int2eV3::int_shiftgcam(int gc1, int gc2, int gc3, int gc4,
                       int tam1, int tam2, int tam3, int tam4, int peAB)
{
  int am1,am2,am3,am4;

  /* Copy the gc{1,2,3,4} into g{1,2,3,4} (static globals). */
  g1 = gc1;
  g2 = gc2;
  g3 = gc3;
  g4 = gc4;

  /* Compute the angular momentum quartet. */
  am1 = tam1;
  am2 = tam2;
  am3 = tam3;
  am4 = tam4;

  // (a0|b0) does need shifting
  if (am2==0 && am4==0) {
    return e0f0_con_ints_array[g1][g2][g3][g4](am1,am3);
    }

  /* Copy the A B equivalency info into a static global variable. */
  eAB = peAB;

  /* Compute the intermediates. */
  AmB[0] =  build.int_v_r10 - build.int_v_r20;
  AmB[1] =  build.int_v_r11 - build.int_v_r21;
  AmB[2] =  build.int_v_r12 - build.int_v_r22;
  CmD[0] =  build.int_v_r30 - build.int_v_r40;
  CmD[1] =  build.int_v_r31 - build.int_v_r41;
  CmD[2] =  build.int_v_r32 - build.int_v_r42;

#if CHECK_INTEGRAL_ALGORITHM > 1
  ExEnv::outn() << "generating ("
       << am1 << "," << am2 << "," << am3 << "," << am4 << ")"
       << ":" << endl;
#endif

  // the (e0|f0) integrals have been initialized
  IntV3Arraydoublep2 &e0f0 = e0f0_con_ints_array[g1][g2][g3][g4];

  // generate (e0|cd) for each needed e
  int e, c, d;
  int off_e = 0;
  int size34 = INT_NCART(am3)*INT_NCART(am4);
  double *buf34_1 = buf34;
  double *buf34_2 = bufshared;
  for (e=am1; e<=am1+am2; e++) {
    int size_e = INT_NCART(e);
    for (d=1; d<=am4; d++) {
      int size_d = INT_NCART(d);
      int size_dm1 = INT_NCART(d-1);
      int off_c_dm1 = 0;
      int off_cp1_dm1 = size_e * INT_NCART(am3) * size_dm1;
      int off_c_d = 0;
      for (c=am3; c<=am3+am4-d; c++) {
        int size_c = INT_NCART(c);
        int size_cp1 = INT_NCART(c+1);
        double *I0001, *I0010, *I0000;
        if (d==am4) {
          I0001 = &buf12[off_e];
          }
        else I0001 = &buf34_1[off_c_d];
        if (d==1) {
          I0010 = e0f0(e,c+1);
          I0000 = e0f0(e,c);
          }
        else {
          I0010 = &buf34_2[off_cp1_dm1];
          I0000 = &buf34_2[off_c_dm1];
          }
        shiftam_34(I0001,I0010,I0000,e,0,c,d);
        off_c_d += size_e * size_c * size_d;
        off_c_dm1 = off_cp1_dm1;
        off_cp1_dm1 += size_e * size_cp1 * size_dm1;
        }
      // swap the buffers.
      double *tmp = buf34_1;
      buf34_1 = buf34_2;
      buf34_2 = tmp;
      }
    off_e += size_e * size34;
    }

  // generate (ab|cd)
  int a, b;
  int size_c_d = size34;
  double *buf12_1 = bufshared;
  double *buf12_2 = buf12;
  for (b=1; b<=am2; b++) {
    int size_b = INT_NCART(b);
    int size_bm1 = INT_NCART(b-1);
    int off_a_b = 0;
    int off_ap1_bm1 = INT_NCART(am1) * size_bm1 * size_c_d;
    int off_a_bm1 = 0;
    for (a=am1; a<=am1+am2-b; a++) {
      int size_a = INT_NCART(a);
      int size_ap1 = INT_NCART(a+1);
      double *I0100 = &buf12_1[off_a_b];
      double *I1000;
      double *I0000;
      if (b==1 && am4 == 0) {
        I1000 = e0f0(a+1,am3);
        if (eAB) I0000 = 0;
        else I0000 = e0f0(a,am3);
        }
      else {
        I1000 = &buf12_2[off_ap1_bm1];
        if (eAB) I0000 = 0;
        else I0000 = &buf12_2[off_a_bm1];
        }
      if (eAB) shiftam_12eAB(I0100,I1000,I0000,a,b,am3,am4);
      else shiftam_12(I0100,I1000,I0000,a,b,am3,am4);
      off_a_b += size_a * size_b * size_c_d;
      off_a_bm1 = off_ap1_bm1;
      off_ap1_bm1 += size_ap1 * size_bm1 * size_c_d;