# ## Basis Set Data
# Note use of curly braces here. Julia assumes that if you have square braces, you want
# things flattened as much as possible (to be as fast as possible, I guess). Curlys 
# preserve the list structure the way I would expect from Python

const sto3g = [
    # H
    [('S',
      [(3.4252509099999999, 0.15432897000000001),
       (0.62391373000000006, 0.53532813999999995),
       (0.16885539999999999, 0.44463454000000002)])],
    # He
    [('S',
      [(6.3624213899999997, 0.15432897000000001),
       (1.1589229999999999, 0.53532813999999995),
       (0.31364978999999998, 0.44463454000000002)])],
    # Li
    [('S',
      [(16.119575000000001, 0.15432897000000001),
       (2.9362007000000001, 0.53532813999999995),
       (0.79465050000000004, 0.44463454000000002)]),
     ('S',
      [(0.63628969999999996, -0.099967230000000004),
       (0.14786009999999999, 0.39951282999999999),
       (0.048088699999999998, 0.70011546999999996)]),
     ('P',
      [(0.63628969999999996, 0.15591627),
       (0.14786009999999999, 0.60768372000000004),
       (0.048088699999999998, 0.39195739000000002)])],
    # Be
    [('S',
      [(30.167871000000002, 0.15432897000000001),
       (5.4951153000000001, 0.53532813999999995),
       (1.4871927, 0.44463454000000002)]),
     ('S',
      [(1.3148331, -0.099967230000000004),
       (0.3055389, 0.39951282999999999),
       (0.099370700000000006, 0.70011546999999996)]),
     ('P',
      [(1.3148331, 0.15591627),
       (0.3055389, 0.60768372000000004),
       (0.099370700000000006, 0.39195739000000002)])],
    # B
    [('S',
      [(48.791113000000003, 0.15432897000000001),
       (8.8873622000000001, 0.53532813999999995),
       (2.4052669999999998, 0.44463454000000002)]),
     ('S',
      [(2.2369561, -0.099967230000000004),
       (0.51982050000000002, 0.39951282999999999),
       (0.16906180000000001, 0.70011546999999996)]),
     ('P',
      [(2.2369561, 0.15591627),
       (0.51982050000000002, 0.60768372000000004),
       (0.16906180000000001, 0.39195739000000002)])],
    # C
    [('S',
      [(71.616837000000004, 0.15432897000000001),
       (13.045095999999999, 0.53532813999999995),
       (3.5305122, 0.44463454000000002)]),
     ('S',
      [(2.9412493999999998, -0.099967230000000004),
       (0.68348310000000001, 0.39951282999999999),
       (0.22228990000000001, 0.70011546999999996)]),
     ('P',
      [(2.9412493999999998, 0.15591627),
       (0.68348310000000001, 0.60768372000000004),
       (0.22228990000000001, 0.39195739000000002)])],
    # N
    [('S',
      [(99.106168999999994, 0.15432897000000001),
       (18.052312000000001, 0.53532813999999995),
       (4.8856602000000002, 0.44463454000000002)]),
     ('S',
      [(3.7804559000000002, -0.099967230000000004),
       (0.87849659999999996, 0.39951282999999999),
       (0.28571439999999998, 0.70011546999999996)]),
     ('P',
      [(3.7804559000000002, 0.15591627),
       (0.87849659999999996, 0.60768372000000004),
       (0.28571439999999998, 0.39195739000000002)])],
    # O
    [('S',
      [(130.70931999999999, 0.15432897000000001),
       (23.808861, 0.53532813999999995),
       (6.4436083000000002, 0.44463454000000002)]),
     ('S',
      [(5.0331513000000001, -0.099967230000000004),
       (1.1695960999999999, 0.39951282999999999),
       (0.38038899999999998, 0.70011546999999996)]),
     ('P',
      [(5.0331513000000001, 0.15591627),
       (1.1695960999999999, 0.60768372000000004),
       (0.38038899999999998, 0.39195739000000002)])],
    # F
    [('S',
      [(166.67912999999999, 0.15432897000000001),
       (30.360811999999999, 0.53532813999999995),
       (8.2168206999999995, 0.44463454000000002)]),
     ('S',
      [(6.4648032000000004, -0.099967230000000004),
       (1.5022812000000001, 0.39951282999999999),
       (0.48858849999999998, 0.70011546999999996)]),
     ('P',
      [(6.4648032000000004, 0.15591627),
       (1.5022812000000001, 0.60768372000000004),
       (0.48858849999999998, 0.39195739000000002)])],
    # Ne
    [('S',
       [(207.01561000000001, 0.15432897000000001),
        (37.708151000000001, 0.53532813999999995),
        (10.205297, 0.44463454000000002)]),
      ('S',
       [(8.2463151000000003, -0.099967230000000004),
        (1.9162661999999999, 0.39951282999999999),
        (0.62322929999999999, 0.70011546999999996)]),
      ('P',
       [(8.2463151000000003, 0.15591627),
        (1.9162661999999999, 0.60768372000000004),
            (0.62322929999999999, 0.39195739000000002)])]
]
basis_set_data = Dict("sto3g" => sto3g)

mutable struct BasisSet # list of CGBFs
    bfs::Array{CGBF,1}
end

basisset() = BasisSet(CGBF[])

function push!(basis::BasisSet,cbf::CGBF)
    Base.push!(basis.bfs,cbf)
end

function build_basis(mol::Molecule,name="sto3g")
    data = basis_set_data[name]
    basis_set = basisset()
    for atom in mol.atomlist
        for btuple in data[atom.atno]
            sym,primlist = btuple
            for (I,J,K) in sym2power[sym]
                cbf = cgbf(atom.x,atom.y,atom.z,I,J,K)
                push!(basis_set,cbf)
                for (expn,coef) in primlist
                    push!(cbf,expn,coef)
                end
            end
        end
    end
    return basis_set
end

const sym2power = Dict(
    'S' => [(0,0,0)],
    'P' => [(1,0,0),(0,1,0),(0,0,1)],
    'D' => [(2,0,0),(0,2,0),(0,0,2),(1,1,0),(1,0,1),(0,1,1)]
    )


function test_geo_basis()
    @testset "test_geo_basis" begin
    @test nuclear_repulsion(h2) ≈ 0.7223600367
    @test nel(h2) == 2
    @test nel(h2o) == 10
    @test length(sto3g)==10
    bfs = build_basis(h2)
    @test length(bfs.bfs)==2
    l,r = bfs.bfs
    @test overlap(l,l)      ≈ 1
    @test overlap(r,r)      ≈ 1
    @test overlap(l,r)      ≈ 0.6647387449282997
    @test kinetic(l,l)      ≈ 0.76003188
    @test kinetic(r,r)      ≈ 0.76003188
    @test kinetic(l,r)      ≈ 0.24141861181119084
    @test coulomb(l,l,l,l)  ≈ 0.7746059439196398
    @test coulomb(r,r,r,r)  ≈ 0.7746059439196398
    @test coulomb(l,l,r,r)  ≈ 0.5727937653511646
    @test coulomb(l,l,l,r)  ≈ 0.4488373301593464
    @test coulomb(l,r,l,r)  ≈ 0.3025451156654606
    bfs = build_basis(h2o)

    s1,s2,px,py,pz,hl,hr = bfs.bfs
    @test coulomb(s1,s2,hl,hr) ≈ 0.03855344493645537
    @test coulomb(s1,pz,hl,hr) ≈ -0.0027720110485359053
    @test coulomb(s1,hl,pz,hr) ≈ -0.010049491284827426
    @test coulomb(s1,py,hl,hr)==0
    @test coulomb(s1,hl,py,hr)==0
    end #testset
end

