LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-08 20:33:26 +01:00
Code Issues Releases Wiki Activity

Avoid extra calculation of S2

Browse Source
This commit is contained in:
Anthony Scemama 2016-10-11 23:42:05 +02:00
parent 7906976c5a
commit c3a6b5ba6b
3 changed files with 29 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
src/Davidson/diagonalization_hs2.irp.f
View File

@ -1,4 +1,4 @@
subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit)
subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit)
use bitmasks
implicit none
BEGIN_DOC
@ -22,7 +22,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nin
integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st)
double precision, intent(out) :: energies(N_st), s2_out(N_st)
double precision, allocatable :: H_jj(:), S2_jj(:)
double precision :: diag_h_mat_elem
@ -46,6 +46,9 @@ subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nin
!$OMP END PARALLEL
call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit)
do i=1,N_st_diag
s2_out(i) = S2_jj(i)
enddo
deallocate (H_jj,S2_jj)
end
@ -79,7 +82,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
END_DOC
integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(in) :: H_jj(sze), S2_jj(sze)
double precision, intent(in) :: H_jj(sze)
double precision, intent(inout) :: S2_jj(sze)
integer, intent(in) :: iunit
double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st_diag)
@ -200,7 +204,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
call normalize(u_in(1,k),sze)
enddo
do while (.not.converged)
do k=1,N_st_diag
@ -239,12 +243,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
! enddo
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), W(1,shift+1), size(W,1), &
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), W(1,shift+1), size(W,1), &
0.d0, h(1,shift+1), size(h,1))
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), S(1,shift+1), size(S,1), &
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), S(1,shift+1), size(S,1), &
0.d0, s_(1,shift+1), size(s_,1))
! Diagonalize h
@ -254,14 +258,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! Compute S2 for each eigenvector
! -------------------------------
call dgemm('N','N',shift2,shift2,shift2, &
1.d0, s_, size(s_,1), y, size(y,1), &
call dgemm('N','N',shift2,shift2,shift2, &
1.d0, s_, size(s_,1), y, size(y,1), &
0.d0, s_tmp, size(s_tmp,1))
call dgemm('T','N',shift2,shift2,shift2, &
1.d0, y, size(y,1), s_tmp, size(s_tmp,1), &
call dgemm('T','N',shift2,shift2,shift2, &
1.d0, y, size(y,1), s_tmp, size(s_tmp,1), &
0.d0, s_, size(s_,1))
do k=1,shift2
s2(k) = s_(k,k) + S_z2_Sz
enddo
@ -326,8 +330,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
do k=1,N_st_diag
do i=1,sze
R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
* (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz)
R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
* (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz)
enddo
if (k <= N_st) then
residual_norm(k) = u_dot_u(R(1,k),sze)
@ -367,10 +371,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
! enddo
!
call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), &
U(1,shift2+k),1,0.d0,c,1)
call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), &
c,1,1.d0,U(1,shift2+k),1)
call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), &
U(1,shift2+k),1,0.d0,c,1)
call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), &
c,1,1.d0,U(1,shift2+k),1)
!
! do l=1,k-1
! c(1) = u_dot_v(U(1,shift2+k),U(1,shift2+l),sze)
@ -379,7 +383,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
! enddo
!
call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), &
call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), &
U(1,shift2+k),1,0.d0,c,1)
call dgemv('N',sze,k-1,-1.d0,U(1,shift2+1),size(U,1), &
c,1,1.d0,U(1,shift2+k),1)

8
src/Davidson/diagonalize_CI.irp.f
View File

@ -59,13 +59,13 @@ END_PROVIDER
! size(CI_eigenvectors,1), &
! N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants)
!
call davidson_diag_HS2(psi_det,CI_eigenvectors, &
! call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,&
! min(N_det,N_states_diag),size(CI_eigenvectors,1))
call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_eigenvectors_s2, &
size(CI_eigenvectors,1),CI_electronic_energy, &
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants)
call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,&
min(N_det,N_states_diag),size(CI_eigenvectors,1))
else if (diag_algorithm == "Lapack") then

2
src/Davidson/u0Hu0.irp.f
View File

@ -252,8 +252,6 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
ave_workload = ave_workload/dble(shortcut(0,1))
print *, 'Ave workload :', ave_workload
do sh=shortcut(0,1),1,-1
workload = shortcut(0,1)+dble(shortcut(sh+1,1) - shortcut(sh,1))**2
do i=sh, shortcut(0,2), shortcut(0,1)