From c23ca860173021f5778020050d1ae522404296f0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 7 Jun 2017 23:05:00 +0200 Subject: [PATCH] Introduced DSYEVR --- plugins/Hartree_Fock/EZFIO.cfg | 2 +- plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f | 5 +- plugins/Hartree_Fock/diagonalize_fock.irp.f | 183 ++++++++++--------- 3 files changed, 99 insertions(+), 91 deletions(-) diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index 4b628ac0..35879330 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -2,7 +2,7 @@ type: Threshold doc: Threshold on the magnitude of the smallest eigenvalues of the overlap matrix in the AO basis interface: ezfio,provider,ocaml -default: 1.e-5 +default: 1.e-6 [max_dim_diis] type: integer diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f index d1706330..032b7cd1 100644 --- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f +++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f @@ -98,10 +98,11 @@ END_DOC double precision :: level_shift_save level_shift_save = level_shift do while (Delta_Energy_SCF > 0.d0) - level_shift = level_shift + 0.05d0 + level_shift = level_shift + 0.1d0 MO_coef = eigenvectors_Fock_matrix_MO TOUCH MO_coef level_shift Delta_Energy_SCF = HF_energy - energy_SCF_previous + energy_SCF = HF_energy if (level_shift-level_shift_save > 0.5d0) exit enddo level_shift = level_shift_save @@ -235,7 +236,7 @@ END_DOC stop 'bug in DIIS' endif - if (rcond > 1.d-10) then + if (rcond > 1.d-12) then ! Compute extrapolated Fock matrix diff --git a/plugins/Hartree_Fock/diagonalize_fock.irp.f b/plugins/Hartree_Fock/diagonalize_fock.irp.f index b303b537..2aabf9bd 100644 --- a/plugins/Hartree_Fock/diagonalize_fock.irp.f +++ b/plugins/Hartree_Fock/diagonalize_fock.irp.f @@ -37,93 +37,100 @@ F(jorb,iorb) = 0.d0 enddo enddo - endif - - - - - ! Insert level shift here - do i = elec_beta_num+1, elec_alpha_num - F(i,i) += 0.5d0*level_shift - enddo - - do i = elec_alpha_num+1, mo_tot_num - F(i,i) += level_shift - enddo - - n = mo_tot_num - lwork = 1+6*n + 2*n*n - liwork = 3 + 5*n - - allocate(work(lwork)) - allocate(iwork(liwork) ) - - lwork = -1 - liwork = -1 - - call dsyevd( 'V', 'U', mo_tot_num, F, & - size(F,1), diagonal_Fock_matrix_mo, & - work, lwork, iwork, liwork, info) - - if (info /= 0) then - print *, irp_here//' DSYEVD failed : ', info - stop 1 - endif - lwork = int(work(1)) - liwork = iwork(1) - deallocate(iwork) - deallocate(work) - - allocate(work(lwork)) - allocate(iwork(liwork) ) - call dsyevd( 'V', 'U', mo_tot_num, F, & - size(F,1), diagonal_Fock_matrix_mo, & - work, lwork, iwork, liwork, info) - deallocate(iwork) - - - if (info /= 0) then - call dsyev( 'V', 'L', mo_tot_num, F, & - size(F,1), diagonal_Fock_matrix_mo, & - work, lwork, info) - - if (info /= 0) then - print *, irp_here//' DSYEV failed : ', info - stop 1 - endif - endif - - call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & - mo_coef, size(mo_coef,1), F, size(F,1), & - 0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1)) - deallocate(work, F) - - - -END_PROVIDER + endif + + + + + ! Insert level shift here + do i = elec_beta_num+1, elec_alpha_num + F(i,i) += 0.5d0*level_shift + enddo + + do i = elec_alpha_num+1, mo_tot_num + F(i,i) += level_shift + enddo + + n = mo_tot_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork)) + allocate(iwork(liwork) ) + + lwork = -1 + liwork = -1 + + + integer :: m, ISUPPZ(mo_tot_num) + call dsyevr('V', 'A', 'U', mo_tot_num, F, size(F,1), & + -100.d0, 100.d0, 1, mo_tot_num, 0.d0, & + m, diagonal_Fock_matrix_mo, & + eigenvectors_Fock_matrix_mo, & + size(eigenvectors_Fock_matrix_mo,1), & + isuppz, work, lwork, iwork, liwork, info) + + + + if (info /= 0) then + print *, irp_here//' DSYEVD failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(iwork) + deallocate(work) + + allocate(work(lwork)) + allocate(iwork(liwork) ) + + call dsyevr('V', 'A', 'U', mo_tot_num, F, size(F,1), & + -100.d0, 100.d0, 1, mo_tot_num, 0.d0, & + m, diagonal_Fock_matrix_mo, & + eigenvectors_Fock_matrix_mo, & + size(eigenvectors_Fock_matrix_mo,1), & + isuppz, work, lwork, iwork, liwork, info) + + deallocate(iwork) + + if (info /= 0) then + print *, irp_here//' DSYEV failed : ', info + stop 1 + endif + + F(1:mo_tot_num,1:mo_tot_num) = eigenvectors_Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num) + + call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & + mo_coef, size(mo_coef,1), F, size(F,1), & + 0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1)) + deallocate(work, F) + + -BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] - implicit none - BEGIN_DOC - ! diagonal element of the fock matrix calculated as the sum over all the interactions - ! with all the electrons in the RHF determinant - ! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij - END_DOC - integer :: i,j - double precision :: accu - do j = 1,elec_alpha_num - accu = 0.d0 - do i = 1, elec_alpha_num - accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) - enddo - diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) - enddo - do j = elec_alpha_num+1,mo_tot_num - accu = 0.d0 - do i = 1, elec_alpha_num - accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) - enddo - diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) - enddo - +END_PROVIDER + +BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] + implicit none + BEGIN_DOC + ! diagonal element of the fock matrix calculated as the sum over all the interactions + ! with all the electrons in the RHF determinant + ! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij + END_DOC + integer :: i,j + double precision :: accu + do j = 1,elec_alpha_num + accu = 0.d0 + do i = 1, elec_alpha_num + accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) + enddo + diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) + enddo + do j = elec_alpha_num+1,mo_tot_num + accu = 0.d0 + do i = 1, elec_alpha_num + accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) + enddo + diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) + enddo + END_PROVIDER