Commit 0948b020 authored by Ivo Filot's avatar Ivo Filot
Browse files

Adding HF scripts

parent fa11d9db
......@@ -116,4 +116,7 @@ dmypy.json
# cache
__pycache__
*.pyc
\ No newline at end of file
*.pyc
# pdf files
*.pdf
from pyquante2 import basisset, molecule, rhf, onee_integrals
# Calculate CO
co = molecule([(6, 0.00000000, 0.00000000, -0.63546711),
(8, 0.00000000, 0.00000000, 0.47832425)],
units='Angstrom',
name="CO")
bfs = basisset(co, 'sto3g') # construct basis functions
solver = rhf(co, bfs) # set up solver object
ens = solver.converge() # perform scf calculation
orbe = solver.orbe # get MO energies
orbs = solver.orbs # get MO coefficients
overlap = onee_integrals(bfs, co).S # get overlap matrix
print('Energies\n', ens)
print('Orbital energies\n', orbe)
print('Orbital coefficients\n', orbs.transpose())
print('Overlap matrix\n', overlap)
\ No newline at end of file
from pyquante2 import basisset, molecule, rhf, onee_integrals, contourplot
# conda install -c rpmuller pyquante2_pure
# Calculate H2
h2 = molecule([(1, 0.00000000, 0.00000000, 0.36628549),
(1, 0.00000000, 0.00000000, -0.36628549)],
units='Angstrom')
bfs = basisset(h2, 'sto3g') # construct basis functions
solver = rhf(h2, bfs) # set up solver object
ens = solver.converge() # perform scf calculation
orbe = solver.orbe # get MO energies
orbs = solver.orbs # get MO coefficients
overlap = onee_integrals(bfs, h2).S # get overlap matrix
print('Energies\n', ens)
print('Orbital energies\n', orbe)
print('Orbital coefficients\n', orbs.transpose())
print('Overlap matrix\n', overlap)
\ No newline at end of file
from pyquante2 import basisset, molecule, rhf, onee_integrals, contourplot
import orbplot.orbplotter as op
# Calculate H2
h2 = molecule([(1, 0.00000000, 0.00000000, 0.36628549),
(1, 0.00000000, 0.00000000, -0.36628549)],
units='Angstrom')
bfs = basisset(h2, 'sto3g') # construct basis functions
solver = rhf(h2, bfs) # set up solver object
ens = solver.converge() # perform scf calculation
orbe = solver.orbe # get MO energies
orbs = solver.orbs # get MO coefficients
overlap = onee_integrals(bfs, h2).S # get overlap matrix
print('Orbital energies\n', orbe)
print('Orbital coefficients\n', orbs.transpose())
op.plot_contour(h2, orbs[:,0], bfs, plane='xz', depth=0,
title='H$_{2}$ $\\sigma_{g}$ orbital')
op.plot_contour(h2, orbs[:,1], bfs, plane='xz', depth=0,
title='H$_{2}$ $\\sigma_{u}$ orbital')
\ No newline at end of file
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
def plot_energies(ens):
plt.figure(dpi=240)
plt.plot(ens, '.-', color='black', markersize=18)
plt.xlabel('Iterations')
plt.ylabel('Energy [Hartree]')
plt.show()
def plot_contour(atoms,orb,bfs,plane='xy',depth=0,npts=200,
title="Contour plot"):
if plane is 'xy':
plot_contour_xy(atoms,orb,bfs,depth,npts,title)
return
if plane is 'yz':
plot_contour_yz(atoms,orb,bfs,depth,npts,title)
return
if plane is 'xz':
plot_contour_xz(atoms,orb,bfs,depth,npts,title)
return
def plot_contour_xy(atoms,orb,bfs,depth=0,npts=200,title='Contour plot xy'):
bbox = atoms.bbox()
x = np.linspace(-5, 5, npts)
y = np.linspace(-5, 5, npts)
f = np.zeros((npts,npts),'d')
X,Y = np.meshgrid(x,y)
for c,bf in zip(orb,bfs):
f += c*bf(X,Y,level)
plt.figure(figsize=(4,4), dpi=240)
ax = plt.gca()
rrange = np.max([np.abs(np.min(f)), np.max(f)])
levels = np.linspace(-rrange, rrange, 32)
contours = plt.contour(X, Y, f, levels=levels, colors='black', zorder=2)
#plt.clabel(contours, inline=True, fontsize=8)
xlim = ax.get_xlim()
ylim = ax.get_ylim()
clr = plt.imshow(f, extent=[xlim[0], xlim[1], ylim[0], ylim[1]],
origin='lower', interpolation='quadric')
#ax.set_aspect(1)
plt.title(title)
plt.xlabel('$x$ in $\\AA$')
plt.ylabel('$y$ in $\\AA$')
plt.xticks(np.linspace(-5,5,5))
plt.yticks(np.linspace(-5,5,5))
for at in atoms:
plt.text(at.r[0], at.r[1], at.tag()[0], color='white')
plt.show()
def plot_contour_yz(atoms,orb,bfs,depth=0,npts=200,title='Contour plot yz'):
bbox = atoms.bbox()
y = np.linspace(-5, 5, npts)
z = np.linspace(-5, 5, npts)
f = np.zeros((npts,npts),'d')
Y,Z = np.meshgrid(y,z)
for c,bf in zip(orb,bfs):
f += c*bf(level,Y,Z)
plt.figure(figsize=(4,4), dpi=240)
ax = plt.gca()
rrange = np.max([np.abs(np.min(f)), np.max(f)])
levels = np.linspace(-rrange, rrange, 32)
contours = plt.contour(Y, Z, f, levels=levels, colors='black', zorder=2)
#plt.clabel(contours, inline=True, fontsize=8)
xlim = ax.get_xlim()
ylim = ax.get_ylim()
clr = plt.imshow(f, extent=[xlim[0], xlim[1], ylim[0], ylim[1]],
origin='lower', interpolation='quadric')
#ax.set_aspect(1)
plt.title(title)
plt.xlabel('$y$ in $\\AA$')
plt.ylabel('$z$ in $\\AA$')
plt.xticks(np.linspace(-5,5,5))
plt.yticks(np.linspace(-5,5,5))
for at in atoms:
plt.text(at.r[1], at.r[2], at.tag()[0], color='white')
plt.show()
def plot_contour_xz(atoms,orb,bfs,depth=0,npts=200,title='Contour plot xz'):
bbox = atoms.bbox()
x = np.linspace(-5, 5, npts)
z = np.linspace(-5, 5, npts)
f = np.zeros((npts,npts),'d')
X,Z = np.meshgrid(x,z)
for c,bf in zip(orb,bfs):
f += c*bf(X,depth,Z)
plt.figure(figsize=(4,4), dpi=240)
ax = plt.gca()
rrange = np.max([np.abs(np.min(f)), np.max(f)])
levels = np.linspace(-rrange, rrange, 32)
contours = plt.contour(X, Z, f, levels=levels, colors='black', zorder=2)
#plt.clabel(contours, inline=True, fontsize=8)
xlim = ax.get_xlim()
ylim = ax.get_ylim()
clr = plt.imshow(f, extent=[xlim[0], xlim[1], ylim[0], ylim[1]],
origin='lower', interpolation='quadric')
#ax.set_aspect(1)
plt.title(title)
plt.xlabel('$x$ in $\\AA$')
plt.ylabel('$z$ in $\\AA$')
plt.xticks(np.linspace(-5,5,5))
plt.yticks(np.linspace(-5,5,5))
for at in atoms:
plt.text(at.r[0], at.r[2], at.tag()[0], color='white')
plt.show()
\ No newline at end of file
import orbplot.orbplotter as op
#
# The command below assumes that you have already performed the SCF calculation
# for the CO molecule. If you have not done so (or cannot remember doing so),
# please open 'co.py' (it resides in the same folder as this file) and run
# that program.
#
for i in range(0,10):
op.plot_contour(co, orbs[:,i], bfs, plane='xz', depth=0,
title='CO MO %i' % i)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment