<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<html><head><meta name="robots" content="noindex">
<meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>Developing Code Using SC</title>
<link href="doxygen.css" rel="stylesheet" type="text/css">
</head><body bgcolor="#ffffff">
<!-- Generated by Doxygen 1.2.5 on Mon Oct 14 14:18:08 2002 -->
<center>
<a class="qindex" href="index.html">Main Page</a> <a class="qindex" href="hierarchy.html">Class Hierarchy</a> <a class="qindex" href="annotated.html">Compound List</a> <a class="qindex" href="files.html">File List</a> <a class="qindex" href="functions.html">Compound Members</a> <a class="qindex" href="pages.html">Related Pages</a> </center>
<hr><a name="develop"><h2>Developing Code Using SC</h2></a>
In addition to the executables, the Scientific Computing toolkit libraries and include files can be installed on your machine. This is described in the <a href="compile.html#compile">Compiling</a> section of this manual.
<p>
The <code>sc-config</code> 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 <code>main</code> subroutine.
<p>
<ul>
<li> <a href="develop.html#scconfig">The sc-config Program</a> <li> <a href="develop.html#scinit">Initializing SC</a> <li> <a href="develop.html#devsampsrc">MP2 Implementation Example: Source</a> <li> <a href="develop.html#devsampmak">MP2 Implementation Example: Makefile</a> <li> <a href="develop.html#devsampinp">MP2 Implementation Example: Input</a> </ul>
<p>
<a name="scconfig"><h3>The sc-config Program</h3></a>
<p>
The sc-config program returns information about how SC was compiled and installed. The following information is available:
<p>
<dl compact>
<dt><code>--prefix</code><dd> The directory where SC is installed. <dt><code>--version</code><dd> The version of SC. <dt><code>--libdir</code><dd> The directory were the libraries are found. <dt><code>--libs</code><dd> The libraries and library paths needed to link. <dt><code>--cppflags</code><dd> The include directories. <dt><code>--cc</code><dd> The C compiler. <dt><code>--cflags</code><dd> The C compiler flags. <dt><code>--cxx</code><dd> The C++ compiler. <dt><code>--cxxflags</code><dd> The C++ compiler flags. <dt><code>--f77</code><dd> The FORTRAN 77 compiler. <dt><code>--f77flags</code><dd> The FORTRAN 77 compiler flags.
</dl>
<p>
To use the sc-config program to link your executable to SC, use a Makefile for GNU make similar to the following:
<p>
<pre>
SCCONFIG = /usr/local/mpqc/current/bin/sc-config
CXX := $(shell $(SCCONFIG) --cxx)
CXXFLAGS := $(shell $(SCCONFIG) --cxxflags)
CPPFLAGS := $(shell $(SCCONFIG) --cppflags)
LIBS := $(shell $(SCCONFIG) --libs)
myprog: myprog.o
$(CXX) $(CXXFLAGS) -o $@ $^ $(LIBS)
</pre>
<p>
<a name="scinit"><h3>Initializing SC</h3></a>
<p>
First the execution environment must be initialized using the ExEnv init member.
<p>
<pre>
ExEnv::init(argc, argv);
</pre>
<p>
By default, all output will go to the console stream, cout. To change this, use the following code:
<p>
<pre>
ostream *outstream = new ofstream(outputfilename);
ExEnv::set_out(outstream);
</pre>
<p>
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.
<p>
<pre>
MPI_Init(&argc, &argv);
</pre>
<p>
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 <code>basename</code> is a <code>const char *</code>:
<p>
<pre>
SCFormIO::set_default_basename(basename);
</pre>
<p>
If your job might run in parallel, then make the following call or the nodes will print redundant information. The <code>myproc</code> argument is the rank of the called node.
<p>
<pre>
SCFormIO::init_mp(myproc);
</pre>
<p>
This segment of code sets up an object object to provide multi-threading:
<p>
<pre>
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();
</pre>
<p>
This segment of code sets up the message passing object:
<p>
<pre>
RefMessageGrp grp = MessageGrp::initial_messagegrp(argc, argv);
if (grp.nonnull())
MessageGrp::set_default_messagegrp(grp);
else
grp = MessageGrp::get_default_messagegrp();
</pre>
<p>
<a name="devsampsrc"><h3>MP2 Implementation Example: Source</h3></a>
<p>
This example code illustrates a complete MP2 energy implementation using the SC 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.
<p>
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 <a href="develop.html#devsampmak">MP2 Implementation Example: Makefile</a> section.
<p>
<div class="fragment"><pre>
<font class="preprocessor">#include <chemistry/qc/wfn/obwfn.h></font>
<font class="preprocessor">#include <chemistry/qc/scf/clhf.h></font>
<font class="keyword">using</font> <font class="keyword">namespace </font>std;
<font class="keyword">using</font> <font class="keyword">namespace </font>sc;
<font class="keyword">class </font>MP2: <font class="keyword">public</font> Wavefunction {
Ref<OneBodyWavefunction> ref_mp2_wfn_;
<font class="keywordtype">double</font> compute_mp2_energy();
<font class="keyword">public</font>:
MP2(<font class="keyword">const</font> Ref<KeyVal> &);
MP2(StateIn &);
<font class="keywordtype">void</font> save_data_state(StateOut &);
<font class="keywordtype">void</font> compute(<font class="keywordtype">void</font>);
<font class="keywordtype">void</font> obsolete(<font class="keywordtype">void</font>);
<font class="keywordtype">int</font> nelectron(<font class="keywordtype">void</font>);
RefSymmSCMatrix density(<font class="keywordtype">void</font>);
<font class="keywordtype">int</font> spin_polarized(<font class="keywordtype">void</font>);
<font class="keywordtype">int</font> value_implemented(<font class="keywordtype">void</font>) <font class="keyword">const</font>;
};
<font class="keyword">static</font> ClassDesc MP2_cd(<font class="keyword">typeid</font>(MP2), <font class="stringliteral">"MP2"</font>, 1, <font class="stringliteral">"public Wavefunction"</font>,
0, create<MP2>, create<MP2>);
MP2::MP2(<font class="keyword">const</font> Ref<KeyVal> &keyval):Wavefunction(keyval) {
ref_mp2_wfn_ << keyval->describedclassvalue(<font class="stringliteral">"reference"</font>);
<font class="keywordflow">if</font>(ref_mp2_wfn_.null()) {
ExEnv::out0() << <font class="stringliteral">"reference is null"</font> << endl;
abort();
}
}
MP2::MP2(StateIn &statein):Wavefunction(statein)
{
ref_mp2_wfn_ << SavableState::restore_state(statein);
}
<font class="keywordtype">void</font>
MP2::save_data_state(StateOut &stateout)<font class="keyword"> </font>{
Wavefunction::save_data_state(stateout);
SavableState::save_state(ref_mp2_wfn_.pointer(),stateout);
}
<font class="keywordtype">void</font>
MP2::compute(<font class="keywordtype">void</font>)<font class="keyword"></font>
<font class="keyword"></font>{
<font class="keywordflow">if</font>(gradient_needed()) {
ExEnv::out0() << <font class="stringliteral">"No gradients yet"</font> << endl;
abort();
}
<font class="keywordtype">double</font> extra_hf_acc = 10.;
ref_mp2_wfn_->set_desired_value_accuracy(<a class="code" href="class_sc__Function.html#a15">desired_value_accuracy</a>()
/ extra_hf_acc);
<font class="keywordtype">double</font> refenergy = ref_mp2_wfn_->energy();
<font class="keywordtype">double</font> mp2energy = compute_mp2_energy();
ExEnv::out0() << indent << <font class="stringliteral">"MP2 Energy = "</font> << mp2energy << endl;
set_value(refenergy + mp2energy);
set_actual_value_accuracy(ref_mp2_wfn_->actual_value_accuracy()
* extra_hf_acc);
}
<font class="keywordtype">void</font>
MP2::obsolete(<font class="keywordtype">void</font>)<font class="keyword"> </font>{
Wavefunction::obsolete();
ref_mp2_wfn_->obsolete();
}
<font class="keywordtype">int</font>
MP2::nelectron(<font class="keywordtype">void</font>)<font class="keyword"> </font>{
<font class="keywordflow">return</font> ref_mp2_wfn_->nelectron();
}
RefSymmSCMatrix
MP2::density(<font class="keywordtype">void</font>)<font class="keyword"> </font>{
ExEnv::out0() << <font class="stringliteral">"No density yet"</font> << endl;
abort();
<font class="keywordflow">return</font> 0;
}
<font class="keywordtype">int</font>
MP2::spin_polarized(<font class="keywordtype">void</font>)<font class="keyword"> </font>{
<font class="keywordflow">return</font> 0;
}
<font class="keywordtype">int</font>
MP2::value_implemented(<font class="keywordtype">void</font>)<font class="keyword"> const </font>{
<font class="keywordflow">return</font> 1;
}
<font class="keywordtype">double</font>
MP2::compute_mp2_energy()<font class="keyword"></font>
<font class="keyword"></font>{
<font class="keywordflow">if</font>(molecule()->point_group()->char_table().order() != 1) {
ExEnv::out0() << <font class="stringliteral">"C1 symmetry only"</font> << endl;
abort();
}
RefSCMatrix vec = ref_mp2_wfn_->eigenvectors();
<font class="keywordtype">int</font> nao = vec.nrow();
<font class="keywordtype">int</font> nmo = vec.ncol();
<font class="keywordtype">int</font> nocc = ref_mp2_wfn_->nelectron()/2;
<font class="keywordtype">int</font> nvir = nmo - nocc;
<font class="keywordtype">double</font> *cvec = <font class="keyword">new</font> <font class="keywordtype">double</font> [vec.nrow() * vec.ncol()];
vec->convert(cvec);
<font class="keywordtype">double</font> *pqrs = <font class="keyword">new</font> <font class="keywordtype">double</font> [nao * nao * nao * nao];
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> n = 0; n < nao*nao*nao*nao; n++) pqrs[n] = 0.0;
Ref<TwoBodyInt> twoint = <a class="code" href="class_sc__Wavefunction.html#a29">integral</a>()->electron_repulsion();
<font class="keyword">const</font> <font class="keywordtype">double</font> *buffer = twoint->buffer();
Ref<GaussianBasisSet> basis = this-><a class="code" href="class_sc__Wavefunction.html#a28">basis</a>();
<font class="keywordtype">int</font> nshell = basis->nshell();
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> P = 0; P < nshell; P++) {
<font class="keywordtype">int</font> nump = basis->shell(P).nfunction();
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> Q = 0; Q < nshell; Q++) {
<font class="keywordtype">int</font> numq = basis->shell(Q).nfunction();
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> R = 0; R < nshell; R++) {
<font class="keywordtype">int</font> numr = basis->shell(R).nfunction();
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> S = 0; S < nshell; S++) {
<font class="keywordtype">int</font> nums = basis->shell(S).nfunction();
twoint->compute_shell(P,Q,R,S);
<font class="keywordtype">int</font> index = 0;
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> p=0; p < nump; p++) {
<font class="keywordtype">int</font> op = basis->shell_to_function(P)+p;
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> q = 0; q < numq; q++) {
<font class="keywordtype">int</font> oq = basis->shell_to_function(Q)+q;
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> r = 0; r < numr; r++) {
<font class="keywordtype">int</font> oor = basis->shell_to_function(R)+r;
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> s = 0; s < nums; s++,index++) {
<font class="keywordtype">int</font> os = basis->shell_to_function(S)+s;
<font class="keywordtype">int</font> ipqrs = (((op*nao+oq)*nao+oor)*nao+os);
pqrs[ipqrs] = buffer[index];
}
}
}
}
}
}
}
}
twoint = 0;
<font class="keywordtype">double</font> *ijkl = <font class="keyword">new</font> <font class="keywordtype">double</font> [nmo * nmo * nmo * nmo];
<font class="keywordtype">int</font> idx = 0;
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> i = 0; i < nmo; i++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> j = 0; j < nmo; j++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> k = 0; k < nmo; k++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> l = 0; l < nmo; l++, idx++) {
ijkl[idx] = 0.0;
<font class="keywordtype">int</font> index = 0;
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> p = 0; p < nao; p++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> q = 0; q < nao; q++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> r = 0; r < nao; r++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> 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];
}
}
}
}
}
}
}
}
<font class="keyword">delete</font> [] pqrs;
<font class="keyword">delete</font> [] cvec;
<font class="keywordtype">double</font> *evals = <font class="keyword">new</font> <font class="keywordtype">double</font> [nmo];
ref_mp2_wfn_->eigenvalues()->convert(evals);
<font class="keywordtype">double</font> energy = 0.0;
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> i=0; i < nocc; i++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> j=0; j < nocc; j++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> a=nocc; a < nmo; a++) {
<font class="keywordflow">for</font>(<font class="keywordtype">int</font> b=nocc; b < nmo; b++) {
<font class="keywordtype">int</font> iajb = (((i*nmo+a)*nmo+j)*nmo+b);
<font class="keywordtype">int</font> ibja = (((i*nmo+b)*nmo+j)*nmo+a);
energy += (2 * ijkl[iajb] - ijkl[ibja]) * ijkl[iajb]/
(evals[i] + evals[j] - evals[a] - evals[b]);
}
}
}
}
<font class="keyword">delete</font> [] ijkl;
<font class="keyword">delete</font> [] evals;
<font class="keywordflow">return</font> energy;
}
</div></pre>
<p>
<a name="devsampmak"><h3>MP2 Implementation Example: Makefile</h3></a>
<p>
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 <a href="develop.html#devsampsrc">MP2 Implementation Example: Source</a> section. The <a href="develop.html#scconfig">sc-config command</a> is used to obtain information about how the SC toolkit was compiled and installed. The library specified with -lmpqc provides the main routine from mpqc.
<p>
<div class="fragment"><pre><font class="preprocessor"># Change this to the path to your installed sc-config script.</font>
<font class="preprocessor"></font>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)
mp2: mp2.o
$(CXX) $(CXXFLAGS) -o $@ $^ -L$(LIBDIR) -lmpqc $(LIBS)
</div></pre>
<p>
<a name="devsampinp"><h3>MP2 Implementation Example: Input</h3></a>
<p>
This input file can be used with the program illustrated in the <a href="develop.html#devsampsrc">MP2 Implementation Example: Source</a> section. It will compute the MP2 energy using the new MP2 class. Note that only the <a href="mpqcoo.html#mpqcoo">object-oriented input format</a> can be used with user provided classes.
<p>
<div class="fragment"><pre>% 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 = <font class="stringliteral">"STO-3G"</font>
molecule = $:molecule
)
mpqc: (
checkpoint = no
savestate = no
% MP2 is the <font class="keyword">new</font> <font class="keyword">class</font>. 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
)
)
)
</div></pre>
<p>
<hr>
<address>
<small>
Generated at Mon Oct 14 14:18:08 2002 for <a
href="http://aros.ca.sandia.gov/~cljanss/mpqc">MPQC</a>
2.1.2 using the documentation package <a
href="http://www.stack.nl/~dimitri/doxygen/index.html">Doxygen</a>
1.2.5.
</small>
</address>
</body>
</html>
distrib
>
Mandriva
>
9.1
>
ppc
>
media
>
contrib
>
by-pkgid
>
263386785cefb9ae5d63b926d214d809
>
files
>
941
