Main Page Class Hierarchy Compound List File List Compound Members Related Pages

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

- The sc-config Program
- Initializing SC
- MP2 Implementation Example: Source
- MP2 Implementation Example: Makefile
- MP2 Implementation Example: Input

### 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 MP2 Implementation Example: Makefile 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 MP2 Implementation Example: Source section. The sc-config command 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)

    mp2: mp2.o
            $(CXX) $(CXXFLAGS) -o $@ $^ -L$(LIBDIR) -lmpqc $(LIBS)

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