comp2e.cc

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

//
// comp2e.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 <stdarg.h>

#include <util/misc/formio.h>
#include <chemistry/qc/intv3/macros.h>
#include <chemistry/qc/intv3/flags.h>
#include <chemistry/qc/intv3/types.h>
#include <chemistry/qc/intv3/int2e.h>
#include <chemistry/qc/intv3/utils.h>
#include <chemistry/qc/intv3/tformv3.h>

using namespace std;
using namespace sc;

#undef DER_TIMING
#undef EREP_TIMING

#if defined(DER_TIMING)||defined(EREP_TIMING)
#  include <util/misc/timer.h>
#endif

static inline void
swtch(GaussianBasisSet* &i,GaussianBasisSet* &j)
{
  GaussianBasisSet *tmp;
  tmp = i;
  i = j;
  j = tmp;
}

static inline void
sswtch(GaussianShell**i,GaussianShell**j)
{
  GaussianShell*tmp;
  tmp = *i;
  *i = *j;
  *j = tmp;
}

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 computes the 2erep integrals for a shell quartet
 * specified by psh1, psh2, psh3, psh4.
 * The routine int_initialize_2erep must be called before
 * any integrals can be computed.
 * This routine may decide to change the shell ordering.
 * The new ordering is placed in *psh1,4 on exit.
 * for the derivatives.
 */
void
Int2eV3::erep(int &psh1, int &psh2, int &psh3, int &psh4)
{
  compute_erep(0,&psh1,&psh2,&psh3,&psh4,0,0,0,0);
  }

/* This is an alternate interface to int_erep.  It takes
 * as arguments the flags, an integer vector of shell numbers
 * and an integer vector which will be filled in with size
 * information, if it is non-NULL. */
void
Int2eV3::erep(int *shells, int  *sizes)
{
  erep(shells[0],shells[1],shells[2],shells[3]);
  if (sizes) {
    sizes[0] = bs1_->shell(shells[0]).nfunction();
    sizes[1] = bs2_->shell(shells[1]).nfunction();
    sizes[2] = bs3_->shell(shells[2]).nfunction();
    sizes[3] = bs4_->shell(shells[3]).nfunction();
    }
  }

/* If we need a computation with adjusted angular momentum, then
 * this lower level routine can be called instead of int_erep.
 * The dam{1,2,3,4} arguments given the amount by which the
 * angular momentum is to adjusted.  This differs from libint version
 * 1 in that it used total angular momentum here.
 */
void
Int2eV3::compute_erep(int flags, int *psh1, int *psh2, int *psh3, int *psh4,
                      int dam1, int dam2, int dam3, int dam4)
{
#ifdef EREP_TIMING
  char section[30];
#endif
  GaussianBasisSet *pbs1=bs1_.pointer();
  GaussianBasisSet *pbs2=bs2_.pointer();
  GaussianBasisSet *pbs3=bs3_.pointer();
  GaussianBasisSet *pbs4=bs4_.pointer();
  int size;
  int ii;
  int size1, size2, size3, size4;
  int tam1,tam2,tam3,tam4;
  int i,j,k,l;
  int ogc1,ogc2,ogc3,ogc4;
  int osh1,osh2,osh3,osh4;
  int sh1,sh2,sh3,sh4;
  int am1,am2,am3,am4,am12, am34;
  int minam1,minam2,minam3,minam4;
  int redundant_index;
  int e12,e13e24,e34;
  int p12,p34,p13p24;
  int eAB;

  /* Compute the offset shell numbers. */
  osh1 = *psh1 + bs1_shell_offset_;
  if (!int_unit2) osh2 = *psh2 + bs2_shell_offset_;
  osh3 = *psh3 + bs3_shell_offset_;
  if (!int_unit4) osh4 = *psh4 + bs4_shell_offset_;

  sh1 = *psh1;
  if (!int_unit2) sh2 = *psh2;
  sh3 = *psh3;
  if (!int_unit4) sh4 = *psh4;

  /* Test the arguments to make sure that they are sensible. */
  if (   sh1 < 0 || sh1 >= bs1_->nbasis()
      ||( !int_unit2 && (sh2 < 0 || sh2 >= bs2_->nbasis()))
      || sh3 < 0 || sh3 >= bs3_->nbasis()
      ||( !int_unit4 && (sh4 < 0 || sh4 >= bs4_->nbasis()))) {
    ExEnv::errn() << scprintf("compute_erep has been incorrectly used\n");
    ExEnv::errn() << scprintf("shells (bounds): %d (%d), %d (%d), %d (%d), %d (%d)\n",
            sh1,bs1_->nbasis()-1,
            sh2,bs2_->nbasis()-1,
            sh3,bs3_->nbasis()-1,
            sh4,bs4_->nbasis()-1);
    fail();
    }

  /* Set up pointers to the current shells. */
  int_shell1 = &bs1_->shell(sh1);
  if (!int_unit2) int_shell2 = &bs2_->shell(sh2);
  else int_shell2 = int_unit_shell;
  int_shell3 = &bs3_->shell(sh3);
  if (!int_unit4) int_shell4 = &bs4_->shell(sh4);
  else int_shell4 = int_unit_shell;


  /* Compute the maximum angular momentum on each centers to
   * determine the most efficient way to invoke the building and shifting
   * routines.  The minimum angular momentum will be computed at the
   * same time. */
  minam1 = int_shell1->min_am();
  minam2 = int_shell2->min_am();
  minam3 = int_shell3->min_am();
  minam4 = int_shell4->min_am();
  am1 = int_shell1->max_am();
  am2 = int_shell2->max_am();
  am3 = int_shell3->max_am();
  am4 = int_shell4->max_am();

  am1 += dam1; minam1 += dam1;
  am2 += dam2; minam2 += dam2;
  am3 += dam3; minam3 += dam3;
  am4 += dam4; minam4 += dam4;
  am12 = am1 + am2;
  am34 = am3 + am4;

  /* if no angular momentum remains on one of the centers return */
  if (am1 < 0 || am2 < 0 || am3 < 0 || am4 < 0) {
    return;
    }

#ifdef EREP_TIMING
  sprintf(section,"erep am=%02d",am12+am34);
  tim_enter(section);
  tim_enter("setup");
#endif

  /* Convert the integral to the most efficient form. */
  p12 = 0;
  p34 = 0;
  p13p24 = 0;

  if (am2 > am1) {
    p12 = 1;
    iswtch(&am1,&am2);iswtch(&sh1,&sh2);iswtch(psh1,psh2);iswtch(&osh1,&osh2);
    iswtch(&dam1,&dam2);
    iswtch(&minam1,&minam2);
    sswtch(&int_shell1,&int_shell2);
    swtch(pbs1,pbs2);
    }
  if (am4 > am3) {
    p34 = 1;
    iswtch(&am3,&am4);iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
    iswtch(&dam3,&dam4);
    iswtch(&minam3,&minam4);
    sswtch(&int_shell3,&int_shell4);
    swtch(pbs3,pbs4);
    }
  if (!(int_unit2||int_unit4) && (osh1 == osh4) && (osh2 == osh3) && (osh1 != osh2)) {
    /* Don't make the permutation unless we won't override what was
     * decided above about p34. */
    if (am4 == am3) {
      p34 = 1;
      iswtch(&am3,&am4);iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
      iswtch(&dam3,&dam4);
      iswtch(&minam3,&minam4);
      sswtch(&int_shell3,&int_shell4);
      swtch(pbs3,pbs4);
      }
    }
  if ((am34 > am12)||((am34 == am12)&&(minam1 > minam3))) {
    p13p24 = 1;
    iswtch(&am1,&am3);iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
    iswtch(&am2,&am4);iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
    iswtch(&int_unit2,&int_unit4);
    iswtch(&am12,&am34);
    iswtch(&dam1,&dam3);
    iswtch(&minam1,&minam3);
    sswtch(&int_shell1,&int_shell3);
    swtch(pbs1,pbs3);
    iswtch(&dam2,&dam4);
    iswtch(&minam2,&minam4);
    sswtch(&int_shell2,&int_shell4);
    swtch(pbs2,pbs4);
    }
  /* This tries to make centers A and B equivalent, if possible. */
  else if (  (am3 == am1)
           &&(am4 == am2)
           && !int_unit2
           && !int_unit4
           &&(minam1 == minam3)
           &&(!(  (bs1_ == bs2_)
                &&(bs1_->shell_to_center(sh1)==bs2_->shell_to_center(sh2))))
           &&(   (bs3_ == bs4_)
               &&(bs3_->shell_to_center(sh3)==bs4_->shell_to_center(sh4)))) {
    p13p24 = 1;
    iswtch(&am1,&am3);iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
    iswtch(&am2,&am4);iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
    iswtch(&am12,&am34);
    iswtch(&dam1,&dam3);
    iswtch(&minam1,&minam3);
    sswtch(&int_shell1,&int_shell3);
    swtch(pbs1,pbs3);
    iswtch(&dam2,&dam4);
    iswtch(&minam2,&minam4);
    sswtch(&int_shell2,&int_shell4);
    swtch(pbs2,pbs4);
    }

  if (  int_unit2
        ||((pbs1 == pbs2)
          &&(pbs1->shell_to_center(sh1)==pbs2->shell_to_center(sh2)))) {
    eAB = 1;
    }
  else {
    eAB = 0;
    }

  /* If the centers were permuted, then the int_expweighted variable may
   * need to be changed. */
  if (p12) {
    iswtch(&int_expweight1,&int_expweight2);
    }
  if (p34) {
    iswtch(&int_expweight3,&int_expweight4);
    }
  if (p13p24) {
    iswtch(&int_expweight1,&int_expweight3);
    iswtch(&int_expweight2,&int_expweight4);
    }

  int nc1 = int_shell1->ncontraction();
  int nc2 = int_shell2->ncontraction();
  int nc3 = int_shell3->ncontraction();
  int nc4 = int_shell4->ncontraction();

  /* Compute the shell sizes. */
  for (ii=size1=0; ii<nc1; ii++)
    size1 += INT_NCART(int_shell1->am(ii)+dam1);
  for (ii=size2=0; ii<nc2; ii++)
    size2 += INT_NCART(int_shell2->am(ii)+dam2);
  for (ii=size3=0; ii<nc3; ii++)
    size3 += INT_NCART(int_shell3->am(ii)+dam3);
  for (ii=size4=0; ii<nc4; ii++)
    size4 += INT_NCART(int_shell4->am(ii)+dam4);
  size = size1*size2*size3*size4;

  if (int_integral_storage) {
#ifdef EREP_TIMING
      tim_change("check storage");
#endif
    if (dam1 || dam2 || dam3 || dam4) {
      ExEnv::errn() << scprintf("cannot use integral storage and dam\n");
      fail();
      }
    if (    !int_unit2
         && !int_unit4
         && int_have_stored_integral(sh1,sh2,sh3,sh4,p12,p34,p13p24)) {
      goto post_computation;
      }
    }

  /* Buildam up on center 1 and 3. */
#ifdef EREP_TIMING
  tim_change("build");
#endif
  int_buildgcam(minam1,minam2,minam3,minam4,
                am1,am2,am3,am4,
                dam1,dam2,dam3,dam4,
                sh1,sh2,sh3,sh4, eAB);
#ifdef EREP_TIMING
  tim_change("cleanup");
#endif

  /* Begin loop over generalized contractions. */
  ogc1 = 0;
  for (i=0; i<nc1; i++) {
    tam1 = int_shell1->am(i) + dam1;
    if (tam1 < 0) continue;
    int tsize1 = INT_NCART_NN(tam1);
    ogc2 = 0;
    for (j=0; j<nc2; j++) {
      tam2 = int_shell2->am(j) + dam2;
      if (tam2 < 0) continue;
      int tsize2 = INT_NCART_NN(tam2);
      ogc3 = 0;
      for (k=0; k<nc3; k++) {
        tam3 = int_shell3->am(k) + dam3;
        if (tam3 < 0) continue;
        int tsize3 = INT_NCART_NN(tam3);
        ogc4 = 0;
        for (l=0; l<nc4; l++) {
          tam4 = int_shell4->am(l) + dam4;
          if (tam4 < 0) continue;
          int tsize4 = INT_NCART_NN(tam4);

#ifdef EREP_TIMING
  tim_change("shift");
#endif
  /* Shift angular momentum from 1 to 2 and from 3 to 4. */
  double *shiftbuffer = int_shiftgcam(i,j,k,l,tam1,tam2,tam3,tam4, eAB);
#ifdef EREP_TIMING
  tim_change("cleanup");
#endif

  /* Place the integrals in the integral buffer. */
  /* If permute_ is not set, then repack the integrals while copying. */
  if ((!permute_)&&(p12||p34||p13p24)) {
    int pam1,pam2,pam3,pam4;
    int psize234,psize34;
    int pogc1,pogc2,pogc3,pogc4;
    int psize1,psize2,psize3,psize4;
    pam1 = tam1;
    pam2 = tam2;
    pam3 = tam3;
    pam4 = tam4;
    pogc1 = ogc1;
    pogc2 = ogc2;
    pogc3 = ogc3;
    pogc4 = ogc4;
    psize1 = size1;
    psize2 = size2;
    psize3 = size3;
    psize4 = size4;
    if (p13p24) {
      iswtch(&pam1,&pam3);
      iswtch(&pam2,&pam4);
      iswtch(&pogc1,&pogc3);
      iswtch(&pogc2,&pogc4);
      iswtch(&psize1,&psize3);
      iswtch(&psize2,&psize4);
      }
    if (p34) {
      iswtch(&pam3,&pam4);
      iswtch(&pogc3,&pogc4);
      iswtch(&psize3,&psize4);
      }
    if (p12) {
      iswtch(&pam1,&pam2);
      iswtch(&pogc1,&pogc2);
      iswtch(&psize1,&psize2);
      }
    psize34 = psize4 * psize3;
    psize234 = psize34 * psize2;
    redundant_index = 0;
    int newindexoffset = pogc1*psize234 + pogc2*psize34 + pogc3*psize4 + pogc4;
    if (p13p24||p34) {
      int stride1=psize234;
      int stride2=psize34;
      int stride3=psize4;
      int stride4=1;
      int tmp;