In [7]:

    

from pylab import*
from scipy.integrate import tplquad
from numpy import *
from sympy import *

init_printing(use_latex=True)


    

In [2]:

    

Ax, Ay, Az = symbols('A_x, A_y, A_z')
Bx, By, Bz = symbols('B_x, B_y, B_z')

R = sqrt((Ax-Bx)**2+ (Ay-By)**2 + (Az-Bz)**2)

V = 1.0/R
diff(R,Ax) *-1/R**2


    

    
Out[2]:





$$- \frac{A_{x} - B_{x}}{\left(\left(A_{x} - B_{x}\right)^{2} + \left(A_{y} - B_{y}\right)^{2} + \left(A_{z} - B_{z}\right)^{2}\right)^{\frac{3}{2}}}$$

In [8]:

    

class gaussian():
    
    def __init__(self,R,alpha):
        self.R = R
        self.alpha = alpha
        self.i, self.j, self.k = 0 ,0,0
        
    def __call__(self,x,y,z):
        i, j, k = self.i, self.j, self.k
        R = self.R; alpha = self.alpha
        
        coff = (x-R[0])**i * (y-R[1])**j * (z-R[2])**k
        Ra = (x-R[0])**2 + (y-R[1])**2 + (z-R[2])**2
        expon = exp(-alpha*Ra)
        
        return coff*expon


    

In [9]:

    

lmax = 2
combo=[]

for iA in range(0,lmax+1):
    for jA in range(0,lmax+1):
        for kA in range(0,lmax+1):
            for iB in range(0,lmax+1):
                for jB in range(0,lmax+1):
                    for kB in range(0,lmax+1):
                        
                        if(iA+jA+kA <= lmax and iB+jB+kB <=lmax):
                            combo.append(array((iA, jA, kA, iB, jB, kB)))

combo = array(combo)
print combo.shape


    

    




(100, 6)

In [16]:

    

upL = -10
lowL =  10


    

In [17]:

    

A = array((1.2,2.3,3.4))
B = array((-1.3,1.4,-2.4))
a = 0.2
b = 0.3
Ga = gaussian(A,a)
Gb = gaussian(B,b)
C = array((2.3,0.9,3.2))

import time

t0 = time.time()
results = []
error = []

#combo = zeros((1,6))
#combo[0,:] = array([1,0,0,0,0,0])


for i in range(0,combo.shape[0]):
    Ga.i, Ga.j, Ga.k = combo[i,:3]
    Gb.i, Gb.j, Gb.k = combo[i,3:6]
    
    I = tplquad(lambda x,y,z: 
                Ga(x,y,z)*Gb(x,y,z)/sqrt((x-C[0])**2 + (y-C[1])**2 + (z-C[2])**2),
                upL, lowL, lambda y: upL, lambda y: lowL, lambda x,y: upL, lambda x,y: lowL)
    
    
    print "primitiveA.setPowers({"  \
          + str(Ga.i) + "," + str(Ga.j) + "," + str(Ga.k) +"});primitiveB.setPowers({" \
          + str(Gb.i) + "," + str(Gb.j) + "," + str(Gb.k) + "});"
    print "integrator.setPrimitiveA(primitiveA);integrator.setPrimitiveB(primitiveB);"
    
    results.append("CHECK_CLOSE(" + str("%.12e" %I[0]) \
                   + ", integrator.nuclearAttractionIntegral(), 1e-5);")
    
    error.append(str("%.12e" %I[0]) + ",  " + str("%.12e" %I[1]) + "   " + str(combo[i,:]))

    print results[-1] + "\n"

dt = time.time()-t0
#print dt/60.0




f = open("coulomb", "w")
f.write("\n".join(map(lambda x: str(x), results)))
f.close()

e = open("coulombError", "w")
e.write("\n".join(map(lambda x: str(x), error)))
e.close()


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-17-c0e7f4a49bc9> in <module>()
     23     I = tplquad(lambda x,y,z: 
     24                 Ga(x,y,z)*Gb(x,y,z)/sqrt((x-C[0])**2 + (y-C[1])**2 + (z-C[2])**2),
---> 25                 upL, lowL, lambda y: upL, lambda y: lowL, lambda x,y: upL, lambda x,y: lowL)
     26 
     27 

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in tplquad(func, a, b, gfun, hfun, qfun, rfun, args, epsabs, epsrel)
    485 
    486     """
--> 487     return dblquad(_infunc2,a,b,gfun,hfun,(func,qfun,rfun,args),epsabs=epsabs,epsrel=epsrel)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in dblquad(func, a, b, gfun, hfun, args, epsabs, epsrel)
    423 
    424     """
--> 425     return quad(_infunc,a,b,(func,gfun,hfun,args),epsabs=epsabs,epsrel=epsrel)
    426 
    427 def _infunc2(y,x,func,qfun,rfun,more_args):

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    245     if type(args) != type(()): args = (args,)
    246     if (weight is None):
--> 247         retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
    248     else:
    249         retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    310     if points is None:
    311         if infbounds == 0:
--> 312             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    313         else:
    314             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _infunc(x, func, gfun, hfun, more_args)
    373     b = hfun(x)
    374     myargs = (x,) + more_args
--> 375     return quad(func,a,b,args=myargs)[0]
    376 
    377 def dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-8, epsrel=1.49e-8):

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    245     if type(args) != type(()): args = (args,)
    246     if (weight is None):
--> 247         retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
    248     else:
    249         retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    310     if points is None:
    311         if infbounds == 0:
--> 312             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    313         else:
    314             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _infunc2(y, x, func, qfun, rfun, more_args)
    429     b2 = rfun(x,y)
    430     myargs = (y,x) + more_args
--> 431     return quad(func,a2,b2,args=myargs)[0]
    432 
    433 def tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-8,

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    245     if type(args) != type(()): args = (args,)
    246     if (weight is None):
--> 247         retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
    248     else:
    249         retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    310     if points is None:
    311         if infbounds == 0:
--> 312             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    313         else:
    314             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

<ipython-input-17-c0e7f4a49bc9> in <lambda>(x, y, z)
     22 
     23     I = tplquad(lambda x,y,z: 
---> 24                 Ga(x,y,z)*Gb(x,y,z)/sqrt((x-C[0])**2 + (y-C[1])**2 + (z-C[2])**2),
     25                 upL, lowL, lambda y: upL, lambda y: lowL, lambda x,y: upL, lambda x,y: lowL)
     26 

/usr/lib/python2.7/dist-packages/sympy/functions/elementary/miscellaneous.pyc in sqrt(arg)
    103     """
    104     # arg = sympify(arg) is handled by Pow
--> 105     return C.Pow(arg, S.Half)
    106 
    107 

/usr/lib/python2.7/dist-packages/sympy/core/cache.pyc in wrapper(*args, **kw_args)
     83         k = tuple(k)
     84         try:
---> 85             return func_cache_it_cache[k]
     86         except KeyError:
     87             pass

/usr/lib/python2.7/dist-packages/sympy/core/numbers.pyc in __hash__(self)
   1229 
   1230     def __hash__(self):
-> 1231         return super(Rational, self).__hash__()
   1232 
   1233     def factors(self, limit=None, use_trial=True,

KeyboardInterrupt: 

In [ ]:

    

print dt/3600


    

In [ ]:

    

Ax, Ay, Az = symbols('A_x, A_y, A_z')
x, y, z = symbols('x, y, z')
A = array([Ax,Ay,Az])



Axn = 1.2; Ayn = 2.3; Azn = 3.4
B = array((-1.3,1.4,-2.4))
a = 0.2
b = 0.3
Ga = gaussian(A,a)
Gb = gaussian(B,b)
C = array((2.3,0.9,3.2))

import time

t0 = time.time()
results = []
error = []

#combo = zeros((1,6))
#combo[0,:] = array([1,0,0,0,0,0])


for i in range(0,combo.shape[0]):
    Ga.i, Ga.j, Ga.k = combo[i,:3]
    Gb.i, Gb.j, Gb.k = combo[i,3:6]
    
    Ga = Ga(x,y,z)
    Gb = Gb(x,y,z)
    S = Ga*Gb

    dSAx = diff(S,Ax)
    dSAy = diff(S,Ay)
    dSAz = diff(S,Az)

    dSAx =dSAx.subs(Ax,  Axn).subs(Ay,  Ayn).subs(Az,  Azn)
    dSAy =dSAy.subs(Ay,  Ayn).subs(Ax,  Axn).subs(Az,  Azn)
    dSAz =dSAz.subs(Az,  Azn).subs(Ay,  Ayn).subs(Ax,  Axn)
    

    
    I = tplquad(lambda x,y,z: 
                Ga(x,y,z)*Gb(x,y,z)/sqrt((x-C[0])**2 + (y-C[1])**2 + (z-C[2])**2),
                upL, lowL, lambda y: upL, lambda y: lowL, lambda x,y: upL, lambda x,y: lowL)
    
    
    results.append("CHECK_CLOSE(" + str("%.12e" %I[0]) + ", integrator.nuclearAttractionIntegral(" \
          + str(Ga.i) + "," + str(Ga.j) + "," + str(Ga.k) + "," \
          + str(Gb.i) + "," + str(Gb.j) + "," + str(Gb.k) + ")" \
            + ", 1e-5);")
    
    error.append(str("%.12e" %I[0]) + ",  " + str("%.12e" %I[1]) + "   " + str(combo[i,:]))

    print results[-1]

dt = time.time()-t0
#print dt/60.0


    

In [ ]: