def test_lih_averaging(self):
     bfs = basisset(lih,'sto-3g')
     hamiltonian = rhf(bfs)
     iterator = AveragingIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -7.860746149768, 6)

 def test_oh(self):
     bfs = basisset(oh,'sto-3g')
     hamiltonian = rohf(bfs)
     iterator = ROSCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -74.359151530162, 5)

def plot_orbs():
    bfs = basisset(h2,'sto3g')
    orbs = np.array([[1.0,1.0],
                     [1.0,-1.0]],'d')
    
    vtk_orbital(h2,orbs,bfs)
    return

 def test_h2o_averaging(self):
     bfs = basisset(h2o,'sto-3g')
     hamiltonian = rhf(bfs)
     iterator = AveragingIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -74.959857776754, 5)

def plot_h2_lineplot():
    bfs = basisset(h2,'6-31G')
    solver = rhf(h2,bfs)
    solver.converge()
    points = [(0,0,z) for z in np.linspace(-5,5)]
    lineplot_orbs(points,solver.orbs[:,:2],bfs,True)
    return

 def test_he_triplet_xlda(self):
     from pyquante2.geo.samples import he
     he.multiplicity = 3
     bfs = basisset(he,'sto-3g')
     solver = dft(he,bfs,'xs')
     ens = solver.converge()
     self.assertAlmostEqual(solver.energy,-1.7819689849)

def pyq2_dft(atomtuples=[(2,0,0,0)],basis = '6-31G**',maxit=10,xcname='svwn'):
    import pyquante2 as pyq2
    print ("pyq2 DFT run")
    geo = pyq2.molecule(atomtuples)
    bfs = pyq2.basisset(geo,name=basis)
    i1 = pyq2.onee_integrals(bfs,geo)
    i2 = pyq2.twoe_integrals(bfs)
    grid = pyq2.grid(geo)
    h = i1.T + i1.V
    orbe,orbs = pyq2.geigh(h,i1.S)
    eold = 0
    grid.setbfamps(bfs)
    E0 = geo.nuclear_repulsion()

    for i in range(maxit):
        D = pyq2.dmat(orbs,geo.nocc())
        E1 = 2*pyq2.trace2(h,D)

        J = i2.get_j(D)
        Ej = 2*pyq2.trace2(J,D)

        Exc,Vxc = pyq2.get_xc(grid,0.5*D,xcname=xcname)

        energy = E0+E1+Ej+Exc
        F = h+2*J+Vxc

        orbe,orbs = pyq2.geigh(F,i1.S)

        print (i,energy,E1,Ej,Exc,E0)
        if np.isclose(energy,eold):
            break
        eold = energy
    return energy

 def test_li(self):
     bfs = basisset(li,'sto-3g')
     hamiltonian = uhf(bfs)
     iterator = USCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -7.315525981280, 6)

 def test_h2_631(self):
     bfs = basisset(h2,'6-31gss')
     hamiltonian = rhf(bfs)
     iterator = SCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -1.131333590574, 7)

 def test_h2(self):
     bfs = basisset(h2,'sto-3g')
     hamiltonian = rhf(bfs)
     iterator = SCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -1.117099435262, 7)

 def test_CH4(self):
     """CH4 symmetry Td"""
     bfs = basisset(ch4,'sto-3g')
     hamiltonian = rhf(bfs, twoe_factory=libint_twoe_integrals)
     iterator = SCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -39.726862723517, 6)

def plot_orb2():
    from pyquante2 import basisset,h2
    bfs = basisset(h2,'sto3g')
    orbs = np.array([[1.0,1.0],
                     [1.0,-1.0]],'d')
    
    vtk_orbs(h2,orbs,bfs)
    return

 def test_CH4(self):
     """CH4 symmetry Td"""
     bfs = basisset(CH4, "sto-3g")
     solver = rhf(CH4, bfs, libint=True)
     ens = solver.converge()
     hp5("Ham.h5", solver.i1.V + solver.i1.T)
     hp5("TwoS.h5", solver.i2._2e_ints)
     self.assertPrecisionEqual(solver.energy, -39.72591203477140)

 def test_h2(self):
     bfs = basisset(h2,'6-31g**')
     solver=rhf(h2,bfs)
     solver.converge()
     nvirt = len(bfs)-h2.nocc()
     emp2 = mp2(solver.i2,solver.orbs,solver.orbe,h2.nocc(),len(bfs)-h2.nocc())
     self.assertAlmostEqual(emp2,-0.02632654197486595)
     return

 def test_H2(self):
     bfs = basisset(h2,'cc-pvdz')
     hamiltonian = rhf(bfs, twoe_factory=libint_twoe_integrals)
     iterator = SCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     nvirt = len(bfs)-h2.nocc()
     eccsd = ccsd(hamiltonian, iterator.orbs, iterator.orbe, h2.nocc(), nvirt)
     self.assertAlmostEqual(eccsd, -0.034544318453406, 8)

 def test_HBr(self):
     """HBr"""
     HBr = read_xyz('./molfiles/HBr.xyz')
     bfs = basisset(HBr,'sto-3g')
     hamiltonian = rhf(bfs, twoe_factory=libint_twoe_integrals)
     iterator = SCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -2545.887434128302, 4)

 def test_C2H2Cl2(self):
     """C2H2Cl2 symmetry C2H"""
     C2H2Cl2 = read_xyz('./molfiles/C2H2Cl2.xyz')
     bfs = basisset(C2H2Cl2,'sto-3g')
     hamiltonian = rhf(bfs, twoe_factory=libint_twoe_integrals)
     iterator = SCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -967.533150337277, 4)

def plot_h2_vtk():
    bfs = basisset(h2,'sto3g')
    solver = rhf(h2,bfs)
    ens = solver.converge()

    # Note: these orbitals are not coming out symmetric. Why not??
    print(solver)
    print(solver.orbs)
    vtk_orbs(h2,solver.orbs,bfs,npts=10)

 def test_CF3(self):
     """CF3 radical"""
     CF3 = read_xyz('./molfiles/CF3.xyz')
     bfs = basisset(CF3,'sto-3g')
     hamiltonian = cuhf(bfs, twoe_factory=libint_twoe_integrals)
     iterator = USCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -331.480688906400, 5)

 def test_H2O_4(self):
     """H2O tethramer symmetry S4"""
     H2O4 = read_xyz('./molfiles/H2O_4.xyz')
     bfs = basisset(H2O4,'sto-3g')
     hamiltonian = rhf(bfs, twoe_factory=libint_twoe_integrals)
     iterator = SCFIterator(hamiltonian)
     iterator.converge()
     self.assertTrue(iterator.converged)
     self.assertAlmostEqual(iterator.energy, -299.909789863537, 5)