FermiQC / Fermi.jl

Fermi quantum chemistry program
MIT License
136 stars 25 forks source link

Using RHF as a guess to new RHF is broken #40

Closed gustavojra closed 3 years ago

gustavojra commented 4 years ago
julia> Fermi.HartreeFock.RHF(rhf)
ERROR: type RHF has no field basis
Stacktrace:
 [1] getproperty(::Fermi.HartreeFock.RHF, ::Symbol) at ./Base.jl:33
 [2] macro expansion at /Users/aroeira/Fermi/Fermi.jl/src/Backend/IO/Output.jl:57 [inlined]
 [3] Fermi.HartreeFock.RHF(::Fermi.HartreeFock.RHF, ::Fermi.Integrals.IntegralHelper{Float64}, ::Fermi.HartreeFock.ConventionalRHF) at /Users/aroeira/Fermi/Fermi.jl/src/Methods/HartreeFock/RHF/RHF.jl:176
 [4] Fermi.HartreeFock.RHF(::Fermi.HartreeFock.RHF) at /Users/aroeira/Fermi/Fermi.jl/src/Methods/HartreeFock/RHF/RHF.jl:175
 [5] top-level scope at REPL[170]:1
gustavojra commented 4 years ago

This is due to the change in the Integrals.jl.

Seems like having only the RHF object is hard to find the basis that was used to create it?

Beforehand we had fields basis and LintsBasis. I don't know where they are now.

gustavojra commented 4 years ago

I created a temporary solution

function RHF(wfn::RHF, aoint::IntegralHelper, Alg::B) where B <: RHFAlgorithm

    # Projection of A→ B done using equations described in Werner 2004 
    # https://doi.org/10.1080/0026897042000274801
    @output "Using {} wave function as initial guess\n" wfn.ints.bname["primary"]

    Ca = Float64[]
    for orb in wfn.ints.orbs["FU"].orbs
        Ca = vcat(Ca, orb.C)
    end

    nbf = Int(√length(Ca))
    Ca = reshape(Ca, (nbf, nbf))

    Sbb = aoint["S"]
    println(typeof(Sbb))
    S = Hermitian(aoint["S"])
    Λ = Array(S^(-1/2))

    open("/tmp/molfile1.xyz","w") do molfile
        natom = length(wfn.molecule.atoms)
        write(molfile,"$natom\n\n")
        write(molfile,Fermi.Geometry.get_xyz(wfn.molecule))
    end

    Lints.libint2_init()
    wfnmol = Lints.Molecule("/tmp/molfile1.xyz")
    wfnbas = Lints.BasisSet(wfn.ints.bname["primary"], wfnmol)

    molecule = Fermi.Geometry.Molecule()
    open("/tmp/molfile2.xyz","w") do molfile
        natom = length(molecule.atoms)
        write(molfile,"$natom\n\n")
        write(molfile,Fermi.Geometry.get_xyz(molecule))
    end

    mol = Lints.Molecule("/tmp/molfile2.xyz")
    bas = Lints.BasisSet(aoint.bname["primary"], wfnmol)

    Sab = Lints.projector(wfnbas, bas)
    T = transpose(Ca)*Sab*(Sbb^-1)*transpose(Sab)*Ca
    Cb = (Sbb^-1)*transpose(Sab)*Ca*T^(-1/2)
    Cb = real.(Cb)
    RHF(molecule, aoint, Cb, Λ, Alg)
end

Obviously this can use a massive simplification once we clean up integrals and orbitals.