PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00
Code Releases Activity

working on rotate

Browse Source
This commit is contained in:
eginer 2022-10-26 18:32:08 +02:00
parent 6436bfc9c8
commit 9716ee2b66
7 changed files with 535 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

75
src/non_hermit_dav/biorthog.irp.f
View File

@ -329,6 +329,7 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
n_good = 0
thr = 1.d-5
do i = 1, n
print*, 'Re(i) + Im(i)', WR(i), WI(i)
if(dabs(WI(i)) .lt. thr) then
n_good += 1
else
@ -392,8 +393,8 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
! -------------------------------------------------------------------------------------
! check bi-orthogonality
thr_d = 1d-8
thr_nd = 1d-8
thr_d = 1d-10 ! -7
thr_nd = 1d-10 ! -7
allocate( S(n_real_eigv,n_real_eigv) )
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .false.)
@ -422,8 +423,8 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
! ---
!call impose_orthog_degen_eigvec(n, eigval, reigvec)
!call impose_orthog_degen_eigvec(n, eigval, leigvec)
! call impose_orthog_degen_eigvec(n, eigval, reigvec)
! call impose_orthog_degen_eigvec(n, eigval, leigvec)
call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
@ -1040,6 +1041,72 @@ subroutine split_matrix_degen(aw,n,shift)
end
subroutine give_degen(a,n,shift,list_degen,n_degen_list)
implicit none
BEGIN_DOC
! returns n_degen_list :: the number of degenerated SET of elements (i.e. with |A(i)-A(i+1)| below shift)
!
! for each of these sets, list_degen(1,i) = first degenerate element of the set i,
!
! list_degen(2,i) = last degenerate element of the set i.
END_DOC
double precision,intent(in) :: A(n)
double precision,intent(in) :: shift
integer, intent(in) :: n
integer, intent(out) :: list_degen(2,n),n_degen_list
integer :: i,j,n_degen,k
logical :: keep_on
double precision,allocatable :: Aw(:)
list_degen = -1
allocate(Aw(n))
Aw = A
i=1
k = 0
do while(i.lt.n)
if(dabs(Aw(i)-Aw(i+1)).lt.shift)then
k+=1
j=1
list_degen(1,k) = i
keep_on = .True.
do while(keep_on)
if(i+j.gt.n)then
keep_on = .False.
exit
endif
if(dabs(Aw(i)-Aw(i+j)).lt.shift)then
j+=1
else
keep_on=.False.
exit
endif
enddo
n_degen = j
list_degen(2,k) = list_degen(1,k)-1 + n_degen
j=0
keep_on = .True.
do while(keep_on)
if(i+j+1.gt.n)then
keep_on = .False.
exit
endif
if(dabs(Aw(i+j)-Aw(i+j+1)).lt.shift)then
Aw(i+j) += (j-n_degen/2) * shift
j+=1
else
keep_on = .False.
exit
endif
enddo
Aw(i+n_degen-1) += (n_degen-1-n_degen/2) * shift
i+=n_degen
else
i+=1
endif
enddo
n_degen_list = k
end
subroutine cancel_small_elmts(aw,n,shift)
implicit none
BEGIN_DOC

212
src/non_hermit_dav/lapack_diag_non_hermit.irp.f
View File

@ -1254,6 +1254,105 @@ end subroutine impose_orthog_svd
! ---
subroutine impose_orthog_svd_overlap(n, m, C,overlap)
implicit none
integer, intent(in) :: n, m
double precision, intent(in ) :: overlap(n,n)
double precision, intent(inout) :: C(n,m)
integer :: i, j, num_linear_dependencies
double precision :: threshold
double precision, allocatable :: S(:,:), tmp(:,:),Stmp(:,:)
double precision, allocatable :: U(:,:), Vt(:,:), D(:)
print *, ' apply SVD to orthogonalize vectors'
! ---
allocate(S(m,m),Stmp(n,m))
! S = C.T x overlap x C
call dgemm( 'N', 'N', n, m, n, 1.d0 &
, overlap, size(overlap, 1), C, size(C, 1) &
, 0.d0, Stmp, size(Stmp, 1) )
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, C, size(C, 1), Stmp, size(Stmp, 1) &
, 0.d0, S, size(S, 1) )
print *, ' eigenvec overlap bef SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
! ---
allocate(U(m,m), Vt(m,m), D(m))
call svd(S, m, U, m, D, Vt, m, m, m)
deallocate(S)
threshold = 1.d-6
num_linear_dependencies = 0
do i = 1, m
if(abs(D(i)) <= threshold) then
D(i) = 0.d0
num_linear_dependencies += 1
else
ASSERT (D(i) > 0.d0)
D(i) = 1.d0 / dsqrt(D(i))
endif
enddo
if(num_linear_dependencies > 0) then
write(*,*) ' linear dependencies = ', num_linear_dependencies
write(*,*) ' m = ', m
stop
endif
! ---
allocate(tmp(n,m))
! tmp <-- C x U
call dgemm( 'N', 'N', n, m, m, 1.d0 &
, C, size(C, 1), U, size(U, 1) &
, 0.d0, tmp, size(tmp, 1) )
deallocate(U, Vt)
! C <-- tmp x sigma^-0.5
do j = 1, m
do i = 1, n
C(i,j) = tmp(i,j) * D(j)
enddo
enddo
deallocate(D, tmp)
! ---
allocate(S(m,m))
! S = C.T x C
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, C, size(C, 1), C, size(C, 1) &
, 0.d0, S, size(S, 1) )
print *, ' eigenvec overlap aft SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
deallocate(S)
! ---
end subroutine impose_orthog_svd
! ---
subroutine impose_orthog_GramSchmidt(n, m, C)
implicit none
@ -2389,3 +2488,116 @@ end subroutine impose_biorthog_svd
! ---
subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
implicit none
integer, intent(in) :: n, m
double precision, intent(in) :: overlap(n,n)
double precision, intent(inout) :: L(n,m), R(n,m)
integer :: i, j, num_linear_dependencies
double precision :: threshold
double precision, allocatable :: S(:,:), tmp(:,:),Stmp(:,:)
double precision, allocatable :: U(:,:), V(:,:), Vt(:,:), D(:)
! ---
allocate(S(m,m),Stmp(n,m))
! S = C.T x overlap x C
call dgemm( 'N', 'N', n, m, n, 1.d0 &
, overlap, size(overlap, 1), R, size(R, 1) &
, 0.d0, Stmp, size(Stmp, 1) )
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, L, size(L, 1), Stmp, size(Stmp, 1) &
, 0.d0, S, size(S, 1) )
print *, ' overlap bef SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
! ---
allocate(U(m,m), Vt(m,m), D(m))
call svd(S, m, U, m, D, Vt, m, m, m)
deallocate(S)
threshold = 1.d-6
num_linear_dependencies = 0
do i = 1, m
if(abs(D(i)) <= threshold) then
D(i) = 0.d0
num_linear_dependencies += 1
else
ASSERT (D(i) > 0.d0)
D(i) = 1.d0 / dsqrt(D(i))
endif
enddo
if(num_linear_dependencies > 0) then
write(*,*) ' linear dependencies = ', num_linear_dependencies
write(*,*) ' m = ', m
stop
endif
allocate(V(m,m))
do i = 1, m
do j = 1, m
V(j,i) = Vt(i,j)
enddo
enddo
deallocate(Vt)
! ---
allocate(tmp(n,m))
! tmp <-- R x V
call dgemm( 'N', 'N', n, m, m, 1.d0 &
, R, size(R, 1), V, size(V, 1) &
, 0.d0, tmp, size(tmp, 1) )
deallocate(V)
! R <-- tmp x sigma^-0.5
do j = 1, m
do i = 1, n
R(i,j) = tmp(i,j) * D(j)
enddo
enddo
! tmp <-- L x U
call dgemm( 'N', 'N', n, m, m, 1.d0 &
, L, size(L, 1), U, size(U, 1) &
, 0.d0, tmp, size(tmp, 1) )
deallocate(U)
! L <-- tmp x sigma^-0.5
do j = 1, m
do i = 1, n
L(i,j) = tmp(i,j) * D(j)
enddo
enddo
deallocate(D, tmp)
! ---
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 aft SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
deallocate(S)
! ---
end subroutine impose_biorthog_svd
! ---

67
src/some_mu_of_r/print_mu_av_tc.irp.f Normal file
View File

@ -0,0 +1,67 @@
program print_mu_av_tc
implicit none
read_wf = .True.
touch read_wf
print*,'average_mu_lda = ',average_mu_lda
! print*,'average_mu_rs = ',average_mu_rs
print*,'average_mu_rs_c = ',average_mu_rs_c
print*,'average_mu_rs_c_lda = ',average_mu_rs_c_lda
call plot_mu_tc
end
subroutine plot_mu_tc
implicit none
integer :: i,nx,m
double precision :: xmin,xmax,dx
double precision :: r(3)
double precision :: sqpi
double precision :: rho_a_hf, rho_b_hf, g0,rho_hf
double precision :: rs,grad_n,grad_n_a(3), grad_n_b(3)
double precision :: g0_UEG_mu_inf
double precision :: cst_rs,alpha_rs,mu_rs, mu_rs_c, mu_lda, mu_grad_n
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
if (no_core_density)then
output=trim(ezfio_filename)//'.tc_mu_fc'
else
output=trim(ezfio_filename)//'.tc_mu'
endif
i_unit_output = getUnitAndOpen(output,'w')
nx = 5000
xmin = -10.d0
xmax = 10.d0
dx = (xmax - xmin)/dble(nx)
r = 0.d0
r(3) = xmin
sqpi = dsqrt(dacos(-1.d0))
cst_rs = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0)
alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi
write(i_unit_output,*)'# r(3) rho_hf mu_rs mu_rs_c mu_lda mu_grad_n'
do i = 1, nx
rho_a_hf = 0.d0
grad_n = 0.d0
call density_and_grad_alpha_beta(r,rho_a_hf,rho_b_hf, grad_n_a, grad_n_b)
do m = 1, 3
grad_n += grad_n_a(m)*grad_n_a(m) + grad_n_b(m)*grad_n_b(m) + 2.d0 * grad_n_a(m) * grad_n_b(m)
enddo
rho_hf = rho_a_hf + rho_b_hf
grad_n = dsqrt(grad_n)
grad_n = grad_n/(4.d0 * rho_hf)
rs = cst_rs * rho_hf**(-1.d0/3.d0)
g0 = g0_UEG_mu_inf(rho_a_hf,rho_b_hf)
mu_rs = 1.d0/rs
mu_rs_c = alpha_rs/dsqrt(rs)
mu_lda = - 1.d0 / (dlog(2.d0 * g0) * sqpi)
mu_grad_n = grad_n
write(i_unit_output,'(100(F16.10,X))')r(3),rho_hf,mu_rs,mu_rs_c,mu_lda,mu_grad_n
r(3) += dx
enddo
end

7
src/some_mu_of_r/some_mu_of_r.irp.f
View File

@ -1,7 +0,0 @@
program some_mu_of_r
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
end

12
src/tc_scf/diago_bi_ort_tcfock.irp.f
View File

@ -18,13 +18,13 @@
PROVIDE Fock_matrix_tc_mo_tot
call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot &
, fock_tc_leigvec_mo, fock_tc_reigvec_mo &
, n_real_tc, eigval_right_tmp )
call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot &
, fock_tc_leigvec_mo, fock_tc_reigvec_mo &
, n_real_tc, eigval_right_tmp )
!if(max_ov_tc_scf)then
! call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &
! , n_real_tc, eigval_right_tmp )
! call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &
! , n_real_tc, eigval_right_tmp )
!else
! call non_hrmt_diag_split_degen_bi_orthog( mo_num, Fock_matrix_tc_mo_tot &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &

175
src/tc_scf/test_rotate_2.irp.f Normal file
View File

@ -0,0 +1,175 @@
program print_angles
implicit none
my_grid_becke = .True.
! my_n_pt_r_grid = 30
! my_n_pt_a_grid = 50
my_n_pt_r_grid = 10 ! small grid for quick debug
my_n_pt_a_grid = 14 ! small grid for quick debug
call routine
end
subroutine routine
implicit none
integer :: i,j,k,n_degen_list,m,n,n_degen,ilast,ifirst
double precision, allocatable :: mo_r_coef_good(:,:),mo_l_coef_good(:,:)
allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
double precision, allocatable :: mo_r_coef_new(:,:)
double precision :: norm
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))
do j = 1, ao_num
mo_r_coef_new(j,i) *= norm
enddo
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))
do i = 1, mo_num
fock_diag(i) = fock_matrix_mo(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_degen_tc,list_degen,n_degen_list)
print*,'fock_matrix_mo'
do i = 1, mo_num
print*,i,fock_diag(i),angle_left_right(i)
enddo
print*,'Overlap '
do i = 1, mo_num
write(*,'(I2,X,100(F8.4,X))')i,s_mat(:,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
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)
enddo
! Orthogonalization of right functions
print*,'Orthogonalization of right functions'
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
! Orthogonalization of left functions
print*,'Orthogonalization of left functions'
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
print*,'Overlap lef-right '
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
do j = 1, n_degen
write(*,'(100(F8.4,X))')stmp(:,j)
enddo
T = 0.d0
Snew = 0.d0
call maxovl(n_degen, n_degen, stmp, T, Snew)
print*,'overlap after'
do j = 1, n_degen
write(*,'(100(F16.10,X))')Snew(:,j)
enddo
! mo_l_coef_new = 0.D0
! do j = 1, n_degen
! do k = 1, n_degen
! do m = 1, ao_num
! mo_l_coef_new(m,j) += T(k,j) * mo_l_coef_tmp(m,k)
! enddo
! enddo
! enddo
call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0 &
, mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) &
, 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) )
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
print*,'Overlap test'
do j = 1, n_degen
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
call impose_biorthog_svd_overlap(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
print*,'LAST OVERLAP '
do j = 1, n_degen
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_l_coef_new,ao_overlap,stmp)
print*,'LEFT OVERLAP '
do j = 1, n_degen
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
print*,'RIGHT OVERLAP '
do j = 1, n_degen
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)
enddo
deallocate(stmp,smat2)
deallocate(mo_r_coef_tmp,mo_l_coef_tmp,mo_l_coef_new)
deallocate(T,Snew)
enddo
allocate(stmp(mo_num, mo_num))
call build_s_matrix(ao_num,n_degen,mo_l_coef_good,mo_r_coef_good,ao_overlap,stmp)
print*,'LEFT/RIGHT OVERLAP '
do j = 1, mo_num
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
call build_s_matrix(ao_num,n_degen,mo_l_coef_good,mo_l_coef_good,ao_overlap,stmp)
print*,'LEFT/LEFT OVERLAP '
do j = 1, mo_num
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
call build_s_matrix(ao_num,n_degen,mo_r_coef_good,mo_r_coef_good,ao_overlap,stmp)
print*,'RIGHT/RIGHT OVERLAP '
do j = 1, mo_num
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
end
subroutine build_s_matrix(m,n,C1,C2,overlap,smat)
implicit none
integer, intent(in) :: m,n
double precision, intent(in) :: C1(m,n),C2(m,n),overlap(m,m)
double precision, intent(out):: smat(n,n)
integer :: i,j,k,l
smat = 0.D0
do i = 1, n
do j = 1, n
do k = 1, m
do l = 1, m
smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j)
enddo
enddo
enddo
enddo
end
subroutine orthog_functions(m,n,coef,overlap)
implicit none
integer, intent(in) :: m,n
double precision, intent(in) :: overlap(m,m)
double precision, intent(inout) :: coef(m,n)
double precision, allocatable :: stmp(:,:)
integer :: j
allocate(stmp(n,n))
call build_s_matrix(m,n,coef,coef,overlap,stmp)
print*,'overlap before'
do j = 1, n
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
call impose_orthog_svd_overlap(m, n, coef,overlap)
call build_s_matrix(m,n,coef,coef,overlap,stmp)
do j = 1, n
coef(1,:m) *= 1.d0/dsqrt(stmp(j,j))
enddo
print*,'overlap after'
call build_s_matrix(m,n,coef,coef,overlap,stmp)
do j = 1, n
write(*,'(100(F16.10,X))')stmp(:,j)
enddo
end

6
src/utils/loc.f
View File

@ -12,8 +12,10 @@ C W: new overlap matrix
C
C
implicit real*8(a-h,o-y),logical*1(z)
parameter (id1=700)
dimension s(id1,id1),t(id1,id1),w(id1,id1)
! parameter (id1=700)
! dimension s(id1,id1),t(id1,id1),w(id1,id1)
double precision, intent(inout) :: s(n,n)
double precision, intent(out) :: t(n,n),w(n,n)
data small/1.d-6/
zprt=.true.