develop(3)

NAME
develop - Developing Code Using SC
 In addition to the executables, the Scientific Computing toolkit libraries and include files can be installed on your machine. This is described in the 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.

  The sc-config Program

The sc-config program returns information about how SC was compiled and installed. The following information is available:

--prefix
     The directory where SC is installed.

--version
     The version of SC.

--libdir
     The directory were the libraries are found.

--libs
     The libraries and library paths needed to link.

--cppflags
     The include directories.

--cc
     The C compiler.

--cflags
     The C compiler flags.

--cxx
     The C++ compiler.

--cxxflags
     The C++ compiler flags.

--f77
     The FORTRAN 77 compiler.

--f77flags
     The FORTRAN 77 compiler flags.

To use the sc-config program to link your executable to SC, use a Makefile for GNU make similar to the following:

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)

  Initializing SC

First the execution environment must be initialized using the ExEnv init member.

  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 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();

  MP2 Implementation Example: Source

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.

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 section.

#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()) {
    ExEnv::out0() << 'reference is null' << endl;
    abort();
  }
}
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()) {
    ExEnv::out0() << 'No gradients yet' << endl;
    abort();
  }
  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) {
  ExEnv::out0() << 'No density yet' << endl;
  abort();
  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) {
    ExEnv::out0() << 'C1 symmetry only' << endl;
    abort();
  }
  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;
  double *cvec = new double [vec.nrow() * vec.ncol()];
  vec->convert(cvec);
  double *pqrs = new double [nao * nao * nao * nao];
  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;
  double *ijkl = new double [nmo * nmo * nmo * nmo];
  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];
                }
              }
            }
          }
        }
      }
    }
  }
  delete [] pqrs;
  delete [] cvec;
  double *evals = new double [nmo];
  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]);
        }
      }
    }
  }
  delete [] ijkl;
  delete [] evals;
  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 section. The 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.

# 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)

  MP2 Implementation Example: Input

This input file can be used with the program illustrated in the section. It will compute the MP2 energy using the new MP2 class. Note that only the 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
    )
  )
)