From c7dd091b632ba8f41b259788d4344728bec514f0 Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 29 Oct 2022 18:53:25 +0200 Subject: [PATCH] split tc_prop.irp.f into tc_natorb.irp.f --- src/tc_bi_ortho/tc_natorb.irp.f | 187 ++++++++++++++++++++++++++++++ src/tc_bi_ortho/tc_prop.irp.f | 186 ----------------------------- src/tc_bi_ortho/test_natorb.irp.f | 2 +- src/tc_scf/routines_rotates.irp.f | 23 ++-- src/utils/block_diag_degen.irp.f | 4 +- 5 files changed, 203 insertions(+), 199 deletions(-) create mode 100644 src/tc_bi_ortho/tc_natorb.irp.f diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f new file mode 100644 index 00000000..71fc6716 --- /dev/null +++ b/src/tc_bi_ortho/tc_natorb.irp.f @@ -0,0 +1,187 @@ + + BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)] + implicit none + BEGIN_DOC + ! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals) + ! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals) + ! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !! + END_DOC + double precision, allocatable :: dm_tmp(:,:) + integer :: i,j,k,n_real + allocate( dm_tmp(mo_num,mo_num)) + dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1) + print*,'dm_tmp' + do i = 1, mo_num + write(*,'(100(F16.10,X))')-dm_tmp(:,i) + enddo + call non_hrmt_bieig( mo_num, dm_tmp& + , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& + , n_real, natorb_tc_eigval ) + double precision :: accu + accu = 0.d0 + do i = 1, n_real + print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i) + accu += -natorb_tc_eigval(i) + enddo + print*,'accu = ',accu + dm_tmp = 0.d0 + do i = 1, n_real + accu = 0.d0 + do k = 1, mo_num + accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i) + enddo + accu = 1.d0/dsqrt(dabs(accu)) + natorb_tc_reigvec_mo(:,i) *= accu + natorb_tc_leigvec_mo(:,i) *= accu + do j = 1, n_real + do k = 1, mo_num + dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j) + enddo + enddo + enddo + double precision :: accu_d, accu_nd + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, mo_num + accu_d += dm_tmp(i,i) + ! write(*,'(100(F16.10,X))')dm_tmp(:,i) + do j = 1, mo_num + if(i==j)cycle + accu_nd += dabs(dm_tmp(j,i)) + enddo + enddo + print*,'Trace of the overlap between TC natural orbitals ',accu_d + print*,'L1 norm of extra diagonal elements of overlap matrix ',accu_nd + + + END_PROVIDER + + BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)] + implicit none + integer ::i,j,k + print*,'Diagonal elements of the Fock matrix before ' + do i = 1, mo_num + write(*,*)i,Fock_matrix_tc_mo_tot(i,i) + enddo + double precision, allocatable :: fock_diag(:) + allocate(fock_diag(mo_num)) + fock_diag = 0.d0 + do i = 1, mo_num + fock_diag(i) = 0.d0 + do j = 1, mo_num + do k = 1, mo_num + fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i) + enddo + enddo + enddo + integer, allocatable :: iorder(:) + allocate(iorder(mo_num)) + do i = 1, mo_num + iorder(i) = i + enddo + call dsort(fock_diag,iorder,mo_num) + print*,'Diagonal elements of the Fock matrix after ' + do i = 1, mo_num + write(*,*)i,fock_diag(i) + enddo + do i = 1, mo_num + fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i)) + do j = 1, mo_num + fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i)) + fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i)) + enddo + enddo + + END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)] + &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)] + &BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ] + + BEGIN_DOC + ! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP + ! + ! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo + END_DOC + + implicit none + integer :: i, j, k, q, p + double precision :: accu, accu_d + double precision, allocatable :: tmp(:,:) + + + ! ! MO_R x R + call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & + , mo_r_coef, size(mo_r_coef, 1) & + , fock_diag_sorted_r_natorb, size(fock_diag_sorted_r_natorb, 1) & + , 0.d0, natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) ) + ! + ! MO_L x L + call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & + , mo_l_coef, size(mo_l_coef, 1) & + , fock_diag_sorted_l_natorb, size(fock_diag_sorted_l_natorb, 1) & + , 0.d0, natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1) ) + + + allocate( tmp(mo_num,ao_num) ) + + ! tmp <-- L.T x S_ao + call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 & + , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + ! S <-- tmp x R + call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 & + , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) & + , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) ) + + deallocate( tmp ) + + ! --- + double precision :: norm + do i = 1, mo_num + norm = 1.d0/dsqrt(dabs(overlap_natorb_tc_eigvec_ao(i,i))) + do j = 1, mo_num + natorb_tc_reigvec_ao(j,i) *= norm + natorb_tc_leigvec_ao(j,i) *= norm + enddo + enddo + + allocate( tmp(mo_num,ao_num) ) + + ! tmp <-- L.T x S_ao + call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 & + , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + ! S <-- tmp x R + call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 & + , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) & + , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) ) + + + + deallocate( tmp ) + + accu_d = 0.d0 + accu = 0.d0 + do i = 1, mo_num + accu_d += overlap_natorb_tc_eigvec_ao(i,i) + do j = 1, mo_num + if(i==j)cycle + accu += dabs(overlap_natorb_tc_eigvec_ao(j,i)) + enddo + enddo + print*,'Trace of the overlap_natorb_tc_eigvec_ao = ',accu_d + print*,'mo_num = ',mo_num + print*,'L1 norm of extra diagonal elements of overlap matrix ',accu + accu = accu / dble(mo_num**2) + + END_PROVIDER + diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f index bc755ada..c7f6c986 100644 --- a/src/tc_bi_ortho/tc_prop.irp.f +++ b/src/tc_bi_ortho/tc_prop.irp.f @@ -49,192 +49,6 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state END_PROVIDER - BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)] - &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)] - &BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)] - implicit none - BEGIN_DOC - ! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals) - ! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals) - ! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !! - END_DOC - double precision, allocatable :: dm_tmp(:,:) - integer :: i,j,k,n_real - allocate( dm_tmp(mo_num,mo_num)) - dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1) - print*,'dm_tmp' - do i = 1, mo_num - write(*,'(100(F16.10,X))')-dm_tmp(:,i) - enddo - call non_hrmt_bieig( mo_num, dm_tmp& - , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& - , n_real, natorb_tc_eigval ) - double precision :: accu - accu = 0.d0 - do i = 1, n_real - print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i) - accu += -natorb_tc_eigval(i) - enddo - print*,'accu = ',accu - dm_tmp = 0.d0 - do i = 1, n_real - accu = 0.d0 - do k = 1, mo_num - accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i) - enddo - accu = 1.d0/dsqrt(dabs(accu)) - natorb_tc_reigvec_mo(:,i) *= accu - natorb_tc_leigvec_mo(:,i) *= accu - do j = 1, n_real - do k = 1, mo_num - dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j) - enddo - enddo - enddo - double precision :: accu_d, accu_nd - accu_d = 0.d0 - accu_nd = 0.d0 - do i = 1, mo_num - accu_d += dm_tmp(i,i) - ! write(*,'(100(F16.10,X))')dm_tmp(:,i) - do j = 1, mo_num - if(i==j)cycle - accu_nd += dabs(dm_tmp(j,i)) - enddo - enddo - print*,'Trace of the overlap between TC natural orbitals ',accu_d - print*,'L1 norm of extra diagonal elements of overlap matrix ',accu_nd - - - END_PROVIDER - - BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)] - &BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)] - &BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)] - implicit none - integer ::i,j,k - print*,'Diagonal elements of the Fock matrix before ' - do i = 1, mo_num - write(*,*)i,Fock_matrix_tc_mo_tot(i,i) - enddo - double precision, allocatable :: fock_diag(:) - allocate(fock_diag(mo_num)) - fock_diag = 0.d0 - do i = 1, mo_num - fock_diag(i) = 0.d0 - do j = 1, mo_num - do k = 1, mo_num - fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i) - enddo - enddo - enddo - integer, allocatable :: iorder(:) - allocate(iorder(mo_num)) - do i = 1, mo_num - iorder(i) = i - enddo - call dsort(fock_diag,iorder,mo_num) - print*,'Diagonal elements of the Fock matrix after ' - do i = 1, mo_num - write(*,*)i,fock_diag(i) - enddo - do i = 1, mo_num - fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i)) - do j = 1, mo_num - fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i)) - fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i)) - enddo - enddo - - END_PROVIDER - - - - BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)] - &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)] - &BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ] - - BEGIN_DOC - ! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP - ! - ! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo - END_DOC - - implicit none - integer :: i, j, k, q, p - double precision :: accu, accu_d - double precision, allocatable :: tmp(:,:) - - - ! ! MO_R x R - call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & - , mo_r_coef, size(mo_r_coef, 1) & - , fock_diag_sorted_r_natorb, size(fock_diag_sorted_r_natorb, 1) & - , 0.d0, natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) ) - ! - ! MO_L x L - call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & - , mo_l_coef, size(mo_l_coef, 1) & - , fock_diag_sorted_l_natorb, size(fock_diag_sorted_l_natorb, 1) & - , 0.d0, natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1) ) - - - allocate( tmp(mo_num,ao_num) ) - - ! tmp <-- L.T x S_ao - call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 & - , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) & - , 0.d0, tmp, size(tmp, 1) ) - - ! S <-- tmp x R - call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 & - , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) & - , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) ) - - deallocate( tmp ) - - ! --- - double precision :: norm - do i = 1, mo_num - norm = 1.d0/dsqrt(dabs(overlap_natorb_tc_eigvec_ao(i,i))) - do j = 1, mo_num - natorb_tc_reigvec_ao(j,i) *= norm - natorb_tc_leigvec_ao(j,i) *= norm - enddo - enddo - - allocate( tmp(mo_num,ao_num) ) - - ! tmp <-- L.T x S_ao - call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 & - , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) & - , 0.d0, tmp, size(tmp, 1) ) - - ! S <-- tmp x R - call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 & - , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) & - , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) ) - - - - deallocate( tmp ) - - accu_d = 0.d0 - accu = 0.d0 - do i = 1, mo_num - accu_d += overlap_natorb_tc_eigvec_ao(i,i) - do j = 1, mo_num - if(i==j)cycle - accu += dabs(overlap_natorb_tc_eigvec_ao(j,i)) - enddo - enddo - print*,'Trace of the overlap_natorb_tc_eigvec_ao = ',accu_d - print*,'mo_num = ',mo_num - print*,'L1 norm of extra diagonal elements of overlap matrix ',accu - accu = accu / dble(mo_num**2) - - END_PROVIDER - BEGIN_PROVIDER [double precision, tc_bi_ortho_dipole, (3,N_states)] implicit none integer :: i,j,istate,m diff --git a/src/tc_bi_ortho/test_natorb.irp.f b/src/tc_bi_ortho/test_natorb.irp.f index f99f3cea..54c9a827 100644 --- a/src/tc_bi_ortho/test_natorb.irp.f +++ b/src/tc_bi_ortho/test_natorb.irp.f @@ -32,7 +32,7 @@ subroutine routine thr_deg = 1.d-6 mat_ref = -one_e_dm_mo print*,'diagonalization by block' - call diagonalize_dm_per_fock_degen(fock_diag,mat_ref,mo_num,thr_deg,leigvec,reigvec,eigval) + call diag_mat_per_fock_degen(fock_diag,mat_ref,mo_num,thr_deg,leigvec,reigvec,eigval) call non_hrmt_bieig( mo_num, mat_ref& , leigvec_ref, reigvec_ref& , n_real, eigval_ref) diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f index f829066e..51c218fa 100644 --- a/src/tc_scf/routines_rotates.irp.f +++ b/src/tc_scf/routines_rotates.irp.f @@ -52,31 +52,32 @@ subroutine routine_save_rotated_mos(thr_deg,good_angles) enddo double precision, allocatable :: fock_diag(:),s_mat(:,:) integer, allocatable :: list_degen(:,:) - allocate(list_degen(2,mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num)) + allocate(list_degen(mo_num,0:mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num)) do i = 1, mo_num fock_diag(i) = Fock_matrix_tc_mo_tot(i,i) enddo ! compute the overlap between the left and rescaled right call build_s_matrix(ao_num,mo_num,mo_r_coef_new,mo_r_coef_new,ao_overlap,s_mat) - call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list) +! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list) + call give_degen_full_list(fock_diag,mo_num,thr_deg,list_degen,n_degen_list) print*,'fock_matrix_mo' do i = 1, mo_num print*,i,fock_diag(i),angle_left_right(i) enddo do i = 1, n_degen_list - ifirst = list_degen(1,i) - ilast = list_degen(2,i) - n_degen = ilast - ifirst +1 - print*,'ifirst,n_degen = ',ifirst,n_degen +! ifirst = list_degen(1,i) +! ilast = list_degen(2,i) +! n_degen = ilast - ifirst +1 + n_degen = list_degen(i,0) double precision, allocatable :: stmp(:,:),T(:,:),Snew(:,:),smat2(:,:) double precision, allocatable :: mo_l_coef_tmp(:,:),mo_r_coef_tmp(:,:),mo_l_coef_new(:,:) allocate(stmp(n_degen,n_degen),smat2(n_degen,n_degen)) allocate(mo_r_coef_tmp(ao_num,n_degen),mo_l_coef_tmp(ao_num,n_degen),mo_l_coef_new(ao_num,n_degen)) allocate(T(n_degen,n_degen),Snew(n_degen,n_degen)) do j = 1, n_degen - mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,j+ifirst-1) - mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,j+ifirst-1) + mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j)) + mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j)) enddo ! Orthogonalization of right functions print*,'Orthogonalization of RIGHT functions' @@ -138,8 +139,10 @@ subroutine routine_save_rotated_mos(thr_deg,good_angles) ! write(*,'(100(F16.10,X))')stmp(:,j) !enddo do j = 1, n_degen - mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j) - mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j) +! mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j) +! mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j) + mo_l_coef_good(1:ao_num,list_degen(i,j)) = mo_l_coef_new(1:ao_num,j) + mo_r_coef_good(1:ao_num,list_degen(i,j)) = mo_r_coef_tmp(1:ao_num,j) enddo deallocate(stmp,smat2) deallocate(mo_r_coef_tmp,mo_l_coef_tmp,mo_l_coef_new) diff --git a/src/utils/block_diag_degen.irp.f b/src/utils/block_diag_degen.irp.f index 65405656..9884bf21 100644 --- a/src/utils/block_diag_degen.irp.f +++ b/src/utils/block_diag_degen.irp.f @@ -1,11 +1,11 @@ -subroutine diagonalize_dm_per_fock_degen(fock_diag,mat_ref,n,thr_deg,leigvec,reigvec,eigval) +subroutine diag_mat_per_fock_degen(fock_diag,mat_ref,n,thr_deg,leigvec,reigvec,eigval) implicit none integer, intent(in) :: n double precision, intent(in) :: fock_diag(n),mat_ref(n,n),thr_deg double precision, intent(out):: leigvec(n,n),reigvec(n,n),eigval(n) BEGIN_DOC -! subroutine that diagonalizes a matrix mat_ref BY BLOCK +! subroutine that diagonalizes a matrix mat_ref BY BLOCK ! ! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag" !