PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 04:13:55 +01:00
Code Releases Activity

Merge pull request #347 from antoine-marie/overlap

Browse Source 
add maximum overlap method in hartree fock
This commit is contained in:
Anthony Scemama 2024-09-20 16:09:11 +02:00 committed by GitHub
commit 360ac7b128
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
3 changed files with 151 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/scf_utils/EZFIO.cfg
View File

@ -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

96
src/scf_utils/reorder_mo_max_overlap.irp.f Normal file
View File

@ -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

84
src/scf_utils/roothaan_hall_scf.irp.f
View File

@ -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
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
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 @@ 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