mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-09 06:53:38 +01:00
added a diagonalization by block, doing some cleaning
This commit is contained in:
eginer
parent
315f2582f6
commit
861f0f9461
177
src/tc_bi_ortho/tc_natorb.irp.f
Normal file
177
src/tc_bi_ortho/tc_natorb.irp.f
Normal file
@ -0,0 +1,177 @@
|
||||
|
||||
subroutine diagonalize_dm_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
|
||||
!
|
||||
! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
|
||||
!
|
||||
! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
|
||||
!
|
||||
! : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
|
||||
!
|
||||
! : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
|
||||
!
|
||||
! etc... the advantage is to guarentee no spurious mixing because of numerical problems.
|
||||
END_DOC
|
||||
double precision, allocatable :: leigvec_unsrtd(:,:),reigvec_unsrtd(:,:),eigval_unsrtd(:)
|
||||
integer, allocatable :: list_degen(:,:),list_same_degen(:)
|
||||
integer, allocatable :: iorder(:),list_degen_sorted(:)
|
||||
integer :: n_degen_list,n_degen,size_mat,i,j,k,icount,m,index_degen
|
||||
integer :: ii,jj,i_good,j_good,n_real
|
||||
double precision, allocatable :: mat_tmp(:,:),eigval_tmp(:),leigvec_tmp(:,:),reigvec_tmp(:,:)
|
||||
|
||||
allocate(leigvec_unsrtd(n,n),reigvec_unsrtd(n,n),eigval_unsrtd(n))
|
||||
leigvec_unsrtd = 0.d0
|
||||
reigvec_unsrtd = 0.d0
|
||||
eigval_unsrtd = 0.d0
|
||||
|
||||
allocate(list_degen(mo_num,0:mo_num))
|
||||
|
||||
! obtain degeneracies
|
||||
call give_degen_full_list(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
|
||||
allocate(iorder(n_degen_list),list_degen_sorted(n_degen_list))
|
||||
do i =1, n_degen_list
|
||||
n_degen = list_degen(i,0)
|
||||
list_degen_sorted(i) = n_degen
|
||||
iorder(i)=i
|
||||
enddo
|
||||
! sort by number of degeneracies
|
||||
call isort(list_degen_sorted,iorder,n_degen_list)
|
||||
integer :: icount_eigval
|
||||
logical, allocatable :: is_ok(:)
|
||||
allocate(is_ok(n_degen_list))
|
||||
is_ok = .True.
|
||||
icount_eigval = 0
|
||||
! loop over degeneracies
|
||||
do i = 1, n_degen_list
|
||||
if(.not.is_ok(i))cycle
|
||||
is_ok(i) = .False.
|
||||
n_degen = list_degen_sorted(i)
|
||||
print*,'diagonalizing for n_degen = ',n_degen
|
||||
k = 1
|
||||
! group all the entries having the same degeneracies
|
||||
! do while (list_degen_sorted(i+k)==n_degen)
|
||||
do m = i+1,n_degen_list
|
||||
if(list_degen_sorted(m)==n_degen)then
|
||||
is_ok(i+k)=.False.
|
||||
k += 1
|
||||
endif
|
||||
enddo
|
||||
print*,'number of identical degeneracies = ',k
|
||||
size_mat = k*n_degen
|
||||
print*,'size_mat = ',size_mat
|
||||
allocate(mat_tmp(size_mat,size_mat),list_same_degen(size_mat))
|
||||
allocate(eigval_tmp(size_mat),leigvec_tmp(size_mat,size_mat),reigvec_tmp(size_mat,size_mat))
|
||||
! group all the elements sharing the same degeneracy
|
||||
icount = 0
|
||||
do j = 1, k ! jth set of degeneracy
|
||||
index_degen = iorder(i+j-1)
|
||||
do m = 1, n_degen
|
||||
icount += 1
|
||||
list_same_degen(icount) = list_degen(index_degen,m)
|
||||
enddo
|
||||
enddo
|
||||
print*,'icount = ',icount
|
||||
print*,'list of elements '
|
||||
do icount = 1, size_mat
|
||||
print*,icount,list_same_degen(icount)
|
||||
enddo
|
||||
! you copy subset of matrix elements having all the same degeneracy in mat_tmp
|
||||
do ii = 1, size_mat
|
||||
i_good = list_same_degen(ii)
|
||||
do jj = 1, size_mat
|
||||
j_good = list_same_degen(jj)
|
||||
mat_tmp(jj,ii) = mat_ref(j_good,i_good)
|
||||
enddo
|
||||
enddo
|
||||
call non_hrmt_bieig( size_mat, mat_tmp&
|
||||
, leigvec_tmp, reigvec_tmp&
|
||||
, n_real, eigval_tmp)
|
||||
do ii = 1, size_mat
|
||||
icount_eigval += 1
|
||||
eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues
|
||||
do jj = 1, size_mat ! copy the eigenvectors
|
||||
j_good = list_same_degen(jj)
|
||||
leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
|
||||
reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
|
||||
enddo
|
||||
enddo
|
||||
deallocate(mat_tmp,list_same_degen)
|
||||
deallocate(eigval_tmp,leigvec_tmp,reigvec_tmp)
|
||||
enddo
|
||||
print*,'icount,n',icount,n
|
||||
|
||||
deallocate(iorder)
|
||||
allocate(iorder(n))
|
||||
print*,'before sort '
|
||||
do i = 1, n
|
||||
print*,'i,eigval(i) = ',i,eigval_unsrtd(i)
|
||||
iorder(i) = i
|
||||
enddo
|
||||
call dsort(eigval_unsrtd,iorder,n)
|
||||
print*,'after sort '
|
||||
do i = 1, n
|
||||
i_good = iorder(i)
|
||||
eigval(i) = eigval_unsrtd(i)
|
||||
print*,'i,eigval(i) = ',i,eigval(i)
|
||||
do j = 1, n
|
||||
leigvec(j,i) = leigvec_unsrtd(j,i_good)
|
||||
reigvec(j,i) = reigvec_unsrtd(j,i_good)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine give_degen_full_list(A,n,thr,list_degen,n_degen_list)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! you enter with an array A(n) and spits out all the elements degenerated up to thr
|
||||
!
|
||||
! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL
|
||||
!
|
||||
! list_degen(i,0) = number of degenerate entries
|
||||
!
|
||||
! list_degen(i,1) = index of the first degenerate entry
|
||||
!
|
||||
! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries
|
||||
!
|
||||
! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
|
||||
END_DOC
|
||||
double precision,intent(in) :: A(n)
|
||||
double precision,intent(in) :: thr
|
||||
integer, intent(in) :: n
|
||||
integer, intent(out) :: list_degen(n,0:n),n_degen_list
|
||||
logical, allocatable :: is_ok(:)
|
||||
allocate(is_ok(n))
|
||||
integer :: i,j,icount,icheck
|
||||
n_degen_list = 0
|
||||
is_ok = .True.
|
||||
do i = 1, n
|
||||
if(.not.is_ok(i))cycle
|
||||
n_degen_list +=1
|
||||
is_ok(i) = .False.
|
||||
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
|
||||
! print*,A(i),A(j)
|
||||
is_ok(j) = .False.
|
||||
icount += 1
|
||||
list_degen(n_degen_list,icount) = j
|
||||
endif
|
||||
enddo
|
||||
list_degen(n_degen_list,0) = icount
|
||||
enddo
|
||||
icheck = 0
|
||||
do i = 1, n_degen_list
|
||||
icheck += list_degen(i,0)
|
||||
enddo
|
||||
if(icheck.ne.n)then
|
||||
print*,'pb ! :: icheck.ne.n'
|
||||
print*,icheck,n
|
||||
stop
|
||||
endif
|
||||
end
|
||||
@ -66,8 +66,6 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state
|
||||
do i = 1, mo_num
|
||||
write(*,'(100(F16.10,X))')-dm_tmp(:,i)
|
||||
enddo
|
||||
! call non_hrmt_diag_split_degen( mo_num, dm_tmp&
|
||||
! call non_hrmt_fock_mat( mo_num, dm_tmp&
|
||||
call non_hrmt_bieig( mo_num, dm_tmp&
|
||||
, natorb_tc_leigvec_mo, natorb_tc_reigvec_mo&
|
||||
, n_real, natorb_tc_eigval )
|
||||
|
||||
51
src/tc_bi_ortho/test_natorb.irp.f
Normal file
51
src/tc_bi_ortho/test_natorb.irp.f
Normal file
@ -0,0 +1,51 @@
|
||||
program test_natorb
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end.
|
||||
END_DOC
|
||||
print *, 'Hello world'
|
||||
my_grid_becke = .True.
|
||||
my_n_pt_r_grid = 30
|
||||
my_n_pt_a_grid = 50
|
||||
read_wf = .True.
|
||||
touch read_wf
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
call routine
|
||||
! call test
|
||||
|
||||
end
|
||||
|
||||
subroutine routine
|
||||
implicit none
|
||||
double precision, allocatable :: fock_diag(:),eigval(:),leigvec(:,:),reigvec(:,:),mat_ref(:,:)
|
||||
allocate(eigval(mo_num),leigvec(mo_num,mo_num),reigvec(mo_num,mo_num),fock_diag(mo_num),mat_ref(mo_num, mo_num))
|
||||
double precision, allocatable :: eigval_ref(:),leigvec_ref(:,:),reigvec_ref(:,:)
|
||||
allocate(eigval_ref(mo_num),leigvec_ref(mo_num,mo_num),reigvec_ref(mo_num,mo_num))
|
||||
|
||||
double precision :: thr_deg
|
||||
integer :: i,n_real,j
|
||||
print*,'fock_matrix'
|
||||
do i = 1, mo_num
|
||||
fock_diag(i) = Fock_matrix_mo(i,i)
|
||||
print*,i,fock_diag(i)
|
||||
enddo
|
||||
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 non_hrmt_bieig( mo_num, mat_ref&
|
||||
, leigvec_ref, reigvec_ref&
|
||||
, n_real, eigval_ref)
|
||||
print*,'TEST ***********************************'
|
||||
double precision :: accu_l, accu_r
|
||||
do i = 1, mo_num
|
||||
accu_l = 0.d0
|
||||
accu_r = 0.d0
|
||||
do j = 1, mo_num
|
||||
accu_r += reigvec_ref(j,i) * reigvec(j,i)
|
||||
accu_l += leigvec_ref(j,i) * leigvec(j,i)
|
||||
enddo
|
||||
print*,i
|
||||
write(*,'(I3,X,100(F16.10,X))')i,eigval(i),eigval_ref(i),accu_l,accu_r
|
||||
enddo
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user