Fix the basis optimization file

Author: Leticia-maria

src/basis-forge.jl

@@ -1,6 +1,6 @@
 ###############################################################################
 #  basis_forge.jl                                                             #
-#  ---------                                                                  #
+#  ------------------------------------------------------------------------   #
 #  Build custom Gaussian‑format basis sets via automatic least‑squares fits   #
 #  to Slater targets.                                                         #
 #                                                                             #
@@ -9,63 +9,61 @@
 ###############################################################################
 
 using LinearAlgebra
-using QuadGK                     # ∫₀^∞ … r² dr integrals
-using Optim                      # derivative‑free optimisation
+using QuadGK          # ∫₀^∞ … r² dr integrals
+using Optim           # derivative‑free optimisation
 
 # ---------- numerical settings ----------------------------------------------
-const ATOL = 1e-12               # quadrature accuracy
+const ATOL = 1e-12    # quadrature accuracy
 
 # ---------- Slater target ----------------------------------------------------
 χ_STO(r, ζ; n, l) = r^(n - 1) * exp(-ζ * r) * r^l
-
 sto_norm(ζ; n, l) = sqrt((2ζ)^(2n + 2l + 1) / factorial(2n + 2l))
 
 # ---------- primitive Gaussian ----------------------------------------------
 χ_G(r, α; l) = r^l * exp(-α * r^2)
-
 radial_integral(f) = QuadGK.quadgk(r -> f(r) * r^2, 0.0, Inf; atol = ATOL)[1]
 
 # ---------- build A, B, C matrices ------------------------------------------
 function build_matrices(alphas, ζ; n, l)
     k      = length(alphas)
-    Nsto   = sto_norm(ζ; n, l)
-    sto(r) = Nsto * χ_STO(r, ζ; n, l)
+    Nsto   = sto_norm(ζ; n = n, l = l)
+    sto(r) = Nsto * χ_STO(r, ζ; n = n, l = l)
 
-    A = 1.0                                             # <STO|STO>
-    B = [ radial_integral(r -> χ_G(r, α; l) * sto(r)) for α in alphas ]
-    C = [ radial_integral(r -> χ_G(r, αi; l) * χ_G(r, αj; l))
-          for αi in alphas, αj in alphas ]
+    A = 1.0                                             # ⟨STO|STO⟩
+    B = [radial_integral(r -> χ_G(r, α; l) * sto(r)) for α in alphas]
+    C = [radial_integral(r -> χ_G(r, αi; l) * χ_G(r, αj; l))
+         for αi in alphas, αj in alphas]
     return A, B, Symmetric(reshape(C, k, k))
 end
 
 # ---------- inner LSQ for coefficients + RMS --------------------------------
 function inner_fit(alphas, ζ; n, l)
-    A, B, C = build_matrices(alphas, ζ; n, l)
+    A, B, C = build_matrices(alphas, ζ; n = n, l = l)
     coeffs  = C \ B
     coeffs ./= sqrt(dot(coeffs, C * coeffs))            # renormalise
-    RMS²    = A - 2dot(coeffs, B) + dot(coeffs, C * coeffs)
-    return coeffs, sqrt(RMS²)
+    rms²    = A - 2dot(coeffs, B) + dot(coeffs, C * coeffs)
+    return coeffs, sqrt(rms²)
 end
 
 # ---------- outer optimisation of β, γ, ζ  (even‑tempered αᵢ) ---------------
-function optimise_basis(k; n, l; β0 = 0.15, γ0 = 3.2, ζ0 = 1.0)
+function optimise_basis(k, n, l, β0 = 0.15, γ0 = 3.2, ζ0 = 1.0)
     obj(p) = begin
         β, γ, ζ = p
-        αs = [β * γ^(i - 1) for i = 1:k]
-        _, rms = inner_fit(αs, ζ; n, l)
+        alphas  = [β * γ^(i - 1) for i = 1:k]
+        _, rms  = inner_fit(alphas, ζ; n = n, l = l)
         rms
     end
 
-    p0     = [β0, γ0, ζ0]
-    lower  = [1e-4, 1.01, 0.1]     # β>0, γ>1, ζ>0
-    upper  = [1e2,  10.0, 5.0]
+    p0    = [β0, γ0, ζ0]
+    lower = [1e-4, 1.01, 0.1]    # β>0, γ>1, ζ>0
+    upper = [1e2,  10.0, 5.0]
 
-    res = optimize(obj, lower, upper, p0,
-                   Fminbox{NelderMead}(); x_tol = 1e-10, f_tol = 1e-12)
+    opts  = Optim.Options(x_tol = 1e-10, f_tol = 1e-12)
+    res   = optimize(obj, lower, upper, p0, Fminbox(NelderMead()), opts)
 
-    β, γ, ζ = Optim.minimizer(res)
-    αs      = [β * γ^(i - 1) for i = 1:k]
-    coeffs, rms = inner_fit(αs, ζ; n, l)
+    β, γ, ζ   = Optim.minimizer(res)
+    alphas    = [β * γ^(i - 1) for i = 1:k]
+    coeffs, rms = inner_fit(alphas, ζ; n = n, l = l)
     return β, γ, ζ, coeffs, rms
 end
 
@@ -76,7 +74,7 @@ function print_gaussian94(io, element, n, l, β, γ, coeffs)
     println(io, " $element   $lchar   $k")
     for (i, c) in enumerate(coeffs)
         α = β * γ^(i - 1)
-        @printf(io, " %14.8f  %14.8f\n", α, c)
+        print(io, " $α $c\n")
     end
 end
 
@@ -88,22 +86,22 @@ shell_table = Dict(
     "H" => [(1, 0, 3)],                          # 1s STO‑3G
     "C" => [(1, 0, 3),                           # 1s
             (2, 0, 3),                           # 2s
-            (2, 1, 3),                           # 2p (shares β,γ with 2s manually)
-            (0, 2, 1)]                           # 3d polarisation (single prim.)
+            (2, 1, 3),                           # 2p (shares β,γ with 2s)
+            (0, 2, 1)]                           # 3d polarisation
 )
 
 open("my-basis.gbs", "w") do io
     for el in elements
         for (idx, (n, l, k)) in enumerate(shell_table[el])
             # share β, γ between 2s and 2p (SP trick) if consecutive in list
             if el == "C" && n == 2 && l == 1 && idx > 1
-                β, γ, _, _, _ = optimise_basis(k; n = 2, l = 0)
-                ζ     = 1.0                           # ζ irrelevant for pure Gaussian fit
+                β, γ, _, _, _ = optimise_basis(k, 2, 0)
+                ζ = 1.0                          # ζ irrelevant for pure Gaussian fit
                 coeffs, _ = inner_fit([β * γ^(i - 1) for i = 1:k],
                                        ζ; n = 2, l = 1)
                 print_gaussian94(io, el, n, l, β, γ, coeffs)
             else
-                β, γ, ζ, coeffs, _ = optimise_basis(k; n, l)
+                β, γ, ζ, coeffs, _ = optimise_basis(k, n, l)
                 print_gaussian94(io, el, n, l, β, γ, coeffs)
             end
         end

src/my-basis.gbs

@@ -0,0 +1,18 @@
+ H   s   3
+ 0.24965746218505253 0.4768772962297831
+ 0.9884344390703059 1.138167295631296
+ 3.913372473585632 0.9768881958678963
+ C   s   3
+ 0.24965746218505253 0.4768772962297831
+ 0.9884344390703059 1.138167295631296
+ 3.913372473585632 0.9768881958678963
+ C   s   3
+ 0.15 0.3508367526991321
+ 0.48 0.8230658425364937
+ 1.5360000000000003 -0.45283881485029753
+ C   p   3
+ 0.15 0.31802869255359717
+ 0.48 -0.3594329396282634
+ 1.5360000000000003 0.4658671835361496
+ C   d   1
+ 0.2536263833250992 0.23652096181090168

src/sto3g-test-2.jl

@@ -16,40 +16,40 @@ import Optim: converged
 const atol = 1e-12   # absolute accuracy for quadgk
 
 # ----------  STO target ------------------------------------------------------
-χ_STO(r, ζ; n=1, l=0) = r^(n-1) * exp(-ζ*r) * r^l      # radial × angular factor
+χ_STO(r, ζ, n=1, l=0) = r^(n-1) * exp(-ζ*r) * r^l      # radial × angular factor
 
 # analytic normalisation  N = √[(2ζ)^{2n+2l+1}/(2n+2l)!]
-function sto_norm(ζ; n, l)
+function sto_norm(ζ, n, l)
     return sqrt( (2ζ)^(2n+2l+1) / factorial(2n+2l) )
 end
 
 # ----------  primitive Gaussian ---------------------------------------------
-χ_G(r, α; l=0) = r^l * exp(-α * r^2)
+χ_G(r, α, l=0) = r^l * exp(-α * r^2)
 
 # ----------  helper: radial integral ----------------------------------------
 radial_integral(f) = QuadGK.quadgk(r -> f(r) * r^2, 0.0, Inf; atol)[1]
 
 # ----------  build matrices for given αs and ζ ------------------------------
-function build_matrices(alphas, ζ; n=1, l=0)
+function build_matrices(alphas, ζ, n=1, l=0)
     k = length(alphas)
-    Nsto = sto_norm(ζ; n, l)
-    sto  = r -> Nsto * χ_STO(r, ζ; n, l)
+    Nsto = sto_norm(ζ, n, l)
+    sto  = r -> Nsto * χ_STO(r, ζ, n, l)
 
     # A = <STO|STO> = 1 (normalised), but keep it explicit for completeness
     A = 1.0
 
     # Bᵢ = <PGFᵢ|STO>
-    B = [ radial_integral(r -> χ_G(r, α; l) * sto(r)) for α in alphas ]
+    B = [ radial_integral(r -> χ_G(r, α, l) * sto(r)) for α in alphas ]
 
     # Cᵢⱼ = <PGFᵢ|PGFⱼ>
-    C = [ radial_integral(r -> χ_G(r, αi; l) * χ_G(r, αj; l))
+    C = [ radial_integral(r -> χ_G(r, αi, l) * χ_G(r, αj, l))
           for αi in alphas, αj in alphas ]
     return A, B, Symmetric(reshape(C, k, k))
 end
 
 # ----------  inner LSQ: coeffs + RMS ----------------------------------------
-function inner_fit(alphas, ζ; n=1, l=0)
-    A, B, C = build_matrices(alphas, ζ; n, l)
+function inner_fit(alphas, ζ, n=1, l=0)
+    A, B, C = build_matrices(alphas, ζ, n, l)
     coeffs  = C \ B                          # normal equations
     coeffs ./= sqrt(dot(coeffs, C * coeffs)) # renormalise contracted AO
     RMS²    = A - 2dot(coeffs,B) + dot(coeffs, C * coeffs)
@@ -62,10 +62,10 @@ end
 
 `p = [β, γ, ζ]`  → returns RMS radial error for k primitives.
 """
-function rms_objective(p, k; n=1, l=0)
+function rms_objective(p, k, n=1, l=0)
     β, γ, ζ = p
     αs = [ β * γ^(i-1) for i = 1:k ]
-    _, rms = inner_fit(αs, ζ; n, l)
+    _, rms = inner_fit(αs, ζ, n, l)
     return rms
 end
 
@@ -75,32 +75,29 @@ end
 
 Returns (β★, γ★, ζ★, coeffs, RMS★).
 """
-function optimise_basis(k; n=1, l=0, β0=0.15, γ0=3.2, ζ0=1.0)
-    # initial parameter vector
-    p0   = [β0, γ0, ζ0]
-
-    # parameter bounds  (β>0, γ>1, ζ>0) – feel free to widen
-    lower = [1e-4, 1.01, 0.1]
-    upper = [1e2 , 10.0,  5.0]
-
-    obj(p) = rms_objective(p, k; n, l)
+function optimise_basis(k, n, l, β0 = 0.15, γ0 = 3.2, ζ0 = 1.0)
+    obj(p) = begin
+        β, γ, ζ = p
+        αs = [β * γ^(i - 1) for i = 1:k]
+        _, rms = inner_fit(αs, ζ, n, l)
+        rms
+    end
 
-    res = optimize(obj, lower, upper, p0,
-                   Fminbox{NelderMead}();
-                   x_tol = 1e-10, f_tol = 1e-12)
+    p0    = [β0, γ0, ζ0]
+    lower = [1e-4, 1.01, 0.1]   # β>0, γ>1, ζ>0
+    upper = [1e2,  10.0, 5.0]
 
-    if !converged(res)
-        error("Optimisation failed!")
-    end
+    opts  = Optim.Options(x_tol = 1e-10, f_tol = 1e-12)
+    res   = optimize(obj, lower, upper, p0, Fminbox(NelderMead()), opts)
 
-    β★, γ★, ζ★ = Optim.minimizer(res)
-    αs  = [ β★ * γ★^(i-1) for i = 1:k ]
-    coeffs, RMS★ = inner_fit(αs, ζ★; n, l)
-    return β★, γ★, ζ★, coeffs, RMS★, res
+    β, γ, ζ = Optim.minimizer(res)
+    αs      = [β * γ^(i - 1) for i = 1:k]
+    coeffs, rms = inner_fit(αs, ζ, n, l)
+    return β, γ, ζ, coeffs, rms
 end
 
 # ----------  example: H 1s STO‑3G (k=3, n=1, l=0) ---------------------------
-β, γ, ζ, c, rms, res = optimise_basis(3; n=1, l=0)
+β, γ, ζ, c, rms, res = optimise_basis(3, 1, 0)
 
 println("=== STO‑3G automatic fit (H 1s) ===")
 println("β       = ", round(β, 8))