In [1]:

    

from pylab import*
from scipy.integrate import tplquad


    

In [2]:

    

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 [3]:

    

lmax = 0
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):
                        for iC in range(0,lmax+1):
                            for jC in range(0,lmax+1):
                                for kC in range(0,lmax+1):
                                    for iD in range(0,lmax+1):
                                        for jD in range(0,lmax+1):
                                            for kD in range(0,lmax+1):
                        
                                                if(iA+jA+kA <= lmax 
                                                   and iB+jB+kB <=lmax 
                                                   and iC+jC+kC <=lmax 
                                                   and iD+jD+kD <=lmax):
                                                    combo.append(array((iA, jA, kA, 
                                                                        iB, jB, kB,
                                                                        iC, jC, kC,
                                                                        iD, jD, kD)))

combo = array(combo)
print combo.shape
print combo


    

    




(1, 12)
[[0 0 0 0 0 0 0 0 0 0 0 0]]

In [4]:

    

def lowL(x, y=0):
    return -1.0
    
def upL(x, y=0):
    return 1.0

def integrand(x1,y1,z1,x2,y2,z2):
    return Ga(x1,y1,z1)*Gb(x1,y1,z1)*Gc(x2,y2,z2)*Gd(x2,y2,z2)/sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
    
def secondIntegrals(x, y, z):
    I = tplquad(integrand, lowL(0), upL(0), lowL, upL, lowL, upL, args=(x, y, z))
    return I[0]


    

In [5]:

    

A = array((1.2,2.3,3.4))
B = array((-1.3,1.4,-2.4))
C = array((2.3,0.9,3.2))
D = array((-2.0,1.9,2.2))
a = 0.2
b = 0.3
c = 0.4
d = 0.5
Ga = gaussian(A,a)
Gb = gaussian(B,b)
Gc = gaussian(C,c)
Gd = gaussian(D,d)

import time

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

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]
    Gc.i, Gc.j, Gc.k = combo[i,6:9]
    Gd.i, Gd.j, Gd.k = combo[i,9:12]
    
    
    I = tplquad(secondIntegrals, lowL(0), upL(0), lowL, upL, lowL, upL)
    
    
    results.append("CHECK_CLOSE(" + str("%.12e" %I[0]) + ", integrator.electronRepulsionIntegral(" \
          + str(Ga.i) + "," + str(Ga.j) + "," + str(Ga.k) + "," \
          + str(Gb.i) + "," + str(Gb.j) + "," + str(Gb.k) + "," \
          + str(Gc.i) + "," + str(Gc.j) + "," + str(Gc.k) + "," \
          + str(Gd.i) + "," + str(Gd.j) + "," + str(Gd.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




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

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


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-5-e0406ea98677> in <module>()
     25 
     26 
---> 27     I = tplquad(secondIntegrals, lowL(0), upL(0), lowL, upL, lowL, upL)
     28 
     29 

/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-4-23414140d954> in secondIntegrals(x, y, z)
      9 
     10 def secondIntegrals(x, y, z):
---> 11     I = tplquad(integrand, lowL(0), upL(0), lowL, upL, lowL, upL, args=(x, y, z))
     12     return I[0]

/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-4-23414140d954> in integrand(x1, y1, z1, x2, y2, z2)
      6 
      7 def integrand(x1,y1,z1,x2,y2,z2):
----> 8     return Ga(x1,y1,z1)*Gb(x1,y1,z1)*Gc(x2,y2,z2)*Gd(x2,y2,z2)/sqrt((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)
      9 
     10 def secondIntegrals(x, y, z):

<ipython-input-2-98dc4a54ccc4> in __call__(self, x, y, z)
     11 
     12         coff = (x-R[0])**i * (y-R[1])**j * (z-R[2])**k
---> 13         Ra = (x-R[0])**2 + (y-R[1])**2 + (z-R[2])**2
     14         expon = exp(-alpha*Ra)
     15 

KeyboardInterrupt: 





    




-c:8: RuntimeWarning: divide by zero encountered in double_scalars

In [ ]:

    

print dt/3600


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    

import numpy
import scipy.integrate
import math

def w(r, theta, phi, alpha, beta, gamma):
    return(-math.log(theta * beta))

def integrand(phi, alpha, gamma, r, theta, beta):
    ww = w(r, theta, phi, alpha, beta, gamma)
    k = 1.
    T = 1.
    return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta)

# limits of integration

def zero(x, y=0):
    return 0.

def one(x, y=0):
    return 1.

def pi(x, y=0):
    return math.pi

def twopi(x, y=0):
    return 2.*math.pi

# integrate over phi [0, Pi), alpha [0, 2 Pi), gamma [0, 2 Pi)
def secondIntegrals(r, theta, beta):
    res, err = scipy.integrate.tplquad(integrand, 0., 2.*math.pi, zero, twopi, zero, pi, args=(r, theta, beta))
    return res

# integrate over r [0, 1), beta [0, 2 Pi), theta [0, 2 Pi)
def integral():
    return scipy.integrate.tplquad(secondIntegrals, 0., 2.*math.pi, zero, twopi, zero, one)

expected = 16*math.pow(math.pi,5)/3.
result, err = integral()
diff = abs(result - expected)

print "Result = ", result, " estimated error = ", err
print "Known result = ", expected, " error = ", diff, " = ", 100.*diff/expected,


    

In [ ]:

    

import mcint
import random
import math

def w(r, theta, phi, alpha, beta, gamma):
    return(-math.log(theta * beta))

def integrand(x):
    r     = x[0]
    theta = x[1]
    alpha = x[2]
    beta  = x[3]
    gamma = x[4]
    phi   = x[5]

    k = 1.
    T = 1.
    ww = w(r, theta, phi, alpha, beta, gamma)
    return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta)

def sampler():
    while True:
        r     = random.uniform(0.,1.)
        theta = random.uniform(0.,2.*math.pi)
        alpha = random.uniform(0.,2.*math.pi)
        beta  = random.uniform(0.,2.*math.pi)
        gamma = random.uniform(0.,2.*math.pi)
        phi   = random.uniform(0.,math.pi)
        yield (r, theta, alpha, beta, gamma, phi)


domainsize = math.pow(2*math.pi,4)*math.pi*1
expected = 16*math.pow(math.pi,5)/3.

for nmc in [1000, 10000, 100000, 1000000, 10000000, 100000000]:
    random.seed(1)
    result, error = mcint.integrate(integrand, sampler(), measure=domainsize, n=nmc)
    diff = abs(result - expected)

    print "Using n = ", nmc
    print "Result = ", result, "estimated error = ", error
    print "Known result = ", expected, " error = ", diff, " = ", 100.*diff/expected, "%"
    print " "