This is described in the Compiling section of this manual.
The sc-config
program can be use to obtain the compilers, compiler options, and libraries needed to use the SC toolkit from your program. This utility is discussed below, along with how the SC toolkit must be initialized in your main
subroutine.
ExEnv::init(argc, argv);
By default, all output will go to the console stream, cout. To change this, use the following code:
ostream *outstream = new ofstream(outputfilename); ExEnv::set_out(outstream);
MPI is allowed wait until MPI_Init is called to fill in argc and argv, so you may have to call MPI_Init before you even know that we ready to construct MPIMessageGrp. So if an MPIMessageGrp is needed, it is up to the developer to call MPI_Init to get the argument list for certain MPI implementations.
MPI_Init(&argc, &argv);
When files are read and written, an extension is added to a basename to construct the file name. The default is "SC". To use another basename, make the following call, where basename
is a const char *
:
SCFormIO::set_default_basename(basename);
If your job might run in parallel, then make the following call or the nodes will print redundant information. The myproc
argument is the rank of the called node.
SCFormIO::init_mp(myproc);
This segment of code sets up an object to provide multi-threading:
RefThreadGrp thread = ThreadGrp::initial_threadgrp(argc, argv); ThreadGrp::set_default_threadgrp(thread); if (thread.nonnull()) ThreadGrp::set_default_threadgrp(thread); else thread = ThreadGrp::get_default_threadgrp();
This segment of code sets up the message passing object:
RefMessageGrp grp = MessageGrp::initial_messagegrp(argc, argv); if (grp.nonnull()) MessageGrp::set_default_messagegrp(grp); else grp = MessageGrp::get_default_messagegrp();
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 { Ref<OneBodyWavefunction> ref_mp2_wfn_; 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); int spin_polarized(void); 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, class_desc()); } } 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(); } RefSymmSCMatrix MP2::density(void) { throw FeatureNotImplemented("no density yet", __FILE__, __LINE__, class_desc()); return 0; } int MP2::spin_polarized(void) { return 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; }
# Change this to the path to your installed sc-config script.
SCCONFIG = /usr/local/mpqc/current/bin/sc-config
CXX := $(shell $(SCCONFIG) --cxx)
CXXFLAGS := $(shell $(SCCONFIG) --cxxflags)
CPPFLAGS := $(shell $(SCCONFIG) --cppflags)
LIBS := $(shell $(SCCONFIG) --libs)
LIBDIR := $(shell $(SCCONFIG) --libdir)
LTLINK := $(shell $(SCCONFIG) --ltlink)
LTLINKBINOPTS := $(shell $(SCCONFIG) --ltlinkbinopts)
mp2: mp2.o
$(LTLINK) $(CXX) $(CXXFLAGS) -o $@ $^ -L$(LIBDIR) -lmpqc $(LIBS) $(LTLINKBINOPTS)
% 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 ) ) )
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.
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 SC 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.
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 }
g(const sc::Ref<sc::RegionTimer> ®tim) { 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> ®tim) { sc::Timer timer(regtim, "f()"); f(); timer.reset(); // or let the destructor exit the region }
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.
$ 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.
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 src/bin/mpqc/validate
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 src/bin/mpqc/validate
directory for examples.
After you have added new inputs and modified the Makefile, change into the src/bin/mpqc/validate
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
.