Added restoration of symmetry

src/determinants/density_matrix.irp.f

@@ -280,6 +280,8 @@ subroutine save_natural_mos
    ! the |MO| basis
    END_DOC
    call set_natural_mos
+   call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10)
+   call orthonormalize_mos
    call save_mos
 end
 

src/hartree_fock/scf.irp.f

@@ -1,34 +1,34 @@
 program scf
   BEGIN_DOC
-! 
+!
 ! The :ref:`scf` program performs *Restricted* Hartree-Fock
 ! calculations (the spatial part of the |MOs| is common for alpha and beta
 ! spinorbitals).
-! 
+!
 ! It performs the following actions:
-! 
+!
 ! #. Compute/Read all the one- and two-electron integrals, and store them
 !    in memory
 ! #. Check in the |EZFIO| database if there is a set of |MOs|.
 !    If there is, it will read them as initial guess. Otherwise, it will
 !    create a guess.
 ! #. Perform the |SCF| iterations
-! 
+!
 ! For the keywords related to the |SCF| procedure, see the ``scf_utils``
 ! directory where you will find all options.
-! 
+!
 ! At each iteration, the |MOs| are saved in the |EZFIO| database. Hence,
 ! if the calculation crashes for any unexpected reason, the calculation
 ! can be restarted by running again the |SCF| with the same |EZFIO|
 ! database.
-! 
+!
 ! To start again a fresh |SCF| calculation, the |MOs| can be reset by
 ! running the :ref:`qp_reset` command.
-! 
-! The `DIIS`_ algorithm is implemented, as well as the `level-shifting`_ 
+!
+! The `DIIS`_ algorithm is implemented, as well as the `level-shifting`_
 ! method. If the |SCF| does not converge, try again with a higher value of
 ! :option:`level_shift`.
-! 
+!
 ! .. _DIIS: https://en.wikipedia.org/w/index.php?title=DIIS
 ! .. _level-shifting: https://doi.org/10.1002/qua.560070407
 !
@@ -55,6 +55,7 @@ subroutine create_guess
           size(mo_one_e_integrals,1),                            &
           size(mo_one_e_integrals,2),                            &
           mo_label,1,.false.)
+      call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10)
       SOFT_TOUCH mo_coef mo_label
     else if (mo_guess_type == "Huckel") then
       call huckel_guess

src/mo_basis/EZFIO.cfg

@@ -36,4 +36,3 @@ size: (mo_basis.mo_num)
 type: character*(32)
 doc: MD5 checksum characterizing the |AO| basis set.
 interface: ezfio
-

src/mo_guess/h_core_guess_routine.irp.f

@@ -8,6 +8,7 @@ subroutine hcore_guess
   call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals,          &
                                      size(mo_one_e_integrals,1),  &
                                      size(mo_one_e_integrals,2),label,1,.false.)
+  call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 )
   call save_mos
   SOFT_TOUCH mo_coef mo_label
 end

src/mo_one_e_ints/orthonormalize.irp.f

@@ -1,11 +1,16 @@
 subroutine orthonormalize_mos
   implicit none
-  integer :: m,p,s
+  integer :: m,p,s,i
   m = size(mo_coef,1)
   p = size(mo_overlap,1)
-  call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num,lin_dep_cutoff)
-  mo_label = 'Orthonormalized'
-  SOFT_TOUCH mo_coef mo_label
+  do i=1,4
+    call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num,lin_dep_cutoff)
+    call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10)
+  enddo
+  if (restore_symm) then
+     call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12)
+  endif
+  SOFT_TOUCH mo_coef
 end
 
 

src/scf_utils/EZFIO.cfg

@@ -51,3 +51,8 @@ doc: If true, leave untouched all the orbitals defined as core and optimize all
 interface: ezfio,provider,ocaml
 default: False
 
+[restore_symm]
+type: logical
+doc: If true, try to find symmetry in the MO coefficient matrices
+interface: ezfio,provider,ocaml
+default: True

src/scf_utils/huckel.irp.f

@@ -25,6 +25,7 @@ subroutine huckel_guess
 
   TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta
   mo_coef = eigenvectors_fock_matrix_mo
+  call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 )
   SOFT_TOUCH mo_coef
   call save_mos
   deallocate(A)

src/scf_utils/roothaan_hall_scf.irp.f

@@ -165,6 +165,8 @@ END_DOC
 
   if(.not.frozen_orb_scf)then
    call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
+   call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10)
+   call orthonormalize_mos
    call save_mos
   endif
 
@@ -303,7 +305,7 @@ END_DOC
       Fock_matrix_AO_(i,j) = 0.d0
     enddo
     do k=1,dim_DIIS
-      if (dabs(X_vector_DIIS(k)) < 1.d-10) cycle
+      if (dabs(X_vector_DIIS(k)) < 1.d-12) cycle
       do i=1,ao_num
         Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) +            &
             X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1)

src/tools/save_ortho_mos.irp.f

@@ -10,5 +10,7 @@ program save_ortho_mos
   ! Thanks to the Lowdin orthonormalization, the new MOs are the most similar to the guess MOs.
   END_DOC
   call orthonormalize_mos
+  mo_label = 'Orthonormalized'
+  SOFT_TOUCH mo_label
   call save_mos
 end

src/utils/linear_algebra.irp.f

@@ -1367,3 +1367,92 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
   deallocate(A,eigenvalues)
 end
 
+subroutine nullify_small_elements(m,n,A,LDA,thresh)
+  implicit none
+  integer, intent(in) :: m,n,LDA
+  double precision, intent(inout) :: A(LDA,n)
+  double precision, intent(in) :: thresh
+  integer :: i,j
+  double precision :: amax
+
+  ! Find max value
+  amax = 0.d0
+  do j=1,n
+    do i=1,m
+      amax = max(dabs(A(i,j)), amax)
+    enddo
+  enddo
+  amax = 1.d0/amax
+
+  ! Remove tiny elements
+  do j=1,n
+    do i=1,m
+      if ( dabs(A(i,j) * amax) < thresh ) then
+         A(i,j) = 0.d0
+      endif
+    enddo
+  enddo
+
+end
+
+subroutine restore_symmetry(m,n,A,LDA,thresh)
+  implicit none
+  integer, intent(in) :: m,n,LDA
+  double precision, intent(inout) :: A(LDA,n)
+  double precision, intent(in) :: thresh
+  integer :: i,j,k,l
+  logical, allocatable :: done(:,:)
+  double precision :: f, g, count, thresh2
+  thresh2 = dsqrt(thresh)
+  call nullify_small_elements(m,n,A,LDA,thresh)
+
+  allocate(done(m,n))
+
+  do j=1,n
+    do i=1,m
+      done(i,j) = A(i,j) == 0.d0
+    enddo
+  enddo
+
+  do j=1,n
+    do i=1,m
+      if ( done(i,j) ) cycle
+      done(i,j) = .True.
+      count = 1.d0
+      f = 1.d0/A(i,j)
+      do l=1,n
+        do k=1,m
+          if ( done(k,l) ) cycle
+          g = f * A(k,l)
+          if ( dabs(dabs(g) - 1.d0) < thresh2 ) then
+            count = count + 1.d0
+            if (g>0.d0) then
+              A(i,j) = A(i,j) + A(k,l)
+            else
+              A(i,j) = A(i,j) - A(k,l)
+            end if
+          endif
+        enddo
+      enddo
+      if (count > 1.d0) then
+        A(i,j) = A(i,j) / count
+        do l=1,n
+          do k=1,m
+            if ( done(k,l) ) cycle
+            g = f * A(k,l)
+            if ( dabs(dabs(g) - 1.d0) < thresh2 ) then
+              done(k,l) = .True.
+              if (g>0.d0) then
+                A(k,l) = A(i,j)
+              else
+                A(k,l) = -A(i,j)
+              end if
+            endif
+          enddo
+        enddo
+      endif
+
+    enddo
+  enddo
+
+end