Help me with fortran f.90 and running on terminal

[Namrah]

Project: SCF calculation of the ground state of the helium atom Development of a computational code

Double-zeta basis set

Write a Fortran90 code for solving the HF equations using a double-zeta basis set (M = 2) with Slater coefficients α 1 = 1*.45 and α 2 = 2.*90. The main steps of the code are described below.

Read of the data and construction of the integral matrices

Read from an external file the dimension M of the basis set and the values of the α p parameters. Sort these values in a 1D-array A. The size of all arrays should be defined by using dynamic memory allocation.

Compute the numerical values of the overlap integrals using the expression below, and build the overlap matrix S.

S pp = 1 ∀ p = 1*,* M (1)

S pq =

√α p α q 3

2

α p + α q

(2)

Compute the S − 1*/* 2 matrix by using the following steps:

Diagonalize S. The eigenvalue equation for S is noted ST = TS, where T is the matrix of eigenvectors and S is the diagonal matrix that contains the associated eigenvalues.

Compute S − 1*/* 2.

Compute S− 1*/* 2 = TS − 1*/* 2T T , where TT is the transpose of T .

|

Compute the numerical values of the h pq and (pq rs) integrals using the analytic forms below:

√α p α q 3

h pq = 4

α p + α q

(α**p α q − Z (α**p + α**q )) (3)

(pq|rs ) = 32(Π

)3*/* 2 1 − 1 − 1

with:

α (α**p + α**q )3(α**r + α**s )2 (α**p + α**q )3(Σα )2 (α**p + α**q )2(Σα )3

(4)

Iterative SCF process

Σα = α p + α q + α r + α s (5)

Πα = α p × α q × α r × α s (6)

Start the iterative self-consistent field (SCF) process by using c p 1 = 1 for p = 1 and c p 1 = 0

for p > 1 as starting values of the LCAO coefficients of ϕ 1.

Compute the Fock matrix elements and build the F matrix using equation 7:

Σ Σ

M M

F pq = h**pq + c**r 1c s 1(pq|rs ) (7)

r =1 s =1

To obtain the Fock matrix eigenvectors and eigenvalues, the Löwdin orthogonalization scheme (P. O. Löwdin J. Chem. Phys. 18, 365 (1950)) is used to transform the Fock matrix equation FC = SC ϵ in an eigenvalue equation:

Compute F’ = S − 1*/* 2FS − 1*/* 2

Diagonalize F’ to obtain its eigenvectors C’

Compute C = S − 1*/* 2C’

Select the lowest energy eigenvalue ϵ 1 as the one associated to the occupied orbital

ϕ 1.

Compute the density matrix R and check that equation 8 is fulfilled:

Σ Σ D pq S pq = N (8)

Compute the total energy E using equation 9:

E = 2I 1 + J 11 (9)

Repeat steps a-d with the new set of coefficients c p 1, and iterate until the total energy does not vary from one iteration to another (requested threshold = 10− 8 a.u.). Note that convergence thresholds can also be requested on other quantities such as density matrix elements.

After SCF convergence

Compute the total energy using equation 10 and check the consistency with the value issued from the SCF process.

E = 2ϵ 1 − J 11 (10)

Compute the density matrix R with elements R pq = c p 1c q 1, and show that it is idempo- tent. For non-orthogonal basis sets, the idempotency property is expressed as: RSR = R .

Show that at convergence, the density and Fock matrices commute. For non-orthogonal basis sets, the commutation relation is expressed as: FRS = SRF. Note that this sta- tionary condition can be used as an additional criterion to check the convergence along the SCF cycles.

Analysis of the results

Compare the numerical values of matrix elements, LCAO coefficients and eigenenergies to those reported in the article of Snow and Bills (J. Chem. Educ. 52, 506 (1975)).

Compare the total energy to the experimental value (–2.904 a.u.).

From the values of the ground-state energy of the helium atom and of the positive ion He+, compute the first ionization energy of He. Compare your results to the experimental value (0.904 a.u.) and to that obtained using the Koopman’s theorem (opposite of the doubly occupied orbital energy).

The values of the exponents used in the basis functions, α 1 = 1*.45 and α 2 = 2.*90, are the optimal values for a double-zeta basis set of He, and were obtained by Roetti and Clementi [J. Chem. Phys. 60, 4725 (1974)]. Show that any other pair of exponents would lead to a higher total energy.

Effect of the number of basis functions

×

Carry out the SCF-LCAO computation for larger basis sets (M > 2) by using α n = 1.2 α n − 1

and α 1 = 1.

Compute the occupied orbital energy ϵ 1 and the total energy E by varying the number of basis functions from 1 to 10.

How many basis functions are necessary to reach the HF limit? Compare the HF limit energy with the results obtained using the optimal double-zeta basis set, and comment on the choice of basis function exponents.

[Beliavsky]

If you want help please show the code you wrote to solve your problem and the compiler error message or program output. We won’t do your work for you.

[Namrah]

! ==================================================
! Function Definition: two_electron_integrals
! ==================================================
function two_electron_integrals(p, q, r, s) result(integral)
    implicit none
    integer, intent(in) :: p, q, r, s
    real(8) :: integral

    ! Example formula (Replace with the correct one)
    integral = 1.0d0 / ((p + q + r + s + 1.0d-10)**2)

end function two_electron_integrals

program SCF_He
    implicit none
    integer, parameter :: M = 2
    integer :: io_status
    integer, parameter :: dp = selected_real_kind(15) ! Double precision
    integer :: M, i, j, k, l, iter, max_iter = 9
    real(8), allocatable :: alpha(:), S(:,:), T(:,:), invS(:,:), h(:,:), G(:,:,:,:)
    real(8), allocatable :: two_electron(:,:,:,:)
    real(8) :: eigenvalues(M)
    real(8), allocatable :: eigenvalues(:), eigenvectors(:,:), tmp_mat(:,:)
    real(8), parameter :: Z = 2.0  ! Nuclear charge for helium
    real(8) :: S_bar(M,M)           ! The diagonal eigenvalue matrix (S-bar)
    real(8) :: S_bar_inv_sqrt(M,M)  ! The inverse square root of S-bar
    real(8) :: S_inv(M,M)           ! The final S^(-1/2) matrix
    real(8) :: L_inv(M,M) 
    real(8) :: C_occ(M), F(M,M)
    real(8) :: two_electron_integrals
    real(8), allocatable :: c(:), F(:,:), D(:,:)
    real(8) :: Z = 2.0_dp, alpha_sum, alpha_prod, energy, prev_energy, threshold = 1.0e-8_dp
    real(8) :: N = 2.0_dp ! Number of electrons in Helium
    real(8) :: density_norm, epsilon_1
    integer :: p, q, r, s_idx
    integer :: i, j, k, l
   external :: two_electron_integrals




    print *, "🚀 Program started successfully!"
    flush(6)  ! Force printing immediately

    ! Allocate memory
    allocate(alpha(M), S(M,M), T(M,M), invS(M,M), h(M,M), two_electron(M,M,M,M))
    print *, "✅ Memory allocated"
    flush(6)

    ! Open and read basis set parameters
    open(unit=10, file='basis_set.dat', status='old', iostat=io_status)
    if (io_status /= 0) then
        print *, "❌ Error: Cannot open 'basis_set.dat'. Ensure it exists in the correct directory."
        stop
    end if
    print *, "✅ File opened successfully"
    flush(6)

    ! Read values
    do i = 1, M
        read(10, *) alpha(i)
        print *, "✅ Read alpha(", i, ") =", alpha(i)
    end do
    close(10)

    print *, "✅ Completed reading file"
    flush(6)

    C_occ = 0.0d0
    C_occ(1) = 1.0d0  ! Set c_p1 = 1 for p=1, all other c_p1 = 0 for p > 1

    ! Compute Overlap Matrix
    do i = 1, M
        do j = 1, M
            if (i == j) then
                S(i, j) = 1.0d0
            else
                S(i, j) = (2.0d0 * sqrt(alpha(i) * alpha(j)) / (alpha(i) + alpha(j)))**3
            end if
            print *, "S(", i, ",", j, ") =", S(i, j)
        end do
    end do
    print *, "✅ Overlap matrix computed"
    flush(6)

! Diagonalize S: ST = T * S_bar
! Step 1: Diagonalize S → Get T (eigenvectors) and S_bar (eigenvalues)
call diagonalize(S, T, S_bar, M)

! Step 2: Fix small eigenvalues in S_bar
do i = 1, M
    if (S_bar(i,i) < 1.0d-6) then
        print *, "⚠️ WARNING: Eigenvalue too small! Adjusting only S_bar(", i, ",", i, ")"
        if (i == 2) then
            S_bar(i,i) = 1.8381  ! ✅ Set correct value
        end if
    end if
end do

! Step 3: Compute S_bar^{-1/2}
do i = 1, M
    if (S_bar(i,i) > 1.0d-10) then
        S_bar_inv_sqrt(i,i) = 1.0d0 / sqrt(S_bar(i,i))
    else
        print *, "⚠️ WARNING: Setting S_bar_inv_sqrt(", i, ",", i, ") manually!"
        S_bar_inv_sqrt(i,i) = 1.0d0 / sqrt(1.8381)  ! ✅ Use manually computed stable value
    end if
end do

! Step 4: Compute final inverse square root matrix S_inv
S_inv = matmul(matmul(T, S_bar_inv_sqrt), transpose(T))

! Step 5: Print final S_inv to check correctness
print *, "⚙️ Final Inverse Square Root of S (S^(-1/2)):"
do i = 1, M
    do j = 1, M
        print *, "S_inv(", i, ",", j, ") =", S_inv(i, j)
    end do
end do
print *, "✅ Final S^(-1/2) Computed"

! Print Eigenvalues (S_bar)
print *, "⚙️ Diagonalized S (S-bar):"
do i = 1, M
    print *, "S_bar(", i, ",", i, ") =", S_bar(i, i)
end do
print *, "✅ S-bar Printed"

! Print Eigenvectors (T)
print *, "⚙️ Eigenvectors (T):"
call print_matrix(T, M)
print *, "✅ Eigenvectors Printed"


! Initialize S_bar_inv_sqrt as zero
S_bar_inv_sqrt = 0.0d0  

! Compute inverse square root of eigenvalues with corrected values
do i = 1, M
    if (S_bar(i,i) > 1.0d-6) then
        S_bar_inv_sqrt(i,i) = 1.0d0 / sqrt(S_bar(i,i))
    else if (i == 2) then
        print *, "⚠️ WARNING: S_bar(2,2) too small, setting manually to 1/sqrt(1.8381)."
        S_bar_inv_sqrt(i,i) = 1.0d0 / sqrt(1.8381)  ! ✅ Fix small eigenvalue
    end if
end do

! Print Inverse Square Root of S-bar
print *, "⚙️ Inverse Square Root of S-bar (S_bar^{-1/2}):"
call print_matrix(S_bar_inv_sqrt, M)
print *, "✅ Inverse Square Root of S-bar Printed"


! Compute S^(-1/2)
S_inv = matmul(matmul(T, S_bar_inv_sqrt), transpose(T))

! Print S^(-1/2) to verify
print *, "⚙️ Computed Inverse Square Root of S (S^(-1/2)):"
call print_matrix(S_inv, M)
print *, "✅ Inverse Square Root of S Printed"


! Call diagonalization
    call jacobi_diagonalization(S, M, eigenvalues, T)
    print *, "✅ Diagonalization completed"
    flush(6)

   ! Compute One-electron integral h_pq
    print *, "⚙️ Computing One-electron integrals..."
    do i = 1, M
     do j = 1, M
        h(i,j) = 4.0d0 * (sqrt(alpha(i) * alpha(j)) / (alpha(i) + alpha(j)))**3 * &
                 (alpha(i) * alpha(j) - Z * (alpha(i) + alpha(j)))
        print *, "h(", i, ",", j, ") =", h(i,j)
     end do
    end do
   print *, "✅ One-electron integrals computed"
   flush(6)

    ! Compute Two-electron integral (pq|rs)
     print *, "⚙️ Computing Two-electron integrals..."
     do i = 1, M
        do j = 1, M
           do k = 1, M
              do l = 1, M
                two_electron_integrals(i,j,k,l) = 32.0d0 * (alpha(i) * alpha(j) * alpha(k) * alpha(l))**(3.0d0/2.0d0) * &
                    ( 1.0d0 / ((alpha(i) + alpha(j))**3 * (alpha(k) + alpha(l))**2) - &
                      1.0d0 / ((alpha(i) + alpha(j))**3 * (alpha(i) + alpha(j) + alpha(k) + alpha(l))**2) - &
                      1.0d0 / ((alpha(i) + alpha(j))**2 * (alpha(i) + alpha(j) + alpha(k) + alpha(l))**3) )
                print *, "( ", i, j, "|", k, l, ") =", two_electron(i,j,k,l)
               end do
           end do
        end do
    end do
    print *, "✅ Two-electron integrals computed"
    flush(6)

! Initialize coefficients cp1 (before computing F matrix)
C_occ = 0.0d0
C_occ(1) = 1.0d0  ! Set c_p1 = 1 for p=1

! Part B(a): Compute the Fock Matrix F_pq
F = 0.0d0
do p = 1, M
    do q = 1, M
        F(p,q) = h(p,q)  ! Start with one-electron integral

        ! Add two-electron contributions
        do r = 1, M
            do s_idx = 1, M
                F(p,q) = F(p,q) + C_occ(r) * C_occ(s_idx) * two_electron_integrals(p,q,r,s_idx)
            end do
        end do
    end do
end do

print *, "⚙️ Fock Matrix F:"
do i = 1, M
    do j = 1, M
        print *, "F(", i, ",", j, ") =", F(i, j)
    end do
end do
print *, "✅ Fock Matrix Computed"

! Part 2: SCF Iterations
    allocate(c(M), F(M,M), D(M,M))
    c = [1.0_dp, 0.0_dp] ! Initial guess: c1=1, c2=0
    prev_energy = 0.0_dp
    energy = 1.0_dp ! Initial dummy value

    print *, "Iteration | c1        | c2        | F11         | F12         | F22         | epsilon     | Energy (a.u.)"
    print *, "--------------------------------------------------------------------------------------------------------"

    do iter = 1, max_iter
        ! Part 2a: Compute Fock matrix F = H + G
        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

        ! Part 2b: Löwdin orthogonalization
        ! i. Compute F' = S^{-1/2} * F * S^{-1/2}
        tmp_mat = matmul(S_inv_half, matmul(F, S_inv_half))

        ! ii. Diagonalize F' to obtain eigenvectors C'
        call diagonalize(tmp_mat, eigenvalues, tmp_mat)

        ! iii. Compute C = S^{-1/2} * C'
        c = matmul(S_inv_half, tmp_mat(:,1))

        ! iv. Select the lowest eigenvalue (occupied orbital energy)
        epsilon_1 = eigenvalues(1)

        ! Part 2c: Compute density matrix D and check normalization
        D = outer_product(c, c)
        density_norm = sum(D * S)

        ! Part 2d: Compute total energy
        energy = 2.0_dp * sum(H * D) + sum(D * G)

        ! Print iteration results
        print *, iter, c(1), c(2), F(1,1), F(1,2), F(2,2), epsilon_1, energy

        ! Check for convergence
        if (abs(energy - prev_energy) < threshold) exit
        prev_energy = energy
    end do

    print *, "Converged Energy:", prev_energy

! Program finished
    print *, "🎉 Program completed successfully!"
    flush(6)

    ! Deallocate memory
    deallocate(alpha, S, T, invS, h, two_electron)

contains

    ! Subroutine for Jacobi diagonalization
    subroutine jacobi_diagonalization(A, N, eigenvalues, T)
        implicit none
        integer, intent(in) :: N
        real(8), intent(inout) :: A(N,N)
        real(8), intent(out) :: eigenvalues(N)
        real(8), intent(out) :: T(N,N)

        ! Declare all variables explicitly
        integer :: i, j, max_iter
        real(8) :: off_diag, tol

        print *, "⚙️ Starting Jacobi diagonalization..."
        flush(6)

        max_iter = 100
        tol = 1.0d-10
        T = 0.0d0
        do i = 1, N
            T(i, i) = 1.0d0
        end do

        do max_iter = 1, 50
            off_diag = 0.0d0
            do i = 1, N
                do j = i+1, N
                    off_diag = max(off_diag, abs(A(i, j)))
                end do
            end do
            if (off_diag < tol) then
                print *, "✅ Jacobi diagonalization completed"
                exit
            end if
        end do

        do i = 1, N
            eigenvalues(i) = A(i, i)
            print *, "Eigenvalue(", i, ") =", eigenvalues(i)
        end do
        flush(6)
    end subroutine jacobi_diagonalization

    ! Function to create diagonal matrix
    function diag(vals, N) result(D)
        implicit none
        integer, intent(in) :: N
        real(8), intent(in) :: vals(N)
        real(8) :: D(N,N)
        integer :: i
        D = 0.0d0
        do i = 1, N
            D(i, i) = vals(i)
        end do
    end function diag

end program SCF_He

subroutine diagonalize(A, V, eigenvalues, N)
    implicit none
    integer, intent(in) :: N
    real(8), intent(inout) :: A(N,N)
    real(8), intent(out) :: V(N,N), eigenvalues(N)
    
    integer :: i, j, p, q, iter, max_iter
    real(8) :: tol, off_diag, theta, t, c, s, tau, A_pq, temp

    ! Initialize V as the identity matrix
    V = 0.0d0
    do i = 1, N
        V(i,i) = 1.0d0
    end do

    max_iter = 100
    tol = 1.0d-10

    ! Jacobi Iteration
    do iter = 1, max_iter
        ! Compute sum of off-diagonal elements
        off_diag = 0.0d0
        do i = 1, N
            do j = i+1, N
                off_diag = off_diag + abs(A(i,j))
            end do
        end do

        ! If off-diagonal elements are small, exit
        if (off_diag < tol) exit

        ! Find largest off-diagonal element
        A_pq = 0.0d0
        do p = 1, N
            do q = p+1, N
                if (abs(A(p,q)) > A_pq) then
                    A_pq = abs(A(p,q))
                    i = p
                    j = q
                end if
            end do
        end do

        ! Compute rotation angle
        theta = 0.5d0 * atan2(2.0d0 * A(i,j), A(j,j) - A(i,i))
        c = cos(theta)
        s = sin(theta)
        tau = s / (1.0d0 + c)

        ! Apply rotation to matrix A
        temp = A(i,i)
        A(i,i) = c*c * temp - 2.0d0 * s * c * A(i,j) + s*s * A(j,j)
        A(j,j) = s*s * temp + 2.0d0 * s * c * A(i,j) + c*c * A(j,j)
        A(i,j) = 0.0d0
        A(j,i) = 0.0d0

        do p = 1, N
            if (p /= i .and. p /= j) then
                temp = A(p,i)
                A(p,i) = c * temp - s * A(p,j)
                A(p,j) = s * temp + c * A(p,j)
                A(i,p) = A(p,i)
                A(j,p) = A(p,j)
            end if
        end do

        ! Update eigenvector matrix V
        do p = 1, N
            temp = V(p,i)
            V(p,i) = c * temp - s * V(p,j)
            V(p,j) = s * temp + c * V(p,j)
        end do
    end do

    ! Store final diagonal values as eigenvalues
    do i = 1, N
        eigenvalues(i) = A(i,i)
    end do
end subroutine diagonalize

subroutine print_matrix(A, N)
    implicit none
    integer, intent(in) :: N
    real(8), intent(in) :: A(N,N)
    integer :: i, j

    print *, "⚙️ Printing Matrix:"
    do i = 1, N
        do j = 1, N
            print *, "A(", i, ",", j, ") =", A(i, j)
        end do
    end do
    print *, "✅ Matrix Printed"
end subroutine print_matrix

    subroutine sort(arr)
        ! Bubble sort to sort the array in ascending order
        real(dp), intent(inout) :: arr(:)
        integer :: n, i, j
        real(dp) :: temp
        n = size(arr)
        do i = 1, n-1
            do j = i+1, n
                if (arr(i) > arr(j)) then
                    temp = arr(i)
                    arr(i) = arr(j)
                    arr(j) = temp
                end if
            end do
        end do
    end subroutine sort

    subroutine diagonalize(matrix, eigenvalues, eigenvectors)
        ! Diagonalize a symmetric matrix using LAPACK dsyev
        real(dp), intent(in) :: matrix(:,:)
        real(dp), intent(out) :: eigenvalues(:), eigenvectors(:,:)
        integer :: n, lwork, info
        real(dp), allocatable :: work(:)
        n = size(matrix, 1)
        lwork = 3*n - 1
        allocate(work(lwork))
        eigenvectors = matrix
        call dsyev('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, info)
        deallocate(work)
    end subroutine diagonalize

    function diag(v) result(d)
        ! Create a diagonal matrix from a vector
        real(dp), intent(in) :: v(:)
        real(dp) :: d(size(v), size(v))
        integer :: i
        d = 0.0_dp
        do i = 1, size(v)
            d(i,i) = v(i)
        end do
    end function diag

    function outer_product(a, b) result(res)
        ! Compute outer product of two vectors
        real(dp), intent(in) :: a(:), b(:)
        real(dp) :: res(size(a), size(b))
        res = spread(a, 2, size(b)) * spread(b, 1, size(a))
    end function outer_product

[kimala]

Please enclose your code in triple backtick to render it correctly and consider rewriting what I suppose should be the math behind the code. You can use LaTeX in Markdown.

Edit: @Namrah Yes sorry, I don’t have permissions to edit your post, so I copy pasted what you wrote into a code block.

[Namrah]

there are more errors in this code then mine

[kimala]

1

@Namrah

program SCF_He
    implicit none
    integer, parameter :: M = 2
    integer :: io_status
    integer, parameter :: dp = selected_real_kind(15) ! Double precision
    integer :: M, i, j, k, l, iter, max_iter = 9
    real(8), allocatable :: alpha(:), S(:,:), T(:,:), invS(:,:), h(:,:), G(:,:,:,:)
    real(8), allocatable :: two_electron(:,:,:,:)
    real(8) :: eigenvalues(M)
    real(8), allocatable :: eigenvalues(:), eigenvectors(:,:), tmp_mat(:,:)
    real(8), parameter :: Z = 2.0  ! Nuclear charge for helium
    real(8) :: S_bar(M,M)           ! The diagonal eigenvalue matrix (S-bar)
    real(8) :: S_bar_inv_sqrt(M,M)  ! The inverse square root of S-bar
    real(8) :: S_inv(M,M)           ! The final S^(-1/2) matrix
    real(8) :: L_inv(M,M) 
    real(8) :: C_occ(M), F(M,M)
    real(8) :: two_electron_integrals
    real(8), allocatable :: c(:), F(:,:), D(:,:)
    real(8) :: Z = 2.0_dp, alpha_sum, alpha_prod, energy, prev_energy, threshold = 1.0e-8_dp
    real(8) :: N = 2.0_dp ! Number of electrons in Helium
    real(8) :: density_norm, epsilon_1
    integer :: p, q, r, s_idx
    integer :: i, j, k, l
    external :: two_electron_integrals

1. Here already I can spot some errors: you have declared the same entities multiple times. Indeed compiler explorer shows me this error multiple times (with gfortran)

/app/example.f90:19:16:

   19 |     integer :: M, i, j, k, l, iter, max_iter = 9
      |                1
Error: Symbol 'm' at (1) already has basic type of INTEGER
/app/example.f90:23:39:

   23 |     real(8), allocatable :: eigenvalues(:), eigenvectors(:,:), tmp_mat(:,:)
      |                                       1
Error: Symbol 'eigenvalues' at (1) already has basic type of REAL
/app/example.f90:31:35:

   31 |     real(8), allocatable :: c(:), F(:,:), D(:,:)
      |                                   1
Error: Symbol 'f' at (1) already has basic type of REAL

This happens because M, eigenvalues and F are declared twice, which is not permitted.

2

sometimes you use real(8) and sometimes real(dp). You should be consistent and define a precision dp and stick to it. A quick way to do is define a precision module as

module precision
    integer, parameter :: dp = selected_real_kind(15) ! Double precision
end module

and import it where you need as

use precision

and then use the suffix _dp for all fixed values, such as 1.0_dp and so on.

3

You have two different subroutines with the same name, diagonalize. this is not allowed, I suggest you to rename one of the two. In the same way you have two routines diag.

4

After some modifications, you have this code, which still does not compile as S_inv_half is not defined anywhere

module precision
    integer, parameter :: dp = selected_real_kind(15) ! Double precision
end module


! ==================================================
! Function Definition: two_electron_integrals
! ==================================================
function two_electron_integrals(p, q, r, s) result(integral)
    use precision
    implicit none
    integer, intent(in) :: p, q, r, s
    real(dp) :: integral

    ! Example formula (Replace with the correct one)
    integral = 1.0_dp  / ((p + q + r + s + 1.0d-10)**2)

end function two_electron_integrals

program SCF_He
    use precision
    implicit none
    integer, parameter :: M = 2
    integer :: io_status
    
    integer :: i, j, k, l, iter, max_iter = 9
    real(dp), allocatable :: alpha(:), S(:,:), T(:,:), invS(:,:), h(:,:), G(:,:,:,:)
    real(dp), allocatable :: two_electron(:,:,:,:)
    real(dp), allocatable :: eigenvalues(:), eigenvectors(:,:), tmp_mat(:,:)
    real(dp), parameter :: Z = 2.0_dp  ! Nuclear charge for helium
    real(dp) :: S_bar(M,M)           ! The diagonal eigenvalue matrix (S-bar)
    real(dp) :: S_bar_inv_sqrt(M,M)  ! The inverse square root of S-bar
    real(dp) :: S_inv(M,M)           ! The final S^(-1/2) matrix
    real(dp) :: L_inv(M,M) 
    real(dp) :: C_occ(M)
    real(dp) :: two_electron_integrals
    real(dp), allocatable :: c(:), D(:,:), F(:,:)
    real(dp) :: alpha_sum, alpha_prod, energy, prev_energy, threshold = 1.0e-8_dp
    real(dp) :: N = 2.0_dp ! Number of electrons in Helium
    real(dp) :: density_norm, epsilon_1
    integer :: p, q, r, s_idx
    external :: two_electron_integrals




    print *, "🚀 Program started successfully!"
    flush(6)  ! Force printing immediately

    ! Allocate memory
    allocate(alpha(M), S(M,M), T(M,M), invS(M,M), h(M,M), two_electron(M,M,M,M))
    print *, "✅ Memory allocated"
    flush(6)

    ! Open and read basis set parameters
    open(unit=10, file='basis_set.dat', status='old', iostat=io_status)
    if (io_status /= 0) then
        print *, "❌ Error: Cannot open 'basis_set.dat'. Ensure it exists in the correct directory."
        stop
    end if
    print *, "✅ File opened successfully"
    flush(6)

    ! Read values
    do i = 1, M
        read(10, *) alpha(i)
        print *, "✅ Read alpha(", i, ") =", alpha(i)
    end do
    close(10)

    print *, "✅ Completed reading file"
    flush(6)

    C_occ = 0.0_dp 
    C_occ(1) = 1.0_dp   ! Set c_p1 = 1 for p=1, all other c_p1 = 0 for p > 1

    ! Compute Overlap Matrix
    do i = 1, M
        do j = 1, M
            if (i == j) then
                S(i, j) = 1.0_dp 
            else
                S(i, j) = (2.0_dp  * sqrt(alpha(i) * alpha(j)) / (alpha(i) + alpha(j)))**3
            end if
            print *, "S(", i, ",", j, ") =", S(i, j)
        end do
    end do
    print *, "✅ Overlap matrix computed"
    flush(6)

! Diagonalize S: ST = T * S_bar
! Step 1: Diagonalize S → Get T (eigenvectors) and S_bar (eigenvalues)
call diagonalize(S, T, S_bar, M)

! Step 2: Fix small eigenvalues in S_bar
do i = 1, M
    if (S_bar(i,i) < 1.0d-6) then
        print *, "⚠️ WARNING: Eigenvalue too small! Adjusting only S_bar(", i, ",", i, ")"
        if (i == 2) then
            S_bar(i,i) = 1.8381  ! ✅ Set correct value
        end if
    end if
end do

! Step 3: Compute S_bar^{-1/2}
do i = 1, M
    if (S_bar(i,i) > 1.0d-10) then
        S_bar_inv_sqrt(i,i) = 1.0_dp  / sqrt(S_bar(i,i))
    else
        print *, "⚠️ WARNING: Setting S_bar_inv_sqrt(", i, ",", i, ") manually!"
        S_bar_inv_sqrt(i,i) = 1.0_dp  / sqrt(1.8381)  ! ✅ Use manually computed stable value
    end if
end do

! Step 4: Compute final inverse square root matrix S_inv
S_inv = matmul(matmul(T, S_bar_inv_sqrt), transpose(T))

! Step 5: Print final S_inv to check correctness
print *, "⚙️ Final Inverse Square Root of S (S^(-1/2)):"
do i = 1, M
    do j = 1, M
        print *, "S_inv(", i, ",", j, ") =", S_inv(i, j)
    end do
end do
print *, "✅ Final S^(-1/2) Computed"

! Print Eigenvalues (S_bar)
print *, "⚙️ Diagonalized S (S-bar):"
do i = 1, M
    print *, "S_bar(", i, ",", i, ") =", S_bar(i, i)
end do
print *, "✅ S-bar Printed"

! Print Eigenvectors (T)
print *, "⚙️ Eigenvectors (T):"
call print_matrix(T, M)
print *, "✅ Eigenvectors Printed"


! Initialize S_bar_inv_sqrt as zero
S_bar_inv_sqrt = 0.0_dp   

! Compute inverse square root of eigenvalues with corrected values
do i = 1, M
    if (S_bar(i,i) > 1.0d-6) then
        S_bar_inv_sqrt(i,i) = 1.0_dp  / sqrt(S_bar(i,i))
    else if (i == 2) then
        print *, "⚠️ WARNING: S_bar(2,2) too small, setting manually to 1/sqrt(1.8381)."
        S_bar_inv_sqrt(i,i) = 1.0_dp  / sqrt(1.8381)  ! ✅ Fix small eigenvalue
    end if
end do

! Print Inverse Square Root of S-bar
print *, "⚙️ Inverse Square Root of S-bar (S_bar^{-1/2}):"
call print_matrix(S_bar_inv_sqrt, M)
print *, "✅ Inverse Square Root of S-bar Printed"


! Compute S^(-1/2)
S_inv = matmul(matmul(T, S_bar_inv_sqrt), transpose(T))

! Print S^(-1/2) to verify
print *, "⚙️ Computed Inverse Square Root of S (S^(-1/2)):"
call print_matrix(S_inv, M)
print *, "✅ Inverse Square Root of S Printed"


! Call diagonalization
    call jacobi_diagonalization(S, M, eigenvalues, T)
    print *, "✅ Diagonalization completed"
    flush(6)

   ! Compute One-electron integral h_pq
    print *, "⚙️ Computing One-electron integrals..."
    do i = 1, M
     do j = 1, M
        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)))
        print *, "h(", i, ",", j, ") =", h(i,j)
     end do
    end do
   print *, "✅ One-electron integrals computed"
   flush(6)

    ! Compute Two-electron integral (pq|rs)
     print *, "⚙️ Computing Two-electron integrals..."
     do i = 1, M
        do j = 1, M
           do k = 1, M
              do l = 1, M
                two_electron_integrals(i,j,k,l) = 32.0_dp  * (alpha(i) * alpha(j) * alpha(k) * alpha(l))**(3.0_dp /2.0_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) )
                print *, "( ", i, j, "|", k, l, ") =", two_electron(i,j,k,l)
               end do
           end do
        end do
    end do
    print *, "✅ Two-electron integrals computed"
    flush(6)

! Initialize coefficients cp1 (before computing F matrix)
C_occ = 0.0_dp 
C_occ(1) = 1.0_dp   ! Set c_p1 = 1 for p=1

! Part B(a): Compute the Fock Matrix F_pq
F = 0.0_dp 
do p = 1, M
    do q = 1, M
        F(p,q) = h(p,q)  ! Start with one-electron integral

        ! Add two-electron contributions
        do r = 1, M
            do s_idx = 1, M
                F(p,q) = F(p,q) + C_occ(r) * C_occ(s_idx) * two_electron_integrals(p,q,r,s_idx)
            end do
        end do
    end do
end do

print *, "⚙️ Fock Matrix F:"
do i = 1, M
    do j = 1, M
        print *, "F(", i, ",", j, ") =", F(i, j)
    end do
end do
print *, "✅ Fock Matrix Computed"

! Part 2: SCF Iterations
    allocate(c(M), F(M,M), D(M,M))
    c = [1.0_dp, 0.0_dp] ! Initial guess: c1=1, c2=0
    prev_energy = 0.0_dp
    energy = 1.0_dp ! Initial dummy value

    print *, "Iteration | c1        | c2        | F11         | F12         | F22         | epsilon     | Energy (a.u.)"
    print *, "--------------------------------------------------------------------------------------------------------"

    do iter = 1, max_iter
        ! Part 2a: Compute Fock matrix F = H + G
        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

        ! Part 2b: Löwdin orthogonalization
        ! i. Compute F' = S^{-1/2} * F * S^{-1/2}
        tmp_mat = matmul(S_inv_half, matmul(F, S_inv_half))

        ! ii. Diagonalize F' to obtain eigenvectors C'
        call diagonalize2(tmp_mat, eigenvalues, tmp_mat)

        ! iii. Compute C = S^{-1/2} * C'
        c = matmul(S_inv_half, tmp_mat(:,1))

        ! iv. Select the lowest eigenvalue (occupied orbital energy)
        epsilon_1 = eigenvalues(1)

        ! Part 2c: Compute density matrix D and check normalization
        D = outer_product(c, c)
        density_norm = sum(D * S)

        ! Part 2d: Compute total energy
        energy = 2.0_dp * sum(H * D) + sum(D * G)

        ! Print iteration results
        print *, iter, c(1), c(2), F(1,1), F(1,2), F(2,2), epsilon_1, energy

        ! Check for convergence
        if (abs(energy - prev_energy) < threshold) exit
        prev_energy = energy
    end do

    print *, "Converged Energy:", prev_energy

! Program finished
    print *, "🎉 Program completed successfully!"
    flush(6)

    ! Deallocate memory
    deallocate(alpha, S, T, invS, h, two_electron)

contains

    ! Subroutine for Jacobi diagonalization
    subroutine jacobi_diagonalization(A, N, eigenvalues, T)
        use precision
        implicit none
        integer, intent(in) :: N
        real(dp), intent(inout) :: A(N,N)
        real(dp), intent(out) :: eigenvalues(N)
        real(dp), intent(out) :: T(N,N)

        ! Declare all variables explicitly
        integer :: i, j, max_iter
        real(dp) :: off_diag, tol

        print *, "⚙️ Starting Jacobi diagonalization..."
        flush(6)

        max_iter = 100
        tol = 1.0d-10
        T = 0.0_dp 
        do i = 1, N
            T(i, i) = 1.0_dp 
        end do

        do max_iter = 1, 50
            off_diag = 0.0_dp 
            do i = 1, N
                do j = i+1, N
                    off_diag = max(off_diag, abs(A(i, j)))
                end do
            end do
            if (off_diag < tol) then
                print *, "✅ Jacobi diagonalization completed"
                exit
            end if
        end do

        do i = 1, N
            eigenvalues(i) = A(i, i)
            print *, "Eigenvalue(", i, ") =", eigenvalues(i)
        end do
        flush(6)
    end subroutine jacobi_diagonalization

    ! Function to create diagonal matrix
    function diag(vals, N) result(D)
        implicit none
        integer, intent(in) :: N
        real(dp), intent(in) :: vals(N)
        real(dp) :: D(N,N)
        integer :: i
        D = 0.0_dp 
        do i = 1, N
            D(i, i) = vals(i)
        end do
    end function diag

    function diag2(v) result(d)
        use precision
        ! Create a diagonal matrix from a vector
        real(dp), intent(in) :: v(:)
        real(dp) :: d(size(v), size(v))
        integer :: i
        d = 0.0_dp
        do i = 1, size(v)
            d(i,i) = v(i)
        end do
    end function diag2


subroutine diagonalize(A, V, eigenvalues, N)
    use precision
    implicit none
    integer, intent(in) :: N
    real(dp), intent(inout) :: A(N,N)
    real(dp), intent(out) :: V(N,N), eigenvalues(N)
    
    integer :: i, j, p, q, iter, max_iter
    real(dp) :: tol, off_diag, theta, t, c, s, tau, A_pq, temp

    ! Initialize V as the identity matrix
    V = 0.0_dp 
    do i = 1, N
        V(i,i) = 1.0_dp 
    end do

    max_iter = 100
    tol = 1.0d-10

    ! Jacobi Iteration
    do iter = 1, max_iter
        ! Compute sum of off-diagonal elements
        off_diag = 0.0_dp 
        do i = 1, N
            do j = i+1, N
                off_diag = off_diag + abs(A(i,j))
            end do
        end do

        ! If off-diagonal elements are small, exit
        if (off_diag < tol) exit

        ! Find largest off-diagonal element
        A_pq = 0.0_dp 
        do p = 1, N
            do q = p+1, N
                if (abs(A(p,q)) > A_pq) then
                    A_pq = abs(A(p,q))
                    i = p
                    j = q
                end if
            end do
        end do

        ! Compute rotation angle
        theta = 0.5d0 * atan2(2.0_dp  * A(i,j), A(j,j) - A(i,i))
        c = cos(theta)
        s = sin(theta)
        tau = s / (1.0_dp  + c)

        ! Apply rotation to matrix A
        temp = A(i,i)
        A(i,i) = c*c * temp - 2.0_dp  * s * c * A(i,j) + s*s * A(j,j)
        A(j,j) = s*s * temp + 2.0_dp  * s * c * A(i,j) + c*c * A(j,j)
        A(i,j) = 0.0_dp 
        A(j,i) = 0.0_dp 

        do p = 1, N
            if (p /= i .and. p /= j) then
                temp = A(p,i)
                A(p,i) = c * temp - s * A(p,j)
                A(p,j) = s * temp + c * A(p,j)
                A(i,p) = A(p,i)
                A(j,p) = A(p,j)
            end if
        end do

        ! Update eigenvector matrix V
        do p = 1, N
            temp = V(p,i)
            V(p,i) = c * temp - s * V(p,j)
            V(p,j) = s * temp + c * V(p,j)
        end do
    end do

    ! Store final diagonal values as eigenvalues
    do i = 1, N
        eigenvalues(i) = A(i,i)
    end do
end subroutine diagonalize

subroutine print_matrix(A, N)
    use precision
    implicit none
    integer, intent(in) :: N
    real(dp), intent(in) :: A(N,N)
    integer :: i, j

    print *, "⚙️ Printing Matrix:"
    do i = 1, N
        do j = 1, N
            print *, "A(", i, ",", j, ") =", A(i, j)
        end do
    end do
    print *, "✅ Matrix Printed"
end subroutine print_matrix

    subroutine sort(arr)
        use precision
        ! Bubble sort to sort the array in ascending order
        real(dp), intent(inout) :: arr(:)
        integer :: n, i, j
        real(dp) :: temp
        n = size(arr)
        do i = 1, n-1
            do j = i+1, n
                if (arr(i) > arr(j)) then
                    temp = arr(i)
                    arr(i) = arr(j)
                    arr(j) = temp
                end if
            end do
        end do
    end subroutine sort

    subroutine diagonalize2(matrix, eigenvalues, eigenvectors)
        use precision
        ! Diagonalize a symmetric matrix using LAPACK dsyev
        real(dp), intent(in) :: matrix(:,:)
        real(dp), intent(out) :: eigenvalues(:), eigenvectors(:,:)
        integer :: n, lwork, info
        real(dp), allocatable :: work(:)
        n = size(matrix, 1)
        lwork = 3*n - 1
        allocate(work(lwork))
        eigenvectors = matrix
        call dsyev('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, info)
        deallocate(work)
    end subroutine diagonalize2



    function outer_product(a, b) result(res)
        use precision
        ! Compute outer product of two vectors
        real(dp), intent(in) :: a(:), b(:)
        real(dp) :: res(size(a), size(b))
        res = spread(a, 2, size(b)) * spread(b, 1, size(a))
    end function outer_product

end program SCF_He

The last errors to check are

/app/example.f90:254:35:

  254 |         tmp_mat = matmul(S_inv_half, matmul(F, S_inv_half))
      |                                   1
Error: Symbol 's_inv_half' at (1) has no IMPLICIT type; did you mean 's_inv'?
/app/example.f90:191:16:

  191 |                 two_electron_integrals(i,j,k,l) = 32.0_dp  * (alpha(i) * alpha(j) * alpha(k) * alpha(l))**(3.0_dp /2.0_dp ) * &
      |                1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.
/app/example.f90:254:35:

  254 |         tmp_mat = matmul(S_inv_half, matmul(F, S_inv_half))
      |                                   1
Error: Symbol 's_inv_half' at (1) has no IMPLICIT type; did you mean 's_inv'?
/app/example.f90:254:35:

  254 |         tmp_mat = matmul(S_inv_half, matmul(F, S_inv_half))
      |                                   1
Error: Symbol 's_inv_half' at (1) has no IMPLICIT type; did you mean 's_inv'?
/app/example.f90:270:44-48:

  270 |         energy = 2.0_dp * sum(H * D) + sum(D * G)
      |                                            1   2
Error: Inconsistent ranks for operator at (1) and (2)

where the one at line 270 is simply because D*G involves very different arrays and you have to do it by yourself.

I have no idea what the following means though, and I would like to hear about it from others, it’s the first time I see this:

/app/example.f90:191:16:

  191 |                 two_electron_integrals(i,j,k,l) = 32.0_dp  * (alpha(i) * alpha(j) * alpha(k) * alpha(l))**(3.0_dp /2.0_dp ) * &
      |                1
Error: The function result on the lhs of the assignment at (1) must have the pointer attribute.

[Namrah]

yes i am doing dear i have now just two errors remaining i am just worried its second hour its not solved yetv:(

[Namrah]

actually i have zero idea of programming

[kimala]

@Namrah

i am just worried its second hour its not solved yetv:(

Are you taking an exam?

[Namrah]

no dear its mine assignment graded and i am not able to do it

[kimala]

@kimala

energy = 2.0_dp * sum(H * D) + sum(D * G)

D and G have different ranks, one is two-dimensional and the other is four-dimensional. This means that D*G is not defined. You will have to compute the energy yourself by writing four do loops, something like

     energy = 0._dp
     do i = 1, M
        do j = 1, M
           do k = 1, M
              do l = 1, M
                energy = energy + 2._dp * H( ) *D( ) + D( )*G( )
               end do
           end do
        end do
    end do

with something inside the parentheses to do the multiplication correctly.

[Namrah]

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

[kimala]

I will not. Please format the code.

[Namrah]

ok np

[jorgeg]

I know this a fortran discourse but clearly your problem is how to structure the logic of the program. I would first do a prototype in Python where you can just use numpy for many things. Once you get correct results in Python start migrating the code to Fortran function by function, checking the correctness of things along the way.