can u see this i did corrections if you just help to remove some errors
! =========================================================
! A Minimal SCF Example for Helium (M = 2 basis functions)
! Program name: namrah
! =========================================================
! - Reads 2 exponents alpha(1..2) from “basis_set.dat”
! - Builds S (overlap), H (one-electron integrals),
! G (two-electron integrals).
! - Performs a basic SCF iteration with Löwdin orthogonalization
! and a Jacobi diagonalization routine.
! - Checks final idempotency and commutator conditions.
!
! NOTE: The integral formulas inside are toy placeholders!
! Replace them with correct expressions for your basis.
!
program namrah
implicit none
! ------------------------------------
! 1) Global Declarations and Constants
! ------------------------------------
integer, parameter :: dp = selected_real_kind(15)
integer, parameter :: M = 2 ! Number of basis functions
real(dp), parameter :: Z = 2.0_dp ! Nuclear charge for Helium
real(dp), parameter :: threshold = 1.0e-8_dp
integer, parameter :: max_iter = 50
! Arrays for exponents, etc.
real(dp), allocatable :: alpha(
! basis exponents
real(dp), allocatable :: S(:,:), H(:,:), G(:,:,:, ! Overlap, One-e, Two-e integrals
real(dp), allocatable :: c(:), D(:, ! MO coeffs (vector) and density
! For diagonalizing S and F’:
real(dp), allocatable :: T(:,:), eS(:), S_bar(:,:), S_bar_inv_sqrt(:,
real(dp), allocatable :: S_inv_half(:,
! Temporary storage for Fock:
real(dp), allocatable :: F(:,:), Fprime(:,:), Vtmp(:,:), eF(
! SCF energies, iteration counters
real(dp) :: energy, prev_energy, diff
integer :: i, j, k, l, iter, io_status
! For post-SCF analysis:
real(dp) :: mo_energy_1, J11, energy_SC
! ------------------------------------
! 2) Allocate Arrays
! ------------------------------------
allocate(alpha(M), S(M,M), H(M,M), G(M,M,M,M), &
c(M), D(M,M), T(M,M), eS(M), &
S_bar(M,M), S_bar_inv_sqrt(M,M), &
S_inv_half(M,M), F(M,M), Fprime(M,M), &
Vtmp(M,M), eF(M))
! ------------------------------------
! 3) Read Basis Exponents alpha
! ------------------------------------
open(unit=10, file=‘basis_set.dat’, status=‘old’, iostat=io_status)
if (io_status /= 0) then
print , “
Error opening ‘basis_set.dat’. Check file path.”
stop
end if
do i=1,M
read(10,) alpha(i)
end do
close(10)
print *, ">>> Read exponents alpha: ", alpha
! ------------------------------------
! 4) Build Overlap Matrix S
! (Replace with correct formula!)
! ------------------------------------
S = 0.0_dp
do i=1,M
do j=1,M
if (i == j) then
! Diagonal elements (very rough placeholder)
S(i,j) = 1.0_dp
else
! Example placeholder formula
S(i,j) = (2.0_dp*sqrt(alpha(i)*alpha(j)) / (alpha(i)+alpha(j)))**3
end if
end do
end do
print *, “>>> Overlap matrix S:”
do i=1,M
print ‘(2F14.6)’, S(i,
end do
! ------------------------------------
! 5) Diagonalize S → get S_bar (diag) and T (evecs)
! We’ll store S eigenvalues in eS(:),
! then build diagonal matrix S_bar from them.
! ------------------------------------
call jacobi_diagonalize(S, T, eS, M)
S_bar = 0.0_dp
do i=1,M
S_bar(i,i) = eS(i)
end do
! Fix any nearly zero eigenvalues
do i=1,M
if (S_bar(i,i) < 1.0e-12_dp) then
print *, “>> WARNING: Overlap eigenvalue too small at i=”, i, " => forcing to 1.0e-6"
S_bar(i,i) = 1.0e-6_dp
end if
end do
! Compute S_bar^(-1/2)
S_bar_inv_sqrt = 0.0_dp
do i=1,M
S_bar_inv_sqrt(i,i) = 1.0_dp / sqrt(S_bar(i,i))
end do
! Build S_inv_half = T * S_bar^(-1/2) * T^T (Löwdin S^(-1/2))
S_inv_half = matmul( matmul(T, S_bar_inv_sqrt), transpose(T) )
! ------------------------------------
! 6) Build the One-Electron Integral Matrix H
! (Replace with correct formula!)
! ------------------------------------
H = 0.0_dp
do i=1,M
do j=1,M
! Example placeholder formula:
! e.g. kinetic + nuclear attraction in minimal form
H(i,j) = 4.0_dp * (sqrt(alpha(i)*alpha(j)) / (alpha(i) + alpha(j)))**3 * &
( alpha(i)alpha(j) - Z(alpha(i) + alpha(j)) )
end do
end do
print *, “>>> One-electron integrals H:”
do i=1,M
print ‘(2F14.6)’, H(i,
end do
! ------------------------------------
! 7) Build Two-Electron Integrals G(p,q,r,s)
! (Replace with correct formula!)
! ------------------------------------
G = 0.0_dp
do i=1,M
do j=1,M
do k=1,M
do l=1,M
! Example placeholder formula:
G(i,j,k,l) = 32.0_dp * (alpha(i)*alpha(j)*alpha(k)*alpha(l))**1.5_dp * &
( 1.0_dp / ((alpha(i)+alpha(j))**3 * (alpha(k)+alpha(l))**2 ) - &
1.0_dp / ((alpha(i)+alpha(j))**3 * (alpha(i)+alpha(j)+alpha(k)+alpha(l))**2 ) - &
1.0_dp / ((alpha(i)+alpha(j))**2 * (alpha(i)+alpha(j)+alpha(k)+alpha(l))**3 ) )
end do
end do
end do
end do
print *, “>>> Selected two-electron integrals G:”
print *, "G(1,1,1,1) = ", G(1,1,1,1)
print *, "G(1,1,2,2) = ", G(1,1,2,2)
print *, "G(1,2,1,2) = ", G(1,2,1,2)
print *, "G(2,2,2,2) = ", G(2,2,2,2)
print *, "G(1,1,1,2) = ", G(1,1,1,2)
print *, "G(1,2,2,2) = ", G(1,2,2,2)
print *, “------------------------------------------------------------”
print *, “>>> Overlap Matrix S (compare with Snow & Bills)”
do i=1,M
print ‘(2F14.6)’, S(i,
end do
print *, “>>> One-electron Integrals H (compare with Snow & Bills)”
do i=1,M
print ‘(2F14.6)’, H(i,
end do
print *, “>>> Two-electron Integrals G (compare with Snow & Bills)”
print *, "G(1,1,1,1) = ", G(1,1,1,1)
print *, "G(1,1,2,2) = ", G(1,1,2,2)
print *, "G(1,2,1,2) = ", G(1,2,1,2)
print *, "G(2,2,2,2) = ", G(2,2,2,2)
! ------------------------------------
! 8) SCF Iteration
! - c(1)=1, c(2)=0 initial guess
! ------------------------------------
c = 0.0_dp
c(1) = 1.0_dp
prev_energy = 0.0_dp
energy = 1.0_dp ! Dummy initial
print *
print *, " Iter | c1 | c2 | F11 | F22 | F12 | eps1 | E"
print *, “----------------------------------------------------------------------------------------------”
print *, “>>> Initial MO Coefficients c:”
print ‘(2F14.6)’, c
do iter=1,max_iter
! (a) Build Fock Matrix F = H + (two-electron part).
! D = c c^T => D(i,j) = c(i)*c(j).
! F(i,j) = H(i,j) + sum_{k,l} D(k,l)*G(i,j,k,l).
F = H
do i=1,M
do j=1,M
do k=1,M
do l=1,M
F(i,j) = F(i,j) + c(k)*c(l)*G(i,j,k,l)
end do
end do
end do
end do
! (b) Löwdin orthogonalization: F' = S^(-1/2) * F * S^(-1/2)
Fprime = matmul( S_inv_half, matmul(F, S_inv_half) )
! (c) Diagonalize F' => eigenvectors Vtmp, eigenvalues eF
call jacobi_diagonalize(Fprime, Vtmp, eF, M)
! (d) The lowest eigenvalue => occupied orbital
! c = S^(-1/2) * eigenvector_of_lowest_eigenvalue
c = matmul(S_inv_half, Vtmp(:,1))
! (e) Normalize c
call normalize_vector(c)
! (f) Compute new density D = c c^T
D = outer_product(c, c)
! (g) Compute SCF energy:
! E = 2 * Tr(D H) + Tr(D D G) in minimal basis
energy = 0.0_dp
do i=1,M
do j=1,M
energy = energy + 2.0_dp * H(i,j)*D(i,j)
! Two-electron part:
do k=1,M
do l=1,M
energy = energy + D(i,j)*D(k,l)*G(i,j,k,l)
end do
end do
end do
end do
print ‘(i3,2x,f10.6,2x,f10.6,2x,f10.6,2x,f10.6,2x,f10.6,2x,f10.6,2x,f10.6)’, &
iter, c(1), c(2), F(1,1), F(2,2), F(1,2), eF(1), energy
! (h) Convergence check
if (abs(energy - prev_energy) < threshold) exit
prev_energy = energy
end do
print *, “------------------------------------------------------------”
print , "Check from formula E = 2epsilon(1) - J(1,1):"
print *, " >> SCF Energy =“, energy
print , " >> 2epsilon(1) =”, 2.0_dp * mo_energy_1
print *, " >> J(1,1) =“, J11
print *, " >> => E (expected) =”, energy_SC
print *, " >> Difference =", abs(energy - energy_SC)
real(dp) :: E_exp = -2.904_dp
print *, “------------------------------------------------------------”
print *, "SCF Total Energy = “, energy, " a.u.”
print *, "Experimental Energy = “, E_exp, " a.u.”
print *, "Difference = “, abs(energy - E_exp), " a.u.”
! ------------------------------------
! 9) Post-SCF Checks
! ------------------------------------
! (A) For He in minimal basis, one might see a formula:
! E ~ 2epsilon(1) - J(1,1) => mo_energy_1 = eF(1), J11=G(1,1,1,1)
mo_energy_1 = eF(1)
J11 = G(1,1,1,1)
energy_SC = 2.0_dpmo_energy_1 - J11
print , “------------------------------------------------------------”
print , "Check from simple formula E = 2epsilon(1) - J(1,1):"
print , " >> 2epsilon(1) =", 2.0_dpmo_energy_1
print *, " >> J(1,1) =“, J11
print *, " >> => E =”, energy_SC
print *, “Difference from SCF Energy =”, abs(energy - energy_SC)
real(dp) :: E_He, E_HePlus, Ionization_energy, Ionization_Koopman
E_He = energy
E_HePlus = -2.0_dp ! Approximate He+ energy from literature
Ionization_energy = E_HePlus - E_He
Ionization_Koopman = -eF(1)
print *, “------------------------------------------------------------”
print *, “Computed First Ionization Energy =”, Ionization_energy, " a.u."
print *, “Experimental Ionization Energy = 0.904 a.u.”
print *, “Ionization Energy using Koopman’s theorem =”, Ionization_Koopman, " a.u."
print *, “Difference from Experimental =”, abs(Ionization_energy - 0.904_dp), " a.u."
print *, “------------------------------------------------------------”
print *, “Testing Alternative Exponents for Basis Set Optimization”
alpha(1) = 1.0_dp
alpha(2) = 3.0_dp
call run_scf(alpha)
print *, “Energy with alpha1=1.0, alpha2=3.0 =”, energy
alpha(1) = 2.0_dp
alpha(2) = 2.5_dp
call run_scf(alpha)
print *, “Energy with alpha1=2.0, alpha2=2.5 =”, energy
print *, “------------------------------------------------------------”
print *, “Conclusion: Any other pair of exponents increases the energy.”
! (B) Idempotency check: D S D ~ D
call check_idempotency(D, S)
! (C) Commutator check: F D S == S D F
call check_commute(F, D, S)
contains
! -------------------------------------------------------
! 1. Jacobi Diagonalization for Symmetric Matrix
! A is overwritten with its diagonal form
! V will store the eigenvectors, evals( the eigenvalues
! -------------------------------------------------------
subroutine jacobi_diagonalize(A, V, evals, n)
implicit none
integer, intent(in) :: n
real(dp), intent(inout):: A(n,n)
real(dp), intent(out) :: V(n,n), evals(n)
integer :: p, q, i, j, iter, max_jacobi
real(dp) :: off_diag, tol, theta, cRot, sRot, tau, tmp
real(dp), parameter :: eps = 1.0e-12_dp
! Initialize V = Identity
V = 0.0_dp
do i=1,n
V(i,i) = 1.0_dp
end do
max_jacobi = 100
tol = 1.0e-12_dp
do iter=1,max_jacobi
! 1) Find largest off-diagonal element
off_diag = 0.0_dp
p = 1; q = 2
do i=1,n
do j=i+1,n
if (abs(A(i,j)) > off_diag) then
off_diag = abs(A(i,j))
p = i
q = j
end if
end do
end do
if (off_diag < tol) exit
! 2) Compute rotation angle
theta = 0.5_dp * atan2(2.0_dp*A(p,q), A(q,q) - A(p,p))
cRot = cos(theta)
sRot = sin(theta)
tau = sRot/(1.0_dp + cRot)
! 3) Update matrix A
tmp = A(p,p)
A(p,p) = tmp - A(p,q)*sRot*2.0_dp
A(q,q) = A(q,q) + A(p,q)*sRot*2.0_dp
A(p,q) = 0.0_dp
A(q,p) = 0.0_dp
! Update off-diagonal elements
do i=1,n
if ( (i /= p) .and. (i /= q) ) then
tmp = A(i,p)
A(i,p) = tmp - sRot*(A(i,q)) &
& + tau*(A(p,p) - tmp)
A(p,i) = A(i,p)
tmp = A(i,q)
A(i,q) = tmp + sRot*(A(i,p)) &
& + tau*(A(q,q) - tmp)
A(q,i) = A(i,q)
end if
end do
! 4) Update eigenvector matrix V
do i=1,n
tmp = V(i,p)
V(i,p) = cRot*tmp - sRot*V(i,q)
V(i,q) = sRot*tmp + cRot*V(i,q)
end do
end do
! Save diagonal => evals
do i=1,n
evals(i) = A(i,i)
end do
end subroutine jacobi_diagonalize
! -------------------------------------------------------
! 2. Normalize a vector
! -------------------------------------------------------
subroutine normalize_vector(vec)
implicit none
real(dp), intent(inout) :: vec(
real(dp) :: norm_val
norm_val = sqrt(sum(vec*vec))
if (norm_val > 1.0e-14_dp) vec = vec / norm_val
end subroutine normalize_vector
! -------------------------------------------------------
! 3. Outer product of two vectors => matrix
! -------------------------------------------------------
function outer_product(a, b) result(res)
implicit none
real(dp), intent(in) :: a(:), b(
real(dp) :: res(size(a), size(b))
integer :: i, j
do i=1,size(a)
do j=1,size(b)
res(i,j) = a(i)*b(j)
end do
end do
end function outer_product
! -------------------------------------------------------
! 4. Check idempotency: R S R ~ R
! -------------------------------------------------------
subroutine check_idempotency(R, S)
implicit none
real(dp), intent(in) :: R(:,:), S(:,
real(dp) :: RSR(size(R,1), size(R,2))
integer :: i,j,n
n = size(R,1)
RSR = matmul(R, matmul(S, R))
do i=1,n
do j=1,n
if (abs(RSR(i,j) - R(i,j)) > 1.0e-6_dp) then
print *, "❌ RSR != R => Not idempotent at (",i,",",j,")"
return
end if
end do
end do
print *, "✅ R is idempotent within tolerance."
end subroutine check_idempotency
! -------------------------------------------------------
! 5. Check Commutator: F R S == S R F
! -------------------------------------------------------
subroutine check_commute(F, R, S)
implicit none
real(dp), intent(in) :: F(:,:), R(:,:), S(:,
real(dp) :: FRS(size(F,1), size(F,2))
real(dp) :: SRF(size(F,1), size(F,2))
integer :: i,j,n
n = size(F,1)
FRS = matmul(F, matmul(R, S))
SRF = matmul(S, matmul(R, F))
do i=1,n
do j=1,n
if (abs(FRS(i,j) - SRF(i,j)) > 1.0e-6_dp) then
print *, "❌ FRS != SRF => SCF might not be fully converged at (",i,",",j,")"
return
end if
end do
end do
print *, "✅ F R S = S R F (commute) within tolerance."
end subroutine check_commute
print *, “------------------------------------------------------------”
print *, “Effect of the Number of Basis Functions (M > 2)”
print *, “Computing SCF energy for M = 1 to 10 basis functions”
print *, “------------------------------------------------------------”
integer :: M_max
real(dp), allocatable :: alpha_var(:), eF_var(:), energy_var(
M_max = 10 ! Maximum basis size
allocate(alpha_var(M_max), eF_var(M_max), energy_var(M_max))
do M = 1, M_max
print *, “Running SCF for M =”, M
! Generate basis set exponents
alpha_var(1) = 1.0_dp
do i = 2, M
alpha_var(i) = 1.2_dp * alpha_var(i-1)
end do
! Reallocate arrays for the new basis size
deallocate(S, H, G, c, D, T, eS, S_bar, S_bar_inv_sqrt, S_inv_half, F, Fprime, Vtmp, eF)
allocate(S(M,M), H(M,M), G(M,M,M,M), c(M), D(M,M), T(M,M), eS(M), &
S_bar(M,M), S_bar_inv_sqrt(M,M), S_inv_half(M,M), F(M,M), Fprime(M,M), &
Vtmp(M,M), eF(M))
! Recompute S, H, G matrices with new basis set
call build_matrices(M, alpha_var, S, H, G)
! Run SCF procedure
call run_scf(M, S, H, G, c, D, T, eS, S_bar, S_bar_inv_sqrt, S_inv_half, F, Fprime, Vtmp, eF, energy)
! Store results
eF_var(M) = eF(1) ! Store lowest orbital energy
energy_var(M) = energy ! Store total energy
end do
! Print results for all basis sizes
print *, “------------------------------------------------------------”
print *, “Basis Set Convergence Results”
print *, “M | Orbital Energy (ϵ1) | Total SCF Energy (E)”
print *, “------------------------------------------------------------”
do M = 1, M_max
print ‘(I3,2X,F15.8,2X,F15.8)’, M, eF_var(M), energy_var(M)
end do
! Compare with the HF limit (-2.8617 a.u.) and double-zeta basis set
real(dp) :: E_HF_limit = -2.8617_dp
print *, “------------------------------------------------------------”
print *, "Hartree-Fock Limit Energy = “, E_HF_limit, " a.u.”
print *, "Best SCF Energy Obtained = “, energy_var(M_max), " a.u.”
print *, "Difference from HF Limit = “, abs(E_HF_limit - energy_var(M_max)), " a.u.”
print *, “------------------------------------------------------------”
print *, “Conclusion: A sufficiently large basis set is needed to reach the HF limit.”
print *, “------------------------------------------------------------”
print *, “------------------------------------------------------------”
print *, “All done!”
end program namrah