From 88cffcb26999f685b9c7ef78d61bb71961cf3d9d Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 5 Apr 2024 17:51:48 +0200
Subject: [PATCH] Force MOs to be on axes. Nice for atoms

---
 src/ao_one_e_ints/ao_ortho_canonical.irp.f |  2 ++
 src/scf_utils/diagonalize_fock.irp.f       |  2 +-
 src/scf_utils/roothaan_hall_scf.irp.f      | 26 +++++++++++++++++++---
 3 files changed, 26 insertions(+), 4 deletions(-)

diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f
index 668b920d..eff7e7be 100644
--- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f
+++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f
@@ -138,6 +138,8 @@ END_PROVIDER
     deallocate(S)
   endif
 
+  FREE ao_overlap
+
 END_PROVIDER
 
 BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]
diff --git a/src/scf_utils/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f
index 5188581a..b9042b29 100644
--- a/src/scf_utils/diagonalize_fock.irp.f
+++ b/src/scf_utils/diagonalize_fock.irp.f
@@ -47,7 +47,7 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
        do j = 1, n_core_orb
          jorb = list_core(j)
          F(iorb,jorb) = 0.d0
-         F(jorb,iorb) = 0.d0                                                                                                                 
+         F(jorb,iorb) = 0.d0
        enddo
      enddo
    endif
diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f
index 730cb496..3f5c8549 100644
--- a/src/scf_utils/roothaan_hall_scf.irp.f
+++ b/src/scf_utils/roothaan_hall_scf.irp.f
@@ -13,9 +13,9 @@ END_DOC
   integer                        :: iteration_SCF,dim_DIIS,index_dim_DIIS
 
   logical                        :: converged
-  integer                        :: i,j
+  integer                        :: i,j,m
   logical, external              :: qp_stop
-  double precision, allocatable :: mo_coef_save(:,:)
+  double precision, allocatable :: mo_coef_save(:,:), S(:,:)
 
   PROVIDE ao_md5 mo_occ level_shift
 
@@ -208,9 +208,29 @@ END_DOC
       size(Fock_matrix_mo,2),mo_label,1,.true.)
    call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10)
    call orthonormalize_mos
-   call save_mos
   endif
 
+
+  ! Identify degenerate MOs and force them on the axes
+  allocate(S(ao_num,ao_num))
+  i=1
+  do while (i<mo_num)
+    j=i
+    m=1
+    do while ( (j<mo_num).and.(fock_matrix_diag_mo(j+1)-fock_matrix_diag_mo(i) < 1.d-8) )
+      j += 1
+      m += 1
+    enddo
+    if (m>1) then
+      call dgemm('N','T',ao_num,ao_num,m,1.d0,mo_coef(1,i),size(mo_coef,1),mo_coef(1,i),size(mo_coef,1),0.d0,S,size(S,1))
+      call pivoted_cholesky( S, m, -1.d0, ao_num, mo_coef(1,i))
+    endif
+    i = j+1
+  enddo
+
+
+  call save_mos
+
   call write_double(6, Energy_SCF, 'SCF energy')
 
   call write_time(6)