Hi,

I tried to use Fermi on two Macs [M1 first, then Intel] both under Julia 1.6 and 1.8.

the failure on the Intel-Mac worries me, it's running Catalina, so that should be stable.

Fermi doesn't compile, because

Precompiling project...
✗ TBLIS
✗ Fermi
0 dependencies successfully precompiled in 4 seconds (103 already precompiled)
2 dependencies errored. To see a full report either run import Pkg; Pkg.precompile() or load the packages
Testing Running tests...
ERROR: LoadError: InitError: could not load library "/Users/marko/.julia/artifacts/bad39b751ed9b819c0e04ac492bb388e5b494858/lib/libtblis.0.dylib"

what am I doing wrong?

thanks.

marko

Currently, the basis field is a map from Atom to BasisFunction objects. Back working on it I found it a little cumbersome and possibly fragile. So I want to add a field Shells to an array of BasisFunctions (possible name shells) and remove the basis field.

libcint: replace the custom cmake build step with libcint_jll (see #111, solved by PR #112)

deps/lib.gz: this seems to be a data dependency so a few things come to mind:

• just extract the files into git repo (instead of shipping a zipped lib.gz)
• use https://github.com/oxinabox/DataDeps.jl
• (maybe it's possible to use / read the data from the lib.gz without unpacking it (similar to what Pkg.jl is doing with the general registry information)?)

@gustavojra where does the lib.gz come from? Do you compile it manually? Or do you download it from an external URL?

Hello,
I am a bit of a QM newbie I was wondering if it was possible to get calculate automatic charges using fermi ?The application was specifically to calculate pKa of some acids e.g. "Using Atomic Charges to Describe the pKa of Carboxylic Acids". Initial was planning to use xtb but I am always looking for a good excuse to improve by julia code.

Adding spin contamination metric (perhaps even as field of UHF) would be great. We might want that to always be displayed at the end.

Currently, integrals are obtained from an integral helper object by using a dictionary-like call

ints["S"]

which involves the getindex function.

I don't love the current system because it relies on several if statements... It is verbose and hard to extent. I want a way to leverage the multiple-dispatch system such that others can easily implement their own integral helper subtypes.

Thus, we should navigate away from

ints["S"]

and start using a functional form, such as

overlap(ints)

We have the feature of allowing the usage of an old RHF object as an initial guess to another RHF (e.g. small basis to larger basis or different geometries).

It is possible to have the same type of projection in CC theory: https://pubs.acs.org/doi/10.1021/ct100636a

This is a highly desirable feature in order to produce thermochemical protocols, such as Focal Point Analysis.

pyscf:

from pyscf import gto
mol_h2o = gto.M(atom = 'N 0 0 0.2; N 0 0 0', basis = 'sto-3g')

from pyscf import scf
rhf_h2o = scf.RHF(mol_h2o)
e_h2o = rhf_h2o.kernel()

prints converged SCF energy = -52.6419222048247

fermi:

using Fermi
@molecule {
    N        0.0 0.0 0.2
    N        0.0 0.0 0.0
}

@set basis "sto-3g"

hfst = @energy rhf

prints -52.51997436116736

Can someone help me understand this energy difference?

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Highly desired features that we must implement asap. This would be the easiest type of gradients since it is independent on the method.

• Serial
• Parallel

For the serial approach we can easily employ the orbitals and amplitudes from the reference structure as initial guess. For parallel we need to figure out how to make this info available (Since parallel in this case means basically generating inputs for individual computations).

A more involved approach for the Parallel gradient would be exploring Distributed.jl.

function RHF(molecule::Molecule, aoint::IntegralHelper, C::Array{Float64,2}, Λ, Alg::ConventionalRHF)
    @output "Computing integrals ..."
    t = @elapsed aoint["μ"]
    @output " done in {:>5.2f} s\n" t
    RHF(molecule,aoint,C,aoint["μ"],Λ)
end
function RHF(molecule::Molecule, aoint::IntegralHelper, C::Array{Float64,2}, Λ, Alg::DFRHF)
    @output "Computing integrals ..."
    t = @elapsed aoint["B"]
    @output " done in {:>5.2f} s\n" t
    RHF(molecule,aoint,C,aoint["B"],Λ)
end

This Λ needs some specification, right now this method will accept anything for Λ. Also some clarification on what is that via docs or comment would be nice.

Basically nothing is up to date in docs

Hi,

I would have a feature request. (Would love to code it myself and PR, but, not being an expert, this would take me several days. Hoping you can help!)

Fermi runs HF with a simple interface. I would like to use the molecular orbitals and ERIs to form the molecular Hamiltonian as a LinearMap, for use in FCI calculations. Right now I’m using a workaround with ‘PyCall’ and ‘PySCF’. A pure Julia version would be more efficient and more elegant.

Thanks!

CCSD Iterations

The vast majority of use-cases involving restricted CCSD(T) would involve a restricted Hartree-Fock reference. In this case, many of the contractions in the energy, T1, T2 functions are trivially zero, and can be removed:

(T) correction

A small amount of time can be saved by pulling redundant computations out of the loops. Here, the occupied delta term can be pre-saved in the i,j,k loop trio above rather than repeatedly computed in the virtual index loops:

A similar idea can be applied to the orbital energy denominator here; the occupied part can be precomputed before the virtual loops, saving two flops:

Something like the following:

                # Compute Energy contribution
                δjk = Int(j == k)
                δ_occ = 2 - δij - δjk
                Dd_occ = fo[i] + fo[j] + fo[k]
                for a in 1:v
                    for b in 1:a
                        δ_vir = 1 + Int(a == b)
                        for c in 1:b
                            Dd = Dd_occ - fv[a] - fv[b] - fv[c]
                            δ_vir = δ_vir + Int(b == c)
                            @inbounds X = (W[a,b,c]*V[a,b,c] + W[a,c,b]*V[a,c,b] + W[b,a,c]*V[b,a,c] + W[b,c,a]*V[b,c,a] + W[c,a,b]*V[c,a,b] + W[c,b,a]*V[c,b,a])
                            @inbounds Y = (V[a,b,c] + V[b,c,a] + V[c,a,b])
                            @inbounds Z = (V[a,c,b] + V[b,a,c] + V[c,b,a])
                            E = (Y - 2*Z)*(W[a,b,c] + W[b,c,a] + W[c,a,b]) + (Z - 2*Y)*(W[a,c,b]+W[b,a,c]+W[c,b,a]) + 3*X
                            Et += E*δ_occ/(Dd*δ_vir)

This pulls about 15% of the flops out of the innermost loop, which will scale to meaningful savings for large basis sets.

If nothing else, these changes are worth a quick local benchmark. If they do not help much, feel free to close this issue.

Hello, I'm getting a OutOfMemoryError error (aug-cc-pvtz, mp2) is there a way to set memory in the same way as in psi4 ?

Hi @gustavojra,

I'm finding the RHF on H2 fails:

using Fermi

@molecule {
  H  0.0 0.0 0.0
  H  0.76 0.0 0.0
}
@set basis sto-3g
wfn = @energy rhf

with the output:

================================================================================
|                                 Hartree-Fock                                 |
|                                  Module  by                                  |
|                         G.J.R Aroeira and M.M. Davis                         |
================================================================================
Collecting necessary integrals...
Done in    0.57259 s
Using GWH Guess
Molecule:

H    0.000000000000    0.000000000000    0.000000000000
H    0.760000000000    0.000000000000    0.000000000000


Charge: 0   Multiplicity: 1   
Nuclear repulsion:    0.6962858038
 Number of AOs:                            2
 Number of Doubly Occupied Orbitals:       1
 Number of Virtual Spatial Orbitals:       1
nt: 6
 Guess Energy    -1.81166646326759

 Iter.            E[RHF]         ΔE       Dᵣₘₛ        t     DIIS     damp
--------------------------------------------------------------------------------
    1     -1.1153806594  -1.115e+00   0.000e+00     0.76    false      NaN
ERROR: LoadError: ArgumentError: matrix contains Infs or NaNs
Stacktrace:
  [1] chkuplofinite(A::Matrix{Float64}, uplo::Char)
    @ LinearAlgebra.LAPACK /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/lapack.jl:109
  [2] syevr!(jobz::Char, range::Char, uplo::Char, A::Matrix{Float64}, vl::Float64, vu::Float64, il::Int64, iu::Int64, abstol::Float64)
    @ LinearAlgebra.LAPACK /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/lapack.jl:5089
  [3] #eigen!#185
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/symmetriceigen.jl:4 [inlined]
  [4] #eigen#186
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/symmetriceigen.jl:9 [inlined]
  [5] #diagonalize#20
    @ ~/.julia/packages/Fermi/phTqJ/src/Backend/Arrays.jl:189 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/Fermi/phTqJ/src/Methods/HartreeFock/RHF/RHFa.jl:138 [inlined]
  [7] macro expansion
    @ ./timing.jl:299 [inlined]
  [8] macro expansion
    @ ~/.julia/packages/Fermi/phTqJ/src/Methods/HartreeFock/RHF/RHFa.jl:133 [inlined]
  [9] macro expansion
    @ ./timing.jl:299 [inlined]
 [10] Fermi.HartreeFock.RHF(ints::Fermi.Integrals.IntegralHelper{Float64, Fermi.Integrals.SparseERI, Fermi.Orbitals.AtomicOrbitals}, C::FermiMDArray{Float64, 2}, Λ::FermiMDArray{Float64, 2}, Alg::Fermi.HartreeFock.RHFa)
    @ Fermi.HartreeFock ~/.julia/packages/Fermi/phTqJ/src/Methods/HartreeFock/RHF/RHFa.jl:132
 [11] Fermi.HartreeFock.RHF(ints::Fermi.Integrals.IntegralHelper{Float64, Fermi.Integrals.SparseERI, Fermi.Orbitals.AtomicOrbitals}, Alg::Fermi.HartreeFock.RHFa)
    @ Fermi.HartreeFock ~/.julia/packages/Fermi/phTqJ/src/Methods/HartreeFock/RHF/RHFa.jl:32
 [12] Fermi.HartreeFock.RHF(Alg::Fermi.HartreeFock.RHFa)
    @ Fermi.HartreeFock ~/.julia/packages/Fermi/phTqJ/src/Methods/HartreeFock/RHF/RHFa.jl:5
 [13] Fermi.HartreeFock.RHF()
    @ Fermi.HartreeFock ~/.julia/packages/Fermi/phTqJ/src/Methods/HartreeFock/RHF/RHF.jl:81
 [14] top-level scope
    @ ~/Dropbox (Simons Foundation)/workdir/Fermi.jl/hartree_fock_H2.jl:8
 [15] include(fname::String)
    @ Base.MainInclude ./client.jl:451
 [16] top-level scope
    @ REPL[1]:1
in expression starting at /home/mfishman/Dropbox (Simons Foundation)/workdir/Fermi.jl/hartree_fock_H2.jl:8

My setup is:

(@v1.7) pkg> st Fermi
      Status `~/.julia/environments/v1.7/Project.toml`
  [9237668d] Fermi v0.4.0

julia> versioninfo()
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) E-2176M  CPU @ 2.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = vim

Maybe I'm running it wrong.

Thanks,
Matt

We should be able to do something like

@set basis cc-pvdz
wf1 = @energy rhf
wf2 = @energy rhf <- wf1

to read in starting wfs.

According to Julia documentation, it is best to avoid fields with abstract containers.

This is exactly what we do with our orbital structures:

struct RHFOrbitals <: AbstractRestrictedOrbitals
    molecule::Molecule
    basis::String
    eps::AbstractArray{Float64,1}
    sd_energy::Float64
    C::AbstractArray{Float64,2}
end

Thus, it would be better (according to Julia documentation) to write

struct RHFOrbitals{A} <: AbstractRestrictedOrbitals
    molecule::Molecule
    basis::String
    eps::AbstractArray{Float64,1}
    sd_energy::Float64
    C::A
end

Thou, I doubt this would change performance since orbitals are very rarely created. (e.g. at the end of a computation)

@energy CCSD(T)

Will not work because characters like () are removed during the Expr -> String convertion.

If you try

julia> @energy uhf

You get an error if the reference is not set. Since we explicitly asking for UHF I think checking the reference keyword is not necessary. For RHF however, we need to check if the number of electrons is even!

I think the following should be done with the macro @energy

uhf maps to UHF no matter what
rhf maps to RHF but checks the number of electrons (done inside the RHF code already)
scf looks up the keyword reference and maps to uhf or rhf

Trying water with cc-pvtz raises

================================================================================
|    Hartree Fock                                                              |
|        Module by M.M. Davis and G.J.R Aroeira                                |
================================================================================
    executing RHF
    Forming initial Fock matrix ... 
    Nuclear Repulsion:     8.8880641698
done in  0.01s
ERROR: MethodError: no method matching TBLIS.TTensor{Complex{Float64}}(::Array{Complex{Float64},2}, ::Float64)
Closest candidates are:
  TBLIS.TTensor{Complex{Float64}}(::Any, ::Any, ::Any, ::Any, ::Any) where T at /home/aroeira/.julia/dev/TBLIS/src/TTensor.jl:12
Stacktrace:
 [1] contract!(::Array{Complex{Float64},2}, ::Array{Complex{Float64},2}, ::Array{Complex{Float64},2}, ::Float64, ::Float64, ::Float64, ::String, ::String, ::String, ::Fermi.Environments.No_IC, ::Fermi.Environments.NoCommunicator, ::Fermi.Environments.NoAccelerator, ::Fermi.Environments.TBLIS) at /home/aroeira/Fermi.jl/src/Backend/contract.jl:77
 [2] contract!(::Array{Complex{Float64},2}, ::Array{Complex{Float64},2}, ::Array{Complex{Float64},2}, ::Float64, ::Float64, ::Float64, ::String, ::String, ::String) at /home/aroeira/Fermi.jl/src/Backend/contract.jl:44
 [3] contract(::Array{Complex{Float64},2}, ::Array{Complex{Float64},2}, ::Float64, ::Float64, ::String, ::String) at /home/aroeira/Fermi.jl/src/Backend/contract.jl:34
 [4] contract(::Array{Complex{Float64},2}, ::Array{Complex{Float64},2}, ::String, ::String) at /home/aroeira/Fermi.jl/src/Backend/contract.jl:29
 [5] Fermi.HartreeFock.RHFWavefunction(::Fermi.Geometry.Molecule, ::Fermi.Integrals.ConventionalAOIntegrals{Float64}, ::Fermi.HartreeFock.ConventionalRHF) at /home/aroeira/Fermi.jl/src/Methods/HartreeFock/RHF/ConventionalRHF.jl:60
 [6] Fermi.HartreeFock.RHFWavefunction(::Fermi.Geometry.Molecule, ::Fermi.Integrals.ConventionalAOIntegrals{Float64}) at /home/aroeira/Fermi.jl/src/Methods/HartreeFock/RHF/ConventionalRHF.jl:29
 [7] RHFWavefunction at /home/aroeira/Fermi.jl/src/Methods/HartreeFock/RHF/ConventionalRHF.jl:19 [inlined]
 [8] Fermi.HartreeFock.RHFWavefunction() at /home/aroeira/Fermi.jl/src/Methods/HartreeFock/RHF/ConventionalRHF.jl:3
 [9] top-level scope at REPL[7]:1

Have you considered using https://github.com/JuliaBinaryWrappers/libcint_jll.jl instead of manually downloading and compiling libcint in deps/builds.jl?

Using a custom deps/build.jl isn't the recommended approach for distributing binary dependencies anymore. Using Julia's artifact system is much cleaner.

👋 Hi there!

As you may be aware, some of us over at JuliaMolSim have been developing AtomsBase, an interface package for specifying atomic geometries. I think Fermi.jl would be a great candidate to support AtomsBase and aid interoperability with other atomistic tools in Julia!

I thought I'd put this on your radar, and would be happy to discuss the best way to go about implementing the interface, and/or contribute to a PR to do so!

It is pointless since a BasisSet is mandatory and it contains molecule already. Moreover, it can be a source of bugs. Currently (v0.3.0). Fermi will ignore molecule inputs given directly to the IntegralHelper. I have that fixed, needs PR.

Currently we have defined CCIntegrals which is a type that is passed into the code so we know which CC formulation is appropriated to handle that type of integrals.

Should this be generalized for all methods?

DF-MP2, DF-HF and DF-CC all depende on a DF structured to be passed on.

I think we should create an AbstractIntegral type, where struct (instances) of it can be DensityFitting or Conventional.

This being defined at higher level can lead to a clean up on the methods. Also, instead of having the integrals being handled at the HF part we should store the integrals inside the concrete Integral type.

Here's how it would work.

RHF is called as

RHF(AOIntegrals::Conventional)

RHF(AOIntegrals::DensityFitting)

etc...

So we need to pass integrals into the RHF module. Of course that we also create a simple function

RHF()
   #Use lints to create integral object, call it X
   RHF(X)
end

Thanks for comming to my tedtalk

Sorry this isnt an issue. I was looking for some place to start discussions and post ideas.

Hi @gustavojra,

Thanks for the nice package!

Is there a simple way to obtain the coefficients of the second quantized Hamiltonian (the one- and two-body integrals) from a Hartree-Fock solution?

Thanks,
Matt

When invalid arguments are passed into a function, e.g. Fermi.MollerPlesset.RMP2 the error message is too long because it prints the actual arguments (that could be huge arrays).

It is best if change it to show only the types. For example:

function RMP2(x...)
    if !any(i-> i isa RMP2Algorithm, x)
        RMP2(x..., get_rmp2_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for RMP2 method: $args"))
    end
end

In the original paper of oda, the criteria for damping value is:
$$\lambda_m = \begin{cases} 1, & c\le-s/2 \ -s/2c, & \text{otherwise} \end{cases}$$.
However, the code seems discrepent with the formula.
Is this a mistake or am I missing something here?

Cluster decomposition is by far the bottleneck for ecCC jobs.

Print functions inside Base.show for custom struct need to receive io object.

Probably very minor. It happened when I ran HF (molecule) with a (10,6) active space

ERROR: BoundsError: attempt to access 36-element Array{Float64,1} at index [37]
Stacktrace:
 [1] getindex(::Array{Float64,1}, ::Int64) at ./array.jl:788
 [2] macro expansion at /home/vulcan/aroeira/Fermi.jl/src/Backend/IO/Output.jl:57 [inlined]
 [3] Fermi.ConfigurationInteraction.CASCI{Float64}(::Fermi.HartreeFock.RHF, ::Array{Float64,2}, ::Array{Float64,4}, ::Int64, ::Int64, ::Int64, ::Fermi.ConfigurationInteraction.SparseHamiltonian) at /home/vulcan/aroeira/Fermi.jl/src/Methods/ConfigurationInteraction/CASCI/SparseHamiltonian.jl:97
 [4] Fermi.ConfigurationInteraction.CASCI{Float64}(::Fermi.ConfigurationInteraction.SparseHamiltonian) at /home/vulcan/aroeira/Fermi.jl/src/Methods/ConfigurationInteraction/CASCI/SparseHamiltonian.jl:48
 [5] Fermi.ConfigurationInteraction.CASCI{Float64}() at /home/vulcan/aroeira/Fermi.jl/src/Methods/ConfigurationInteraction/CASCI/CASCI.jl:69
 [6] Fermi.ConfigurationInteraction.CASCI() at /home/vulcan/aroeira/Fermi.jl/src/Methods/ConfigurationInteraction/CASCI/CASCI.jl:58
 [7] top-level scope at REPL[12]:1

I was browsing your repo and noticed that the Atom typed defined in Molecules.jl is parametric in {I,F}: but here, this parametric dependence is not accounted for.
As a result, any access of the elements of the atoms field's elements will be type-unstable.

For C2 TZ

   • CAS Decomposition started:
Active Occupied Orbitals: [3, 4, 5, 6]
Active Virtual Orbitals:  [7, 8, 9, 10]
Getting T1 and T2...Done.
Computing ecT1, ecT2
3 out of 6
ERROR: BoundsError: attempt to access 4×4×4×4 Array{Float64,4} at index [6, 6, 2, 2]
Stacktrace:
 [1] setindex!(::Array{Float64,4}, ::Float64, ::Int64, ::Int64, ::Int64, ::Vararg{Int64,N} where N) at ./array.jl:828
 [2] get_casT4αβ!(::Array{Float64,4}, ::Int64, ::Int64, ::Int64, ::Int64, ::Array{Float64,1}, ::Array{Fermi.ConfigurationInteraction.Determinant,1}, ::Array{Float64,1}, ::Array{Fermi.ConfigurationInteraction.Determinant,1}, ::Fermi.ConfigurationInteraction.Determinant, ::Int64, ::Int64, ::Array{Float64,2}, ::Array{Float64,4}, ::Array{Float64,4}, ::Array{Float64,4}) at /home/mmd01986/.julia/dev/Fermi.jl/src/Methods/CoupledCluster/ecRCCSD/CTF.jl:509
 [3] cas_decomposition(::Tuple{Fermi.ConfigurationInteraction.Determinant,Array{Float64,1},Array{Fermi.ConfigurationInteraction.Determinant,1},Array{Float64,1},Array{Fermi.ConfigurationInteraction.Determinant,1},Array{Float64,1},Array{Fermi.ConfigurationInteraction.Determinant,1}}, ::Int64, ::Int64, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float64,2}, ::Array{Float64,4}, ::Array{Float64,4}, ::Array{Float64,4}) at /home/mmd01986/.julia/dev/Fermi.jl/src/Methods/CoupledCluster/ecRCCSD/CTF.jl:925
 [4] macro expansion at ./util.jl:234 [inlined]
 [5] Fermi.CoupledCluster.ecRCCSD{Float64}(::Fermi.ConfigurationInteraction.CASCI{Float64}, ::Fermi.CoupledCluster.CTF) at /home/mmd01986/.julia/dev/Fermi.jl/src/Methods/CoupledCluster/ecRCCSD/CTF.jl:51
 [6] top-level scope at REPL[117]:1
 [7] eval(::Module, ::Any) at ./boot.jl:331
 [8] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [9] run_backend(::REPL.REPLBackend) at /home/mmd01986/.julia/packages/Revise/BqeJF/src/Revise.jl:1184
 [10] top-level scope at REPL[1]:0