scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2025-01-03 01:55:52 +01:00
Code Issues Releases Wiki Activity

Hdiag

Browse Source
This commit is contained in:
Anthony Scemama 2021-05-20 23:46:00 +02:00
parent 4372cc8403
commit 693b81265e
1 changed files with 97 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
devel/svdwf/buildpsi_diagSVDit_modif_v2.irp.f
View File

@ -35,7 +35,7 @@ subroutine run
double precision :: norm_coeff_psi, inv_sqrt_norm_coeff_psi
double precision :: overlapU, overlapU_mat, overlapV, overlapV_mat, overlap_psi
double precision, allocatable :: Hdiag(:), Hkl(:,:), H0(:,:), H(:,:,:,:)
double precision, allocatable :: H_diag(:,:), Hkl(:,:), H0(:,:), H(:,:,:,:)
double precision, allocatable :: psi_postsvd(:,:), coeff_psi_perturb(:)
integer :: n_FSVD, n_selected, n_toselect, n_tmp, it_svd, it_svd_max
@ -121,6 +121,7 @@ subroutine run
E_prev = 0.d0
allocate(H(n_det_alpha_unique,n_det_beta_unique,n_det_alpha_unique,n_det_beta_unique))
allocate(H_diag(n_det_alpha_unique,n_det_beta_unique))
allocate(psi_postsvd(n_det_alpha_unique,n_det_beta_unique))
do while( ( it_svd .lt. it_svd_max) .and. ( tol_energy .gt. 1d-6 ) )
@ -137,7 +138,7 @@ subroutine run
Dref = Dref / dsqrt(norm)
! print *, '-- Compute H --'
call const_H_uv(Uref, Vref, H)
call const_H_uv(Uref, Vref, H, H_diag, n_selected)
! H0(i,j) = < u_i v_j | H | u_i v_j >
! E0 = < psi_0 | H | psi_0 >
@ -148,7 +149,7 @@ subroutine run
enddo
enddo
E0 = E0 + nuclear_repulsion
! print *,' E0 =', E0
print *,' E0 =', E0
double precision, allocatable :: eigval0(:)
double precision, allocatable :: eigvec0(:,:,:)
@ -169,6 +170,7 @@ subroutine run
! print *, ' --- Diag post-SVD --- '
call lapack_diag(eigval0, eigvec0, H_tmp, n_selected**2, n_selected**2)
print *, 'eig =', eigval0(1) + nuclear_repulsion
deallocate(H_tmp, eigval0)
! print *, ' --- SVD --- '
@ -177,7 +179,7 @@ subroutine run
deallocate(eigvec0)
! print *, ' --- Compute H --- '
call const_H_uv(Uref, Vref, H)
call const_H_uv(Uref, Vref, H, H_diag, n_selected)
! H0(i,j) = < u_i v_j | H | u_i v_j >
! E0 = < psi_0 | H | psi_0 >
@ -208,8 +210,8 @@ subroutine run
ctmp = ctmp + H(k,l,i,j) * psi_postsvd(k,l)
enddo
enddo
psi_postsvd(i,j) = ctmp / (E0 - (H(i,j,i,j)+nuclear_repulsion) )
Ept2 += ctmp*ctmp / (E0 - (H(i,j,i,j)+nuclear_repulsion) )
psi_postsvd(i,j) = ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
Ept2 += ctmp*ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
enddo
enddo
@ -221,8 +223,8 @@ subroutine run
ctmp = ctmp + H(k,l,i,j) * psi_postsvd(k,l)
enddo
enddo
psi_postsvd(i,j) = ctmp / (E0 - (H(i,j,i,j)+nuclear_repulsion) )
Ept2 += ctmp*ctmp / (E0 - (H(i,j,i,j)+nuclear_repulsion) )
psi_postsvd(i,j) = ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
Ept2 += ctmp*ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
enddo
enddo
! do j=n_selected+1,n_det_beta_unique
@ -326,15 +328,17 @@ end subroutine perform_newpostSVD
subroutine const_H_uv(Uref, Vref, H)
subroutine const_H_uv(Uref, Vref, H, H_diag, n_selected)
USE OMP_LIB
implicit none
integer, intent(in) :: n_selected
double precision, intent(in) :: Uref(n_det_alpha_unique,n_det_beta_unique)
double precision, intent(in) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(out) :: H(n_det_alpha_unique,n_det_beta_unique, n_det_alpha_unique,n_det_beta_unique)
double precision, intent(out) :: H_diag(n_det_alpha_unique,n_det_beta_unique)
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: i, j, k, l, degree
@ -350,7 +354,6 @@ subroutine const_H_uv(Uref, Vref, H)
nb = n_det_beta_unique
allocate( H0(na,nb,na,nb) )
allocate( H1(nb,na,nb,na) )
det1(:,1) = psi_det_alpha_unique(:,1)
det2(:,1) = psi_det_alpha_unique(:,1)
@ -387,29 +390,108 @@ subroutine const_H_uv(Uref, Vref, H)
!$OMP END DO
!$OMP END PARALLEL
call wall_time(t1)
double precision :: H0_d(n_det_alpha_unique,n_det_beta_unique)
double precision :: H1_d(n_det_alpha_unique,n_det_beta_unique)
double precision :: tmp3(n_det_alpha_unique,n_det_beta_unique,n_det_alpha_unique)
double precision, allocatable :: tmp1(:,:), tmp0(:,:)
tmp3 = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,k,l,m,det1,det2,degree,h12,tmp1,tmp0)&
!$OMP SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique,&
!$OMP N_int,tmp3,Uref,Vref,H_diag)
allocate(tmp1(na,na), tmp0(na,na))
do i=1,na
do m=1,na
tmp1(m,i) = Uref(i,m)
enddo
enddo
!$OMP DO
do l = 1, nb
det2(:,2) = psi_det_beta_unique(:,l)
do j = 1, nb
det1(:,2) = psi_det_beta_unique(:,j)
call get_excitation_degree_spin(det1(1,2),det2(1,2),degree,N_int)
if (degree > 2) cycle
do k = 1, na
det2(:,1) = psi_det_alpha_unique(:,k)
do i = 1, na
det1(:,1) = psi_det_alpha_unique(:,i)
call get_excitation_degree(det1,det2,degree,N_int)
if ( degree > 2) cycle
call i_H_j(det1, det2, N_int, h12)
do m=1,nb
tmp3(m,j,l) = tmp3(m,j,l) + h12 * tmp1(m,i) * tmp1(m,k)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do m=1,nb
do l=1,nb
do j=1,nb
tmp1(j,l) = tmp3(m,j,l)
enddo
enddo
call DGEMM('N','N',nb,nb,nb,1.d0, &
tmp1, size(tmp1,1), &
Vref, size(Vref,1), &
0.d0, tmp0, size(tmp0,1))
do n=1,na
H_diag(m,n) = 0.d0
do j=1,nb
H_diag(m,n) = H_diag(m,n) + tmp0(j,n) * Vref(j,n)
enddo
enddo
enddo
!$OMP END DO
deallocate(tmp1, tmp0)
!$OMP END PARALLEL
! (i,j,k,l) -> (j,k,l,p)
allocate( H1(nb,na,nb,na) )
call DGEMM('T','N', nb * na * nb, na, na, &
1.d0, H0, size(H0,1), Uref, size(Uref,1), 0.d0, H1, size(H1,1)*size(H1,2)*size(H1,3))
deallocate( H0 )
! (j,k,l,p) -> (k,l,p,q)
allocate( H0(na,nb,na,nb) )
call DGEMM('T','N', na * nb * na, nb, nb, &
1.d0, H1, size(H1,1), Vref, size(Vref,1), 0.d0, H0, size(H0,1)*size(H0,2)*size(H0,3))
deallocate( H1 )
! (k,l,p,q) -> (l,p,q,r)
allocate( H1(nb,na,nb,na) )
call DGEMM('T','N', nb * na * nb, na, na, &
1.d0, H0, size(H0,1), Uref, size(Uref,1), 0.d0, H1, size(H1,1)*size(H1,2)*size(H1,3))
deallocate( H0 )
! (l,p,q,r) -> (p,q,r,s)
call DGEMM('T','N', na * nb * na, nb, nb, &
1.d0, H1, size(H1,1), Vref, size(Vref,1), 0.d0, H, size(H,1)*size(H,2)*size(H,3))
do j=1,n_selected
do i=1,n_selected
print *, H_diag(i,j), H(i,j,i,j)
enddo
enddo
deallocate(H1)
call wall_time(t2)
print *, t1-t0, t2-t1
print *, 't=', t1-t0, t2-t1
double precision :: t0, t1, t2
deallocate(H1,H0)
stop
end