From a64d02ab427a2b886e19390a33c4fa395f69ed9c Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 11 Oct 2023 15:45:51 +0200 Subject: [PATCH] trying to work on natorb --- bin/qp_reset | 2 + external/ezfio | 2 +- src/non_hermit_dav/biorthog.irp.f | 6 +- .../lapack_diag_non_hermit.irp.f | 90 ++++++++++++++----- src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f | 8 +- src/tc_scf/routines_rotates.irp.f | 1 + 6 files changed, 81 insertions(+), 28 deletions(-) diff --git a/bin/qp_reset b/bin/qp_reset index d94ab24c..b144c4ce 100755 --- a/bin/qp_reset +++ b/bin/qp_reset @@ -97,6 +97,8 @@ if [[ $dets -eq 1 ]] ; then rm --force -- ${ezfio}/determinants/psi_{det,coef}.gz rm --force -- ${ezfio}/determinants/n_det_qp_edit rm --force -- ${ezfio}/determinants/psi_{det,coef}_qp_edit.gz + rm --force -- ${ezfio}/tc_bi_ortho/psi_{l,r}_coef_bi_ortho.gz + fi if [[ $mos -eq 1 ]] ; then diff --git a/external/ezfio b/external/ezfio index d5805497..dba01c4f 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93 +Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3 diff --git a/src/non_hermit_dav/biorthog.irp.f b/src/non_hermit_dav/biorthog.irp.f index 78fddf54..da33f75a 100644 --- a/src/non_hermit_dav/biorthog.irp.f +++ b/src/non_hermit_dav/biorthog.irp.f @@ -331,7 +331,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei !thr = 100d0 thr = Im_thresh_tcscf do i = 1, n - !print*, 'Re(i) + Im(i)', WR(i), WI(i) + print*, 'Re(i) + Im(i)', WR(i), WI(i) if(dabs(WI(i)) .lt. thr) then n_good += 1 else @@ -405,7 +405,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then - !print *, ' lapack vectors are normalized and bi-orthogonalized' + print *, ' lapack vectors are normalized and bi-orthogonalized' deallocate(S) return @@ -422,7 +422,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei else - !print *, ' lapack vectors are not normalized neither bi-orthogonalized' + print *, ' lapack vectors are not normalized neither bi-orthogonalized' ! --- diff --git a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f index 0d652af4..6e5719c1 100644 --- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -1857,7 +1857,7 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ integer :: i, j double precision, allocatable :: SS(:,:) - !print *, ' check bi-orthogonality' + print *, ' check bi-orthogonality' ! --- @@ -1865,10 +1865,10 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ , Vl, size(Vl, 1), Vr, size(Vr, 1) & , 0.d0, S, size(S, 1) ) - !print *, ' overlap matrix:' - !do i = 1, m - ! write(*,'(1000(F16.10,X))') S(i,:) - !enddo + print *, ' overlap matrix:' + do i = 1, m + write(*,'(1000(F16.10,X))') S(i,:) + enddo accu_d = 0.d0 accu_nd = 0.d0 @@ -1883,8 +1883,8 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ enddo accu_nd = dsqrt(accu_nd) / dble(m) - !print *, ' accu_nd = ', accu_nd - !print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + print *, ' accu_nd = ', accu_nd + print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) ! --- @@ -1987,11 +1987,11 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) enddo enddo - !do i = 1, n - ! if(deg_num(i) .gt. 1) then - ! print *, ' degen on', i, deg_num(i), e0(i) - ! endif - !enddo + do i = 1, n + if(deg_num(i) .gt. 1) then + print *, ' degen on', i, deg_num(i), e0(i) + endif + enddo ! --- @@ -2010,7 +2010,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) ! --- - call impose_orthog_svd(n, m, L) +! call impose_orthog_svd(n, m, L) call impose_orthog_svd(n, m, R) !call impose_orthog_GramSchmidt(n, m, L) !call impose_orthog_GramSchmidt(n, m, R) @@ -2030,7 +2030,8 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) !call bi_ortho_s_inv_half(m, L, R, S_inv_half) !deallocate(S, S_inv_half) - call impose_biorthog_svd(n, m, L, R) +! call impose_biorthog_svd(n, m, L, R) +! call impose_biorthog_inverse(n, m, L, R) !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R) @@ -2045,6 +2046,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) endif enddo + call impose_biorthog_inverse(n, n, L0, R0) end subroutine impose_biorthog_degen_eigvec @@ -2420,10 +2422,10 @@ subroutine impose_biorthog_svd(n, m, L, R) , L, size(L, 1), R, size(R, 1) & , 0.d0, S, size(S, 1) ) - !print *, ' overlap bef SVD: ' - !do i = 1, m - ! write(*, '(1000(F16.10,X))') S(i,:) - !enddo + print *, ' overlap bef SVD: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo ! --- @@ -2495,10 +2497,10 @@ subroutine impose_biorthog_svd(n, m, L, R) , L, size(L, 1), R, size(R, 1) & , 0.d0, S, size(S, 1) ) - !print *, ' overlap aft SVD: ' - !do i = 1, m - ! write(*, '(1000(F16.10,X))') S(i,:) - !enddo + print *, ' overlap aft SVD: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo deallocate(S) @@ -2506,6 +2508,50 @@ subroutine impose_biorthog_svd(n, m, L, R) end subroutine impose_biorthog_svd +subroutine impose_biorthog_inverse(n, m, L, R) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(inout) :: L(n,m) + double precision, intent(in) :: R(n,m) + double precision, allocatable :: Lt(:,:),S(:,:) + integer :: i,j + allocate(Lt(m,n)) + allocate(S(m,m)) + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + print *, ' overlap bef SVD: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo + + call get_pseudo_inverse(R,n,n,m,Lt,m,1.d-6) + do i = 1, m + do j = 1, n + L(j,i) = Lt(i,j) + enddo + enddo + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + print *, ' overlap aft SVD: ' + do i = 1, m + write(*, '(1000(F16.10,X))') S(i,:) + enddo + + deallocate(S,Lt) + + +end subroutine impose_biorthog_svd + + ! --- subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr) diff --git a/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f index 47ade8df..6b239cfc 100644 --- a/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f +++ b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f @@ -22,6 +22,7 @@ program tc_natorb_bi_ortho call print_energy_and_mos() call save_tc_natorb() + call print_angles_tc() !call minimize_tc_orb_angles() end @@ -35,9 +36,12 @@ subroutine save_tc_natorb() print*,'Saving the natorbs ' provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao + mo_l_coef = natorb_tc_leigvec_ao + mo_r_coef = natorb_tc_reigvec_ao + touch mo_l_coef mo_r_coef - call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao) - call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao) + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) call save_ref_determinant_nstates_1() call ezfio_set_determinants_read_wf(.False.) diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f index 588382b5..cc825429 100644 --- a/src/tc_scf/routines_rotates.irp.f +++ b/src/tc_scf/routines_rotates.irp.f @@ -402,6 +402,7 @@ subroutine print_energy_and_mos(good_angles) print *, ' TC energy = ', TC_HF_energy print *, ' TC SCF energy gradient = ', grad_non_hermit print *, ' Max angle Left/right = ', max_angle_left_right + call print_angles_tc() if(max_angle_left_right .lt. thresh_lr_angle) then print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'