3. Basis Set
MoHa supports Cartesian Gaussian type orbtial (GTO), it is an atomic orbital used in linear combinations forming molecular orbitals. Cartesian GTOs are defiend by an angular par which is homogeneous polynomial in the compoents x, y, and z of the position vector \(\mathbf{r}\). That is, \begin{equation} G_{ijk}(r_A,\alpha) = x^i_A y^j_A z^k_A e^{-\alpha r^2_A} \end{equation}
here \(r_A = r - A\) is the atom-centered function \(\alpha > 0\) is the real orbital exponent and \(i+j+k = n\).
MoHa use the same basis set input file as NWChem, all the available basis set file can be found in the moha/data/examples directory. If you need basis set more than MoHa offered, you can downloaded it from the EMSL webpage (https://bse.pnl.gov/bse/portal).
STO-3G EMSL basis set of hydrogen and oxygen:
# STO-3G EMSL Basis Set Exchange Library 8/21/18 3:24 AM
# Elements References
# -------- ----------
# H - Ne: W.J. Hehre, R.F. Stewart and J.A. Pople, J. Chem. Phys. 2657 (1969).
# Na - Ar: W.J. Hehre, R. Ditchfield, R.F. Stewart, J.A. Pople,
# J. Chem. Phys. 2769 (1970).
# K,Ca - : W.J. Pietro, B.A. Levy, W.J. Hehre and R.F. Stewart,
# Ga - Kr: J. Am. Chem. Soc. 19, 2225 (1980).
# Sc - Zn: W.J. Pietro and W.J. Hehre, J. Comp. Chem. 4, 241 (1983) + Gaussian.
# Y - Cd: W.J. Pietro and W.J. Hehre, J. Comp. Chem. 4, 241 (1983). + Gaussian
#
BASIS "ao basis" PRINT
#BASIS SET: (3s) -> [1s]
H S
3.42525091 0.15432897
0.62391373 0.53532814
0.16885540 0.44463454
#BASIS SET: (6s,3p) -> [2s,1p]
O S
130.7093200 0.15432897
23.8088610 0.53532814
6.4436083 0.44463454
O SP
5.0331513 -0.09996723 0.15591627
1.1695961 0.39951283 0.60768372
0.3803890 0.70011547 0.39195739
END
Atom |
number |
name |
symbol |
coordinate |
mass |
|---|---|---|---|---|---|
mol[0] |
8 |
Oxygen |
O |
(0.000000,-0.143225,0.000000) |
15.9994 |
mol[1] |
1 |
Hydrogen |
H |
(1.638036,1.136548,-0.000000) |
1.007975 |
mol[2] |
1 |
Hydrogen |
H |
(-1.638036,1.136548,-0.000000) |
1.007975 |
[1]:
from moha.molecule import Molecule
from moha.basis import BasisSet
geo = [[8, 0.000000000000, -0.143225816552, 0.000000000000],
['h', 1.638036840407, 1.136548822547, -0.000000000000],
["H", -1.638036840407, 1.136548822547, -0.000000000000]]
mol = Molecule.build(geo,pg=False)
orb = BasisSet.build(mol,'sto-3g.nwchem')
[2]:
print(dir(orb[0]))
print(-orb[0])
print(orb[0]+orb[1])
print(2*orb[0])
print(orb[0]/2)
['__add__', '__call__', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__mul__', '__ne__', '__neg__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__rmul__', '__setattr__', '__sizeof__', '__str__', '__sub__', '__subclasshook__', '__truediv__', '__weakref__', 'coefficients', 'expansion', 'exponents', 'norm', 'norm_memo', 'origin', 'parallel', 'shell']
<moha.basis.gauss.ContractedGaussian object at 0x7fc492f0d2a0>
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Input In [2], in <cell line: 3>()
1 print(dir(orb[0]))
2 print(-orb[0])
----> 3 print(orb[0]+orb[1])
4 print(2*orb[0])
5 print(orb[0]/2)
TypeError: unsupported operand type(s) for +: 'ContractedGaussian' and 'ContractedGaussian'
Atom |
shell |
Angular momentum |
Spherical part |
|---|---|---|---|
O |
8 |
Oxygen |
O |
3.1. Primitive Gaussian Type Orbital
Most quantum chemistry package supports Contracted Gaussian type orbtial(CGTO), it is an atomic orbital forming in linear combinations of primitive Gaussians type orbital(PGTO). Here in this note, only the integral over primitive Gaussian type orbital is to be considered which are defiend by an angular part which is homogeneous polynomial in the compoents x, y, and z of the position vector \(\mathbf{r}\). That is,
\begin{equation} G_{ijk}(r_A,a) = x^i_A y^j_A z^k_A e^{-a r^2_A} \end{equation}
a>0 is the orbital exponent.
\(r_A = r − A\) is the electronic coordinate.
i ≥ 0, j ≥ 0, k ≥ 0 is the quantum number.
Total angular-momentum quantum number l = i + j + k ≥ 0
The Cartesian Gaussians may be factorized in the three Cartesian directons \begin{equation} G_{ijk}(r_A,a) = G_i(r_A,a)G_j(r_A,a)G_k(r_A,a) \end{equation} where \begin{equation} G_i(r_A,a)= x^i_A e^{-a x^2_A} \end{equation}
Gaussians satisfy a simple recurrence relation:
\begin{equation} x_A G_i(a, x_A ) = G_{i+1}(a, x_A ) \end{equation}
The differentiation of a Gaussian yields a linear combination of two Gaussians \begin{equation} \frac{\partial G_i(a, x_A )}{\partial A_x} = -\frac{\partial G_i(a, x_A)}{\partial x} = 2aG_{i+1}(a, x_A) − iG_{i−1}(a, x_A) \end{equation}
An important property of Gaussians is the Gaussian product rule. The product of two Gaussians is another Gaussian centered somewhere on the line connecting the original Gaussians; its exponent is the sum of the original exponents:
\begin{equation} e^{-a x^2_A} e^{-b x^2_B} = e^{-\mu X_{AB}^2} e^{-p x^2_P} \end{equation}
while the centre-of-charge \(P_x\) and the relative coordinate or Gaussian separation \(X_{AB}\) are given by
\begin{equation} \begin{aligned} p &= a+b\\ \mu &= \frac{ab}{a+b}\\ P_x &= \frac{aA_x+bB_x}{a+b}\\ X_{AB} &= A_x - B_x \end{aligned} \end{equation}
[ ]:
import numpy as np
class PrimitiveGaussian(object):
"""Primitive Gaussian functions.
Attributes
----------
origin : 3-list
Coordinate of the nuclei.
n : int
Principal quantum number.
shell : 3-tuple
Angular momentum.
coef : float
Coefficent of Primitive Gaussian function.
exp : float
Primitive Gaussian exponent.
Methods
-------
__init__(self,type,origin,n,shell,coef,exp)
Initialize the instance.
"""
def __init__(self,coefficient,origin,shell,exponent):
"""Initialize the instance.
Parameters
----------
origin : list
Coordinate of the nuclei.
n : int
Principal quantum number.
shell : list
Angular momentum.
coef : float
Primitive Gaussian coefficients.
exp : float
Primitive Gaussian exponent.
"""
self.coefficient = coefficient
self.origin = origin
self.shell = shell
self.exponent = exponent
if __name__ == '__main__':
# Molecular coordinates
H2O = [[0., 1.43233673, -0.96104039],
[0., -1.43233673, -0.96104039],
[0., 0., 0.24026010]]
# Orbital exponents
OrbCoeff = np.array([[3.425250914, 0.6239137298, 0.168855404],
[3.425250914, 0.6239137298, 0.168855404],
[130.7093214, 23.80886605, 6.443608313],
[5.033151319, 1.169596125, 0.38038896],
[5.033151319, 1.169596125, 0.38038896],
[5.033151319, 1.169596125, 0.38038896],
[5.033151319, 1.169596125, 0.38038896]])
# H1s, H2s, O1s, O2s, O2px , O2py, O2p
FCenter = [H2O[0], H2O[1], H2O[2], H2O[2], H2O[2], H2O[2], H2O[2]]
CartAng = [[0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0],
[1, 0, 0], [0, 1, 0], [0, 0, 1]]
pga = PrimitiveGaussian(1.0,FCenter[0],CartAng[0],OrbCoeff[0,0])
pgb = PrimitiveGaussian(1.0,FCenter[6],CartAng[6],OrbCoeff[6,0])