From 2e9f539827ad7aa208f2ce0dd52291a18ed44c98 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 20 Jun 2022 19:55:44 +0200 Subject: [PATCH] Fixed some bugs for S>0. --- src/csf/configuration_CI_sigma_helpers.irp.f | 22 ++++++++++++++++++-- src/csf/obtain_I_foralpha.irp.f | 22 ++++++++++++++------ src/csf/sigma_vector.irp.f | 7 ++++--- src/davidson/diagonalization_hcfg.irp.f | 5 ++++- 4 files changed, 44 insertions(+), 12 deletions(-) diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index 8d447818..cea7640c 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -44,10 +44,13 @@ use bitmasks logical :: pqAlreadyGenQ logical :: pqExistsQ logical :: ppExistsQ + integer*8 :: MS double precision :: t0, t1 call wall_time(t0) + MS = elec_alpha_num-elec_beta_num + allocate(tableUniqueAlphas(mo_num,mo_num)) NalphaIcfg_list = 0 @@ -128,8 +131,13 @@ use bitmasks Jsomo = IBCLR(Isomo,p-1) Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBSET(Idomo,q-1) - kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) - kend = idxI-1 + ! Check for Minimal alpha electrons (MS) + if(POPCNT(Jsomo).ge.MS)then + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) + kend = idxI-1 + else + cycle + endif else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then ! DOMO -> VMO Jsomo = IBSET(Isomo,p-1) @@ -148,6 +156,10 @@ use bitmasks else print*,"Something went wrong in obtain_associated_alphaI" endif + ! Check for Minimal alpha electrons (MS) + if(POPCNT(Jsomo).lt.MS)then + cycle + endif ! Again, we don't have to search from 1 ! we just use seniority to find the @@ -211,6 +223,12 @@ use bitmasks Jsomo = IBCLR(Isomo,p-1) Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBSET(Idomo,q-1) + if(POPCNT(Jsomo).ge.MS)then + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) + kend = idxI-1 + else + cycle + endif else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then ! DOMO -> VMO Jsomo = IBSET(Isomo,p-1) diff --git a/src/csf/obtain_I_foralpha.irp.f b/src/csf/obtain_I_foralpha.irp.f index 3926b908..7d7ae09b 100644 --- a/src/csf/obtain_I_foralpha.irp.f +++ b/src/csf/obtain_I_foralpha.irp.f @@ -205,12 +205,14 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: xordiffSOMODOMO integer :: ndiffSOMO integer :: ndiffDOMO - integer :: nxordiffSOMODOMO - integer :: iii,ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes - integer :: listholes(mo_num) - integer :: holetype(mo_num) - integer :: end_index - integer :: Nsomo_alpha + integer :: nxordiffSOMODOMO + integer :: iii,ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes + integer :: listholes(mo_num) + integer :: holetype(mo_num) + integer :: end_index + integer :: Nsomo_alpha + integer*8 :: MS + MS = elec_alpha_num-elec_beta_num nconnectedI = 0 end_index = N_configuration @@ -233,6 +235,10 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Idomo = Ialpha(1,2) Jsomo = psi_configuration(1,1,i) Jdomo = psi_configuration(1,2,i) + ! Check for Minimal alpha electrons (MS) + if(POPCNT(Isomo).lt.MS)then + cycle + endif diffSOMO = IEOR(Isomo,Jsomo) ndiffSOMO = POPCNT(diffSOMO) !if(idxI.eq.1)then @@ -299,6 +305,10 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !IRP_ELSE p = TRAILZ(Isomo) + 1 !IRP_ENDIF + ! Check for Minimal alpha electrons (MS) + !if(POPCNT(Isomo).lt.MS)then + ! cycle + !endif end if case (2) ! DOMO -> SOMO diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index c9299cae..1c2357b2 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1578,7 +1578,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze deallocate(excitationIds_single) deallocate(excitationTypes_single) - !print *," singles part psi(1,177)=",psi_out(1,177) + !print *," singles part psi(1,5)=",psi_out(1,5) allocate(listconnectedJ(N_INT,2,max(sze,10000))) allocate(alphas_Icfg(N_INT,2,max(sze,10000))) @@ -1622,6 +1622,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze call obtain_connected_J_givenI(i, Icfg, listconnectedJ, idslistconnectedJ, nconnectedJ, ntotJ) ! TODO : remove doubly excited for return + !print *,"I=",i," isomo=",psi_configuration(1,1,i)," idomo=",psi_configuration(1,2,i), " psiout=",psi_out(1,5) do k = 1,Nalphas_Icfg ! Now generate all singly excited with respect to a given alpha CFG @@ -1640,7 +1641,6 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze totcolsTKI = 0 rowsTKI = -1 NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k)) - !print *,"alphas_Icfg=",alphas_Icfg(1,1,k) do j = 1,nconnectedI NSOMOI = getNSOMO(connectedI_alpha(:,:,j)) p = excitationIds(1,j) @@ -1690,6 +1690,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze moj = excitationIds(2,l) ! s mol = excitationIds(1,l) ! r diagfac = diagfactors_0 * diagfactors(l)* mo_two_e_integral(mok,mol,moi,moj)! g(pq,sr) = + !print *,"p=",mok,"q=",mol,"r=",moi,"s=",moj do m = 1,colsikpq ! = (ik|jl) GIJpqrs(totcolsTKI+m,l) = diagfac @@ -1749,7 +1750,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze deallocate(excitationTypes) deallocate(diagfactors) - !print *," psi(1,177)=",psi_out(1,177) + !print *," psi(1,823)=",psi_out(1,823), " g(1 8, 3 15)=",mo_two_e_integral(1,8,3,15), " ncore=",n_core_orb ! Add the diagonal contribution !$OMP DO diff --git a/src/davidson/diagonalization_hcfg.irp.f b/src/davidson/diagonalization_hcfg.irp.f index 92097a2f..7bd5dd8d 100644 --- a/src/davidson/diagonalization_hcfg.irp.f +++ b/src/davidson/diagonalization_hcfg.irp.f @@ -339,7 +339,10 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N !call H_u_0_nstates_openmp(u_in,U2,N_st_diag,sze) !call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),W_csf2(1,1)) !do i=1,sze_csf - ! print *,"I=",i," qp=",W_csf2(i,1)," my=",W_csf(i,1), " diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1)) + ! print *,"I=",i," qp=",W_csf2(i,1)," my=",W_csf(i,1)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1)) + ! if(dabs(dabs(W_csf2(i,1))-dabs(W_csf(i,1))) .gt. 1.0e-10)then + ! print *,"somo=",psi_configuration(1,1,i)," domo=",psi_configuration(1,2,i)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1)) + ! endif !end do !stop deallocate(tmpW)