In [1]:

    

from Molecule import Molecule
from SCF import SCF
import numpy as np
from mayavi import mlab
import timeit
%load_ext cythonmagic


    

    




---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-83f6146086c1> in <module>()
      2 from SCF import SCF
      3 import numpy as np
----> 4 from mayavi import mlab
      5 import timeit
      6 get_ipython().magic(u'load_ext cythonmagic')

ImportError: No module named mayavi

In [2]:

    

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


    

In [4]:

    

res = SCF(water,basis='6-31G')


    

In [5]:

    

X,Y,Z = np.mgrid[-10:10:100j,-10:10:100j,-10:10:100j]


    

In [6]:

    

%%cython --compile-args=-O3,-ffast-math
cimport cython
from libc.math cimport pow , exp

@cython.nonecheck(False)
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef double cProb(double x, double y, double z, object orbList,double [:] C, int nOrb):
    cdef int i, j, xnum, ynum, znum, nPrim
    cdef object orb
    cdef double phi = 0
    cdef double xN, yN, zN, R2, Const, ang
    cdef double [:] cent, a, d
    for i in xrange(nOrb):
        orb = orbList[i]
        xnum, ynum, znum = orb.qnums
        cent = orb.Center
        xN, yN, zN = x - cent[0], y - cent[1], z - cent[2]
        ang = pow(xN, xnum) * pow(yN, ynum) * pow(zN, znum)
        Const = C[i] * ang
        R_2 = xN * xN + yN * yN + zN * zN
        nPrim = orb.nPrim
        d, a = orb.d, orb.a
        for j in xrange(nPrim):
            phi += Const * d[j] * exp(-a[j] * R_2)
    return phi


    

In [7]:

    

def Prob(x,y,z):
    return cProb(x,y,z,res.orbList,res.AOcoefs[:,12],res.nOrb)


    

In [9]:

    

v_Prob = np.vectorize(Prob)


    

In [10]:

    

P = v_Prob(X,Y,Z)


    

In [ ]:

    

mlab.contour3d(X,Y,Z,P,contours=16,transparent=True)


    

    
Out[ ]:





<mayavi.modules.iso_surface.IsoSurface at 0x8d5be30>

In [ ]:

    

mlab.show()


    

In [ ]: