from pyquante2 import rhf,h2,h2o,c6h6,basisset

h2o

c6h6

print c6h6

h2

bfs = basisset(h2)
bfs

solver = rhf(h2,bfs)
solver

solver.converge()
solver

print solver

from pyquante2 import lineplot_orbs,contourplot,line

points = line((0,0,-5),(0,0,5))
lineplot_orbs(points,solver.orbs[:,:2],bfs,
    title="Plots of bonding and antibonding orbitals of H2")

contourplot('yz',h2,solver.orbs[:,1],bfs,
    title="Contours of H2 antibonding orbital")

solver = rhf(h2,basisset(h2))

solver

solver.converge(maxiters=1)

solver.energies

import PyQuante, pyquante2

print map(PyQuante.pyints.fact2,[0,1,3,8,-1]), map(pyquante2.utils.fact2,[0,1,3,8,-1])
for a,b in [(5,2),(10,5)]:
    print PyQuante.pyints.binomial(a,b),pyquante2.utils.binomial(a,b)
print PyQuante.pyints.Fgamma(0,0),pyquante2.utils.Fgamma(0,0)

s1 = PyQuante.PGBF.PGBF(1,(0.,0.,0.))
s2 = pyquante2.pgbf(1)
print s1.overlap(s1),pyquante2.S(s2,s2)
print s1.kinetic(s1),pyquante2.T(s2,s2)
print s1.nuclear(s1,(0,0,0)),pyquante2.V(s2,s2,(0,0,0))

from pyquante2.utils import *
from math import gamma
from scipy.special import gammainc
for x in [0.1,2,3,4]:
    print gamma(x),gammainc(x,1e-10),gammainc(x,1e10),gamm_inc(x,1e-10),gamm_inc(x,1e10)

from PyQuante.CGBF import CGBF
from pyquante2 import cgbf,S,T
c1 = CGBF((0,0,0),(0,0,0))
c1b = CGBF((0,0,1.0),(0,0,0))
exps,coefs = [],[]
for ex,co in [(3.4252509099999999, 0.15432897000000001),
           (0.62391373000000006, 0.53532813999999995),
           (0.16885539999999999, 0.44463454000000002)]:
    exps.append(ex)
    coefs.append(co)
    c1.add_primitive(ex,co)
    c1b.add_primitive(ex,co)
c1.normalize()
c1b.normalize()

c2 = cgbf(exps=exps,coefs=coefs)
c2b = cgbf((0,0,1.0),exps=exps,coefs=coefs)

print "Overlaps"
print c1.overlap(c1), S(c2,c2)
print c1.overlap(c1b), S(c2,c2b)

print "Kinetics"
print c1.kinetic(c1), T(c2,c2)
print c1.kinetic(c1b), T(c2,c2b)

s = pyquante2.pgbf(1)
from pyquante2.ints import two,hgp
%timeit two.ERI(s,s,s,s)

%timeit hgp.ERI_hgp(s,s,s,s)

from pyquante2 import ctwo
%timeit ctwo.ERI(s,s,s,s)

%timeit ctwo.ERI_hgp(s,s,s,s)

from pyquante2.ctwo import ERI,ERI_hgp,vrr,vrr_recursive,vrr_nonrecursive
from pyquante2.ints.hgp import vrr as pyvrr

zero = array([0,0,0],'d')
%timeit pyvrr(zero,1.,(0,0,0),1.,zero,1.,1., zero,1.,(0,0,0),1.,zero,1.,1.,0)

%timeit vrr(0,0,0,1.,0,0,0,1.,0,0,0,1.,1.,0,0,0,1.,0,0,0,1.,0,0,0,1.,1.,0)

%timeit vrr_recursive(0,0,0,1.,0,0,0,1.,0,0,0,1.,1.,0,0,0,1.,0,0,0,1.,0,0,0,1.,1.,0)

%timeit vrr_nonrecursive(0,0,0,1.,0,0,0,1.,0,0,0,1.,1.,0,0,0,1.,0,0,0,1.,0,0,0,1.,1.,0)

from PyQuante.PGBF import coulomb
from pyquante2.ctwo import ERI_hgp, ERI

s1 = PyQuante.PGBF.PGBF(1,(0.,0.,0.))
s2 = pyquante2.pgbf(1)

s1b = PyQuante.PGBF.PGBF(1,(0.,0.,1.))
s2b = pyquante2.pgbf(1,(0,0,1))

%timeit coulomb(s1,s1,s1,s1)
%timeit ERI_hgp(s2,s2,s2,s2)
print coulomb(s1,s1,s1,s1), ERI_hgp(s2,s2,s2,s2)
print coulomb(s1,s1,s1b,s1b), ERI_hgp(s2,s2,s2b,s2b), ERI(s2,s2,s2b,s2b)
print coulomb(s1,s1b,s1,s1b), ERI_hgp(s2,s2b,s2,s2b), ERI(s2,s2b,s2,s2b)

from PyQuante.CGBF import CGBF,coulomb
from pyquante2 import cgbf
from pyquante2.ctwo import ERI_hgp

c1 = CGBF((0,0,0),(0,0,0))
c1b = CGBF((0,0,1.0),(0,0,0))
exps,coefs = [],[]
for ex,co in [(3.4252509099999999, 0.15432897000000001),
           (0.62391373000000006, 0.53532813999999995),
           (0.16885539999999999, 0.44463454000000002)]:
    exps.append(ex)
    coefs.append(co)
    c1.add_primitive(ex,co)
    c1b.add_primitive(ex,co)
c1.normalize()
c1b.normalize()

c2 = cgbf(exps=exps,coefs=coefs)
c2b = cgbf((0,0,1.0),exps=exps,coefs=coefs)

%timeit coulomb(c1,c1,c1,c1)
%timeit ERI_hgp(c2,c2,c2,c2)
%timeit ERI(c2,c2,c2,c2)
print coulomb(c1,c1,c1,c1), ERI_hgp(c2,c2,c2,c2)
print coulomb(c1,c1,c1b,c1b), ERI_hgp(c2,c2,c2b,c2b)
print coulomb(c1,c1b,c1,c1b), ERI_hgp(c2,c2b,c2,c2b)

from PyQuante.Ints import getbasis,getints,getT,getV
from PyQuante.hartree_fock import rhf
from PyQuante.Molecule import Molecule

LiH = Molecule('lih',
                 [(3,( .0000000000, .0000000000, .0000000000)),
                  (1,( .0000000000, .0000000000,1.629912))],
                 units='Angstroms')
bfs = getbasis(LiH,'sto-3g')
nbf = len(bfs)
nocc,nopen = LiH.get_closedopen()
assert nopen==0
S,h,Ints = getints(bfs,LiH)
en,orbe,orbs = rhf(LiH,integrals=(S,h,Ints),verbose=True)
print "SCF completed, E = ",en 
print "S = \n",S
print "h = \n",h
print "I2 = \n",Ints[:5]

from PyQuante.TestMolecules import h2o
bfs = getbasis(h2o,'sto-3g')
nbf = len(bfs)
S,h,Ints = getints(bfs,h2o)
en,orbe,orbs = rhf(h2o,integrals=(S,h,Ints),verbose=True)
print "SCF completed, E = ",en 
print "S = \n",S
print "h = \n",h
print "I2 = \n",Ints[:5]


from pyquante2 import molecule,lih,basisset,rhf
from pyquante2.ints.integrals import onee_integrals,twoe_integrals
from pyquante2.utils import dmat,trace2,geigh

lih = molecule([
    (3,.0000000000, .0000000000, .0000000000),
    (1, .0000000000, .0000000000,1.629912)],
    units='Angstroms')

bfs = basisset(lih,'sto3g')

i1 = onee_integrals(bfs,lih)
i2 = twoe_integrals(bfs)
print "S = \n",i1.S
print "h = \n",i1.T + i1.V
print "I2 = \n",i2._2e_ints[:5]
s = rhf(lih,bfs)
s.converge()

from pyquante2 import molecule,h2o,basisset,rhf
from pyquante2.ints.integrals import onee_integrals,twoe_integrals
from pyquante2.utils import dmat,trace2,geigh

bfs = basisset(h2o,'sto3g')

i1 = onee_integrals(bfs,h2o)
i2 = twoe_integrals(bfs)
print "S = \n",i1.S
print "h = \n",i1.T + i1.V
print "I2 = \n",i2._2e_ints[:5]
s = rhf(h2o,bfs)
s.converge()

from pyquante2 import h2o,basisset,rhf
from pyquante2.utils import dmat
from pyquante2.ints.integrals import twoe_integrals,twoe_integrals_compressed
bfs = basisset(h2o,'sto3g')
solver = rhf(h2o,bfs)
solver.converge()
D = dmat(solver.orbs,h2o.nocc())
i2 = twoe_integrals(bfs)
i2c = twoe_integrals_compressed(bfs)

%timeit twoe_integrals(bfs)

%timeit twoe_integrals_compressed(bfs)

%timeit i2.get_j(D)

%timeit i2c.get_j(D)

(i2.get_j(D)-i2c.get_j(D)).max()

from pyquante2 import h2o,basisset,rhf
from pyquante2.utils import dmat
from pyquante2.ints.integrals import twoe_integrals,twoe_integrals_compressed
bfs = basisset(h2o,'6-31g**')
solver = rhf(h2o,bfs)
solver.converge()
D = dmat(solver.orbs,h2o.nocc())
i2 = twoe_integrals(bfs)
i2c = twoe_integrals_compressed(bfs)