obint.cc

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

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

#ifdef __GNUC__
#pragma implementation
#endif

#include <stdexcept>

#include <util/misc/formio.h>

#include <math/scmat/block.h>
#include <math/scmat/blkiter.h>
#include <math/scmat/offset.h>

#include <chemistry/qc/basis/obint.h>
#include <chemistry/qc/basis/basis.h>

using namespace sc;

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

void
EfieldDotVectorData::set_position(double*p)
{
  position[0] = p[0];
  position[1] = p[1];
  position[2] = p[2];
}

void
EfieldDotVectorData::set_vector(double*v)
{
  vector[0] = v[0];
  vector[1] = v[1];
  vector[2] = v[2];
}

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

void
DipoleData::set_origin(double*o)
{
  origin[0] = o[0];
  origin[1] = o[1];
  origin[2] = o[2];
}

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

PointChargeData::PointChargeData(int ncharges,
                                 const double *const*positions,
                                 const double *charges,
                                 int copy_data)
{
  ncharges_ = ncharges;
  if (copy_data) {
    alloced_positions_ = new double*[ncharges];
    alloced_charges_ = new double[ncharges];
    memcpy(alloced_charges_, charges, sizeof(double)*ncharges);
    double *tmp = new double[ncharges*3];
    for (int i=0; i<ncharges; i++) {
      alloced_positions_[i] = tmp;
      for (int j=0; j<3; j++) {
        *tmp++ = positions[i][j];
      }
    }
    positions_ = alloced_positions_;
    charges_ = alloced_charges_;
  }
  else {
    charges_ = charges;
    alloced_charges_ = 0;
    alloced_positions_ = 0;
    charges_ = charges;
    positions_ = positions;
  }
}

PointChargeData::~PointChargeData()
{
  if (alloced_positions_) {
    delete[] alloced_positions_[0];
    delete[] alloced_positions_;
  }
  delete[] alloced_charges_;
}

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

OneBodyInt::OneBodyInt(Integral* integral,
                       const Ref<GaussianBasisSet>&bs1,
                       const Ref<GaussianBasisSet>&bs2) :
  integral_(integral),
  bs1_(bs1), bs2_(bs2)
{
  buffer_ = 0;
}

OneBodyInt::~OneBodyInt()
{
}

int
OneBodyInt::nbasis() const
{
  return bs1_->nbasis();
}

int
OneBodyInt::nbasis1() const
{
  return bs1_->nbasis();
}

int
OneBodyInt::nbasis2() const
{
  return bs2_->nbasis();
}

int
OneBodyInt::nshell() const
{
  return bs1_->nshell();
}

int
OneBodyInt::nshell1() const
{
  return bs1_->nshell();
}

int
OneBodyInt::nshell2() const
{
  return bs2_->nshell();
}

Ref<GaussianBasisSet>
OneBodyInt::basis()
{
  return bs1_;
}

Ref<GaussianBasisSet>
OneBodyInt::basis1()
{
  return bs1_;
}

Ref<GaussianBasisSet>
OneBodyInt::basis2()
{
  return bs2_;
}

const double *
OneBodyInt::buffer() const
{
  return buffer_;
}

void
OneBodyInt::reinitialize()
{
}

bool
OneBodyInt::cloneable()
{
  return false;
}

Ref<OneBodyInt>
OneBodyInt::clone()
{
  throw std::runtime_error("OneBodyInt::clone() not implemented");
}

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

ShellPairIter::ShellPairIter()
{
}

ShellPairIter::~ShellPairIter()
{
}

void
ShellPairIter::init(const double * b, int ishell, int jshell,
                    int fi, int fj, int ni, int nj,
                    int red, double scl)
{
  e12 = ((ishell==jshell) && red);
  
  ioffset=fi;
  joffset=fj;

  iend=ni;
  jend=nj;

  buf=b;
  scale_=scl;
}

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

OneBodyIntIter::OneBodyIntIter()
{
}

OneBodyIntIter::OneBodyIntIter(const Ref<OneBodyInt>& o) :
  obi(o)
{
}

OneBodyIntIter::~OneBodyIntIter()
{
}

void
OneBodyIntIter::start(int ist, int jst, int ien, int jen)
{
  istart=ist;
  jstart=jst;
  iend=ien;
  jend=jen;
  
  icur=istart;
  jcur=jstart;

  if (!iend) {
    iend=obi->nshell1();
    jend=obi->nshell2();
  }

  ij = (icur*(icur+1)>>1) + jcur;
}

static inline int
min(int i, int j)
{
  return (i<j) ? i : j;
}

void
OneBodyIntIter::next()
{
  int jlast = (redund) ? min(icur,jend-1) : jend-1;
  
  if (jcur < jlast) {
    jcur++;
    ij++;
    return;
  }

  jcur=jstart;
  icur++;

  ij = (icur*(icur+1)>>1) + jcur;
}

double
OneBodyIntIter::scale() const
{
  return 1.0;
}

ShellPairIter&
OneBodyIntIter::current_pair()
{
  obi->compute_shell(icur,jcur);
  spi.init(obi->buffer(), icur, jcur,
           obi->basis1()->shell_to_function(icur),
           obi->basis2()->shell_to_function(jcur),
           obi->basis1()->operator()(icur).nfunction(),
           obi->basis2()->operator()(jcur).nfunction(),
           redund, scale()
           );

  return spi;
}

bool
OneBodyIntIter::cloneable()
{
  return obi->cloneable();
}

Ref<OneBodyIntIter>
OneBodyIntIter::clone()
{
  return new OneBodyIntIter(obi->clone());
}

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

OneBodyIntOp::OneBodyIntOp(const Ref<OneBodyInt>& it)
{
  iter = new OneBodyIntIter(it);
}

OneBodyIntOp::OneBodyIntOp(const Ref<OneBodyIntIter>& it) :
  iter(it)
{
}

OneBodyIntOp::~OneBodyIntOp()
{
}

bool
OneBodyIntOp::cloneable()
{
  return iter->cloneable();
}

Ref<SCElementOp>
OneBodyIntOp::clone()
{
  return new OneBodyIntOp(iter->clone());
}

void
OneBodyIntOp::process(SCMatrixBlockIter& b)
{
  ExEnv::err0() << indent
       << "OneBodyIntOp::process: cannot handle generic case\n";
  abort();
}

void
OneBodyIntOp::process_spec_rect(SCMatrixRectBlock* b)
{
  Ref<GaussianBasisSet> bs1 = iter->one_body_int()->basis1();
  Ref<GaussianBasisSet> bs2 = iter->one_body_int()->basis2();
  
  // convert basis function indices into shell indices
  int ishstart = bs1->function_to_shell(b->istart);
  int jshstart = bs2->function_to_shell(b->jstart);

  int b1end = b->iend;
  int ishend = (b1end?bs1->function_to_shell(b1end-1) + 1 : 0);

  int b2end = b->jend;
  int jshend = (b2end?bs2->function_to_shell(b2end-1) + 1 : 0);

  int njdata = b->jend - b->jstart;

  iter->set_redundant(0);

  for (iter->start(ishstart,jshstart,ishend,jshend);
       iter->ready(); iter->next()) {
    ShellPairIter& spi = iter->current_pair();

    for (spi.start(); spi.ready(); spi.next()) {
      int ifn = spi.i();
      int jfn = spi.j();