Introduction

This manual documents the different functions in the PyQuante suite of Quantum Chemistry programs. See also the information on the PyQuante Home Page.

Molecule Objects

PyQuante programs use the Molecule object to contain the information about the molecule - the atoms, the charge, and the multiplicity. The docstring information for the Molecule object is:

Molecule(name,atomlist,**opts) - Object to hold PyQuante molecular data

name       The name of the molecule
atomlist   A list of (atno,(x,y,z)) information for the molecule

Options:      Value   Description
--------      -----   -----------
units         Bohr    Units for the coordinates of the molecule
charge        0       The molecular charge
multiplicity  1       The spin multiplicity of the molecule

Here's an example for constructing a molecule object for water:

h2o=Molecule('h2o',
             atomlist = [(8,(0,0,0)),
                         (1,(1.0,0,0)),
                         (1,(0,1.0,0))],
             units = 'Angstrom')

(of course the bond-angle is 90 degrees here, and is thus completely wrong, but this is only an example). Here's an example for the hydroxide ion that shows the use of the charge field:

oh = Molecule('OH-',
              atomlist = [(8,(0,0,0)),(1,(0.96,0,0))],
              units = 'Angstrom',
              charge=-1)

Here's an example for the NO molecule that shows the use of the multiplicity field

no = Molecule('NO',
              atomlist = [(7,(0,0,0)),(8,(2.12955,0,0))],
              multiplicity=2)

As of version 1.5.1, you may construct molecules using the atomic symbol instead of the atomic number, e.g.

h2o=Molecule('h2o',
             atomlist = [('O',(0,0,0)),
                         ('H',(1.0,0,0)),
                         ('H',(0,1.0,0))],
             units = 'Angstrom')

Currently, the semiempirical code uses an extended verion of the Molecule object that adds a variety of additional features. Upcoming releases will hopefully unify the use of the Molecule between the HF, DFT, and semiempirical codes.

Construction of a Gaussian Basis Set

Basis functions are constructed using the CGBF (contracted Gaussian basis function) object, which, in turn, uses the PGBF (primitive Gaussian basis function) object.

Basis sets are simply lists of CGBF's. In the Ints module there is a convenience function getbasis that constructs basis sets for different molecules. Here is the docstring information for the getbasis function:

bfs = getbasis(atoms,basis_data=None)

Given a Molecule object and a basis library, form a basis set
constructed as a list of CGBF basis functions objects.

The basis data can be input from a number of data files in the PyQuante suite. Here are some of the more commonly used basis sets:

For, example, to construct a basis set for the h2o object for water created above, we would call

from basis_631ss import basis_data
bfs = getbasis(h2o,basis_data)

If the basis_data argument is omitted, the program will default to 6-31G**.

As of version 1.6, you may now specify a string for the basis_data argument, e.g. "6-31G**" so you may do things like

bfs = getbasis(h2o,'6-31G**')
bfs2 = getbasis(h2o,'sto-3g')

and so on. Use PyQuante.Basis.Tools.basis_map.keys() for a list of the supported basis strings.

Computation of One-Electron Integrals

The one-electron integrals consist of the overlap matrix S, the kinetic energy matrix t, and the nuclear attraction matrix vn. The latter two are often combined into the one-electron Hamiltonian h.

There are a number of helper functions in the Ints module:

These functions actually call instance functions of the CGBF objects, which can themselves be used individually. For example, the getS function is no more than the following code

def getS(bfs):
    "Form the overlap matrix"
    nbf = len(bfs)
    S = zeros((nbf,nbf),Float)

    for i in range(nbf):
        bfi = bfs[i]
        for j in range(nbf):
            bfj = bfs[j]
            S[i,j] = bfi.overlap(bfj)
    return S

Some of the instance functions in the CGBF module use functions in the pyints and cints modules.

Computation of Two-Electron Integrals

The two-electron integrals consist of the electron-electron Coulomb repulsion interactions. The easiest way to construct these is to use the functions in the Ints module

The Ints object (not to be confused with the Ints module, which is just a collection of helper functions) consists of a list of the two-electron integrals.

There are actually three different methods to computing the two-electron integrals, and the helper functions in the Ints module default to one of these functions.

There are python versions of these methods implemented in the modules pyints, rys, and hgp, respectively. For speed, there are also C-versions of these modules in cints, crys, and chgp. The program defaults to the Coulomb repulsion routines in crys, since these are the fastest (although recent improvements to chgp make it not much slower).

The coulomb_repulsion function is the same in all six routine:

def coulomb_repulsion((xa,ya,za),norma,(la,ma,na),alphaa,
                      (xb,yb,zb),normb,(lb,mb,nb),alphab,
                      (xc,yc,zc),normc,(lc,mc,nc),alphac,
                      (xd,yd,zd),normd,(ld,md,nd),alphad):

This routine computes the Coulomb repulsion integral between four basis functions. The terms xi, yi, zi are the origins of the different basis functions. The terms normi are the normalization constants for the basis function. The terms li, mi, ni are the exponents of the Cartesian powers for the basis function. And the terms alphai are the Gaussian exponents.

Hartree-Fock Calculations

Here is the docstring information for the rhf function that details the basic calling options:

rhf(atoms,**opts) - Closed-shell HF driving routine
atoms       A Molecule object containing the molecule

Options:      Value   Description
--------      -----   -----------
verbose       False   Print out extra information during DFT calc
ConvCriteria  1e-4    Convergence Criteria
MaxIter       20      Maximum SCF iterations
DoAveraging   True    Use DIIS averaging for convergence acceleration
ETemp         False   Use ETemp value for finite temperature DFT
                      If not False, set to temperature (float)
bfs           None    The basis functions to use. List of CGBF's
basis_data    None    The basis data to use to construct bfs
integrals     None    The one- and two-electron integrals to use
                      If not None, S,h,Ints
orbs          None    If not None, the guess orbitals

See the PyQuante Cookbook below for examples of this function's use.

Density Functional Theory Calculations

Here is the docstring information for the dft function that details the basic calling options:

dft(atoms,**opts) - DFT driving routine

atoms       A Molecule object containing the molecule

Options:      Value   Description
--------      -----   -----------
verbose       False   Print out extra information during DFT calc
ConvCriteria  1e-4    Convergence Criteria
MaxIter       20      Maximum SCF iterations
DoAveraging   True    Use DIIS averaging for convergence acceleration
ETemp         False   Use ETemp value for finite temperature DFT
                      If not False, set to temperature (float)
bfs           None    The basis functions to use. List of CGBF's
basis_data    None    The basis data to use to construct bfs
integrals     None    The one- and two-electron integrals to use
                      If not None, S,h,Ints
orbs          None    If not none, the guess orbitals
functional    SVWN    DFT functional (SVWN, BLYP, S0)
grid_nrad     32      Number of radial shells per atom
grid_fineness 1       Radial shell fineness. 0->coarse, 1->medium, 2->fine

See the PyQuante Cookbook below for examples of this function's use.

MINDO/3 (semiempirical) Calculations

To be completed.

PyQuante Cookbook

A Simple RHF Calculation on Hydrogen

Here's an example of running a simple restricted Hartree-Fock (RHF) calculation on the hydrogen molecule

from PyQuante.hartree_fock import rhf
from PyQuante.Molecule import Molecule
h2 = Molecule('h2',atomlist=[(1,(0,0,0)),(1,(1.4,0,0))])
en,orbe,orbs = rhf(h2)

Since no basis set is used, the program defaults to 6-31G**. At this geometry (R=1.4 bohr) this should produce an energy of -1.1313 hartrees.

A Verbose RHF Calculation on Hydrogen

Here's an example of running the same RHF calculation on H2, but outputting a bit more information by using the verbose tag:

en,orbe,orbs = rhf(h2,verbose=True)

Here is typical output from this run (note, the exact numbers may differ a bit because of floating point differences):

Nbf =  10
Nclosed =  1
Optimization of HF orbitals
0 -1.06918535085
1 -1.13004381838
2 -1.13125922261
3 -1.1312837961

Running a DFT Calculation on Hydrogen

Let's compare the results of this calculation to the results from DFT. We can run a DFT calculation on the same molecule by running the commands

from PyQuante.dft import dft
en,orbe,orbs = dft(h2)

This will produce an energy of -1.1353 hartrees. Again, the 6-31G** basis set is used by default. In DFT calculations, the functional defaults to SVWN (LDA). To use a different functional, you can type

en,orbe,orbs = dft(h2,functional='BLYP')

which will produce an energy of -1.1665 hartrees.

Running Multiple Calculations on the Same Molecule

Suppose we want to construct a molecule and run several different calculations (HF, LDA, BLYP) on the system. Normally the routines dft and rhf compute the integrals, but we can also pass them in as arguments. For example:

h2o = Molecule('H2O',
               atomlist = [(8,(0,0,0)),
                           (1,(0.959,0,0)),
                           (1,(-.230,0.930,0))],
               units = 'Angstrom')
bfs = getbasis(h2o)
S,h,Ints = getints(bfs,h2o)
enhf,orbehf,orbshf = rhf(h2o,integrals=(S,h,Ints))
print "HF  total energy = ",enhf
enlda,orbelda,orbslda = dft(h2o,
                           integrals=(S,h,Ints),
                           bfs = bfs)
print "LDA total energy = ",enlda
enblyp,orbeblyp,orbsblyp = dft(h2o,
                           integrals=(S,h,Ints),
                           bfs = bfs,
                           functional='BLYP')
print "BLYP total energy = ",enblyp