tformv3.cc

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

//
// tformv3.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.
//

#if defined(__GNUC__)
#pragma implementation
#endif

#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <util/misc/formio.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/intv3/macros.h>
#include <chemistry/qc/intv3/tformv3.h>
#include <chemistry/qc/intv3/utils.h>

using namespace sc;

////////////////////////////////////////////////////////////////////////////

#define PRINT 0

void
Int2eV3::transform_init()
{
  source = 0;
  nsourcemax = 0;
}

void
Int2eV3::transform_done()
{
  delete[] source;
}

void
Int1eV3::transform_init()
{
  source = 0;
  nsourcemax = 0;
}

void
Int1eV3::transform_done()
{
  delete[] source;
}

static void
do_copy1(double *source, double *target, int chunk,
         int n1, int s1, int offset1,
         int n2, int s2, int offset2)
{
  int i1, i2;

  for (i1=0; i1<n1; i1++) {
      int off = ((offset1 + i1)*s2 + offset2)*chunk;
      for (i2=0; i2<n2*chunk; i2++, off++) {
          target[off] = source[off];
        }
    }
}

static void
do_copy2(double *source, double *target,
         int n1, int s1, int offset1,
         int n2, int s2, int offset2,
         int n3, int s3, int offset3,
         int n4, int s4, int offset4)
{
  int i1, i2, i3, i4;

  for (i1=0; i1<n1; i1++) {
      for (i2=0; i2<n2; i2++) {
          for (i3=0; i3<n3; i3++) {
              int off = (((offset1 + i1)*s2 + offset2 + i2)
                         *s3 + offset3 + i3)*s4 + offset4;
              for (i4=0; i4<n4; i4++, off++) {
                  target[off] = source[off];
                }
            }
        }
    }
}

static void
do_sparse_transform11(double *source, double *target, int chunk,
                      SphericalTransformIter& trans,
                      int offsetcart1,
                      int offsetpure1,
                      int n2, int s2, int offset2)
{
  int i2;

  for (trans.begin(); trans.ready(); trans.next()) {
      double coef = trans.coef();
      int pure = trans.pureindex();
      int cart = trans.cartindex();
      int offtarget = ((offsetpure1 + pure)*s2 + offset2)*chunk;
      int offsource = ((offsetcart1 + cart)*s2 + offset2)*chunk;
      for (i2=0; i2<n2*chunk; i2++) {
          target[offtarget++] += coef * source[offsource++];
        }
    }
}

static void
do_sparse_transform12(double *source, double *target, int chunk,
                      SphericalTransformIter& trans,
                      int n1, int offset1,
                      int s2cart, int offsetcart2,
                      int s2pure, int offsetpure2)
{
  int i1, ichunk;

  for (trans.begin(); trans.ready(); trans.next()) {
      double coef = trans.coef();
      int pure = trans.pureindex();
      int cart = trans.cartindex();
      for (i1=0; i1<n1; i1++) {
          int offtarget = ((offset1 + i1)*s2pure + offsetpure2 + pure)*chunk;
          int offsource = ((offset1 + i1)*s2cart + offsetcart2 + cart)*chunk; 
          for (ichunk=0; ichunk<chunk; ichunk++) {
              target[offtarget++] += coef * source[offsource++];
            }
        }
    }
}

static void
do_sparse_transform2_1(double *source, double *target,
                       SphericalTransformIter& trans,
                       int stcart, int stpure,
                       int ogctcart, int ogctpure,
                       int n2, int s2, int ogc2,
                       int n3, int s3, int ogc3,
                       int n4, int s4, int ogc4)
{
  int i2, i3, i4;

  for (trans.begin(); trans.ready(); trans.next()) {
      double coef = trans.coef();
      int pure = trans.pureindex();
      int cart = trans.cartindex();
      int offtarget1 = ogctpure + pure;
      int offsource1 = ogctcart + cart;
      int offtarget2 = offtarget1*s2 + ogc2;
      int offsource2 = offsource1*s2 + ogc2;
      for (i2=0; i2<n2; i2++,offtarget2++,offsource2++) {
          int offtarget3 = offtarget2*s3 + ogc3;
          int offsource3 = offsource2*s3 + ogc3;
          for (i3=0; i3<n3; i3++,offtarget3++,offsource3++) {
              int offtarget4 = offtarget3*s4 + ogc4;
              int offsource4 = offsource3*s4 + ogc4;
              for (i4=0; i4<n4; i4++,offtarget4++,offsource4++) {
                  target[offtarget4] += coef * source[offsource4];
                }
            }
        }
    }
}

static void
do_sparse_transform2_2(double *source, double *target,
                       SphericalTransformIter& trans,
                       int stcart, int stpure,
                       int ogctcart, int ogctpure,
                       int n1, int s1, int ogc1,
                       int n3, int s3, int ogc3,
                       int n4, int s4, int ogc4)
{
  int i1, i3, i4;

  for (trans.begin(); trans.ready(); trans.next()) {
      double coef = trans.coef();
      int pure = trans.pureindex();
      int cart = trans.cartindex();
      int offtarget1 = ogc1;
      int offsource1 = ogc1;
      for (i1=0; i1<n1; i1++,offtarget1++,offsource1++) {
          int offtarget2 = offtarget1*stpure + ogctpure + pure;
          int offsource2 = offsource1*stcart + ogctcart + cart;
          int offtarget3 = offtarget2*s3 + ogc3;
          int offsource3 = offsource2*s3 + ogc3;
          for (i3=0; i3<n3; i3++,offtarget3++,offsource3++) {
              int offtarget4 = offtarget3*s4 + ogc4;
              int offsource4 = offsource3*s4 + ogc4;
              for (i4=0; i4<n4; i4++,offtarget4++,offsource4++) {
                  target[offtarget4] += coef * source[offsource4];
                }
            }
        }
    }
}

static void
do_sparse_transform2_3(double *source, double *target,
                       SphericalTransformIter& trans,
                       int stcart, int stpure,
                       int ogctcart, int ogctpure,
                       int n1, int s1, int ogc1,
                       int n2, int s2, int ogc2,
                       int n4, int s4, int ogc4)
{
  int i1, i2, i4;

  for (trans.begin(); trans.ready(); trans.next()) {
      double coef = trans.coef();
      int pure = trans.pureindex();
      int cart = trans.cartindex();
      int offtarget1 = ogc1;
      int offsource1 = ogc1;
      for (i1=0; i1<n1; i1++,offtarget1++,offsource1++) {
          int offtarget2 = offtarget1*s2 + ogc2;
          int offsource2 = offsource1*s2 + ogc2;
          for (i2=0; i2<n2; i2++,offtarget2++,offsource2++) {
              int offtarget3 = offtarget2*stpure + ogctpure + pure;
              int offsource3 = offsource2*stcart + ogctcart + cart;
              int offtarget4 = offtarget3*s4 + ogc4;
              int offsource4 = offsource3*s4 + ogc4;
              for (i4=0; i4<n4; i4++,offtarget4++,offsource4++) {
                  target[offtarget4] += coef * source[offsource4];
                }
            }
        }
    }
}

static void
do_sparse_transform2_4(double *source, double *target,
                       SphericalTransformIter& trans,
                       int stcart, int stpure,
                       int ogctcart, int ogctpure,
                       int n1, int s1, int ogc1,
                       int n2, int s2, int ogc2,
                       int n3, int s3, int ogc3)
{
  int i1, i2, i3;

  for (trans.begin(); trans.ready(); trans.next()) {
      double coef = trans.coef();
      int pure = trans.pureindex();
      int cart = trans.cartindex();
      int offtarget1 = ogc1;
      int offsource1 = ogc1;
      for (i1=0; i1<n1; i1++,offtarget1++,offsource1++) {
          int offtarget2 = offtarget1*s2 + ogc2;
          int offsource2 = offsource1*s2 + ogc2;
          for (i2=0; i2<n2; i2++,offtarget2++,offsource2++) {
              int offtarget3 = offtarget2*s3 + ogc3;
              int offsource3 = offsource2*s3 + ogc3;
              int offtarget4 = offtarget3*stpure + ogctpure + pure;
              int offsource4 = offsource3*stcart + ogctcart + cart;
              for (i3=0; i3<n3; i3++,offtarget4+=stpure,offsource4+=stcart) {
                  //for (i3=0; i3<n3; i3++,offtarget3++,offsource3++) {
                  //int offtarget4 = offtarget3*stpure + ogctpure + pure;
                  //int offsource4 = offsource3*stcart + ogctcart + cart;
                  target[offtarget4] += coef * source[offsource4];
                }
            }
        }
    }
}

static void
do_sparse_transform2(double *source, double *target,
                     int index, SphericalTransformIter& trans,
                     int stcart, int stpure,
                     int ogctcart, int ogctpure,
                     int n1, int s1, int ogc1,
                     int n2, int s2, int ogc2,
                     int n3, int s3, int ogc3,
                     int n4, int s4, int ogc4)
{
  switch (index) {
  case 0:
      do_sparse_transform2_1(source, target, trans,
                             stcart,stpure,
                             ogctcart, ogctpure,
                             n2, s2, ogc2,
                             n3, s3, ogc3,
                             n4, s4, ogc4);
      break;
  case 1:
      do_sparse_transform2_2(source, target, trans,
                             stcart,stpure,
                             ogctcart, ogctpure,
                             n1, s1, ogc1,
                             n3, s3, ogc3,
                             n4, s4, ogc4);
      break;
  case 2:
      do_sparse_transform2_3(source, target, trans,
                             stcart,stpure,
                             ogctcart, ogctpure,
                             n1, s1, ogc1,
                             n2, s2, ogc2,
                             n4, s4, ogc4);
      break;
  case 3:
      do_sparse_transform2_4(source, target, trans,
                             stcart,stpure,
                             ogctcart, ogctpure,
                             n1, s1, ogc1,
                             n2, s2, ogc2,
                             n3, s3, ogc3);
      break;
    }
}

/* make sure enough space exists for the source integrals */
void
Int2eV3::source_space(int nsource)
{
  if (nsourcemax < nsource) {
      delete[] source;
      source = new double[nsource*3];
      target = &source[nsource];
      scratch = &source[nsource*2];
      nsourcemax = nsource;
    }
}

void
Int2eV3::copy_to_source(double *integrals, int nsource)
{
  int i;
  double *tmp, *tmp2;

  /* Allocate more temporary space if needed. */
  source_space(nsource);

  tmp = source;
  tmp2 = integrals;
  for (i=0; i<nsource; i++) *tmp++ = *tmp2++;
}

/* make sure enough space exists for the source integrals */
void
Int1eV3::source_space(int nsource)
{
  if (nsourcemax < nsource) {
      delete[] source;
      source = new double[nsource];
      nsourcemax = nsource;
    }
}

void
Int1eV3::copy_to_source(double *integrals, int nsource)
{