User Tools

Site Tools


crawdad:programming:project7

Programming Project #7: Interfacing Your Code to PSI4

PSI4: Creating a Plugin

Project 7 provides instructions for creating a plugin in PSI4. If, for historical reasons, you'd like to view the PSI3 version, click here.

The purpose of this project is to help you to connect your Hartree-Fock self-consistent field program to the PSI4 program package by creating a plugin. Instead of reading the requisite integrals and other data from text files, you will be computing them using PSI4. This will allow you to try your own molecules as test cases.

A prerequisite to this project is that you must either have a working copy of PSI4 or access to one, including the corresponding libraries and header files. Furthermore, the copy of PSI4 you use must be compiled with the --with-plugins flag (see the link below for more information about flags when compiling). Hokies can access a stable copy on cerebro by adding /cerebro/home/software/psi4/bin/ to your path. Instructions for obtaining and compiling PSI4, along with much more helpful information, can be found at the PSI4 Documentation page. General PSI4 information, with a link to documentation, can be found at psicode.org.

Note that, in order to make use of PSI4's libraries, your program must be compatible with a C++ compiler. Hence, some of the material in this section involves C++ directives and syntax. However, your existing C programs can be easily adjusted for this purpose with no significant structural changes. Furthermore, hopefully this project will ease you into familiarity with aspects of C++ object-oriented programming.

An Introduction

The Structure of a PSI4 Plugin

The following is a generic/blank plugin that allows for easy interfacing with PSI4. In addition to a skeleton for your code, creating a plugin also provides you with a Makefile and sample input.

To execute the plugin, you need only type make and then run psi4.

#include <libplugin/plugin.h>
#include <psi4-dec.h>
#include <libparallel/parallel.h>
#include <liboptions/liboptions.h>
#include <libmints/mints.h>
#include <libpsio/psio.hpp>

INIT_PLUGIN

using namespace boost;

namespace psi{ namespace some_plugin {

extern "C"
int read_options(std::string name, Options& options)
{
    if (name == "SOME_PLUGIN"|| options.read_globals()) {
        /*- The amount of information printed to the output file -*/
        options.add_int("PRINT", 1);
        /*- Add your own options here -*/
    }

    return true;
}

extern "C"
PsiReturnType plain_plugin(Options& options)
{
    int print = options.get_int("PRINT");

    /* Your code goes here */

    return Success;
}

}} // End namespaces

                                      

Points to note about the above code:

  • The list of “#include” directives at the top of the code includes some major PSI4 library headers, but you may need to add more depending on the scope of your projects in the future. (You shouldn't need to add any for your scf plugin.)
  • PSI4 is a single executable when it is run. Fortunately, the Makefile provided with the plugin allows you to recompile your code separately from PSI4, but still be able to link against its capabilities.
  • The read options routine allows you to create options and keywords that can be set in the input file. The print example in the skeleton above demonstrates first adding an option (and setting its default value), and then accessing the option with the options.get_int function. Of course, you can add different types of options (like integers, double, booleans, strings). Learn more by checking out PSI4 page on handling Options.

The plugin's Makefile

Below is the top snippet of the automatically generated Makefile from a PSI4 plugin. You shouldn't need to edit this file, but do pay attention to the top_srcdir and top_objdir variables. You may potentially need to edit these if you move your plugin directory around.

#
# Plugin Makefile generated by Psi4.
#
# You shouldn't need to modify anything in this file.
#
 
# Location of your PSI4 source
top_srcdir = /path/to/psi4/
# Location of your PSI4 install, by default as listed
top_objdir = /path/to/psi4/objdir
 
# Start by figuring out whether we're on Linux or Mac (sorry, Mr. Gates)
UNAME := $(shell uname)
 
include $(top_objdir)/src/bin/MakeVars

Points to note about the above Makefile:

  • The file is tailor-made just for you and your plugin - the path/to/psi4 lines should be automatically set to the appropriate path.
  • A general introduction to Makefiles can be found here.

PSI4 Input File

As mentioned before, creating a plugin provides you with a sample input so you can run your plugin immediately. The input structure has been a major focus in the development of PSI4 and, in the venture to make it as user friendly as possibly, has become a Python script. General details on the input file can be found elsewhere. As the input suggests, the plugin is loaded at runtime. The first set keyword allows input of global variables, such as the basis set or convergence criteria. set some_plugin allows you to input desired values for the options specific to your plugin. Options for other “modules” can be specified in a similar manner.

For a more straightforward transition of your SCF code into the plugin, it is perhaps simplest if we ignore symmetry to begin with. To accomplish this, simply specify symmetry c1 in your input file (see below). Also notice that the pound sign “#” is used for commenting out in the input.

# PYTHONPATH must include directory above plugin directory.
#     Define either externally or here, then import plugin.
sys.path.insert(0, './..')
import some_plugin

molecule {
O
H 1 R
H 1 R 2 A

R = .9
A = 104.5

symmetry c1
}

set {
  basis sto-3g
}

set some_plugin {
  print 1
}

energy('some_plugin')

some_plugin.exampleFN()                             

The top two lines involve appropriately handling your plugin so that Python knows how to deal with it. The line with energy('some_plugin') reveals that your plugin has been hooked to the main PSI4 driver function energy(). The final line of the input illustrates how you would execute Python functions you wrote for your plugin (that you would define in pymodule.py. However, it can certainly be ignored for now. More “plugin relevant” input file information can be found here.

Other Plugin Files

Python Files

Python (the scripting language) has become an integral component of PSI4, especially in regards to the input file (which actually is a Python script!). The plugin you created will also come with two Python files: __init__.py and pymodule.py. In short, these files are necessary for enabling the Python components of the plugin, and connecting the plugin “pythonically” to PSI4. A more detailed explanation can be found here.

Documentation

The doc.rst file is incorporated in the plugin looking to the future - that is, perhaps your plugin involves some new functionality that you'd like to add to PSI4. The doc.rst file helps incorporate the documentation with the new code you have to offer. For more information, go here.

Creating a Plugin

Now that you've been briefly introduced to PSI4 plugins, it's time for you to make your own.

AO Integrals Template

Type the following command to create a plugin called myscf_plugin that computes the AO integrals using PSI4:

psi4 --new-plugin myscf_plugin +aointegrals

A few notes:

  • +aointegrals provides you with the integral objects you'll need to write your SCF code. FYI, there are other templates available, such as MO integrals (mointegrals).
  • The above command should have created a directory called myscf_plugin. There should be six files in the new directory:
    1. input.dat
    2. Makefile
    3. myscf_plugin.cc
    4. __init__.py
    5. pymodule.py
    6. doc.rst
  • You can create your plugin wherever you like (it does not have to be in PSI4 source directory; in fact, outside of the PSI4 source directory is safest).
  • If the above command failed, you may need to recompile your copy of PSI4 with the --with-plugins option.

Don't forget to set symmetry to C1 in the input file:

molecule {
  O
  H 1 R
  H 1 R 2 A

  R = .9
  A = 104.5

  symmetry c1
}

Objects, Classes, and Member Functions (oh my)

Many of the common quantities and procedures in PSI3 have been implemented in an object-oriented manner in PSI4. Creating an object, which belongs to a particular class, will have member functions associated with it by virtue of its class. A simple example is the matrix object (or shared matrix object). Typical matrix manipulations, such as printing, diagonalizing, adding, etc., can be conveniently handled through the Matrix Class's member functions. Assuming you've included the <libmints/matrix.h> in your code, here are some of the things you can do:

/* Here we declare a Matrix object X (similar to int x),
 * and it takes arguments for its name and dimensions
 */
Matrix X("A matrix", dim1, dim2);
SharedMatrix Y(new Matrix("Another matrix", dim1, dim2);

/* "Instantiating" the Matrix object conveniently zeroes the elements,
 * and so we should be able to print a zero matrix no problem
 * using the print member function
 */
X.print();   //Access object member functions with "."
Y->print();  //Access shared object member functions with "->"

Now you've seen a basic example of objects and member functions in action. The wavefunction and molecules objects in particular can tell you lots of helpful information about your system (number of basis functions, molecular orbitals, etc.). Details on objects, the classes to which they belong, and their member functions can be found in the appropriate library in the source code. If you have no idea where to look for an object, try starting with libmints. Once in the right library, the header files (molecule.h for example) are most helpful for quickly finding an object, what arguments it takes, and necessary details on member functions.

The integral and matrix factories are helpful tools that take care of a lot of work you would normally have to do yourself.

Grabbing Basic Info

As mentioned before, certain objects allow you access through them to useful information. As an example, let's get PSI4 to tell us the number of molecular orbitals and the number of doubly occupied orbitals using the wavefunction object.

...
/* Making this wavefunction object allows us to access Wavefunction member objects
 * using "wfn" rather than Process::environment.wavefunction() everytime
 * we need it.
 */
boost::shared_ptr<Wavefunction> wfn = Process::environment.wavefunction();
...

//Basic Info
int nmo = wfn->nmo();  /* Less convenient, but equivalent would be Process::environment.wavefunction()->nmo(); */
int *doccpi = wfn->doccpi();
int docc = doccpi[0];

Here we declare a shared pointer (object) of type Wavefunction called wfn, and then we put into wfn the reference wavefunction. The syntax is directly analagous to what we do when we set int nmo = wfn->nmo.

One last note is that the function doccpi() returns an array of the doubly occupied orbitals per irrep. Since we are not yet taking advantage of symmetry (remember to include symmetry c1 in the input file), there should only be one irrep. Taking the zeroth element of the array we created (doccpi[0]) gives us the total number of doubly occupied orbitals.

Reading One- and Two-Electron Integrals

One-Electron Integrals

The overlap, kinetic-energy, and nuclear-attraction integrals are symmetric with respect to the basis-function indices, e.g.

$$S_{\mu \nu} = S_{\nu\mu}
$$

To save disk space, only the unique integrals are computed (the lower triangle, as it were). That is, for N basis functions, only N(N+1)/2 one-electron integrals of each type are computed. This is what is happening with OneBodyAOInt objects (sOBI, tOBI, and vOBI). Utilizing the compute function, the integrals are stored in convenient, full N x N matrix form (sMat, tMat, and vMat).

The Object Logic to Get the Overlap Matrix

Let's walk through the generation of the overlap matrix (sMat) to get a feel for what's going on with the integral objects, matrix factory objects, and matrix objects.

// The matrix factory can create matrices of the correct dimensions...
boost::shared_ptr<MatrixFactory> factory(new MatrixFactory);
factory->init_with(1, nbf, nbf);

// Form the Overlap Matrix
SharedMatrix sMat(factory->create_matrix("Overlap"));

The matrix factory (which here is named factory) is not necessary, but only a tool of convenience. As we already know, we have several matrices (sMat, tMat, kMat) that will be the same size (AO by AO). Rather than setting them all up explicitly, we are using the factory here to keep up with the nitty-gritty details, and just use the factory to create the desired sized matrix. All we have to be sure of is that the matrix factory factory will set up a matrix with the right dimensions. An alternative, equally valid approach for creating the overlap matrix (sMat) would be:

SharedMatrix sMat(new Matrix("Overlap", nbf, nbf));

So far, we've only declared and allocated memory for our overlap matrix. We want to actually put the integral values into it now. First, we need to compute the integrals. Enter the one body or one electron integral object.

// The integral factory oversees the creation of integral objects
boost::shared_ptr<IntegralFactory> integral(new IntegralFactory
        (aoBasis, aoBasis, aoBasis, aoBasis));
// Form the one-electron integral objects from the integral factory
boost::shared_ptr<OneBodyAOInt> sOBI(integral->ao_overlap());

Notice, once again we've used a factory object (in this case an integral factory) to construct our overlap OneBodyAOInt object sOBI. At this point, though we have an overlap integral object, we still haven't computed the integrals yet. To do this, we use the OneBodyAOInt object's compute function:

// Compute the one electron integrals, telling each object where to store the result
sOBI->compute(sMat);

Basic Matrix Manipulations

Now that we have the one electron integrals in matrix form, it might be useful to know how to access and manipulate their individual elements. In your previous, C code, things looked like this:

// Matrix allocation
double **X;
X= (double **) malloc(dim * sizeof(double *));
for(int i=0; i < dim; ++i)  {
  X[i] = (double *) malloc(dim * sizeof(double));
}

// Matrix multiplication
for(i=0; i<dim; i++)  {
  for(j=0; j<dim; j++)  {
    X[i][j] = 0.0;
    for(k=0; k<dim; k++)  {
      X[i][j] += Y[i][k] * Z[k][j];

Within PSI4, C++ style, we would do this:

// Shared Matrices
SharedMatrix X(new Matrix(dim,dim));
for(i=0; i<dim; i++)  {
  for(j=0; j<dim; j++)  {
    for(k=0; k<dim; k++)  {
      X->add(i,j, Y->get(i,k) * Z->get(k,j));
      
// Matrices
Matrix A(dim,dim);
for(i=0; i<dim; i++)  {
  for(j=0; j<dim; j++)  {
    for(k=0; k<dim; k++)  {
      A.add(i,j, B.get(i,k) * C.get(k,j));

Alternatively, we could have used the Matrix class gemm function:

X->gemm(1, 0, 1, Y, Z, 0);

where looking in libmints/matrix.h tells me about gemm - what arguments it takes and what they mean.

/** General matrix multiply, saves result to this
     * \param transa Transpose the left matrix
     * \param transb Transpose the right matrix
     * \param alpha Prefactor for the matrix multiplication
     * \param a Left matrix
     * \param b Right matrix
     * \param beta Prefactor for the resulting matrix
     */
    void gemm(bool transa, bool transb, double alpha, const Matrix* const a, const Matrix* const b, double beta);
    void gemm(bool transa, bool transb, double alpha, const SharedMatrix& a, const SharedMatrix& b, double beta);
    void gemm(bool transa, bool transb, double alpha, const SharedMatrix& a, const Matrix& b, double beta);
    void gemm(bool transa, bool transb, double alpha, const Matrix& a, const SharedMatrix& b, double beta);
*/

The four different listings of void gemm are an example of overloading a function. As a result, gemm can easily accomodate two Shared Matrix objects, two Matrix objects, or one of each without a problem.

Two-Electron Integrals

The two-electron integrals are stored using Mulliken-ordered indices with full, eight-fold permutational symmetry. The plugin currently computes and prints them for you. One course of action is to grab them as they are computed and put them into a vector. The code below illustrates how to form a properly indexed TEI array with minimal work:

    SharedVector TEI(new Vector("TEI array", ntri*(ntri+1)/2));
      
    // Now, the two-electron integrals
    boost::shared_ptr<TwoBodyAOInt> eri(integral->eri());
    // The buffer will hold the integrals for each shell, as they're computed
    const double *buffer = eri->buffer();
    // The iterator conveniently lets us iterate over functions within shells
    AOShellCombinationsIterator shellIter = integral->shells_iterator();
    int count=0;
    for (shellIter.first(); shellIter.is_done() == false; shellIter.next()) {
        // Compute quartet
        eri->compute_shell(shellIter);
        // From the quartet get all the integrals
        AOIntegralsIterator intIter = shellIter.integrals_iterator();
        for (intIter.first(); intIter.is_done() == false; intIter.next()) {
            int p = intIter.i();
            int q = intIter.j();
            int r = intIter.k();
            int s = intIter.l();
                
            fprintf(outfile, "\t(%2d %2d | %2d %2d) = %20.15f\n",
                p, q, r, s, buffer[intIter.index()]);
            ++count;
                
            //GRABBING TEI's HERE
            /* The buffer[intIter.index()] value is the unique pqrs integral
             * so we just have to put it in our array properly
             */                
            int pq = INDEX(p,q);
            int rs = INDEX(r,s);
            int pqrs = INDEX(pq,rs);
            TEI->set(pqrs, buffer[intIter.index()]); 
          
        }
    }
    fprintf(outfile, "\n\tThere are %d unique integrals\n\n", count);

Other Useful PSI Features

The PSI libraries contain a number of other useful objects/functions that you may need in your programs. This section, along with further demonstrating the handling and use of objects/member functions, highlights a few that are particularly useful for your SCF program accessible through molecule, wavefunction, and matrix objects.

The Molecule Object

You'll definitely need the nuclear repulsion energy. Here's some simple examples of how to obtain it.

Without a molecule object:

// Access the nuclear repulsion energy directly from Process::environment.molecule()
double nucrep = Process::environment.molecule()->nuclear_repulsion_energy();

With a molecule object:

// Create a molecule object
// (Think of this as creating a boost shared pointer of type Molecule, called molecule)
boost::shared_ptr<Molecule> molecule;

// Set your new object to the global molecule
molecule = Process::environment.molecule();

// Use the molecule member function (accessed with an "arrow" 
// since molecule is a "shared" object)
double nucrep = molecule->nuclear_repulsion_energy();

Or equivalently, in two lines:

boost::shared_ptr<Molecule> molecule = Process::environment.molecule();
double nucrep = molecule->nuclear_repulsion_energy();

A few other, noteworthy functions you may or may not need for an scf plugin (taken directly from libmints/molecule.h):

/// Returns the geometry in a Matrix
Matrix geometry() const;
/// Number of atoms
int natom() const { return atoms_.size(); }
/// Print the molecule
void print() const;

As further illustration of handling objects and member functions, let's say we want to print the geometry of the molecule. The above member functions of the molecule class suggest three ways of doing so:

Perhaps the simplest approach is

// Use the molecule object's print function
molecule->print();


Or, since molecule->geometry() returns a Matrix, we can utilize the Matrix object's print function on top of the it.

// Print the geometry using the Matrix object's print function
molecule->geometry().print();


Lastly, we could put the geometry into our own Matrix object, rather than printing the molecule object's member geometry.

// Put the geometry into a matrix
Matrix my_geom = molecule->geometry();
// Print my_geom using the matrix object's print function
my_geom.print();

The Wavefunction Object

The wavefunction class is another area of good, general information. See libmints/wavefunction.h for more on the class.

// Instantiate a wavefunction object, and assign it the global wavefunction.
boost::shared_ptr<Wavefunction> wfn = Process::environment.wavefunction();

NOTE: The above snippet assumes the wavefunction already exists (e.g. SCF has been run already, in which case the wavefunction would be SCF). This is actually the case for your new myscf_plugin. Type make and the run psi4 in your new plugin directory. Look at the output file and you'll see that a Hartree Fock calculaiton has been run prior to the formation and printing of sMat, tMat, and vMat.

Some useful functions (as presented in libmints/wavefunction.h):

/// Returns the basis set object that pertains to this wavefunction.
boost::shared_ptr<BasisSet> basisset() const;
/// Returns the DOCC per irrep array.
const Dimension& doccpi() const { return doccpi_; }
/// Returns the number of SOs
int nso() const { return nso_; }
/// Returns the number of MOs
int nmo() const { return nmo_; }
/// Returns the reference energy
double reference_energy () const { return energy_; }
/// Returns the overlap matrix
SharedMatrix S() const { return S_; }
/// Returns the alpha electrons MO coefficients
SharedMatrix Ca() const;
/// Returns the (SO basis) alpha Fock matrix
SharedMatrix Fa() const;

Process::environment

Something we've ignored so far is the Process::environment business. The short explanation is to think of Process::environment as a special object that gives you access to essential items (like the wavefunction and molecule) through its member functions.

It also provides a means to grab global quantities. For instance, say we are working on some post-Hartree Fock method, and need the HF energy. We can grab it like so:

// Grab the reference SCF total energy
double ref_scf = Process::environment.globals["SCF TOTAL ENERGY"];

For a listing of all such variables accessible through Process::environment, click here.

The Matrix Object

You've already seen the Matrix/SharedMatrix in action above. libmints/matrix.h is the source for Matrix functions, arguments, etc. Rather than go into detail, a simple list of expected functions is given:

  • set
  • zero
  • add
  • copy
  • gemm
  • diagonlize
crawdad/programming/project7.txt · Last modified: 2013/06/26 16:25 by harley