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:
-
basis_631ss The Pople 6-31G** basis set
-
basis_sto3g The Pople STO-3G basis set
-
basis_321 The Pople 3-21G basis set
-
basis_ccpvtz The Dunning cc-pVTZ basis set
-
basis_ccpvtzmf The Dunning cc-pVTZ(-f) basis set (cc-pVTZ without
f-functions)
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**.
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:
-
getT Form the kinetic energy matrix t
-
getS Form the overlap matrix S
-
getV Form the nuclear attraction matrix Vn
-
get1ints Form and return S,h
-
getints Form and return S,h,Ints, where Ints are the
two-electron integrals (see below).
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
-
get2ints Form and return Ints, where Ints are the
two-electron integrals (see below).
-
getints Form and return S,h,Ints, where Ints are the
two-electron integrals.
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
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