LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-26 21:33:30 +01:00
Code Releases Activity

Debugging sigma-vector.

Browse Source
This commit is contained in:
v1j4y 2022-06-06 21:40:53 +02:00
parent 1c41bfdd17
commit a3200652ce
3 changed files with 105 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/csf/cfgCI_utils.c
View File

@ -341,7 +341,6 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
Get BFtoDeterminant Matrix Get BFtoDeterminant Matrix
************************************/ ************************************/
printf("In convertcsftodet\n");
convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI); convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
int rowsI = 0; int rowsI = 0;

121
src/csf/sigma_vector.irp.f
View File

@ -40,49 +40,68 @@
detDimperBF = 0 detDimperBF = 0
MS = elec_alpha_num-elec_beta_num MS = elec_alpha_num-elec_beta_num
! number of cfgs = number of dets for 0 somos ! number of cfgs = number of dets for 0 somos
n_CSF = 0 n_CSF = cfg_seniority_index(NSOMOMin)-1
ncfgprev = cfg_seniority_index(0) print *,"start=",n_CSF
ncfgpersomo = ncfgprev ncfgprev = cfg_seniority_index(NSOMOMin)
do i = iand(MS,1), NSOMOMax-2,2 !do i = 0-iand(MS,1)+2, NSOMOMax,2
if(cfg_seniority_index(i) .EQ. -1) then do i = NSOMOMin+2, NSOMOMax,2
cycle if(cfg_seniority_index(i) .EQ. -1)then
endif ncfgpersomo = N_configuration + 1
if(cfg_seniority_index(i+2) .EQ. -1) then else
ncfgpersomo = N_configuration + 1 ncfgpersomo = cfg_seniority_index(i)
else endif
if(cfg_seniority_index(i+2) > ncfgpersomo) then ncfg = ncfgpersomo - ncfgprev
ncfgpersomo = cfg_seniority_index(i+2) !detDimperBF = max(1,nint((binom(i,(i+1)/2))))
else dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1))))
k = 0 n_CSF += ncfg * dimcsfpercfg
do while(cfg_seniority_index(i+2+k) < ncfgpersomo) print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", n_CSF
k = k + 2 !if(cfg_seniority_index(i+2) == -1) EXIT
ncfgpersomo = cfg_seniority_index(i+2+k) !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF
enddo ncfgprev = cfg_seniority_index(i)
endif
endif
ncfg = ncfgpersomo - ncfgprev
if(i .EQ. 0 .OR. i .EQ. 1) then
dimcsfpercfg = 1
elseif( i .EQ. 3) then
dimcsfpercfg = 2
else
if(iand(MS,1) .EQ. 0) then
dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1))))
else
dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2))))
endif
endif
n_CSF += ncfg * dimcsfpercfg
if(cfg_seniority_index(i+2) > ncfgprev) then
ncfgprev = cfg_seniority_index(i+2)
else
k = 0
do while(cfg_seniority_index(i+2+k) < ncfgprev)
k = k + 2
ncfgprev = cfg_seniority_index(i+2+k)
enddo
endif
enddo enddo
!n_CSF = 0
!ncfgprev = cfg_seniority_index(0)
!ncfgpersomo = ncfgprev
!do i = iand(MS,1), NSOMOMax-2,2
! if(cfg_seniority_index(i) .EQ. -1) then
! cycle
! endif
! if(cfg_seniority_index(i+2) .EQ. -1) then
! ncfgpersomo = N_configuration + 1
! else
! if(cfg_seniority_index(i+2) > ncfgpersomo) then
! ncfgpersomo = cfg_seniority_index(i+2)
! else
! k = 0
! do while(cfg_seniority_index(i+2+k) < ncfgpersomo)
! k = k + 2
! ncfgpersomo = cfg_seniority_index(i+2+k)
! enddo
! endif
! endif
! ncfg = ncfgpersomo - ncfgprev
! if(i .EQ. 0 .OR. i .EQ. 1) then
! dimcsfpercfg = 1
! elseif( i .EQ. 3) then
! dimcsfpercfg = 2
! else
! if(iand(MS,1) .EQ. 0) then
! dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1))))
! else
! dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2))))
! endif
! endif
! n_CSF += ncfg * dimcsfpercfg
! if(cfg_seniority_index(i+2) > ncfgprev) then
! ncfgprev = cfg_seniority_index(i+2)
! else
! k = 0
! do while(cfg_seniority_index(i+2+k) < ncfgprev)
! k = k + 2
! ncfgprev = cfg_seniority_index(i+2+k)
! enddo
! endif
!enddo
END_PROVIDER END_PROVIDER
@ -1477,14 +1496,14 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
deallocate(excitationIds_single) deallocate(excitationIds_single)
deallocate(excitationTypes_single) deallocate(excitationTypes_single)
allocate(listconnectedJ(N_INT,2,max(sze,100))) allocate(listconnectedJ(N_INT,2,max(sze,1000)))
allocate(alphas_Icfg(N_INT,2,max(sze,100))) allocate(alphas_Icfg(N_INT,2,max(sze,1000)))
allocate(connectedI_alpha(N_INT,2,max(sze,100))) allocate(connectedI_alpha(N_INT,2,max(sze,1000)))
allocate(idxs_connectedI_alpha(max(sze,100))) allocate(idxs_connectedI_alpha(max(sze,1000)))
allocate(excitationIds(2,max(sze,100))) allocate(excitationIds(2,max(sze,1000)))
allocate(excitationTypes(max(sze,100))) allocate(excitationTypes(max(sze,1000)))
allocate(diagfactors(max(sze,100))) allocate(diagfactors(max(sze,1000)))
allocate(idslistconnectedJ(max(sze,100))) allocate(idslistconnectedJ(max(sze,1000)))
allocate(CCmattmp(n_st,NBFmax)) allocate(CCmattmp(n_st,NBFmax))
! Loop over all selected configurations ! Loop over all selected configurations
@ -1643,11 +1662,13 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
do i = 1,n_CSF do i = 1,n_CSF
do kk=1,n_st do kk=1,n_st
psi_out(kk,i) += diag_energies(i)*psi_in(kk,i) psi_out(kk,i) += diag_energies(i)*psi_in(kk,i)
print *,psi_in(kk,i)," | ",psi_out(kk,i)
enddo enddo
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
stop
call omp_set_max_active_levels(4) call omp_set_max_active_levels(4)
deallocate(diag_energies) deallocate(diag_energies)

55
src/davidson/diagonalization_hcsf_dressed.irp.f
View File

@ -267,8 +267,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
do k=N_st+1,N_st_diag do k=N_st+1,N_st_diag
do i=1,sze do i=1,sze
call random_number(r1) !call random_number(r1)
call random_number(r2) !call random_number(r2)
r1 = 0.5
r2 = 0.5
r1 = dsqrt(-2.d0*dlog(r1)) r1 = dsqrt(-2.d0*dlog(r1))
r2 = dtwo_pi*r2 r2 = dtwo_pi*r2
u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
@ -286,8 +288,8 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
enddo enddo
! Make random verctors eigenstates of S2 ! Make random verctors eigenstates of S2
call convertWFfromDETtoCSF(N_st_diag,U(1,1),U_csf(1,1)) call convertWFfromDETtoCSF(N_st_diag,U,U_csf)
call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1)) !call convertWFfromCSFtoDET(N_st_diag,U_csf,U)
do while (.not.converged) do while (.not.converged)
itertot = itertot+1 itertot = itertot+1
@ -330,8 +332,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
do kk=1,N_st_diag do kk=1,N_st_diag
do ii=1,sze_csf do ii=1,sze_csf
tmpU(kk,ii) = U_csf(ii,shift+kk) tmpU(kk,ii) = U_csf(ii,shift+kk)
tmpU(kk,ii) = 0.0d0
enddo enddo
enddo enddo
tmpU(1,1) = 1.0d0
call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1) call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1)
do kk=1,N_st_diag do kk=1,N_st_diag
do ii=1,sze_csf do ii=1,sze_csf
@ -479,24 +483,24 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
! Compute residual vector and davidson step ! Compute residual vector and davidson step
! ----------------------------------------- ! -----------------------------------------
if (without_diagonal) then !if (without_diagonal) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
do k=1,N_st_diag do k=1,N_st_diag
do i=1,sze do i=1,sze
U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) & U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) &
/max(H_jj(i) - lambda (k),1.d-2) /max(H_jj(i) - lambda (k),1.d-2)
enddo
enddo enddo
!$OMP END PARALLEL DO enddo
else !$OMP END PARALLEL DO
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) !else
do k=1,N_st_diag ! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
do i=1,sze ! do k=1,N_st_diag
U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) ! do i=1,sze
enddo ! U(i,k) = (lambda(k) * U(i,k) - W(i,k) )
enddo ! enddo
!$OMP END PARALLEL DO ! enddo
endif ! !$OMP END PARALLEL DO
!endif
do k=1,N_st do k=1,N_st
residual_norm(k) = u_dot_u(U(1,k),sze) residual_norm(k) = u_dot_u(U(1,k),sze)
@ -543,6 +547,15 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
! Re-contract U ! Re-contract U
! ------------- ! -------------
call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0, &
W_csf, size(W_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
do k=1,N_st_diag
do i=1,sze_csf
W_csf(i,k) = u_in(i,k)
enddo
enddo
call convertWFfromCSFtoDET(N_st_diag,W_csf,W)
call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0, & call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0, &
U_csf, size(U_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) U_csf, size(U_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
do k=1,N_st_diag do k=1,N_st_diag