# Removal of error from fortran code f.90

**URL:** <https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431>\
**Category:** Homework\
**Created:** [March 21, 2025, 2:58pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431 "2025-03-21T14:58:34Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Namrah](https://avatars.discourse-cdn.com/v4/letter/n/35a633/32.png) [@Namrah](https://fortran-lang.discourse.group/u/Namrah)\
**Post date:** [March 21, 2025, 2:58pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/1 "2025-03-21T14:58:34Z")

</div>

_Can anyone help me to remove just errors from this code. I was trying to solve the problem but due to errors its not running. Its fortran f.90_  
! =========================================================  
! 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 = 2_epsilon(1) - J(1,1):"  
print \*, " \>\> SCF Energy =“, energy  
print _, " \>\> 2_epsilon(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 ~ 2_epsilon(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\_dp_mo\_energy\_1 - J11

print _, “------------------------------------------------------------”  
print , "Check from simple formula E = 2epsilon(1) - J(1,1):"  
print , " \>\> 2epsilon(1) =", 2.0\_dp_mo\_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

---

<div class="post-metadata">

**Author:** ![themos](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/themos/32/521_2.png) [@themos](https://fortran-lang.discourse.group/u/themos)\
**Post date:** [March 21, 2025, 3:26pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/2 "2025-03-21T15:26:29Z")

</div>

Please upload the source on [https://godbolt.org/](https://godbolt.org/) and share the link to it.

---

<div class="post-metadata">

**Author:** ![kimala](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/kimala/32/2707_2.png) [@kimala](https://fortran-lang.discourse.group/u/kimala)\
**Post date:** [March 21, 2025, 3:28pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/3 "2025-03-21T15:28:48Z")

</div>

Is this the same code as here?

> [@Help me with fortran f.90 and running on terminal](https://fortran-lang.discourse.group/t/help-me-with-fortran-f-90-and-running-on-terminal/9419/13):
>
> 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 ! an…

Please do not duplicate topics

---

<div class="post-metadata">

**Author:** ![Namrah](https://avatars.discourse-cdn.com/v4/letter/n/35a633/32.png) [@Namrah](https://fortran-lang.discourse.group/u/Namrah)\
**Post date:** [March 21, 2025, 3:42pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/4 "2025-03-21T15:42:38Z")

</div>

> **[Compiler Explorer - Fortran (x86-64 gfortran (trunk))](https://godbolt.org/z/bbYx94P6s)**
>
> ! =========================================================
> ! A Minimal SCF Example for Helium (M = 2 basis functions)
> ! Program name: namrah
> ! =========================================================
> ! - Reads 2 exponents alpha(1..2) from...

---

<div class="post-metadata">

**Author:** ![Namrah](https://avatars.discourse-cdn.com/v4/letter/n/35a633/32.png) [@Namrah](https://fortran-lang.discourse.group/u/Namrah)\
**Post date:** [March 21, 2025, 3:43pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/5 "2025-03-21T15:43:20Z")

</div>

ok its my first time on this website so please ignore

---

<div class="post-metadata">

**Author:** ![themos](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/themos/32/521_2.png) [@themos](https://fortran-lang.discourse.group/u/themos)\
**Post date:** [March 21, 2025, 4:03pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/6 "2025-03-21T16:03:56Z")

</div>

Ok, so now go through the compiler errors one by one and understand why the compiler is complaining. First one I see is at line 268, “Unexpected data declaration statement”. In Fortran, you are not allowed to sprinkle data declaration statements among the executable statements. You may have meant to write only the assignment at line 268, so do that, and add a declaration for this variable near the top with all the other declarations of real variables.

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [March 21, 2025, 5:52pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/7 "2025-03-21T17:52:44Z")

</div>

> [@Namrah](#):
>
> ```auto
> ! (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)
> 
> ```

Since this is a homework problem, these steps are perhaps intended to be part of the education. Indeed, I remember back in the 1970s when I wrote my first jacobi diagonalization routine, a light bulb lit up in my head and I began to understand what matrix diagonalization really was. But with that in mind, I might point out that those steps are done here in order to solve what is called the generalized symmetric eigenproblem. Nowadays, there are standard routines that solve that problem directly. The LAPACK routine that does so is DSYGV which is described here: [https://www.netlib.org/lapack/explore-html-3.6.1/d2/d8a/group\_\_double\_s\_yeigen\_gaafa856f208a7f658918a11e134b2c2d2.html](https://www.netlib.org/lapack/explore-html-3.6.1/d2/d8a/group__double_s_yeigen_gaafa856f208a7f658918a11e134b2c2d2.html)

That routine actually solves three closely related problems that basically amount to whether the metric matrix or its inverse is given as the argument. If you dig down into the lapack code you will see that it also uses a slightly different sequence of steps to solve the problem because the symmetric orthonormalization and jacobi steps are not very efficient, or accurate for that matter. However, unless the student works through the steps doing it the wrong way, it doesn’t mean as much when he finally begins to do it the right way, so I’ll just point out that there are better ways to solve that linear algebra problem that the student should learn eventually, and leave it there.

---

<div class="post-metadata">

**Author:** ![Namrah](https://avatars.discourse-cdn.com/v4/letter/n/35a633/32.png) [@Namrah](https://fortran-lang.discourse.group/u/Namrah)\
**Post date:** [March 21, 2025, 7:51pm UTC](https://fortran-lang.discourse.group/t/removal-of-error-from-fortran-code-f-90/9431/8 "2025-03-21T19:51:55Z")

</div>

Thank you for your help !