General combinations of functions with shared parameters

[JamesB-1qbit]

I would like to create a basis set with the following 5 functions and unique sets of parameters alpha, beta, gamma, and position a.

1) 1S type orbital GTO(alpha, [0,0,0]) 2) 2S type orbital GTO(beta, [0,0,0]) 2P type delocalized orbitals 3) GTO(gamma, [a, 0, 0]) - GTO(gamma, [-a, 0, 0]) 4) GTO(gamma, [0, a, 0]) - GTO(gamma, [0, -a, 0]) 5) GTO(gamma, [0, 0, a]) - GTO(gamma, [0, 0, -a])

Where GTO are Gaussian Type Orbitals and alpha/beta/gamma are the GTO parameters and "a" is the delocalization.

[frankwswang]

You can write the following code to realize your goals. I took some liberty to assume some information about the basis functions while constructing the basis set, as you didn't provide enough information to completely describe them. For example, I wasn't sure how you specifically define alpha, beta, and gamma in terms of the coefficients of GTOs, so I assumed they are exponent coefficients. You can modify the code to meet your specific configuration.

using Quiqbox # Quiqbox version: >= 0.5.6, < 0.6

# Constructing the four primitive (input) parameters: α, β, γ, a.
##          |in_val |out_var_sym |in_var_sym |optimize w.r.t in_var
α = ParamBox(1.0,    :ζ_α,        :α,         canDiff=true)
β = ParamBox(1.0,    :ζ_β,        :β,         canDiff=true)
γ = ParamBox(1.0,    :ζ_γ,        :γ,         canDiff=true)
a = ParamBox(1.0,    :c_a,        :a,         canDiff=true)

# Constructing the contraction coefficients from input variables named as "dx" 
# for the three GTOs that will be subtracted, in case you want to keep the 
# subtraction even after alternating their (input) values.
d1 = ParamBox(1.0, :d, :dx, canDiff=true)
d2 = ParamBox(1.0, :d, :dx, canDiff=true)
d3 = ParamBox(1.0, :d, :dx, canDiff=true)

# Constructing `-a` in the GTO center coordinates.
a2 = changeMapping(a, -)

# Specify your angular momentums for the P subshell.
p1 = (1, 0, 0)
p2 = (0, 1, 0)
p3 = (0, 0, 1)

# Generating the basis functions.
bf1s  = genBasisFunc([0., 0., 0.], (α, 1.0)) # S orbital is set in default.
bf2s  = genBasisFunc([0., 0., 0.], (β, 1.0), (0, 0, 0))
bf2p1 = genBasisFunc([a , 0., 0.], (γ, 1.0), p1)  +  (-1) * 
        genBasisFunc([a2, 0., 0.], (γ, d1 ), p1)
bf2p2 = genBasisFunc([0., a , 0.], (γ, 1.0), p2)  +  (-1) * 
        genBasisFunc([0., a2, 0.], (γ, d2 ), p2)
bf2p3 = genBasisFunc([0., 0., a ], (γ, 1.0), p3)  +  (-1) * 
        genBasisFunc([0., 0., a2], (γ, d3 ), p3)

# Constructing the basis set.
bs = [bf1s, bf2s, bf2p1, bf2p2, bf2p3]

# Checking all the differentiable parameters with unique in_var (α, β, γ, a, dx).
filter(isDiffParam, markParams!(bs, true))

[frankwswang]

Hi @JamesB-1qbit! Thank you for opening the issue! If you don't have further questions, I'll close this issue soon.

[JamesB-1qbit]

Thank you for the help. However, it appears that Quiqbox does not appear to like mixing floats and ParamBox in general so I had to do something like.

using Quiqbox # Quiqbox version: >= 0.5.6, < 0.6

# Constructing the four primitive (input) parameters: α, β, γ, a.
##          |in_val |out_var_sym |in_var_sym |optimize w.r.t in_var
α = ParamBox(1.0,    :ζ_α,        :α,         canDiff=true)
β = ParamBox(4.0,    :ζ_β,        :β,         canDiff=true)
γ = ParamBox(1.0,    :ζ_γ,        :γ,         canDiff=true)
a = ParamBox(1.0,    :c_a,        :a,         canDiff=true)
z = ParamBox(0.0,    :c_z,        :z,         canDiff=false)

# Constructing the contraction coefficients from input variables named as "dx"
# for the three GTOs that will be subtracted, in case you want to keep the
# subtraction even after alternating their (input) values.
d1 = ParamBox(1.0, :d, :dx, canDiff=true)
d2 = ParamBox(1.0, :d, :dx, canDiff=true)
d3 = ParamBox(1.0, :d, :dx, canDiff=true)

# Constructing `-a` in the GTO center coordinates.
a2 = changeMapping(a, -)

# Specify your angular momentums for the P subshell.
p1 = (1, 0, 0)
p2 = (0, 1, 0)
p3 = (0, 0, 1)

genSpatialPoint([a, z, z])
# Generating the basis functions.
bf1s  = genBasisFunc([0., 0., 0.], GaussFunc(α, 1.0), [(0,0,0)]) # S orbital is set in default.
bf2s  = genBasisFunc([0., 0., 0.], GaussFunc(β, 1.0))
bf2p1 = genBasisFunc([a , z, z], GaussFunc(γ, d1), p1)  +  (-1) *
        genBasisFunc([a2, z, z], GaussFunc(γ, d1), p1)
bf2p2 = genBasisFunc([z, a , z], GaussFunc(γ, d2), p2)  +  (-1) *
        genBasisFunc([z, a2, z], GaussFunc(γ, d2), p2)
bf2p3 = genBasisFunc([z, z, a ], GaussFunc(γ, d3), p3)  +  (-1) *
        genBasisFunc([z, z, a2], GaussFunc(γ, d3), p3)

# Constructing the basis set.
bs = [bf1s, bf2s, bf2p1, bf2p2, bf2p3]

# Checking all the differentiable parameters with unique in_var (α, β, γ, a, dx).
filter(isDiffParam, markParams!(bs, true))

[frankwswang]

When you use a Number instead of a ParamBox, either directly or indirectly (by encapsulating it into an array/tuple) as the input argument to construct any parameterized container, what the constructor function (e.g., genBasisFunc) does is constructing a ParamBox on the fly with it as the input value and itself as the mapping function.

Thus, in my example, all the 0. in the center coordinates of 2p orbitals became independent parameters with their input value being 0.

When you replace them with the same ParamBox, z, genBasisFunc assumes you want to maintain the reference to the original z when using it as multiple basis function parameters at the same time. This is how you reduce the overall number of basis set parameters, and it is one of the ways to correlate GTO parameters.

To make a simplified analogy of what happened under the hood, you can think of the parameters of a basis function as a Vector{Array{T, 0}} that holds all the "ParamBox" as if they are Array{T, 0} (T is a concrete subtype of Number). There are two methods of a parameter generator genPar for the input argument:

genPar(v::Number) = fill(v)

genPar(v::Array{<:Any, 0}) = v

such that

genBfPars(inputPars::Vector) = genPar.(inputPars)

The difference between using a Number and a "ParamBox" is then similar to:

cz = fill(3.0)
cen1 = [1.0, 2.0, cz]
cen2 = [1.0, 2.0, cz]

bf1Pars = genBfPars(cen1)

bf2Pars = genBfPars(cen2)

bf1Pars[1][] == bf2Pars[1][] # Should return `true`

bf1Pars[1][] = 1.1

bf1Pars[1][] == bf2Pars[1][] # Should return `false`

bf1Pars[3][] == bf2Pars[3][] # Should return `true`

bf1Pars[3][] = 3.3

bf1Pars[3][] == bf2Pars[3][] # Should still return `true`

Overall, this is a design choice to give the user the freedom to control the balance between flexibility and efficiency of tuning basis set parameters.

[JamesB-1qbit]

Thank you for the explanation of your theory. It makes sense. I suppose what I meant to say is that the code you provide threw argument type errors.

ERROR: LoadError: MethodError: no method matching genBasisFunc(::SpatialPoint{Float64, 3, Tuple{ParamBox{Float64, :X, typeof(itself)}, ParamBox{Float64, :Y, typeof(itself)}, ParamBox{Float64, :Z, typeof(itself)}}}, ::Tuple{ParamBox{Float64, :ζ_α, typeof(itself)}, Float64})
Closest candidates are:
  genBasisFunc(::SpatialPoint{T, D}, ::Tuple, ::String, ::Tuple{Vararg{Bool}}; normalizeGTO) where {T, D} at ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:595
  genBasisFunc(::Union{Tuple, AbstractVector}, ::Any...; kws...) at ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:617
  genBasisFunc(::SpatialPoint{T}, ::GaussFunc{T}, ::Any...; kws...) where T at ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:614
  ...
Stacktrace:
 [1] genBasisFunc(coord::Vector{Float64}, args::Tuple{ParamBox{Float64, :ζ_α, typeof(itself)}, Float64}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Quiqbox ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:617
 [2] genBasisFunc(coord::Vector{Float64}, args::Tuple{ParamBox{Float64, :ζ_α, typeof(itself)}, Float64})
   @ Quiqbox ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:617
 [3] top-level scope
   @ ~/Documents/Coding/Quiqbox.jl/examples/github.jl:26
in expression starting at /Users/jamesbrown/Documents/Coding/Quiqbox.jl/examples/github.jl:26

Adding in GaussFunc to the arguments fixed this one. i.e. changing line 26 to

bf1s  = genBasisFunc([0., 0., 0.], GaussFunc(α, 1.0), [(0,0,0)]) # S orbital is set in default.

Then the combination of ParamBox with a float was not a valid argument for genSpatialPoint threw an error on line 28.

ERROR: LoadError: MethodError: no method matching genSpatialPoint(::Tuple{ParamBox{Float64, :c_a, typeof(itself)}, Float64, Float64})
Closest candidates are:
  genSpatialPoint(::Tuple{Vararg{Union{Array{T, 0}, T}, D}}) where {D, T<:AbstractFloat} at ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:201
  genSpatialPoint(::Pair{Array{T, 0}, Symbol}, ::Int64) where T<:AbstractFloat at ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:252
  genSpatialPoint(::Pair{Array{T, 0}, Symbol}, ::Int64, ::Function; canDiff) where T<:AbstractFloat at ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:252
  ...
Stacktrace:
 [1] genSpatialPoint(v::Vector{Any})
   @ Quiqbox ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:200
 [2] genBasisFunc(::Vector{Any}, ::GaussFunc{Float64, typeof(itself), typeof(itself)}, ::Vararg{Any}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Quiqbox ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:617
 [3] genBasisFunc(::Vector{Any}, ::GaussFunc{Float64, typeof(itself), typeof(itself)}, ::Tuple{Int64, Int64, Int64})
   @ Quiqbox ~/.julia/packages/Quiqbox/0Acru/src/Basis.jl:617
 [4] top-level scope
   @ ~/Documents/Coding/Quiqbox.jl/examples/github.jl:28
in expression starting at /Users/jamesbrown/Documents/Coding/Quiqbox.jl/examples/github.jl:28

[frankwswang]

Have you updated Quiqbox to version 0.5.6? Mixing the use of Number and ParamBox in arguments is a new syntax sugar added in 0.5.6.

[JamesB-1qbit]

I updated and it now works. Thanks!