diff --git a/plugins/local/bi_ortho_mos/overlap.irp.f b/plugins/local/bi_ortho_mos/overlap.irp.f index ff5d5c84..7f07929f 100644 --- a/plugins/local/bi_ortho_mos/overlap.irp.f +++ b/plugins/local/bi_ortho_mos/overlap.irp.f @@ -56,10 +56,10 @@ print*,'Average trace of overlap_bi_ortho is different from 1 by ', dabs(accu_d-1.d0) print*,'And bi orthogonality is off by an average of ',accu_nd print*,'****************' - print*,'Overlap matrix betwee mo_l_coef and mo_r_coef ' - do i = 1, mo_num - write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:) - enddo + !print*,'Overlap matrix betwee mo_l_coef and mo_r_coef ' + !do i = 1, mo_num + ! write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:) + !enddo endif print*,'Average trace of overlap_bi_ortho (should be 1.)' print*,'accu_d = ',accu_d diff --git a/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f b/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f index cb38347e..4d4bc047 100644 --- a/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -2144,6 +2144,7 @@ subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0) enddo !print*,' accu_nd after = ', accu_nd if(accu_nd .gt. 1d-12) then + print*, ' accu_nd =', accu_nd print*, ' your strategy for degenerates orbitals failed !' print*, m, 'deg on', i stop diff --git a/plugins/local/tc_scf/minimize_tc_angles.irp.f b/plugins/local/tc_scf/minimize_tc_angles.irp.f index c7752930..e5f6cf87 100644 --- a/plugins/local/tc_scf/minimize_tc_angles.irp.f +++ b/plugins/local/tc_scf/minimize_tc_angles.irp.f @@ -20,7 +20,7 @@ program minimize_tc_angles ! TODO ! check if rotations of orbitals affect the TC energy ! and refuse the step - call minimize_tc_orb_angles + call minimize_tc_orb_angles() end diff --git a/plugins/local/tc_scf/routines_rotates.irp.f b/plugins/local/tc_scf/routines_rotates.irp.f index c42e846e..2c5510f2 100644 --- a/plugins/local/tc_scf/routines_rotates.irp.f +++ b/plugins/local/tc_scf/routines_rotates.irp.f @@ -40,9 +40,6 @@ subroutine LTxSxR(n, m, L, S, R, C) end subroutine LTxR -! --- - - ! --- subroutine minimize_tc_orb_angles() @@ -103,7 +100,10 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:) double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:) - E_thr = 1d-04 + PROVIDE TC_HF_energy + PROVIDE mo_r_coef mo_l_coef + + E_thr = 1d-07 E_old = TC_HF_energy allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num)) mo_r_coef_old = mo_r_coef @@ -111,7 +111,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) good_angles = .False. - allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num)) + allocate(mo_l_coef_good(ao_num,mo_num), mo_r_coef_good(ao_num,mo_num)) print *, ' ***************************************' print *, ' ***************************************' @@ -123,7 +123,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) mo_r_coef_good = mo_r_coef mo_l_coef_good = mo_l_coef - allocate(mo_r_coef_new(ao_num, mo_num)) + allocate(mo_r_coef_new(ao_num,mo_num)) mo_r_coef_new = mo_r_coef do i = 1, mo_num norm = 1.d0/dsqrt(overlap_mo_r(i,i)) @@ -141,10 +141,11 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) 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) if(n_core_orb.ne.0)then - call give_degen_full_listcore(fock_diag, mo_num, list_core, n_core_orb, thr_deg, list_degen, n_degen_list) + call give_degen_full_listcore(fock_diag, mo_num, list_core, n_core_orb, thr_deg, list_degen, n_degen_list) else - call give_degen_full_list(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) endif + print *, ' fock_matrix_mo' do i = 1, mo_num print *, i, fock_diag(i), angle_left_right(i) @@ -156,50 +157,52 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) ! n_degen = ilast - ifirst +1 n_degen = list_degen(i,0) - if(n_degen .ge. 1000)n_degen = 1 ! convention for core orbitals + if(n_degen .ge. 1000) n_degen = 1 ! convention for core orbitals if(n_degen .eq. 1) cycle + print*, ' working on orbital', i + print*, ' multiplicity =', n_degen 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)) print*,'Right orbitals before' - do j = 1, n_degen - write(*,'(100(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j)) - enddo + do j = 1, n_degen + write(*,'(1000(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j)) + enddo print*,'Left orbitals before' - do j = 1, n_degen - write(*,'(100(F16.10,X))')mo_l_coef(1:ao_num,list_degen(i,j)) - enddo + do j = 1, n_degen + write(*,'(1000(F16.10,X))') mo_l_coef(1:ao_num,list_degen(i,j)) + enddo if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then - integer :: i_list, j_list - i_list = list_degen(i,1) - j_list = list_degen(i,2) - print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(list_degen(i,2)) - print*,'i_list = ',i_list - print*,'i_list = ',j_list - print*,'Swapping left/right orbitals' - call print_strong_overlap(i_list, j_list) - mo_r_coef_tmp(1:ao_num,1) = mo_r_coef_new(1:ao_num,i_list) - mo_r_coef_tmp(1:ao_num,2) = mo_l_coef(1:ao_num,i_list) - mo_l_coef_tmp(1:ao_num,1) = mo_l_coef(1:ao_num,j_list) - mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list) + integer :: i_list, j_list + i_list = list_degen(i,1) + j_list = list_degen(i,2) + print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(list_degen(i,2)) + print*,'i_list = ',i_list + print*,'i_list = ',j_list + print*,'Swapping left/right orbitals' + call print_strong_overlap(i_list, j_list) + mo_r_coef_tmp(1:ao_num,1) = mo_r_coef_new(1:ao_num,i_list) + mo_r_coef_tmp(1:ao_num,2) = mo_l_coef(1:ao_num,i_list) + mo_l_coef_tmp(1:ao_num,1) = mo_l_coef(1:ao_num,j_list) + mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list) else - do j = 1, n_degen - print*,'i_list = ',list_degen(i,j) - 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 + do j = 1, n_degen + print*,'i_list = ',list_degen(i,j) + 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 endif print*,'Right orbitals ' - do j = 1, n_degen - write(*,'(100(F16.10,X))')mo_r_coef_tmp(1:ao_num,j) - enddo + do j = 1, n_degen + write(*,'(1000(F16.10,X))') mo_r_coef_tmp(1:ao_num,j) + enddo print*,'Left orbitals ' - do j = 1, n_degen - write(*,'(100(F16.10,X))')mo_l_coef_tmp(1:ao_num,j) - enddo + do j = 1, n_degen + write(*,'(100(F16.10,X))') mo_l_coef_tmp(1:ao_num,j) + enddo ! Orthogonalization of right functions print *, ' Orthogonalization of RIGHT functions' print *, ' ------------------------------------' diff --git a/src/tools/print_detweights.irp.f b/src/tools/print_detweights.irp.f index d5b0f2c9..5e5f2bb1 100644 --- a/src/tools/print_detweights.irp.f +++ b/src/tools/print_detweights.irp.f @@ -5,7 +5,8 @@ program print_detweights read_wf = .True. touch read_wf - call main() + call print_exc() + !call main() end @@ -41,6 +42,7 @@ subroutine main() do i = 1, N_det deg_sorted(i) = deg(ii(i)) + print *, deg_sorted(i), c(i) enddo print *, ' saving psi' @@ -52,7 +54,7 @@ subroutine main() print *, ' Error opening file!' stop endif - + write(10) N_det write(10) deg_sorted write(10) c @@ -63,4 +65,33 @@ subroutine main() end +! --- + +subroutine print_exc() + + implicit none + + integer :: i + integer, allocatable :: deg(:) + + PROVIDE N_int + PROVIDE N_det + PROVIDE psi_det + + allocate(deg(N_det)) + + do i = 1, N_det + call get_excitation_degree(psi_det(1,1,1), psi_det(1,1,i), deg(i), N_int) + enddo + + open(unit=10, file="exc.dat", action="write") + write(10,*) N_det + write(10,*) deg + close(10) + + deallocate(deg) + +end + + diff --git a/src/utils/block_diag_degen.irp.f b/src/utils/block_diag_degen.irp.f index 188bfa58..1a9ca8d6 100644 --- a/src/utils/block_diag_degen.irp.f +++ b/src/utils/block_diag_degen.irp.f @@ -191,7 +191,7 @@ subroutine give_degen_full_list(A, n, thr, list_degen, n_degen_list) list_degen(n_degen_list,1) = i icount = 1 do j = i+1, n - if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then + if(dabs(A(i)-A(j)).lt.thr .and. is_ok(j)) then is_ok(j) = .False. icount += 1 list_degen(n_degen_list,icount) = j