MPQC  3.0.0-alpha
Developing Code Using MPQC

In addition to the executables, the MPQC toolkit libraries and include files can be installed on your machine. This is described in the Compiling section of this manual.

The mpqc-config program can be use to obtain the compilers, compiler options, and libraries needed to use the MPQC toolkit from your program. This utility is discussed below, along with how the MPQC toolkit must be initialized in your main subroutine.

The mpqc-config Program

The mpqc-config program returns information about how MPQC was compiled and installed. See The mpqc-config Program for more information.

Initializing MPQC

A variety of initializations must be done when a program using the MPQC classes is started. This includes finding out what sc::MessageGrp to be used for communication, what sc::ThreadGrp is to be used for multi-threading, what default sc::Integral factory should be used, reading the input file in a way that is efficient on a parallel machine, setting up the output, and so on. The sc::MPQCInit class simplifies this setup. It is used by first initializing an sc::GetLongOpt object. This is used to retrieve flags specified on the command line. The programmer enrolls any desired options into the sc::GetLongOpt object and then passes this, along with argc and argv, to the constructor for sc::MPQCInit. The sc::MPQCInit class will enroll the options it requires: –messagegrp, –threadgrp, –memory, –resources, and –integral. Then the programmer can parse the options and retrieve the value of any flag or obtain anything that was specified on the command line that was not associated with a flag.

Below is a code fragment illustrating the use of sc::MPQCInit.

  sc::GetLongOpt opt;
  opt.usage("[options] [filename]");
  opt.enroll("verbose", GetLongOpt::NoValue, "enable extra printing", 0);
  sc::MPQCInit init(opt,argc,argv);
  int optind = opt.parse(argc, argv);
  std::string inputfile;
  if (argc - optind == 0) {
    inputfile = "input.in";
  }
  else if (argc - optind == 1) {
    inputfile = argv[optind];
  }
  else {
    opt.usage(std::cout);
    return 1;
  }
  sc::Ref<sc::KeyVal> keyval = init.init(inputfile);

Using MPQC's main: MP2 Implementation Example

This section illustrates how to add a new quantum chemistry method implementation to MPQC.

MP2 Implementation Example: Source

This example code illustrates a complete MP2 energy implementation using the MPQC Toolkit. First an MP2 class is declared and the necesary base class member functions are provided. Next a ClassDesc is defined. Finally, the member functions are defined.

Note that no main routine is provided. This is because this file is designed to be used to extend the functionality of the mpqc executable. To generate a new mpqc executable with the new class available for use, see the MP2 Implementation Example: Makefile section.

#include <stddef.h>
#include <util/misc/autovec.h>
#include <util/misc/scexception.h>
#include <chemistry/qc/wfn/obwfn.h>
#include <chemistry/qc/scf/clhf.h>
using namespace std;
using namespace sc;
class MP2: public Wavefunction {
double compute_mp2_energy();
public:
MP2(const Ref<KeyVal> &);
MP2(StateIn &);
void save_data_state(StateOut &);
void compute(void);
void obsolete(void);
int nelectron(void);
RefSymmSCMatrix density(void);
double magnetic_moment() const;
int value_implemented(void) const;
};
static ClassDesc MP2_cd(typeid(MP2), "MP2", 1, "public Wavefunction",
0, create<MP2>, create<MP2>);
MP2::MP2(const Ref<KeyVal> &keyval):Wavefunction(keyval) {
ref_mp2_wfn_ << keyval->describedclassvalue("reference");
if(ref_mp2_wfn_.null()) {
throw InputError("require a OneBodyWavefunction object",
__FILE__, __LINE__, "reference", 0,
}
}
MP2::MP2(StateIn &statein):Wavefunction(statein)
{
ref_mp2_wfn_ << SavableState::restore_state(statein);
}
void
MP2::save_data_state(StateOut &stateout) {
Wavefunction::save_data_state(stateout);
SavableState::save_state(ref_mp2_wfn_.pointer(),stateout);
}
void
MP2::compute(void)
{
if(gradient_needed()) {
throw FeatureNotImplemented("no gradients yet",
__FILE__, __LINE__, class_desc());
}
double extra_hf_acc = 10.;
ref_mp2_wfn_->set_desired_value_accuracy(desired_value_accuracy()
/ extra_hf_acc);
double refenergy = ref_mp2_wfn_->energy();
double mp2energy = compute_mp2_energy();
ExEnv::out0() << indent << "MP2 Energy = " << mp2energy << endl;
set_value(refenergy + mp2energy);
set_actual_value_accuracy(ref_mp2_wfn_->actual_value_accuracy()
* extra_hf_acc);
}
void
MP2::obsolete(void) {
Wavefunction::obsolete();
ref_mp2_wfn_->obsolete();
}
int
MP2::nelectron(void) {
return ref_mp2_wfn_->nelectron();
}
MP2::density(void) {
throw FeatureNotImplemented("no density yet",
__FILE__, __LINE__, class_desc());
return 0;
}
double
MP2::magnetic_moment() const {
return 0.0;
}
int
MP2::value_implemented(void) const {
return 1;
}
double
MP2::compute_mp2_energy()
{
if(molecule()->point_group()->char_table().order() != 1) {
throw FeatureNotImplemented("C1 symmetry only",
__FILE__, __LINE__, class_desc());
}
RefSCMatrix vec = ref_mp2_wfn_->eigenvectors();
int nao = vec.nrow();
int nmo = vec.ncol();
int nocc = ref_mp2_wfn_->nelectron()/2;
int nvir = nmo - nocc;
auto_vec<double> cvec_av(new double [vec.nrow() * vec.ncol()]);
double *cvec = cvec_av.get();
vec->convert(cvec);
auto_vec<double> pqrs_av(new double [nao * nao * nao * nao]);
double *pqrs = pqrs_av.get();
for(int n = 0; n < nao*nao*nao*nao; n++) pqrs[n] = 0.0;
Ref<TwoBodyInt> twoint = integral()->electron_repulsion();
const double *buffer = twoint->buffer();
Ref<GaussianBasisSet> basis = this->basis();
int nshell = basis->nshell();
for(int P = 0; P < nshell; P++) {
int nump = basis->shell(P).nfunction();
for(int Q = 0; Q < nshell; Q++) {
int numq = basis->shell(Q).nfunction();
for(int R = 0; R < nshell; R++) {
int numr = basis->shell(R).nfunction();
for(int S = 0; S < nshell; S++) {
int nums = basis->shell(S).nfunction();
twoint->compute_shell(P,Q,R,S);
int index = 0;
for(int p=0; p < nump; p++) {
int op = basis->shell_to_function(P)+p;
for(int q = 0; q < numq; q++) {
int oq = basis->shell_to_function(Q)+q;
for(int r = 0; r < numr; r++) {
int oor = basis->shell_to_function(R)+r;
for(int s = 0; s < nums; s++,index++) {
int os = basis->shell_to_function(S)+s;
int ipqrs = (((op*nao+oq)*nao+oor)*nao+os);
pqrs[ipqrs] = buffer[index];
}
}
}
}
}
}
}
}
twoint = 0;
auto_vec<double> ijkl_av(new double [nmo * nmo * nmo * nmo]);
double *ijkl = ijkl_av.get();
int idx = 0;
for(int i = 0; i < nmo; i++) {
for(int j = 0; j < nmo; j++) {
for(int k = 0; k < nmo; k++) {
for(int l = 0; l < nmo; l++, idx++) {
ijkl[idx] = 0.0;
int index = 0;
for(int p = 0; p < nao; p++) {
for(int q = 0; q < nao; q++) {
for(int r = 0; r < nao; r++) {
for(int s = 0; s < nao; s++,index++) {
ijkl[idx] += cvec[p*nmo + i] * cvec[q*nmo +j]
* cvec[r*nmo + k] * cvec[s*nmo + l]
* pqrs[index];
}
}
}
}
}
}
}
}
pqrs_av.release(); pqrs = 0;
cvec_av.release(); cvec = 0;
auto_vec<double> evals_av(new double [nmo]);
double *evals = evals_av.get();
ref_mp2_wfn_->eigenvalues()->convert(evals);
double energy = 0.0;
for(int i=0; i < nocc; i++) {
for(int j=0; j < nocc; j++) {
for(int a=nocc; a < nmo; a++) {
for(int b=nocc; b < nmo; b++) {
int iajb = (((i*nmo+a)*nmo+j)*nmo+b);
int ibja = (((i*nmo+b)*nmo+j)*nmo+a);
energy += (2 * ijkl[iajb] - ijkl[ibja]) * ijkl[iajb]/
(evals[i] + evals[j] - evals[a] - evals[b]);
}
}
}
}
ijkl_av.release(); ijkl = 0;
evals_av.release(); evals = 0;
return energy;
}

MP2 Implementation Example: Makefile

This example Makefile demonstrates how to link in a new class to form a new mpqc executable, here named mp2. The code is given in the MP2 Implementation Example: Source section. The mpqc-config command is used to obtain information about how the MPQC toolkit was compiled and installed. The library specified with -lmpqc provides the main routine from mpqc.

# Change this to the (absolute) path to your installed mpqc-config script, if needed
MPQCCONFIG = ../../../../../bin/mpqc-config
CXX := $(shell $(MPQCCONFIG) --cxx)
CXXFLAGS := $(shell $(MPQCCONFIG) --cxxflags)
CC := $(shell $(MPQCCONFIG) --cc)
CCFLAGS := $(shell $(MPQCCONFIG) --cflags)
F77 := $(shell $(MPQCCONFIG) --f77)
F90 := $(F77)
FFLAGS := $(shell $(MPQCCONFIG) --f77flags)
F90FLAGS := $(FFLAGS)
CPPFLAGS := $(shell $(MPQCCONFIG) --cppflags)
LIBS := $(shell $(MPQCCONFIG) --libs)
LIBDIR := $(shell $(MPQCCONFIG) --libdir)
LDFLAGS := $(shell $(MPQCCONFIG) --ldflags)
TARGET = mp2
$(TARGET): $(TARGET).o
$(CXX) $(CXXFLAGS) -o $@ $^ -L$(LIBDIR) $(LDFLAGS) -lmpqcmain $(LIBS)
clean:
-rm -f $(TARGET) $(TARGET).o

MP2 Implementation Example: Input

This input file can be used with the program illustrated in the MP2 Implementation Example: Source section. It will compute the MP2 energy using the new MP2 class. Note that only the object-oriented input format can be used with user provided classes.

% emacs should use -*- KeyVal -*- mode
molecule<Molecule>: (
symmetry = C1
unit = angstrom
{ atoms geometry } = {
O [ 0.00000000 0.00000000 0.37000000 ]
H [ 0.78000000 0.00000000 -0.18000000 ]
H [ -0.78000000 0.00000000 -0.18000000 ]
}
)
basis<GaussianBasisSet>: (
name = "STO-3G"
molecule = $:molecule
)
mpqc: (
checkpoint = no
savestate = no
% MP2 is the new class. Change MP2 to MBPT2 to test
% against the standard MP2 code
mole<MP2>: (
molecule = $:molecule
basis = $:basis
reference<CLHF>: (
molecule = $:molecule
basis = $:basis
memory = 16000000
)
)
)

A user-provided main: HF Implementation Example

The MP2 sample program above used the MPQC main routine. This does not always provide sufficient flexibility, and a user-provide main routine is sometimes more appropriate. However, quite a few steps are need to prepare for a run, most of which have to do with initializing various libraries to prepare for a parallel run. The HF example illustrates how an sc::MPQCInit object can be used to prepare for a computation. It implements a complete closed-shell Hartree-Fock program using various MPQC functionality, without provide a new MPQC sc::Wavefunction implementation as done in the previous example (Using MPQC's main: MP2 Implementation Example).

This example program can be found in the doc/devsamp/hf subdirectory of the MPQC source distribution and is shown here in its entirety.

HF Implementation Example: Source Code

#include <chemistry/molecule/molecule.h>
#include <chemistry/qc/basis/integral.h>
#include <chemistry/qc/basis/symmint.h>
#include <chemistry/qc/basis/orthog.h>
#include <chemistry/qc/intv3/intv3.h>
#include <chemistry/qc/lcao/fockbuild.h>
#include <chemistry/qc/lcao/fockdist.h>
#include <chemistry/qc/lcao/clhfcontrib.h>
#include <math/optimize/scextrap.h>
#include <math/optimize/diis.h>
#include <math/optimize/scextrapmat.h>
#include <util/group/mstate.h>
#include <util/group/pregtime.h>
#include <util/group/memory.h>
#include <util/group/thread.h>
#include <exception>
#include <mpqcinit.h>
using namespace sc;
using std::cout;
using std::endl;
int
try_main(int argc, char **argv)
{
opt.usage("[options] [filename]");
opt.enroll("verbose", GetLongOpt::NoValue, "enable extra printing", 0);
MPQCInit init(opt,argc,argv);
int optind = opt.parse(argc, argv);
std::string inputfile;
if (argc - optind == 0) {
inputfile = "hf.in";
}
else if (argc - optind == 1) {
inputfile = argv[optind];
}
else {
opt.usage(std::cout);
return 1;
}
Ref<sc::KeyVal> keyval = init.init(inputfile);
if (opt.retrieve("verbose")) {
sc::ExEnv::out0() << indent
<< "nthread = " << thr->nthread() << std::endl;
sc::ExEnv::out0() << indent
<< "nnode = " << msg->n() << std::endl;
}
basis << keyval->describedclassvalue("basis");
// make an integral factory and get the petite list
Ref<Integral> integral = new IntegralV3(basis);
Ref<PetiteList> pl = integral->petite_list();
integral->set_storage(32000000);
// construct the overlap matrix, s
RefSymmSCMatrix s_skel(basis->basisdim(), basis->matrixkit());
= new OneBodyIntOp(new SymmOneBodyIntIter(integral->overlap(), pl));
s_skel.assign(0.0);
s_skel.element_op(s_op);
s_op=0;
RefSymmSCMatrix s_SO(pl->SO_basisdim(), basis->so_matrixkit());
pl->symmetrize(s_skel,s_SO);
// construct the basis set orthogonalizer
= new OverlapOrthog(OverlapOrthog::Symmetric,
s_SO,
basis->so_matrixkit(),
1.0e-6 /*lindep tolerance*/,
0 /* debug */
);
// construct the core hamiltonian (nuclear attraction + kinetic energy)
RefSymmSCMatrix hcore_skel(basis->basisdim(), basis->matrixkit());
= new OneBodyIntOp(new SymmOneBodyIntIter(integral->hcore(), pl));
hcore_skel.assign(0.0);
hcore_skel.element_op(hcore_op);
hcore_op=0;
RefSymmSCMatrix hcore_SO(pl->SO_basisdim(), basis->so_matrixkit());
pl->symmetrize(hcore_skel,hcore_SO);
RefSymmSCMatrix density_AO(pl->AO_basisdim(), basis->so_matrixkit());
density_AO.assign(0.0);
RefSCMatrix vector_OSO(pl->SO_basisdim(), orthog->orthog_dim(),
basis->so_matrixkit());
vector_OSO.assign(0.0);
RefDiagSCMatrix evals(orthog->orthog_dim(), basis->so_matrixkit());
RefSymmSCMatrix density_SO(pl->SO_basisdim(), basis->so_matrixkit());
extrap->set_tolerance(1.0e-6);
// read the occupation numbers
RefSymmSCMatrix occ(orthog->orthog_dim(), basis->so_matrixkit());
occ.assign(0.0);
for (int i=0; i<orthog->orthog_dim(); i++) {
occ(i,i) = keyval->doublevalue("occ",i);
}
// construct the initial guess to vector_OSO and density_AO
RefSymmSCMatrix hcore_OSO(orthog->orthog_dim(), basis->so_matrixkit());
hcore_OSO.assign(0.0);
hcore_OSO.accumulate_transform(orthog->basis_to_orthog_basis(),hcore_SO);
hcore_OSO.diagonalize(evals, vector_OSO);
density_SO.assign(0.0);
RefSCMatrix vector_SO = orthog->basis_to_orthog_basis().t() * vector_OSO;
density_SO.accumulate_transform(vector_SO, occ);
density_AO = pl->to_AO_basis(density_SO);
while (!extrap->converged()) {
//while (!extrap->converged()) {
// construct the G matrix
RefSymmSCMatrix g_skel(basis->basisdim(), basis->matrixkit());
g_skel.assign(0.0);
= new CLHFContribution(basis,basis,basis,"replicated");
g_contrib->set_fmat(0, g_skel);
g_contrib->set_pmat(0, density_AO);
Ref<FockBuild> fb = new FockBuild(fockdist, g_contrib, false,
basis, basis, basis);
fb->set_accuracy(1e-12);
fb->build();
g_skel.scale(1.0/(double)pl->order());
RefSymmSCMatrix g_SO(pl->SO_basisdim(), basis->so_matrixkit());
pl->symmetrize(g_skel,g_SO);
// compute the exchange/correlation potential
double exc = 0.0;
// construct the fock matrix
RefSymmSCMatrix f_SO = hcore_SO + g_SO;
// transform the fock matrix into the orthogonal basis
RefSymmSCMatrix f_OSO(orthog->orthog_dim(), basis->so_matrixkit());
f_OSO.assign(0.0);
f_OSO.accumulate_transform(orthog->basis_to_orthog_basis(),f_SO);
// Extrapolate the fock matrix
RefSymmSCMatrix error_MO(orthog->orthog_dim(), basis->so_matrixkit());
error_MO.assign(0.0);
error_MO.accumulate_transform(vector_OSO.t(),f_OSO);
error_MO->scale_diagonal(0.0);
RefSymmSCMatrix error_SO(pl->SO_basisdim(), basis->so_matrixkit());
error_SO.assign(0.0);
error_SO.accumulate_transform(orthog->basis_to_orthog_basis().t()
*vector_OSO, error_MO);
extrap->extrapolate(data,error);
// Diagonalize the fock matrix
f_OSO.diagonalize(evals, vector_OSO);
// Compute the density matrix
density_SO.assign(0.0);
vector_SO = orthog->basis_to_orthog_basis().t() * vector_OSO;
density_SO.accumulate_transform(vector_SO, occ);
density_AO = pl->to_AO_basis(density_SO);
double energy =exc + basis->molecule()->nuclear_repulsion_energy()
+(density_SO * 0.5 * (hcore_SO + f_SO)).trace();
ExEnv::out0() << indent
<< scprintf("energy = %14.8f, error = %10.8f",
energy, extrap->error()) << std::endl;
}
return 0;
}
int
main(int argc, char *argv[])
{
try {
try_main(argc, argv);
}
catch (std::bad_alloc &e) {
cout << argv[0] << ": ERROR: MEMORY ALLOCATION FAILED:" << endl
<< e.what()
<< endl;
exit(1);
}
catch (std::exception &e) {
cout << argv[0] << ": ERROR: EXCEPTION RAISED:" << endl
<< e.what()
<< endl;
exit(1);
}
catch (...) {
cout << argv[0] << ": ERROR: UNKNOWN EXCEPTION RAISED" << endl;
exit(1);
}
return 0;
}

HF Implementation Example: Makefile

This example Makefile will build an executable named hf. The mpqc-config command is used to obtain information about how the MPQC toolkit was compiled and installed. The library specified with -lmpqcinit provides the initialization routines from mpqc.

# Change this to the (absolute) path to your installed mpqc-config script, if needed
MPQCCONFIG = ../../../../../bin/mpqc-config
CXX := $(shell $(MPQCCONFIG) --cxx)
CXXFLAGS := $(shell $(MPQCCONFIG) --cxxflags)
CC := $(shell $(MPQCCONFIG) --cc)
CCFLAGS := $(shell $(MPQCCONFIG) --cflags)
F77 := $(shell $(MPQCCONFIG) --f77)
F90 := $(F77)
FFLAGS := $(shell $(MPQCCONFIG) --f77flags)
F90FLAGS := $(FFLAGS)
CPPFLAGS := $(shell $(MPQCCONFIG) --cppflags)
LIBS := $(shell $(MPQCCONFIG) --libs)
LIBDIR := $(shell $(MPQCCONFIG) --libdir)
LDFLAGS := $(shell $(MPQCCONFIG) --ldflags)
TARGET = hf
$(TARGET): $(TARGET).o
$(CXX) $(CXXFLAGS) -o $@ $^ -L$(LIBDIR) $(LDFLAGS) $(LIBS)
clean:
-rm -f $(TARGET) $(TARGET).o

HF Implementation Example: Input

This input will run a Hartree-Fock calculation on a water molecule. The section beginning mpqc: provides the input needed to run the mpqc program with the same basis set and geometry as a check. The hf program only reads the occ and basis keywords at the beginning of the sample input.

% occupation vector elements must be 0 or 2
occ = [ 2 2 2 2 2 ]
basis<GaussianBasisSet>: (
name = "STO-3G"
puream = 0
molecule<Molecule>: (
symmetry = c1
{ atoms geometry } = {
O [ 0 0 0 ]
H [ 0 -1 1]
H [ 0 1 1]
}
)
)
% This is not needed by the hf.cc example. It is here to allow this
% input file to be used with MPQC for testing purposes.
mpqc: (
mole<CLHF>: (
molecule = $:basis:molecule
basis = $:basis
density_reset_frequency = 1
guess_wavefunction<HCoreWfn>: (
molecule = $:basis:molecule
basis = $:basis
)
)
)

Exception Handling in MPQC

The development of MPQC began before exception handling was available in C++. A retrofit of the code to use exceptions is in progress. It is difficult to retrofit a code, especially a parallel code, to do exception handling. There will be some limitations: exception handling will not work well for parallel jobs, objects whose members throw might be left in a questionable state, etc. However, it is intended that MPQC objects will be usable in an interactive environment. It is also planned that exceptions be used internally to facilitate recover from certain problems.

All new code should use exceptions instead of exit or abort and allocate resources in such a way that, if an exception occurs, all resources such as memory or locks are released. A hierarchy of exception classes has been created that maps better to scientific computing than the standard exceptions. More information is below, as well as in the documentation for the SCException class and its derivatives.

Exceptions and Memory Allocation

Consider the following code fragment:

Object *obj = new Object;
double *array = new double[n];
f(obj, array, mol);
delete obj;
delete[] array;

If an exception is thrown in the function f(), then storage for array and obj will not be released. The standard C++ library provides a class, auto_ptr, to deal with obj, and the MPQC toolkit provides a class, auto_vec, to deal with array.

The include files for these two classes are:

#include <memory>
#include <util/misc/autovec.h>

the code would be modified as follows:

std::auto_ptr<Object> obj(new Object);
sc::auto_vec<double> array(new double[n]);
f(obj.get(), array.get());
obj.release();  // or just let the destructor release it
array.release();  // or just let the destructor release it

Note that when sc::Ref is used to store pointers, the storage will automatically be released when necessary. No special treatment is needed to deal with exceptions.

Exceptions and Locks

Consider the following code fragment:

g(const sc::Ref<sc::ThreadLock> &lock)
{
  lock->lock();
  f();
  lock->unlock();
}

If f() throws, then the lock is never released. The ThreadLock lock() and unlock() members should not be used anymore. Now do the following:

g(const sc::Ref<sc::ThreadLock> &lock)
{
  sc::ThreadLockHolder lockholder(lock);
  f();
  lockholder->unlock(); // or let the destructor unlock it
}

Exceptions and Region Timers

Consider the following code fragment:

g(const sc::Ref<sc::RegionTimer> &regtim)
{
  regtim->enter("f()");
  f();
  regtim->exit();
}

If f() throws, then the "f()" timing region is never exited. Instead use the following:

g(const sc::Ref<sc::RegionTimer> &regtim)
{
  sc::Timer timer(regtim, "f()");
  f();
  timer.reset(); // or let the destructor exit the region
}

Using the MPQC Exception Classes

The MPQC exceptions provide information that can be used into two ways. First, text information is provided so that if the exception is not caught at a lower level, then the mpqc executable will catch it and write information about the problem to the terminal or an output file. Second, information about the nature of the problem is provided, to permit developers to catch the exception and deal with it in some way. The documentation for sc::SCException and all of its derivatives gives more information about the exceptions that are available. As an example, consider the following loop, where a maximum number of iterations is permitted:

XYZ::update()
{
  for (int i=0; i<maxiter; i++) {
    // ... compute xyz update ...
    if (residual < threshold) return;
  }
  throw MaxIterExceeded("too many iterations xyz computation",
                        __FILE__, __LINE__, maxiter, class_desc());
}

The first argument to the exception class is a brief description of the error. Additional information can be provided, see SCException::elaborate() description below. The next two arguments are the filename and line number. The C preprocessor provides these for you with the FILE and LINE macros. The next argument is specific to the MaxIterExceeded exception; it is the maximum number of iterations. Finally, a ClassDesc* can be given, which will be used to print out the class name of the object that failed. All of these arguments are optional; however, the first three should always be given.

It is possible to provide additional information using the SCException::elaborate() member. This will return a ostream, and the additional information can be written to this stream. However, if for some reason it is not possible to write to this stream (say, there wasn't enough memory to allocate it), then an exception will be thrown. For this reason, the string description given as the first argument should be informative since the additional information might not be available, and attempts to use elaborate() should be in a try block. So, for example, the elaborate() member could be used in the above example as follows:

XYZ::update()
{
  for (int i=0; i<maxiter; i++) {
    // ... compute xyz update ...
    if (residual < threshold) return;
  }
  MaxIterExceeded ex("too many iterations in xyz computation",
                     __FILE__, __LINE__, maxiter, class_desc());
  try {
    ex.elaborate() << "this can happen when the stepsize is too small"
                   << std::endl
                   << "the stepsize is " << stepsize
                   << std::endl;
  }
  catch (...) {}
  throw ex;
}

Note that writing to stream returned by elaborate() won't necessarily cause anything to get written to the terminal or an output file. The information will be available when the what() member is called, if writing to the stream succeeds. If the exception is caught by the mpqc main routine, then it will be printed for the user to see. If the program catches the exception and determines that it is possible to proceed in a different way, then the user will never see the text.

Debugging Code with Exceptions

Usually, exceptions are not the desired behaviour in a program, and it is necessary to debug a program that throws an exception. This was easy when abort was called, because abort would raise a signal that was caught by the debugger and the code is stopped at the appropriate place. With exceptions the matter is more complex, because the stack is unwound when an exception is thrown and most debugging information is lost. To work around this problem, a breakpoint can be set in code that will be reached only in an exception, and will be run before the stack unwind begins. A useful place to do this when GCC is used as the compiler is in the routine __cxa_allocate_exception(). So, in gdb, the following could be done:

$ gdb ./scextest
(gdb) b main
(gdb) run
Breakpoint 1, main () at /home/cljanss/src/SC/src/lib/util/misc/scextest.cc:172
172           f();
(gdb) b __cxa_allocate_exception
(gdb) cont
Breakpoint 2, 0x40582d46 in __cxa_allocate_exception ()
   from /usr/lib/gcc-lib/i686-pc-linux-gnu/3.3.5/libstdc++.so.5
(gdb) where
#0  0x40582d46 in __cxa_allocate_exception ()
   from /usr/lib/gcc-lib/i686-pc-linux-gnu/3.3.5/libstdc++.so.5
#1  0x0804b3f7 in f () at /home/cljanss/src/SC/src/lib/util/misc/scextest.cc:60
#2  0x0804b9e9 in main ()
    at /home/cljanss/src/SC/src/lib/util/misc/scextest.cc:172

Giving gdb "b main" followed by "run" was required before gdb could find the __cxa_allocate_exception symbol.

Adding Test Cases to the Verification Suite

There are two ways to test an MPQC build. The testbuild and testrun make targets can be used to run test programs in various library directories, and the check and related make targets can be used to run MPQC on sets of input files. See Validating MPQC for more information about how to run the tests.

Test programs can be added to the library directories by providing a source file with a main routine. The set of test programs that is to be built and run by testbuild and testrun, respectively, is given by the TESTPROGS variable in the library's Makefile. It may be necessary for an explicit rule to be given for building the test program to ensure that necessary libraries are linked in. If a file named after the test program with a .out suffix is found in the source directory, then testrun fail if the command's output differs from that file. Care must be taken to ensure that the output is architecture independent in this case. Otherwise, testrun will fail only if running the command results in a nonzero return code.

Additional MPQC test inputs can be added in the test/ref directory. These inputs can be provided in one of two ways. An input which is used to automatically generate multiple test cases can be written (with a .qci suffix), or a subdirectory with each input can be made. See Makefile, basis1.qci, and input in the test/ref directory for examples.

After you have added new inputs and modified the Makefile, change into the test/ref subdirectory of your object directory (where you compiled MPQC) and type make inputs. This will create a input subdirectory containing MPQC input files with a .in suffix. Files ending with a .qci suffix will also be placed in the input directory. These contain a description of the calculation that is used by the utility program that checks the results of the validation suite. Both the .in and .qci files for the new test cases must be copied into the ref directory in the source tree. Note that inputs that are not useful in your build environment are not created by make inputs.

sc::MessageGrp::get_default_messagegrp
static MessageGrp * get_default_messagegrp()
Returns the default message group.
sc::PetiteList::AO_basisdim
RefSCDimension AO_basisdim()
blocked AO dimension (if symmetry = c1, equivalent to SO_basisdim(); otherwise number of blocks = 1,...
sc::MPQCInit
This helper class simplifies initialization of MPQC.
Definition: mpqcinit.h:21
sc::GetLongOpt::enroll
int enroll(const char *const opt, const OptType t, const char *const desc, const char *const val)
Enroll an option.
sc::auto_vec::release
T * release() MPQC__NOEXCEPT
Release ownership.
Definition: autovec.h:91
sc::auto_vec::get
T * get() const MPQC__NOEXCEPT
Returns the pointer.
Definition: autovec.h:82
sc::GaussianBasisSet::nshell
unsigned int nshell() const
Return the number of shells.
Definition: gaussbas.h:500
sc::RefSymmSCMatrix
The RefSymmSCMatrix class is a smart pointer to an SCSymmSCMatrix specialization.
Definition: matrix.h:265
sc::Integral::overlap
virtual Ref< OneBodyInt > overlap()=0
Return a OneBodyInt that computes the overlap.
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition: matrix.h:135
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::GetLongOpt::retrieve
const char * retrieve(const char *const opt) const
Retrieve an option.
sc::Wavefunction
A Wavefunction is a MolecularEnergy that utilizies a GaussianBasisSet.
Definition: wfn.h:52
sc::auto_vec< double >
sc::OverlapOrthog::basis_to_orthog_basis
RefSCMatrix basis_to_orthog_basis()
Returns a matrix which does the requested transform from a basis to an orthogonal basis.
sc::OverlapOrthog
This class computes the orthogonalizing transform for a basis set.
Definition: orthog.h:38
sc::GetLongOpt::usage
void usage(std::ostream &outfile=std::cout) const
Print usage information.
sc::RefDiagSCMatrix
The RefDiagSCMatrix class is a smart pointer to an DiagSCMatrix specialization.
Definition: matrix.h:389
sc::GetLongOpt::parse
int parse(int argc, char *const *argv)
Parse command line options.
sc::PetiteList::SO_basisdim
RefSCDimension SO_basisdim()
blocked SO dimension (number of blocks = order of the point group, each subdimension has 1 block)
sc::TwoBodyInt::compute_shell
virtual void compute_shell(int, int, int, int)=0
Given four shell indices, integrals will be computed and placed in the buffer.
sc::InputError
This is thrown when invalid input is provided.
Definition: scexception.h:153
sc::FeatureNotImplemented
This is thrown when an attempt is made to use a feature that is not yet implemented.
Definition: scexception.h:120
sc::ThreadGrp::nthread
int nthread() const
The number of threads that will be run by start_thread.
Definition: thread.h:107
sc::SymmOneBodyIntIter
Iterator over symmetry unique shell pairs.
Definition: symmint.h:41
sc::KeyVal::describedclassvalue
Ref< DescribedClass > describedclassvalue(const char *key=0, const KeyValValue &def=KeyValValueRefDescribedClass())
Returns a reference to an object of type DescribedClass.
sc::Integral::petite_list
Ref< PetiteList > petite_list()
Return the PetiteList object.
sc::Integral::set_default_integral
static void set_default_integral(const Ref< Integral > &)
Specifies a new default Integral factory.
sc::GaussianBasisSet::basisdim
RefSCDimension basisdim()
Returns the SCDimension object for the dimension.
Definition: gaussbas.h:495
sc::StateIn
Definition: statein.h:79
sc::ClassDesc
This class is used to contain information about classes.
Definition: class.h:147
sc::PetiteList::to_AO_basis
RefSymmSCMatrix to_AO_basis(const RefSymmSCMatrix &O_so)
converts an operator matrix from SO to AO basis.
sc::MessageGrp::n
int n()
Returns the number of processors.
Definition: message.h:192
sc::FockDistribution
FockDistribution is a factory for constructing the desired FockDist specialization.
Definition: fockdist.h:187
sc::IntegralV3
IntegralV3 computes integrals between Gaussian basis functions.
Definition: intv3.h:41
sc::SymmSCMatrixSCExtrapData
Definition: scextrapmat.h:36
sc::Integral::hcore
virtual Ref< OneBodyInt > hcore()=0
Return a OneBodyInt that computes the core Hamiltonian integrals.
sc::SymmSCMatrixSCExtrapError
Definition: scextrapmat.h:98
sc::GaussianBasisSet::molecule
Ref< Molecule > molecule() const
Return the Molecule object.
Definition: gaussbas.h:489
sc::class_desc
ClassDesc * class_desc()
Return the ClassDesc corresponding to template argument.
Definition: class.h:257
sc::StateOut
Definition: stateout.h:71
sc::TwoBodyInt::buffer
virtual const double * buffer(TwoBodyOper::type type=TwoBodyOper::eri) const
The computed shell integrals will be put in the buffer returned by this member.
sc::RefSCMatrix::t
RefSCMatrix t() const
Return the transpose of this.
sc::FockBuild
The FockBuild class works with the FockBuildThread class to generate Fock matrices for both closed sh...
Definition: fockbuild.h:805
sc::FockBuild::build
void build()
Contruct the Fock matrices.
sc::GaussianBasisSet::matrixkit
Ref< SCMatrixKit > matrixkit()
Returns the SCMatrixKit that is to be used for AO bases.
Definition: gaussbas.h:491
sc::GaussianShell::nfunction
unsigned int nfunction() const
The number of basis functions.
Definition: gaussshell.h:191
sc::CLHFContribution
Computes components of the Fock matrix necessary for closed-shell calculations (i....
Definition: clhfcontrib.h:17
sc::OneBodyIntOp
Definition: obint.h:300
sc::Integral::set_storage
virtual void set_storage(size_t i)
Sets the total amount of storage, in bytes, that is available.
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc::GaussianBasisSet::so_matrixkit
Ref< SCMatrixKit > so_matrixkit()
Returns the SCMatrixKit that is to be used for SO bases.
Definition: gaussbas.h:493
sc::ThreadGrp::get_default_threadgrp
static ThreadGrp * get_default_threadgrp()
Returns the default ThreadGrp.
sc::KeyVal::doublevalue
double doublevalue(const char *key=0, const KeyValValue &def=KeyValValuedouble())
Returns the double value of key.
sc::GaussianBasisSet::shell_to_function
int shell_to_function(int i) const
Return the number of the first function in the given shell.
Definition: gaussbas.h:540
sc::DIIS
The DIIS class provides DIIS extrapolation.
Definition: diis.h:36
sc::GaussianBasisSet::shell
const Shell & shell(int i) const
Return a reference to Shell number i.
Definition: gaussbas.h:558
sc::GetLongOpt
Parse command line options.
Definition: GetLongOpt.h:16
sc::scprintf
This class allows printf-like output to be sent to an ostream.
Definition: formio.h:97
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14

Generated at Sun Jan 26 2020 23:24:02 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.