comp1e.cc

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

//
// comp1e.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 <math.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/int1e.h>
#include <chemistry/qc/intv3/tformv3.h>

using namespace std;
using namespace sc;

#define IN(i,j) ((i)==(j)?1:0)
#define SELECT(x1,x2,x3,s) (((s)==0)?x1:(((s)==1)?(x2):(x3)))

#define DEBUG_NUC_SHELL_DER 0
#define DEBUG_NUC_PRIM 0
#define DEBUG_NUC_SHELL 0

#define DEBUG_EFIELD_PRIM 0

/* ------------ Initialization of 1e routines. ------------------- */
/* This routine returns a buffer large enough to hold a shell doublet
 * of integrals (if order == 0) or derivative integrals (if order == 1).
 */
void
Int1eV3::int_initialize_1e(int flags, int order)
{
  int jmax1,jmax2,jmax;
  int scratchsize=0,nshell2;

  /* The efield routines look like derivatives so bump up order if
   * it is zero to allow efield integrals to be computed.
   */
  if (order == 0) order = 1;

  jmax1 = bs1_->max_angular_momentum();
  jmax2 = bs2_->max_angular_momentum();
  jmax = jmax1 + jmax2;

  fjt_ = new FJT(jmax + 2*order);

  nshell2 = bs1_->max_ncartesian_in_shell()*bs2_->max_ncartesian_in_shell();

  if (order == 0) {
    init_order = 0;
    scratchsize = nshell2;
    }
  else if (order == 1) {
    init_order = 1;
    scratchsize = nshell2*3;
    }
  else {
    ExEnv::errn() << scprintf("int_initialize_1e: invalid order: %d\n",order);
    exit(1);
    }

  buff = new double[scratchsize];
  cartesianbuffer = new double[scratchsize];
  cartesianbuffer_scratch = new double[scratchsize];

  inter.set_dim(jmax1+order+1,jmax2+order+1,jmax+2*order+1);
  efield_inter.set_dim(jmax1+order+1,jmax2+order+1,jmax+2*order+1);
  int i,j,m;
  for (i=0; i<=jmax1+order; i++) {
    int sizei = INT_NCART_NN(i);
    for (j=0; j<=jmax2+order; j++) {
      int sizej = INT_NCART_NN(j);
      for (m=0; m<=jmax+2*order-i-j; m++) {
        inter(i,j,m) = new double[sizei*sizej];
        efield_inter(i,j,m) = new double[sizei*sizej*3];
        }
      for (; m<=jmax+2*order; m++) {
        inter(i,j,m) = 0;
        efield_inter(i,j,m) = 0;
        }
      }
    }

  }

void
Int1eV3::int_done_1e()
{
  init_order = -1;
  delete[] buff;
  delete[] cartesianbuffer;
  delete[] cartesianbuffer_scratch;
  buff = 0;
  cartesianbuffer = 0;
  inter.delete_data();
  efield_inter.delete_data();
}


/* --------------------------------------------------------------- */
/* ------------- Routines for the overlap matrix ----------------- */
/* --------------------------------------------------------------- */

/* This computes the overlap integrals between functions in two shells.
 * The result is placed in the buffer.
 */
void
Int1eV3::overlap(int ish, int jsh)
{
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int gc1,gc2;
  int index,index1,index2;

  c1 = bs1_->shell_to_center(ish);
  c2 = bs2_->shell_to_center(jsh);
  for (int xyz=0; xyz<3; xyz++) {
      A[xyz] = bs1_->r(c1,xyz);
      B[xyz] = bs2_->r(c2,xyz);
    }
  gshell1 = &bs1_->shell(ish);
  gshell2 = &bs2_->shell(jsh);
  index = 0;
  FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
      cartesianbuffer[index] = comp_shell_overlap(gc1,i1,j1,k1,gc2,i2,j2,k2);
      index++;
      END_FOR_GCCART_GS(index2)
    END_FOR_GCCART_GS(index1)

  transform_1e(integral_, cartesianbuffer, buff, gshell1, gshell2);
  }

/* This computes the overlap ints between functions in two shells.
 * The result is placed in the buffer.
 */
void
Int1eV3::overlap_1der(int ish, int jsh,
                      int idercs, int centernum)
{
  int i;
  int c1,c2;
  int ni,nj;

  if (!(init_order >= 0)) {
    ExEnv::errn() << scprintf("int_shell_overlap: one electron routines are not init'ed\n");
    exit(1);
    }

  Ref<GaussianBasisSet> dercs;
  if (idercs == 0) dercs = bs1_;
  else dercs = bs2_;

  c1 = bs1_->shell_to_center(ish);
  c2 = bs2_->shell_to_center(jsh);
  for (int xyz=0; xyz<3; xyz++) {
      A[xyz] = bs1_->r(c1,xyz);
      B[xyz] = bs2_->r(c2,xyz);
    }
  gshell1 = &bs1_->shell(ish);
  gshell2 = &bs2_->shell(jsh);

  ni = gshell1->nfunction();
  nj = gshell2->nfunction();

#if 0
  ExEnv::outn() << scprintf("zeroing %d*%d*3 elements of buff\n",ni,nj);
#endif
  for (i=0; i<ni*nj*3; i++) {
    buff[i] = 0.0;
    }

  int_accum_shell_overlap_1der(ish,jsh,dercs,centernum);
  }

/* This computes the overlap derivative integrals between functions
 * in two shells.  The result is accumulated in the buffer which is ordered
 * atom 0 x, y, z, atom 1, ... .
 */
void
Int1eV3::int_accum_shell_overlap_1der(int ish, int jsh,
                                      Ref<GaussianBasisSet> dercs, int centernum)
{
  accum_shell_1der(buff,ish,jsh,dercs,centernum,&Int1eV3::comp_shell_overlap);
  }

/* Compute the overlap for the shell.  The arguments are the
 * cartesian exponents for centers 1 and 2.  The shell1 and shell2
 * globals are used. */
double
Int1eV3::comp_shell_overlap(int gc1, int i1, int j1, int k1,
                            int gc2, int i2, int j2, int k2)
{
  double exp1,exp2;
  int i,j,xyz;
  double result;
  double Pi;
  double oozeta;
  double AmB,AmB2;
  double tmp;

  if ((i1<0)||(j1<0)||(k1<0)||(i2<0)||(j2<0)||(k2<0)) return 0.0;

  /* Loop over the primitives in the shells. */
  result = 0.0;
  for (i=0; i<gshell1->nprimitive(); i++) {
    for (j=0; j<gshell2->nprimitive(); j++) {

      /* Compute the intermediates. */
      exp1 = gshell1->exponent(i);
      exp2 = gshell2->exponent(j);
      oozeta = 1.0/(exp1 + exp2);
      oo2zeta = 0.5*oozeta;
      AmB2 = 0.0;
      for (xyz=0; xyz<3; xyz++) {
        Pi = oozeta*(exp1 * A[xyz] + exp2 * B[xyz]);
        PmA[xyz] = Pi - A[xyz];
        PmB[xyz] = Pi - B[xyz];
        AmB = A[xyz] - B[xyz];
        AmB2 += AmB*AmB;
        }
      ss =   pow(3.141592653589793/(exp1+exp2),1.5)
           * exp(- oozeta * exp1 * exp2 * AmB2);
      tmp     =  gshell1->coefficient_unnorm(gc1,i)
               * gshell2->coefficient_unnorm(gc2,j)
               * comp_prim_overlap(i1,j1,k1,i2,j2,k2);
      if (exponent_weighted == 0) tmp *= exp1;
      else if (exponent_weighted == 1) tmp *= exp2;
      result += tmp;
      }
    }

  return result;
  }

/* Compute the overlap between two primitive functions. */
#if 0
double
Int1eV3::int_prim_overlap(shell_t *pshell1, shell_t *pshell2,
                          double *pA, double *pB,
                          int prim1, int prim2,
                          int i1, int j1, int k1,
                          int i2, int j2, int k2)
{
  int xyz;
  double Pi;
  double oozeta;
  double AmB,AmB2;

  /* Compute the intermediates. */
  oozeta = 1.0/(gshell1->exponent(prim1) + gshell2->exponent(prim2));
  oo2zeta = 0.5*oozeta;
  AmB2 = 0.0;
  for (xyz=0; xyz<3; xyz++) {
    Pi = oozeta*(gshell1->exponent(prim1) * A[xyz]
                 + gshell2->exponent(prim2) * B[xyz]);
    PmA[xyz] = Pi - A[xyz];
    PmB[xyz] = Pi - B[xyz];
    AmB = A[xyz] - B[xyz];
    AmB2 += AmB*AmB;
    }
  ss =   pow(3.141592653589793/(gshell1->exponent(prim1)
                                +gshell2->exponent(prim2)),1.5)
       * exp(- oozeta * gshell1->exponent(prim1)
             * gshell2->exponent(prim2) * AmB2);
  return comp_prim_overlap(i1,j1,k1,i2,j2,k2);
  }
#endif

double
Int1eV3::comp_prim_overlap(int i1, int j1, int k1,
                         int i2, int j2, int k2)
{
  double result;

  if (i1) {
    result = PmA[0] * comp_prim_overlap(i1-1,j1,k1,i2,j2,k2);
    if (i1>1) result += oo2zeta*(i1-1) * comp_prim_overlap(i1-2,j1,k1,i2,j2,k2);
    if (i2>0) result += oo2zeta*i2 * comp_prim_overlap(i1-1,j1,k1,i2-1,j2,k2);
    return result;
    }
  if (j1) {
    result = PmA[1] * comp_prim_overlap(i1,j1-1,k1,i2,j2,k2);
    if (j1>1) result += oo2zeta*(j1-1) * comp_prim_overlap(i1,j1-2,k1,i2,j2,k2);
    if (j2>0) result += oo2zeta*j2 * comp_prim_overlap(i1,j1-1,k1,i2,j2-1,k2);
    return result;
    }
  if (k1) {
    result = PmA[2] * comp_prim_overlap(i1,j1,k1-1,i2,j2,k2);
    if (k1>1) result += oo2zeta*(k1-1) * comp_prim_overlap(i1,j1,k1-2,i2,j2,k2);
    if (k2>0) result += oo2zeta*k2 * comp_prim_overlap(i1,j1,k1-1,i2,j2,k2-1);
    return result;
    }

  if (i2) {
    result = PmB[0] * comp_prim_overlap(i1,j1,k1,i2-1,j2,k2);
    if (i2>1) result += oo2zeta*(i2-1) * comp_prim_overlap(i1,j1,k1,i2-2,j2,k2);
    if (i1>0) result += oo2zeta*i1 * comp_prim_overlap(i1-1,j1,k1,i2-1,j2,k2);
    return result;
    }
  if (j2) {
    result = PmB[1] * comp_prim_overlap(i1,j1,k1,i2,j2-1,k2);
    if (j2>1) result += oo2zeta*(j2-1) * comp_prim_overlap(i1,j1,k1,i2,j2-2,k2);
    if (j1>0) result += oo2zeta*j1 * comp_prim_overlap(i1,j1-1,k1,i2,j2-1,k2);
    return result;
    }
  if (k2) {
    result = PmB[2] * comp_prim_overlap(i1,j1,k1,i2,j2,k2-1);
    if (k2>1) result += oo2zeta*(k2-1) * comp_prim_overlap(i1,j1,k1,i2,j2,k2-2);
    if (k1>0) result += oo2zeta*k1 * comp_prim_overlap(i1,j1,k1-1,i2,j2,k2-1);
    return result;
    }

  return ss;
  }

/* --------------------------------------------------------------- */
/* ------------- Routines for the kinetic energy ----------------- */
/* --------------------------------------------------------------- */

/* This computes the kinetic energy integrals between functions in two shells.
 * The result is placed in the buffer.
 */
void
Int1eV3::kinetic(int ish, int jsh)
{
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int cart1,cart2;
  int index;
  int gc1,gc2;

  c1 = bs1_->shell_to_center(ish);
  c2 = bs2_->shell_to_center(jsh);
  for (int xyz=0; xyz<3; xyz++) {
      A[xyz] = bs1_->r(c1,xyz);
      B[xyz] = bs2_->r(c2,xyz);
    }
  gshell1 = &bs1_->shell(ish);
  gshell2 = &bs2_->shell(jsh);
  index = 0;
  FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
      cartesianbuffer[index] = comp_shell_kinetic(gc1,i1,j1,k1,gc2,i2,j2,k2);
      index++;
      END_FOR_GCCART_GS(cart2)
    END_FOR_GCCART_GS(cart1)

  transform_1e(integral_, cartesianbuffer, buff, gshell1, gshell2);
  }

void
Int1eV3::int_accum_shell_kinetic(int ish, int jsh)
{
  int c1,i1,j1,k1,c2,i2,j2,k2;
  int cart1,cart2;
  int index;
  int gc1,gc2;

  c1 = bs1_->shell_to_center(ish);
  c2 = bs2_->shell_to_center(jsh);
  for (int xyz=0; xyz<3; xyz++) {
      A[xyz] = bs1_->r(c1,xyz);
      B[xyz] = bs2_->r(c2,xyz);
    }
  gshell1 = &bs1_->shell(ish);
  gshell2 = &bs2_->shell(jsh);
  index = 0;

  FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
    FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
      cartesianbuffer[index] = comp_shell_kinetic(gc1,i1,j1,k1,gc2,i2,j2,k2);
      index++;
      END_FOR_GCCART_GS(cart2)
    END_FOR_GCCART_GS(cart1)
  accum_transform_1e(integral_, cartesianbuffer, buff, gshell1, gshell2);
  }

/* This computes the kinetic energy derivative integrals between functions
 * in two shells.  The result is accumulated in the buffer which is ordered
 * atom 0 x, y, z, atom 1, ... .
 */
void
Int1eV3::int_accum_shell_kinetic_1der(int ish, int jsh,
                                      Ref<GaussianBasisSet> dercs, int centernum)
{
  accum_shell_1der(buff,ish,jsh,dercs,centernum,&Int1eV3::comp_shell_kinetic);
  }

/* This computes the basis function part of 
 * generic derivative integrals between functions
 * in two shells.  The result is accumulated in the buffer which is ordered
 * atom 0 x, y, z, atom 1, ... .
 * The function used to compute the nonderivative integrals is shell_function.
 */
void
Int1eV3::accum_shell_1der(double *buff, int ish, int jsh,
                          Ref<GaussianBasisSet> dercs, int centernum,
                          double (Int1eV3::*shell_function)
                          (int,int,int,int,int,int,int,int))