From b68fe9feb24413c4e85a1458430b44ad4b3440fa Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Wed, 18 Sep 2024 09:43:11 +0200 Subject: [PATCH] add maximum overlap method in hartree fock --- src/scf_utils/EZFIO.cfg | 6 ++ src/scf_utils/reorder_mo_max_overlap.irp.f | 88 +++++++++++++++++++++ src/scf_utils/roothaan_hall_scf.irp.f | 90 +++++++++++++--------- 3 files changed, 146 insertions(+), 38 deletions(-) create mode 100644 src/scf_utils/reorder_mo_max_overlap.irp.f diff --git a/src/scf_utils/EZFIO.cfg b/src/scf_utils/EZFIO.cfg index 7b950b14..1ca97217 100644 --- a/src/scf_utils/EZFIO.cfg +++ b/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 diff --git a/src/scf_utils/reorder_mo_max_overlap.irp.f b/src/scf_utils/reorder_mo_max_overlap.irp.f new file mode 100644 index 00000000..bdd0e252 --- /dev/null +++ b/src/scf_utils/reorder_mo_max_overlap.irp.f @@ -0,0 +1,88 @@ + 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(:,:) + allocate(overlap(mo_num,mo_num),proj(mo_num),iorder(mo_num),mo_coef_tmp(ao_num,mo_num)) + + overlap(:,:) = 0d0 + mo_coef_tmp(:,:) = 0d0 + proj(:) = 0d0 + iorder(:) = 0d0 + + ! this loop compute the overlap between the initial and the current MOS + do i = 1, mo_num ! old mo + do j = 1, mo_num ! curent mo + do k = 1, ao_num + do l = 1, ao_num + overlap(i,j) += mo_coef_begin_iteration(k,i)* ao_overlap(k,l) * mo_coef(l,j) + enddo + enddo + enddo + enddo + + ! 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 + proj_alpha(:) = 0d0 + iorder_alpha(:) = 0d0 + do i = 1, mo_num ! old mo + do j = 1, elec_alpha_num ! curent mo + do k = 1, ao_num + do l = 1, ao_num + overlap_alpha(i,j) += mo_coef_begin_iteration(k,i) * ao_overlap(k,l) * mo_coef(l,j) + enddo + enddo + enddo + enddo + 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 + endif + + deallocate(overlap, proj, iorder, mo_coef_tmp) + + + end + diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index e0fe5319..bbf77bea 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -51,6 +51,11 @@ END_DOC ! 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 @@ END_DOC 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 endif TOUCH MO_coef - ! Calculate error vectors max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO)) @@ -106,41 +112,46 @@ END_DOC energy_SCF = SCF_energy Delta_Energy_SCF = energy_SCF - energy_SCF_previous - if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) 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 + ! if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) 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 @@ 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.) + 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 @@ END_DOC i = j+1 enddo + if(do_mom)then + call reorder_mo_max_overlap + endif call save_mos