diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c index d55da188..6094f5bb 100644 --- a/src/csf/cfgCI_utils.c +++ b/src/csf/cfgCI_utils.c @@ -32,8 +32,8 @@ void getncsfs1(int *inpnsomo, int *inpms, int *outncsfs){ } void getncsfs(int NSOMO, int MS, int *outncsfs){ - int nparcoupl = (NSOMO + MS)/2; - int nparcouplp1 = ((NSOMO + MS)/2)+1; + int nparcoupl = (NSOMO + MS)/2; // n_alpha + int nparcouplp1 = ((NSOMO + MS)/2)+1; // n_alpha + 1 double tmpndets=0.0; if(NSOMO == 0){ (*outncsfs) = 1; diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index a5bc845e..8d447818 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -107,7 +107,7 @@ use bitmasks if(Nsomo_I .EQ. 0) then kstart = 1 else - kstart = cfg_seniority_index(max(0,Nsomo_I-2)) + kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)) endif kend = idxI-1 @@ -121,14 +121,14 @@ use bitmasks Jsomo = IBCLR(Isomo,p-1) Jsomo = IBSET(Jsomo,q-1) Jdomo = Idomo - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) kend = idxI-1 else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then ! SOMO -> SOMO Jsomo = IBCLR(Isomo,p-1) Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBSET(Idomo,q-1) - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) kend = idxI-1 else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then ! DOMO -> VMO @@ -143,7 +143,7 @@ use bitmasks Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBCLR(Idomo,p-1) Jdomo = IBSET(Jdomo,q-1) - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) kend = idxI-1 else print*,"Something went wrong in obtain_associated_alphaI" @@ -227,11 +227,13 @@ use bitmasks endif ! SOMO - NalphaIcfg += 1 - !print *,i,j,"|",NalphaIcfg - alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo - alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) - NalphaIcfg_list(idxI) = NalphaIcfg + !print *,i,j,"|",NalphaIcfg, Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + if(POPCNT(Jsomo) .ge. NSOMOMin) then + NalphaIcfg += 1 + alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo + alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + NalphaIcfg_list(idxI) = NalphaIcfg + endif endif end do end do @@ -240,7 +242,7 @@ use bitmasks ppExistsQ = .False. Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) do k = kstart, idxI-1 diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k))) ndiffSOMO = POPCNT(diffSOMO) @@ -257,10 +259,12 @@ use bitmasks ! Diagonal part (pp,qq) if(nholes > 0 .AND. (.NOT. ppExistsQ))then ! SOMO - NalphaIcfg += 1 - alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1) - alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2) - NalphaIcfg_list(idxI) = NalphaIcfg + if(POPCNT(Jsomo) .ge. NSOMOMin) then + NalphaIcfg += 1 + alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1) + alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2) + NalphaIcfg_list(idxI) = NalphaIcfg + endif endif NalphaIcfg = 0 @@ -373,7 +377,7 @@ END_PROVIDER if(Nsomo_I .EQ. 0) then kstart = 1 else - kstart = cfg_seniority_index(max(0,Nsomo_I-2)) + kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)) endif kend = idxI-1 !print *,"Isomo" @@ -398,14 +402,14 @@ END_PROVIDER Jsomo = IBCLR(Isomo,p-1) Jsomo = IBSET(Jsomo,q-1) Jdomo = Idomo - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) kend = idxI-1 else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then ! SOMO -> SOMO Jsomo = IBCLR(Isomo,p-1) Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBSET(Idomo,q-1) - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) kend = idxI-1 else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then ! DOMO -> VMO @@ -420,7 +424,7 @@ END_PROVIDER Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBCLR(Idomo,p-1) Jdomo = IBSET(Jdomo,q-1) - kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) kend = idxI-1 else print*,"Something went wrong in obtain_associated_alphaI" diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index 9d6e5219..9ddd4fee 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -73,7 +73,6 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) ! perhaps blocking with CFGs of same seniority ! can be more efficient - print *,"bfIcfg=",bfIcfg, " ndetI=",ndetI allocate(tempBuffer(bfIcfg,ndetI)) tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 60c4631e..c9299cae 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -234,7 +234,7 @@ end subroutine get_phase_qp_to_cfg ! initialization psi_coef_config = 0.d0 DetToCSFTransformationMatrix(0,:,:) = 1.d0 - do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2 + do i = 2-iand(MS,1), NSOMOMax,2 Isomo = IBSET(0_8, i) - 1_8 ! rows = Ncsfs ! cols = Ndets @@ -287,7 +287,9 @@ end subroutine get_phase_qp_to_cfg if (psi_configuration(k,1,i) == 0_bit_kind) cycle s = s + popcnt(psi_configuration(k,1,i)) enddo - bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) + salpha = (s+MS)/2 + bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1)))) + !bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) ! perhaps blocking with CFGs of same seniority ! can be more efficient @@ -1374,7 +1376,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze integer(omp_lock_kind), allocatable :: lock(:) call omp_set_max_active_levels(1) - print *," sze = ",sze + !print *," sze = ",sze allocate(lock(sze)) do i=1,sze call omp_init_lock(lock(i)) @@ -1383,7 +1385,6 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze ! print *,"i=",i," psi_cfg_data_1=",psi_config_data(i,1)," psi_cfg_data_2=",psi_config_data(i,2) !end do - !print *," sze = ",sze allocate(diag_energies(n_CSF)) call calculate_preconditioner_cfg(diag_energies) @@ -1639,6 +1640,7 @@ 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) @@ -1648,6 +1650,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze ! for E_pp E_rs and E_ppE_rr case rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1) colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) + !print *,"j=",j," Nsomo=",NSOMOalpha," rowsikpq=",rowsikpq," colsikpq=",colsikpq, " p=",pmodel," q=",qmodel, " extyp=",extype totcolsTKI += colsikpq rowsTKI = rowsikpq enddo diff --git a/src/davidson/diagonalization_hcfg.irp.f b/src/davidson/diagonalization_hcfg.irp.f index bde2bdd9..92097a2f 100644 --- a/src/davidson/diagonalization_hcfg.irp.f +++ b/src/davidson/diagonalization_hcfg.irp.f @@ -324,6 +324,7 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N tmpU(kk,ii) = U_csf(ii,shift+kk) 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) do kk=1,N_st_diag do ii=1,sze_csf @@ -331,6 +332,16 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N enddo enddo + !U_csf = 0.0d0 + !U_csf(1,1) = 1.0d0 + !u_in = 0.0d0 + !call convertWFfromCSFtoDET(N_st_diag,tmpU,U2) + !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)) + !end do + !stop deallocate(tmpW) deallocate(tmpU) endif diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index 42059e98..78478c57 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -71,10 +71,8 @@ END_PROVIDER if (diag_algorithm == "Davidson") then - !if (do_csf) then - if (.true.) then - !if (sigma_vector_algorithm == 'det') then - if (.false.) then + if (do_csf) then + if (sigma_vector_algorithm == 'det') then call davidson_diag_H_csf(psi_det,CI_eigenvectors, & size(CI_eigenvectors,1),CI_electronic_energy, & N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)