btest.cc

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

//
// btest.cc --- test program
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Edward Seidl <[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 <scconfig.h>
#ifdef HAVE_SSTREAM
#  include <sstream>
#else
#  include <strstream.h>
#endif

#include <util/misc/formio.h>
#include <util/keyval/keyval.h>
#include <util/state/stateio.h>
#include <util/state/state_text.h>
#include <util/state/state_bin.h>
#include <chemistry/qc/basis/basis.h>
#include <chemistry/qc/basis/files.h>
#include <chemistry/qc/basis/petite.h>
#include <chemistry/qc/basis/gpetite.h>
#include <chemistry/qc/basis/symmint.h>
#include <chemistry/qc/intv3/intv3.h>
#include <chemistry/qc/basis/sobasis.h>
#include <chemistry/qc/basis/sointegral.h>
#include <chemistry/qc/basis/extent.h>
#include <chemistry/qc/basis/orthog.h>

using namespace std;
using namespace sc;

static void
do_so_shell_test(const Ref<SOBasis>& sobas, const Ref<TwoBodySOInt> &soer,
                 int i, int j, int k, int l)
{
  if (i>=soer->basis1()->nshell()
      ||j>=soer->basis2()->nshell()
      ||k>=soer->basis3()->nshell()
      ||l>=soer->basis4()->nshell()) return;

  int p,q,r,s;
  soer->compute_shell(i,j,k,l);
  const double *buf = soer->buffer();
  int np = sobas->nfunction(i);
  int nq = sobas->nfunction(j);
  int nr = sobas->nfunction(k);
  int ns = sobas->nfunction(l);
  int off = 0;
  cout << "SHELL ("<<i<<j<<"|"<<k<<l<<"):" << endl;
  for (p=0; p<np; p++) {
      for (q=0; q<nq; q++) {
          for (r=0; r<nr; r++) {
              cout << "      ("<<p<<q<<"|"<<r<<"*) =";
              for (s=0; s<ns; s++, off++) {
                  cout << scprintf(" % 10.6f",buf[off]);
                }
              cout << endl;
            }
        }
    }
}

static void
do_so_shell_test(const Ref<SOBasis>& sobas, const Ref<OneBodySOInt> &soov,
                 int i, int j)
{
  if (i>=soov->basis1()->nshell()
      ||j>=soov->basis2()->nshell()) return;

  int p,q;
  soov->compute_shell(i,j);
  const double *buf = soov->buffer();
  int np = sobas->nfunction(i);
  int nq = sobas->nfunction(j);
  int off = 0;
  cout << "SHELL ("<<i<<"|"<<j<<"):" << endl;
  for (p=0; p<np; p++) {
      cout << "      ("<<p<<"|"<<"*) =";
      for (q=0; q<nq; q++, off++) {
          cout << scprintf(" % 10.6f",buf[off]);
        }
      cout << endl;
    }
}

static void
do_so_test(const Ref<KeyVal> &keyval,
           const Ref<Integral>& intgrl, const Ref<GaussianBasisSet> &gbs)
{
  intgrl->set_basis(gbs);

  Ref<SOBasis> sobas = new SOBasis(gbs, intgrl);
  sobas->print(cout);

  Ref<TwoBodyInt> aoer = intgrl->electron_repulsion();
  Ref<TwoBodySOInt> soer = new TwoBodySOInt(aoer);

  Ref<OneBodyInt> aoov = intgrl->overlap();
  Ref<OneBodySOInt> soov = new OneBodySOInt(aoov);

  sobas = soer->basis();
  sobas->print(cout);

  if (keyval->exists(":shell")) {
    do_so_shell_test(sobas, soer,
                     keyval->intvalue(":shell",0),
                     keyval->intvalue(":shell",1),
                     keyval->intvalue(":shell",2),
                     keyval->intvalue(":shell",3));
    do_so_shell_test(sobas, soov,
                     keyval->intvalue(":shell",0),
                     keyval->intvalue(":shell",1));
    }
  else {
      int i,j,k,l;
      cout << "SO Electron Repulsion:" << endl;
      for (i=0; i<sobas->nshell(); i++) {
          for (j=0; j<sobas->nshell(); j++) {
              for (k=0; k<sobas->nshell(); k++) {
                  for (l=0; l<sobas->nshell(); l++) {
                      do_so_shell_test(sobas, soer, i, j, k, l);
                    }
                }
            }
        }
      cout << "SO Overlap:" << endl;
      for (i=0; i<sobas->nshell(); i++) {
          for (j=0; j<sobas->nshell(); j++) {
              do_so_shell_test(sobas, soov, i, j);
            }
        }
    }
}

static void
test_overlap(const Ref<GaussianBasisSet>& gbs, const Ref<GaussianBasisSet>& gbs2,
             const Ref<Integral>& intgrl)
{
  intgrl->set_basis(gbs);

  // first form AO basis overlap
  RefSymmSCMatrix s(gbs->basisdim(), gbs->matrixkit());
  Ref<SCElementOp> ov
      = new OneBodyIntOp(new OneBodyIntIter(intgrl->overlap()));
  s.assign(0.0);
  s.element_op(ov);
  ov=0;
  s.print("overlap");
      
  // now transform s to SO basis
  Ref<PetiteList> pl = intgrl->petite_list();
  RefSymmSCMatrix sb = pl->to_SO_basis(s);
  sb.print("blocked s");
      
  // and back to AO basis
  s = pl->to_AO_basis(sb);
  s.print("reconstituted s");

  // form skeleton overlap
  ov = new OneBodyIntOp(new SymmOneBodyIntIter(intgrl->overlap(),pl));
  s.assign(0.0);
  s.element_op(ov);
  ov=0;
  s.print("overlap");

  // and symmetrize to get blocked overlap again
  sb.assign(0.0);
  pl->symmetrize(s,sb);
  sb.print("blocked again");

  s=0; sb=0;

  // now try overlap between two basis sets
  RefSCMatrix ssq(gbs2->basisdim(),gbs->basisdim(),gbs2->matrixkit());
  intgrl->set_basis(gbs2,gbs);

  ov = new OneBodyIntOp(new OneBodyIntIter(intgrl->overlap()));
  ssq.assign(0.0);
  ssq.element_op(ov);
  ssq.print("overlap sq");
  ov=0;

  Ref<PetiteList> pl2 = intgrl->petite_list(gbs2);
  RefSCMatrix ssqb(pl2->AO_basisdim(), pl->AO_basisdim(), gbs->so_matrixkit());
  ssqb->convert(ssq);

  RefSCMatrix syms2 = pl2->aotoso().t() * ssqb * pl->aotoso();
  syms2.print("symm S2");
}

static void
test_aoorthog(const Ref<GaussianBasisSet>& gbs,
             const Ref<Integral>& intgrl)
{
  cout << "Beginning AO Orthog test" << endl;

  intgrl->set_basis(gbs);

  // first form AO basis overlap
  RefSymmSCMatrix s(gbs->basisdim(), gbs->matrixkit());
  Ref<SCElementOp> ov
      = new OneBodyIntOp(new OneBodyIntIter(intgrl->overlap()));
  s.assign(0.0);
  s.element_op(ov);
  ov=0;

  Ref<OverlapOrthog> orthog;
  orthog = new OverlapOrthog(OverlapOrthog::Symmetric,
                             s,
                             s.kit(),
                             0.0001,
                             1);
  orthog->basis_to_orthog_basis().print("basis to orthog basis");

  s = 0;
  orthog = 0;

  cout << "Ending AO Orthog test" << endl;
}

static void
test_eigvals(const Ref<GaussianBasisSet>& gbs, const Ref<Integral>& intgrl)
{
  intgrl->set_basis(gbs);
  Ref<PetiteList> pl = intgrl->petite_list();

  // form AO Hcore and evecs
  RefSymmSCMatrix hcore_ao(gbs->basisdim(), gbs->matrixkit());
  RefSCMatrix ao_evecs(gbs->basisdim(), gbs->basisdim(), gbs->matrixkit());
  RefDiagSCMatrix ao_evals(gbs->basisdim(), gbs->matrixkit());
  
  hcore_ao.assign(0.0);

  Ref<SCElementOp> op
      = new OneBodyIntOp(new OneBodyIntIter(intgrl->kinetic()));
  hcore_ao.element_op(op);
  op=0;

  Ref<OneBodyInt> nuc = intgrl->nuclear();
  nuc->reinitialize();
  op = new OneBodyIntOp(nuc);
  hcore_ao.element_op(op);
  op=0;
  
  hcore_ao.print("Hcore (AO)");
  
  hcore_ao.diagonalize(ao_evals, ao_evecs);
  ao_evecs.print("AO Evecs");
  ao_evals.print("AO Evals");

  // form SO Hcore and evecs
  RefSymmSCMatrix hcore_so(pl->SO_basisdim(), gbs->so_matrixkit());
  RefSCMatrix so_evecs(pl->SO_basisdim(), pl->SO_basisdim(),
                       gbs->so_matrixkit());
  RefDiagSCMatrix so_evals(pl->SO_basisdim(), gbs->so_matrixkit());
  
  // reuse hcore_ao to get skeleton Hcore
  hcore_ao.assign(0.0);

  op = new OneBodyIntOp(new SymmOneBodyIntIter(intgrl->kinetic(),pl));
  hcore_ao.element_op(op);
  op=0;

  nuc = intgrl->nuclear();
  nuc->reinitialize();
  op = new OneBodyIntOp(new SymmOneBodyIntIter(nuc,pl));
  hcore_ao.element_op(op);
  op=0;
  
  pl->symmetrize(hcore_ao, hcore_so);

  hcore_so.print("Hcore (SO)");
  
  hcore_so.diagonalize(so_evals, so_evecs);
  so_evecs.print("SO Evecs");
  so_evals.print("SO Evals");

  RefSCMatrix new_ao_evecs = pl->evecs_to_AO_basis(so_evecs);
  new_ao_evecs.print("AO Evecs again");

  //RefSCMatrix new_so_evecs = pl->evecs_to_SO_basis(ao_evecs);
  //new_so_evecs.print("SO Evecs again");

  pl->to_AO_basis(hcore_so).print("Hcore (AO) again");
}

void
checkerror(const char *name, int shell, int func,
           double numerical, double check)
{
  double mag = fabs(check);
  double err = fabs(numerical - check);
  cout << scprintf("%2s %2d %2d %12.8f %12.8f er = %6.4f",
                   name, shell, func,
                   numerical, check, err/mag) << endl;
  if (mag > 0.001) {
      if (err/mag > 0.05) {
          cout << scprintf("ERROR %2s %2d %2d %12.8f %12.8f er = %6.4f",
                           name, shell, func,
                           numerical, check, err/mag) << endl;
        }
    }
  else if (err > 0.02) {
      cout << scprintf("ERROR %2s %2d %2d %12.8f %12.8f ea = %16.14f",
                       name, shell, func, numerical, check, err) << endl;
    }
}

void
test_func_values(const Ref<GaussianBasisSet> &gbs,
                 const Ref<Integral> &integral)
{
  cout << "testing basis function value gradient and hessian numerically"
       << endl;

  int nbasis = gbs->nbasis();
  double *b_val = new double[nbasis];
  double *b_val_plsx = new double[nbasis];
  double *b_val_mnsx = new double[nbasis];
  double *b_val_plsy = new double[nbasis];
  double *b_val_mnsy = new double[nbasis];
  double *b_val_plsz = new double[nbasis];
  double *b_val_mnsz = new double[nbasis];
  double *b_val_plsyx = new double[nbasis];
  double *b_val_mnsyx = new double[nbasis];
  double *b_val_plszy = new double[nbasis];
  double *b_val_mnszy = new double[nbasis];