From 61bcd8c10328bfc0ab94909dc2be5b3dfb3f805c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 26 Nov 2015 15:00:02 +0100 Subject: [PATCH] Eigenvalues are now printed with the correct sign --- plugins/Hartree_Fock/damping_SCF.irp.f | 2 +- plugins/Hartree_Fock/huckel.irp.f | 2 +- src/Determinants/density_matrix.irp.f | 4 +-- src/MOGuess/H_CORE_guess.irp.f | 2 +- src/MOGuess/guess_overlap.irp.f | 4 +-- src/MOGuess/h_core_guess_routine.irp.f | 2 +- src/MO_Basis/utils.irp.f | 37 +++++++++++++++++++------- 7 files changed, 36 insertions(+), 17 deletions(-) diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f index c7c005c9..d77c91c5 100644 --- a/plugins/Hartree_Fock/damping_SCF.irp.f +++ b/plugins/Hartree_Fock/damping_SCF.irp.f @@ -119,7 +119,7 @@ subroutine damping_SCF write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '====' write(output_hartree_fock,*) - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label) + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') call ezfio_set_hartree_fock_energy(E_min) diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f index 6139ac46..8f61f0cf 100644 --- a/plugins/Hartree_Fock/huckel.irp.f +++ b/plugins/Hartree_Fock/huckel.irp.f @@ -13,7 +13,7 @@ subroutine huckel_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label) + size(mo_mono_elec_integral,2),label,1) TOUCH mo_coef c = 0.5d0 * 1.75d0 diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 5f087642..9aeb658e 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -185,9 +185,9 @@ subroutine set_natural_mos allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2))) ! Negation to have the occupied MOs first after the diagonalization - tmp = -one_body_dm_mo + tmp = one_body_dm_mo label = "Natural" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,-1) deallocate(tmp) end diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f index 1893c08b..b65fe07d 100644 --- a/src/MOGuess/H_CORE_guess.irp.f +++ b/src/MOGuess/H_CORE_guess.irp.f @@ -10,7 +10,7 @@ program H_CORE_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label) + size(mo_mono_elec_integral,2),label,1) print *, 'save mos' call save_mos diff --git a/src/MOGuess/guess_overlap.irp.f b/src/MOGuess/guess_overlap.irp.f index 768724c6..cf222507 100644 --- a/src/MOGuess/guess_overlap.irp.f +++ b/src/MOGuess/guess_overlap.irp.f @@ -7,8 +7,8 @@ program guess_mimi mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef label = "Guess" - call mo_as_eigvectors_of_mo_matrix(-ao_overlap, & + call mo_as_eigvectors_of_mo_matrix(ao_overlap, & size(ao_overlap,1), & - size(ao_overlap,2),label) + size(ao_overlap,2),label,-1) call save_mos end diff --git a/src/MOGuess/h_core_guess_routine.irp.f b/src/MOGuess/h_core_guess_routine.irp.f index 566592ba..605c7c8a 100644 --- a/src/MOGuess/h_core_guess_routine.irp.f +++ b/src/MOGuess/h_core_guess_routine.irp.f @@ -9,7 +9,7 @@ subroutine hcore_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label) + size(mo_mono_elec_integral,2),label,1) print *, 'save mos' call save_mos SOFT_TOUCH mo_coef mo_label diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 6f96ab93..3f4a260a 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -44,13 +44,14 @@ subroutine save_mos_truncated(n) end -subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label) +subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) implicit none - integer,intent(in) :: n,m + integer,intent(in) :: n,m, sign character*(64), intent(in) :: label double precision, intent(in) :: matrix(n,m) - double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:) + integer :: i,j + double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R call write_time(output_mo_basis) @@ -58,30 +59,48 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label) print *, irp_here, ': Error : m/= mo_tot_num' stop 1 endif - allocate(R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m)) + allocate(A(n,m),R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m)) + if (sign == -1) then + do j=1,m + do i=1,n + A(i,j) = -matrix(i,j) + enddo + enddo + else + do j=1,m + do i=1,n + A(i,j) = matrix(i,j) + enddo + enddo + endif mo_coef_new = mo_coef - call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2)) - integer :: i + call lapack_diag(eigvalues,R,A,n,m) write (output_mo_basis,'(A)'), 'MOs are now **'//trim(label)//'**' write (output_mo_basis,'(A)'), '' write (output_mo_basis,'(A)'), 'Eigenvalues' write (output_mo_basis,'(A)'), '-----------' write (output_mo_basis,'(A)'), '' write (output_mo_basis,'(A)'), '======== ================' - do i = 1, m - write (output_mo_basis,'(I8,X,F16.10)'), i,eigvalues(i) + if (sign == -1) then + do i=1,m + eigvalues(i) = -eigvalues(i) + enddo + endif + do i=1,m + write (output_mo_basis,'(I8,X,F16.10)'), i,eigvalues(i) enddo write (output_mo_basis,'(A)'), '======== ================' write (output_mo_basis,'(A)'), '' call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1)) - deallocate(mo_coef_new,R,eigvalues) + deallocate(A,mo_coef_new,R,eigvalues) call write_time(output_mo_basis) mo_label = label SOFT_TOUCH mo_coef mo_label end + subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label) implicit none integer,intent(in) :: n,m