10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-03 20:53:54 +01:00

saving olympe2 modif

This commit is contained in:
AbdAmmar 2024-04-07 00:29:40 +02:00
parent 83ed57312d
commit d872d60e70
6 changed files with 81 additions and 46 deletions

View File

@ -56,10 +56,10 @@
print*,'Average trace of overlap_bi_ortho is different from 1 by ', dabs(accu_d-1.d0) 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*,'And bi orthogonality is off by an average of ',accu_nd
print*,'****************' print*,'****************'
print*,'Overlap matrix betwee mo_l_coef and mo_r_coef ' !print*,'Overlap matrix betwee mo_l_coef and mo_r_coef '
do i = 1, mo_num !do i = 1, mo_num
write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:) ! write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:)
enddo !enddo
endif endif
print*,'Average trace of overlap_bi_ortho (should be 1.)' print*,'Average trace of overlap_bi_ortho (should be 1.)'
print*,'accu_d = ',accu_d print*,'accu_d = ',accu_d

View File

@ -2144,6 +2144,7 @@ subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0)
enddo enddo
!print*,' accu_nd after = ', accu_nd !print*,' accu_nd after = ', accu_nd
if(accu_nd .gt. 1d-12) then if(accu_nd .gt. 1d-12) then
print*, ' accu_nd =', accu_nd
print*, ' your strategy for degenerates orbitals failed !' print*, ' your strategy for degenerates orbitals failed !'
print*, m, 'deg on', i print*, m, 'deg on', i
stop stop

View File

@ -20,7 +20,7 @@ program minimize_tc_angles
! TODO ! TODO
! check if rotations of orbitals affect the TC energy ! check if rotations of orbitals affect the TC energy
! and refuse the step ! and refuse the step
call minimize_tc_orb_angles call minimize_tc_orb_angles()
end end

View File

@ -40,9 +40,6 @@ subroutine LTxSxR(n, m, L, S, R, C)
end subroutine LTxR end subroutine LTxR
! ---
! --- ! ---
subroutine minimize_tc_orb_angles() 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 :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:) 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 E_old = TC_HF_energy
allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num)) allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num))
mo_r_coef_old = mo_r_coef mo_r_coef_old = mo_r_coef
@ -111,7 +111,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
good_angles = .False. 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 *, ' ***************************************'
print *, ' ***************************************' print *, ' ***************************************'
@ -123,7 +123,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
mo_r_coef_good = mo_r_coef mo_r_coef_good = mo_r_coef
mo_l_coef_good = mo_l_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 mo_r_coef_new = mo_r_coef
do i = 1, mo_num do i = 1, mo_num
norm = 1.d0/dsqrt(overlap_mo_r(i,i)) 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 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)
if(n_core_orb.ne.0)then 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 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 endif
print *, ' fock_matrix_mo' print *, ' fock_matrix_mo'
do i = 1, mo_num do i = 1, mo_num
print *, i, fock_diag(i), angle_left_right(i) 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 = ilast - ifirst +1
n_degen = list_degen(i,0) 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 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(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(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)) allocate(T(n_degen,n_degen), Snew(n_degen,n_degen))
print*,'Right orbitals before' print*,'Right orbitals before'
do j = 1, n_degen do j = 1, n_degen
write(*,'(100(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j)) write(*,'(1000(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j))
enddo enddo
print*,'Left orbitals before' print*,'Left orbitals before'
do j = 1, n_degen do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_l_coef(1:ao_num,list_degen(i,j)) write(*,'(1000(F16.10,X))') mo_l_coef(1:ao_num,list_degen(i,j))
enddo enddo
if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then
integer :: i_list, j_list integer :: i_list, j_list
i_list = list_degen(i,1) i_list = list_degen(i,1)
j_list = list_degen(i,2) j_list = list_degen(i,2)
print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(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 = ',i_list
print*,'i_list = ',j_list print*,'i_list = ',j_list
print*,'Swapping left/right orbitals' print*,'Swapping left/right orbitals'
call print_strong_overlap(i_list, j_list) 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,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_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,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) mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list)
else else
do j = 1, n_degen do j = 1, n_degen
print*,'i_list = ',list_degen(i,j) 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_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)) mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
enddo enddo
endif endif
print*,'Right orbitals ' print*,'Right orbitals '
do j = 1, n_degen do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_r_coef_tmp(1:ao_num,j) write(*,'(1000(F16.10,X))') mo_r_coef_tmp(1:ao_num,j)
enddo enddo
print*,'Left orbitals ' print*,'Left orbitals '
do j = 1, n_degen do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_l_coef_tmp(1:ao_num,j) write(*,'(100(F16.10,X))') mo_l_coef_tmp(1:ao_num,j)
enddo enddo
! Orthogonalization of right functions ! Orthogonalization of right functions
print *, ' Orthogonalization of RIGHT functions' print *, ' Orthogonalization of RIGHT functions'
print *, ' ------------------------------------' print *, ' ------------------------------------'

View File

@ -5,7 +5,8 @@ program print_detweights
read_wf = .True. read_wf = .True.
touch read_wf touch read_wf
call main() call print_exc()
!call main()
end end
@ -41,6 +42,7 @@ subroutine main()
do i = 1, N_det do i = 1, N_det
deg_sorted(i) = deg(ii(i)) deg_sorted(i) = deg(ii(i))
print *, deg_sorted(i), c(i)
enddo enddo
print *, ' saving psi' print *, ' saving psi'
@ -52,7 +54,7 @@ subroutine main()
print *, ' Error opening file!' print *, ' Error opening file!'
stop stop
endif endif
write(10) N_det write(10) N_det
write(10) deg_sorted write(10) deg_sorted
write(10) c write(10) c
@ -63,4 +65,33 @@ subroutine main()
end 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

View File

@ -191,7 +191,7 @@ subroutine give_degen_full_list(A, n, thr, list_degen, n_degen_list)
list_degen(n_degen_list,1) = i list_degen(n_degen_list,1) = i
icount = 1 icount = 1
do j = i+1, n 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. is_ok(j) = .False.
icount += 1 icount += 1
list_degen(n_degen_list,icount) = j list_degen(n_degen_list,icount) = j