add maximum overlap method in hartree fock

src/scf_utils/EZFIO.cfg

@@ -45,6 +45,12 @@ type: double precision
 doc: Calculated HF energy
 interface: ezfio
 
+[do_mom]
+type: logical
+doc: If true, this will run a MOM calculation. The overlap will be computed at each step with respect to the initial MOs. After an initial Hartree-Fock calculation, the guess can be created by swapping molecular orbitals through the qp run swap_mos command. 
+interface: ezfio,provider,ocaml
+default: False
+
 [frozen_orb_scf]
 type: logical
 doc: If true, leave untouched all the orbitals defined as core and optimize all the orbitals defined as active with qp_set_mo_class

src/scf_utils/reorder_mo_max_overlap.irp.f

@@ -0,0 +1,96 @@
+ subroutine reorder_mo_max_overlap
+  implicit none
+   BEGIN_DOC
+   ! routines that compute the projection of each MO of the current `mo_coef` on the space spanned by the occupied orbitals of `mo_coef_begin_iteration`
+   END_DOC
+   integer :: i,j,k,l
+   double precision, allocatable :: overlap(:,:)
+   double precision, allocatable :: proj(:)
+   integer, allocatable :: iorder(:)
+   double precision, allocatable :: mo_coef_tmp(:,:)
+   double precision, allocatable :: tmp(:,:)
+   allocate(overlap(mo_num,mo_num),proj(mo_num),iorder(mo_num),mo_coef_tmp(ao_num,mo_num),tmp(mo_num,ao_num))
+
+   overlap(:,:) = 0d0
+   mo_coef_tmp(:,:) = 0d0
+   proj(:) = 0d0
+   iorder(:) = 0d0
+   tmp(:,:) = 0d0
+
+   ! These matrix products compute the overlap bewteen the initial and the current MOs
+   call dgemm('T','N', mo_num, ao_num, ao_num, 1.d0, &
+        mo_coef_begin_iteration, size(mo_coef_begin_iteration,1), &
+        ao_overlap, size(ao_overlap,1), 0.d0, &
+        tmp, size(tmp,1))
+
+   call dgemm('N','N', mo_num, mo_num, ao_num, 1.d0, &
+        tmp, size(tmp,1), &
+        mo_coef, size(mo_coef, 1), 0.d0, &
+        overlap, size(overlap,1) )
+
+
+   ! for each orbital compute the best overlap
+   do i = 1, mo_num
+      iorder(i) = i ! initialize the iorder list as we need it to sort later
+      do j = 1, elec_alpha_num
+         proj(i) += overlap(j,i)*overlap(j,i) ! compute the projection of current orbital i on the occupied space of the initial orbitals
+      enddo
+      proj(i) = dsqrt(proj(i))
+   enddo
+   ! sort the list of projection to find the mos with the largest overlap
+   call dsort(proj(:),iorder(:),mo_num)
+   ! reorder orbitals according to projection
+   do i=1,mo_num
+      mo_coef_tmp(:,i) = mo_coef(:,iorder(mo_num+1-i)) 
+   enddo
+
+   ! update the orbitals
+   mo_coef(:,:) =  mo_coef_tmp(:,:)
+
+   ! if the determinant is open-shell we need to make sure that the singly occupied orbital correspond to the initial ones
+   if (elec_alpha_num > elec_beta_num) then
+      double precision, allocatable :: overlap_alpha(:,:)
+      double precision, allocatable :: proj_alpha(:)
+      integer, allocatable :: iorder_alpha(:)
+      allocate(overlap_alpha(mo_num,elec_alpha_num),proj_alpha(elec_alpha_num),iorder_alpha(elec_alpha_num))
+      overlap_alpha(:,:) = 0d0
+      mo_coef_tmp(:,:) = 0d0
+      proj_alpha(:) = 0d0
+      iorder_alpha(:) = 0d0
+      tmp(:,:) = 0d0 
+      ! These matrix products compute the overlap bewteen the initial and the current MOs
+      call dgemm('T','N', mo_num, ao_num, ao_num, 1.d0, &
+           mo_coef_begin_iteration, size(mo_coef_begin_iteration,1), &
+           ao_overlap, size(ao_overlap,1), 0.d0, &
+           tmp, size(tmp,1))
+
+      call dgemm('N','N', mo_num, elec_alpha_num, ao_num, 1.d0, &
+           tmp, size(tmp,1), &
+           mo_coef, size(mo_coef, 1), 0.d0, &
+           overlap_alpha, size(overlap_alpha,1) )
+
+      do i = 1, elec_alpha_num
+         iorder_alpha(i) = i ! initialize the iorder list as we need it to sort later
+         do j = 1, elec_beta_num
+            proj_alpha(i) += overlap_alpha(j,i)*overlap_alpha(j,i) ! compute the projection of current orbital i on the beta occupied space of the initial orbitals
+         enddo
+         proj_alpha(i) = dsqrt(proj_alpha(i))
+      enddo
+
+      ! sort the list of projection to find the mos with the largest overlap
+      call dsort(proj_alpha(:),iorder_alpha(:),elec_alpha_num)
+      ! reorder orbitals according to projection
+      do i=1,elec_alpha_num
+         mo_coef_tmp(:,i) = mo_coef(:,iorder_alpha(elec_alpha_num+1-i)) 
+      enddo
+      do i=1,elec_alpha_num
+         mo_coef(:,i) = mo_coef_tmp(:,i) 
+      enddo
+
+      deallocate(overlap_alpha, proj_alpha, iorder_alpha)
+   endif
+
+   deallocate(overlap, proj, iorder, mo_coef_tmp, tmp)
+
+ end
+

src/scf_utils/roothaan_hall_scf.irp.f

@@ -51,6 +51,11 @@ subroutine Roothaan_Hall_SCF
 !
   PROVIDE FPS_SPF_matrix_AO Fock_matrix_AO 
 
+  ! Initialize MO to run IMOM
+  if(do_mom)then
+     call initialize_mo_coef_begin_iteration
+  endif
+
   converged = .False.
   do while ( .not.converged .and. (iteration_SCF < n_it_SCF_max) )
 
@@ -88,16 +93,17 @@ subroutine Roothaan_Hall_SCF
       Fock_matrix_AO_beta  = Fock_matrix_AO*0.5d0
       TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta
 
-    endif
-
+   endif
     MO_coef = eigenvectors_Fock_matrix_MO
+    if(do_mom)then
+       call reorder_mo_max_overlap
+    endif
     if(frozen_orb_scf)then
-     call reorder_core_orb
-     call initialize_mo_coef_begin_iteration
+      call reorder_core_orb
+      call initialize_mo_coef_begin_iteration
     endif
 
     TOUCH MO_coef
-
 !   Calculate error vectors
 
     max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO))
@@ -106,41 +112,46 @@ subroutine Roothaan_Hall_SCF
 
     energy_SCF = SCF_energy
     Delta_Energy_SCF = energy_SCF - energy_SCF_previous
-    if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then
+    if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0).and.(.not.do_mom) ) then
       Fock_matrix_AO(1:ao_num,1:ao_num) = Fock_matrix_DIIS (1:ao_num,1:ao_num,index_dim_DIIS)
       Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0
       Fock_matrix_AO_beta  = Fock_matrix_AO*0.5d0
       TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta
     endif
 
-    double precision :: level_shift_save
-    level_shift_save = level_shift
-    mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num)
-    do while (Delta_energy_SCF > 0.d0)
-      mo_coef(1:ao_num,1:mo_num) = mo_coef_save
-      if (level_shift <= .1d0) then
-        level_shift = 1.d0
-      else
-        level_shift = level_shift * 3.0d0
-      endif
-      TOUCH mo_coef level_shift
-      mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_num)
-      if(frozen_orb_scf)then
-        call reorder_core_orb
-        call initialize_mo_coef_begin_iteration
-      endif
-      TOUCH mo_coef
-      Delta_Energy_SCF = SCF_energy - energy_SCF_previous
-      energy_SCF = SCF_energy
-      if (level_shift-level_shift_save > 40.d0) then
-        level_shift = level_shift_save * 4.d0
-        SOFT_TOUCH level_shift
-        exit
-      endif
-      dim_DIIS=0
-    enddo
-    level_shift = level_shift * 0.5d0
-    SOFT_TOUCH level_shift
+    if (.not.do_mom) then
+       double precision :: level_shift_save
+       level_shift_save = level_shift
+       mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num)
+       do while (Delta_energy_SCF > 0.d0)
+          mo_coef(1:ao_num,1:mo_num) = mo_coef_save
+          if (level_shift <= .1d0) then
+             level_shift = 1.d0
+          else
+             level_shift = level_shift * 3.0d0
+          endif
+          TOUCH mo_coef level_shift
+          mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_num)
+          if(do_mom)then
+             call reorder_mo_max_overlap
+          endif
+          if(frozen_orb_scf)then
+             call reorder_core_orb
+             call initialize_mo_coef_begin_iteration
+          endif
+          TOUCH mo_coef
+          Delta_Energy_SCF = SCF_energy - energy_SCF_previous
+          energy_SCF = SCF_energy
+          if (level_shift-level_shift_save > 40.d0) then
+             level_shift = level_shift_save * 4.d0
+             SOFT_TOUCH level_shift
+             exit
+          endif
+          dim_DIIS=0
+       enddo
+       level_shift = level_shift * 0.5d0
+       SOFT_TOUCH level_shift
+    endif
     energy_SCF_previous = energy_SCF
 
     converged = ( (max_error_DIIS <= threshold_DIIS_nonzero) .and. &
@@ -205,7 +216,7 @@ subroutine Roothaan_Hall_SCF
 
   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.)
+        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
   endif
@@ -228,6 +239,9 @@ subroutine Roothaan_Hall_SCF
     i = j+1
   enddo
 
+  if(do_mom)then
+     call reorder_mo_max_overlap
+  endif
 
   call save_mos