In [1]:

    

from Molecule import Molecule
from SCF import SCF
import scipy as sp
from renderMolecule import renderMolecule
import matplotlib.pyplot as plt
%matplotlib inline


    

In [2]:

    

water = Molecule("Molecules/wat.mol")


    

In [3]:

    

theta = sp.linspace(20,160,400)


    

In [4]:

    

def coordinates(theta):
    return sp.array([[0,0,0],[1.88972612*sp.cos(theta*sp.pi/180.),1.88972612*sp.sin(theta*sp.pi/180.),0],
                             [1.88972612*sp.cos(theta*sp.pi/180.), -1.88972612*sp.sin(theta*sp.pi/180.),0]])


    

In [5]:

    

coordinates(90)


    

    
Out[5]:





array([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  1.15712352e-16,   1.88972612e+00,   0.00000000e+00],
       [  1.15712352e-16,  -1.88972612e+00,   0.00000000e+00]])

In [6]:

    

v_coordinates = sp.vectorize(coordinates,otypes=[sp.ndarray])
coordArray = v_coordinates(theta)


    

In [7]:

    

def EnergyVal(coord):
    res = SCF(water,MP2=True,basis="6-31++G",cartMatrix=coord)
    return res.TotalEnergy


    

In [8]:

    

v_EnergyVal = sp.vectorize(EnergyVal,otypes=[sp.float64])
Energies = v_EnergyVal(coordArray)


    

In [9]:

    

plt.xkcd()
plt.ylabel("RHF Energy (Hartrees)")
plt.xlabel("H-O-H Bond Angle ($^\circ$)")
plt.plot(2 * theta,Energies)


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0xb1b0d10>]






    

In [10]:

    

print "Optimal H-O-H angle is equal to: " +str(2*theta[sp.argmax(Energies[theta.size/2-150:theta.size/2+150])])


    

    




Optimal H-O-H angle is equal to: 145.263157895

In [ ]: