From 5622b9790da599e15880d7aa0362db13b49e8ed3 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Thu, 15 Dec 2022 17:15:30 +0100 Subject: [PATCH] Fixed some bugs. Diagonal energy is OK. Some bugs still present. --- src/csf/configuration_CI_sigma_helpers.irp.f | 55 ++-- src/csf/conversion.irp.f | 1 + src/csf/obtain_I_foralpha.irp.f | 72 +++-- src/csf/sigma_vector.irp.f | 282 +++++++++++++++---- src/davidson/diagonalization_hcfg.irp.f | 35 ++- 5 files changed, 329 insertions(+), 116 deletions(-) diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index 8fb8383e..19533bd5 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -524,7 +524,7 @@ use bitmasks !!! !!!END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*(mo_num))] + BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*12)] &BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ] implicit none !use bitmasks @@ -549,7 +549,8 @@ use bitmasks integer(bit_kind) :: Isomo(N_int), Isomop(N_int), Isomoq(N_int) integer(bit_kind) :: Jdomo(N_int), Jdomop(N_int), Jdomoq(N_int) integer(bit_kind) :: Jsomo(N_int), Jsomop(N_int), Jsomoq(N_int) - integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int) + !integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int) + integer(bit_kind) :: diffDOMO, xordiffSOMODOMO, diffSOMO integer :: ndiffSOMO integer :: ndiffDOMO integer :: nxordiffSOMODOMO @@ -674,13 +675,15 @@ use bitmasks Jsomo(jint) = IBCLR(Jsomo(jint),jpos) endif + Nsomo_J=0 do ii=1, N_int Jcfg(ii,1) = Jsomo(ii) Jcfg(ii,2) = Jdomo(ii) + Nsomo_J += POPCNT(Jsomo(ii)) enddo - call bitstring_to_list(Jcfg,listall,nelall,N_int) - Nsomo_J = nelall + !call bitstring_to_list(Jcfg(1,1),listall,nelall,N_int) + !Nsomo_J = nelall ! Check for Minimal alpha electrons (MS) if(Nsomo_J.lt.MS)then @@ -705,15 +708,13 @@ use bitmasks ndiffDOMO = 0 nxordiffSOMODOMO = 0 do ii = 1, N_int - Jsomo = Jcfg(ii,1) - Jdomo = Jcfg(ii,2) - diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k))) - ndiffSOMO += POPCNT(diffSOMO(ii)) - diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k))) + diffSOMO = IEOR(Jcfg(ii,1),iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k))) + ndiffSOMO += POPCNT(diffSOMO) + diffDOMO = IEOR(Jcfg(ii,2),iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k))) xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO += POPCNT(diffDOMO(ii)) - nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii)) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) end do if((ndiffSOMO .ne. 0) .and. (ndiffSOMO .ne. 2)) cycle @@ -818,6 +819,7 @@ use bitmasks if(Nsomo_J .ge. NSOMOMin) then !print *," Idx = ",idxI, "p = ",pp, " q = ",qq," Jsomo=",Jsomo(1), " Jdomo=",IOR(Jdomo(1),ISHFT(1_8,n_core_orb)-1) NalphaIcfg += 1 + !print *," Idx = ",idxI, " Nalpha=",NalphaIcfg alphasIcfg_list(:,1,idxI,NalphaIcfg) = Jcfg(:,1) if(n_core_orb .le. 63)then alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jcfg(1,2),ISHFT(1_8,n_core_orb)-1) @@ -864,7 +866,7 @@ use bitmasks ndiffSOMO = 0 do ii=1,N_int diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k))) - ndiffSOMO += POPCNT(diffSOMO(ii)) + ndiffSOMO += POPCNT(diffSOMO) end do ! ndiffSOMO cannot be 0 (I /= k) ! if ndiffSOMO /= 2 then it has to be greater than 2 and hense @@ -875,8 +877,8 @@ use bitmasks do ii=1,N_int diffDOMO = IEOR(Icfg(ii,2),iand(act_bitmask(ii,1),psi_configuration(ii,2,k))) xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO += POPCNT(diffDOMO(ii)) - nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii)) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) end do if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then ppExistsQ = .TRUE. @@ -926,7 +928,8 @@ END_PROVIDER integer(bit_kind) :: Isomo(N_int), Isomop(N_int), Isomoq(N_int) integer(bit_kind) :: Jdomo(N_int), Jdomop(N_int), Jdomoq(N_int) integer(bit_kind) :: Jsomo(N_int), Jsomop(N_int), Jsomoq(N_int) - integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int) + !integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int) + integer(bit_kind) :: diffDOMO, xordiffSOMODOMO, diffSOMO integer :: ndiffSOMO integer :: ndiffDOMO integer :: nxordiffSOMODOMO @@ -1084,15 +1087,13 @@ END_PROVIDER ndiffDOMO = 0 nxordiffSOMODOMO = 0 do ii = 1, N_int - Jsomo = Jcfg(ii,1) - Jdomo = Jcfg(ii,2) - diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k))) - ndiffSOMO += POPCNT(diffSOMO(ii)) - diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k))) + diffSOMO = IEOR(Jcfg(ii,1),iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k))) + ndiffSOMO += POPCNT(diffSOMO) + diffDOMO = IEOR(Jcfg(ii,2),iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k))) xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO += POPCNT(diffDOMO(ii)) - nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii)) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) end do if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then @@ -1211,7 +1212,7 @@ END_PROVIDER do k = 1, idxI-1 do ii=1,N_int diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k))) - ndiffSOMO += POPCNT(diffSOMO(ii)) + ndiffSOMO += POPCNT(diffSOMO) end do ! ndiffSOMO cannot be 0 (I /= k) ! if ndiffSOMO /= 2 then it has to be greater than 2 and hense @@ -1222,8 +1223,8 @@ END_PROVIDER do ii=1,N_int diffDOMO = IEOR(Icfg(ii,2),iand(act_bitmask(ii,1),psi_configuration(ii,2,k))) xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO += POPCNT(diffDOMO(ii)) - nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii)) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) end do if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then ppExistsQ = .TRUE. diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index 494c3bfa..7c6c8363 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -114,6 +114,7 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) integer :: idx integer MS MS = elec_alpha_num-elec_beta_num + print *,"size=",size(tmp_psi_coef_det,1)," ",size(tmp_psi_coef_det,2) countcsf = 0 diff --git a/src/csf/obtain_I_foralpha.irp.f b/src/csf/obtain_I_foralpha.irp.f index 8606f556..1d4e81fc 100644 --- a/src/csf/obtain_I_foralpha.irp.f +++ b/src/csf/obtain_I_foralpha.irp.f @@ -102,7 +102,7 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI, xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) ndiffDOMO += POPCNT(diffDOMO) nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) end do if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then @@ -243,13 +243,16 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer :: listholes(mo_num) integer :: holetype(mo_num) integer :: end_index, ishift - integer :: Nsomo_alpha, pp,qq, nperm + integer :: Nsomo_alpha, pp,qq, nperm, iint, ipos integer*8 :: MS integer :: exc(0:2,2,2), tz, m, n, high, low integer :: listall(N_int*bit_kind_size), nelall + integer :: nconnectedExtradiag, nconnectedDiag integer(bit_kind) :: hole, particle, tmp MS = elec_alpha_num-elec_beta_num + nconnectedExtradiag=0 + nconnectedDiag=0 nconnectedI = 0 end_index = N_configuration @@ -260,10 +263,12 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !Nsomo_alpha = POPCNT(Isomo) Icfg = Ialpha Nsomo_alpha = 0 + !print *," Ialpha=" do i=1,N_int Isomo = Ialpha(i,1) Idomo = Ialpha(i,2) Nsomo_alpha += POPCNT(Isomo) + !print *,Isomo, Idomo, "Nsomo=",Nsomo_alpha end do end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_alpha+4,elec_num))-1) if(end_index .LT. 0) end_index= N_configuration @@ -293,20 +298,25 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI ndiffSOMO = 0 ndiffDOMO = 0 nxordiffSOMODOMO = 0 + nsomoJ=0 + nsomoalpha=0 do ii=1,N_int Isomo = Ialpha(ii,1) Idomo = Ialpha(ii,2) Jsomo = psi_configuration(ii,1,i) Jdomo = psi_configuration(ii,2,i) + nsomoJ += POPCNT(Jsomo) + nsomoalpha += POPCNT(Isomo) diffSOMO = IEOR(Isomo,Jsomo) ndiffSOMO += POPCNT(diffSOMO) diffDOMO = IEOR(Idomo,Jdomo) xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) ndiffDOMO += POPCNT(diffDOMO) nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) end do - Jcfg = psi_configuration(:,:,i) + !Jcfg = psi_configuration(:,:,i) + !print *,"nxordiffSOMODOMO(4)=",nxordiffSOMODOMO, " ndiffSOMO(2)=",ndiffSOMO if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then select case(ndiffDOMO) @@ -328,7 +338,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Jsomo = psi_configuration(ii,1,i) IJsomo = IEOR(Isomo, Jsomo) if(popcnt(IAND(Isomo,IJsomo)) > 0)then - p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + ii * bit_kind_size + p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + (ii-1) * bit_kind_size EXIT endif end do @@ -337,22 +347,24 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Isomo = Ialpha(ii,1) Jsomo = psi_configuration(ii,1,i) IJsomo = IEOR(Isomo, Jsomo) - IJsomo = IBCLR(IJsomo,p-1) + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + IJsomo = IBCLR(IJsomo,ipos-1) + endif if(popcnt(IJsomo) > 0)then - q = TRAILZ(IJsomo) + 1 + ii * bit_kind_size + q = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size EXIT endif enddo endif !assert ( p == pp) !assert ( q == qq) - !print *," --- p=",p," q=",q + !print *," 1--- p=",p," q=",q case (1) ! DOMO -> VMO ! or ! SOMO -> SOMO - nsomoJ = POPCNT(Jsomo) - nsomoalpha = POPCNT(Isomo) if(nsomoJ .GT. nsomoalpha) then ! DOMO -> VMO !print *,"obt DOMO -> VMO" @@ -371,7 +383,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Idomo = Ialpha(ii,2) Jdomo = psi_configuration(ii,2,i) if(popcnt(IEOR(Idomo,Jdomo)) > 0)then - p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + ii * bit_kind_size + p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + (ii-1) * bit_kind_size EXIT endif end do @@ -380,9 +392,13 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Isomo = Ialpha(ii,1) Jsomo = psi_configuration(ii,1,i) Isomo = IEOR(Isomo, Jsomo) - Isomo = IBCLR(Isomo,p-1) + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + Isomo = IBCLR(Isomo,ipos-1) + endif if(popcnt(Isomo) > 0)then - q = TRAILZ(Isomo) + 1 + ii * bit_kind_size + q = TRAILZ(Isomo) + 1 + (ii-1) * bit_kind_size EXIT endif end do @@ -404,13 +420,16 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !endif else ! Find p + !print *,"Ialpha somo=",Ialpha(1,1), Ialpha(2,1)," Ialpha domo=",Ialpha(1,2), Ialpha(2,2) + !print *,"J somo=",psi_configuration(1,1,i), psi_configuration(2,1,i)," J domo=",psi_configuration(1,2,i),& + !psi_configuration(2,2,i) do ii=1,N_int Isomo = Ialpha(ii,1) Jsomo = psi_configuration(ii,1,i) Idomo = Ialpha(ii,2) Jdomo = psi_configuration(ii,2,i) if(popcnt(IEOR(Idomo,Jdomo)) > 0)then - q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + ii * bit_kind_size + q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + (ii-1) * bit_kind_size EXIT endif enddo @@ -419,9 +438,14 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Isomo = Ialpha(ii,1) Jsomo = psi_configuration(ii,1,i) Isomo = IEOR(Isomo, Jsomo) - Isomo = IBCLR(Isomo,q-1) + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + Isomo = IBCLR(Isomo,ipos-1) + endif + !print *,"ii=",ii," Isomo=",Isomo if(popcnt(Isomo) > 0)then - p = TRAILZ(Isomo) + 1 + ii * bit_kind_size + p = TRAILZ(Isomo) + 1 + (ii-1) * bit_kind_size EXIT endif enddo @@ -429,6 +453,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !assert ( p == pp) !assert ( q == qq) endif + !print *," 2--- p=",p," q=",q case (2) ! DOMO -> SOMO !print *,"obt DOMO -> SOMO" @@ -447,7 +472,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Jdomo = psi_configuration(ii,2,i) IJsomo = IEOR(Isomo, Jsomo) if(popcnt(IAND(Jsomo,IJsomo)) > 0)then - p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + ii * bit_kind_size + p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + (ii-1) * bit_kind_size EXIT endif enddo @@ -456,20 +481,26 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI Isomo = Ialpha(ii,1) Jsomo = psi_configuration(ii,1,i) IJsomo = IEOR(Isomo, Jsomo) - IJsomo = IBCLR(IJsomo,p-1) + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + IJsomo = IBCLR(IJsomo,ipos-1) + endif if(popcnt(IJsomo) > 0)then - q = TRAILZ(IJsomo) + 1 + ii * bit_kind_size + q = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size EXIT endif enddo endif !assert ( p == pp) !assert ( q == qq) + !print *," 3--- p=",p," q=",q case default print *,"something went wront in get connectedI" end select starti = psi_config_data(i,1) endi = psi_config_data(i,2) + nconnectedExtradiag+=1 nconnectedI += 1 do k=1,N_int connectedI(k,1,nconnectedI) = psi_configuration(k,1,i) @@ -534,6 +565,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI if(holetype(k) .EQ. 1) then starti = psi_config_data(i,1) endi = psi_config_data(i,2) + nconnectedDiag+=1 nconnectedI += 1 connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) idxs_connectedI(nconnectedI)=starti @@ -544,6 +576,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI else starti = psi_config_data(i,1) endi = psi_config_data(i,2) + nconnectedDiag+=1 nconnectedI += 1 connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) idxs_connectedI(nconnectedI)=starti @@ -556,5 +589,6 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI enddo endif end do + !print *,"nconnectedExtradiag=",nconnectedExtradiag," nconnectedDiad=",nconnectedDiag end subroutine obtain_connected_I_foralpha diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 76f9bfc3..9fe81fe9 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -885,7 +885,7 @@ subroutine calculate_preconditioner_cfg(diag_energies) Idomo = psi_configuration(1,2,i) Icfg(1,1) = psi_configuration(1,1,i) Icfg(1,2) = psi_configuration(1,2,i) - NSOMOI = getNSOMO(psi_configuration(:,:,i)) + !NSOMOI = getNSOMO(psi_configuration(:,:,i)) starti = psi_config_data(i,1) endi = psi_config_data(i,2) @@ -1239,27 +1239,34 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod integer,intent(in) :: p,q integer,intent(in) :: extype integer,intent(out) :: pmodel,qmodel - !integer(bit_kind) :: Isomo(N_int) - !integer(bit_kind) :: Idomo(N_int) - !integer(bit_kind) :: Jsomo(N_int) - !integer(bit_kind) :: Jdomo(N_int) - integer*8 :: Isomo - integer*8 :: Idomo - integer*8 :: Jsomo - integer*8 :: Jdomo + integer(bit_kind) :: Isomo(N_int) + integer(bit_kind) :: Idomo(N_int) + integer(bit_kind) :: Jsomo(N_int) + integer(bit_kind) :: Jdomo(N_int) + !integer*8 :: Isomo + !integer*8 :: Idomo + !integer*8 :: Jsomo + !integer*8 :: Jdomo integer*8 :: mask - integer :: iint, ipos + integer :: iint, ipos, ii !integer(bit_kind) :: Isomotmp(N_int) !integer(bit_kind) :: Jsomotmp(N_int) integer*8 :: Isomotmp integer*8 :: Jsomotmp integer :: pos0,pos0prev + integer :: tmpp, tmpq ! TODO Flag (print) when model space indices is > 64 - Isomo = Ialpha(1,1) - Idomo = Ialpha(1,2) - Jsomo = Jcfg(1,1) - Jdomo = Jcfg(1,2) + do ii=1,N_int + !Isomo = Ialpha(ii,1) + !Idomo = Ialpha(ii,2) + !Jsomo = Jcfg(ii,1) + !Jdomo = Jcfg(ii,2) + Isomo(ii) = Ialpha(ii,1) + Idomo(ii) = Ialpha(ii,2) + Jsomo(ii) = Jcfg(ii,1) + Jdomo(ii) = Jcfg(ii,2) + end do pos0prev = 0 pmodel = p qmodel = q @@ -1273,40 +1280,139 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod ! SOMO -> SOMO ! remove all domos !print *,"type -> SOMO -> SOMO" - mask = ISHFT(1_8,p) - 1 - Isomotmp = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Isomotmp = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + !mask = ISHFT(1_8,p) - 1 + !Isomotmp = IAND(Isomo,mask) + !pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + !mask = ISHFT(1_8,q) - 1 + !Isomotmp = IAND(Isomo,mask) + !qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + !print *,"iint=",iint, " p=",p + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Isomotmp = IAND(Isomo(ii),mask) + tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + !print *,"iint=",iint, " ipos=",ipos,"pmodel=",pmodel, XOR(Isomotmp,mask),Isomo(iint) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Isomotmp = IAND(Isomo(ii),mask) + tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + !print *,"iint=",iint, " ipos=",ipos,"qmodel=",qmodel case (2) ! DOMO -> VMO ! remove all domos except one at p !print *,"type -> DOMO -> VMO" - mask = ISHFT(1_8,p) - 1 - Jsomotmp = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomotmp = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + !mask = ISHFT(1_8,p) - 1 + !Jsomotmp = IAND(Jsomo,mask) + !pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + !mask = ISHFT(1_8,q) - 1 + !Jsomotmp = IAND(Jsomo,mask) + !qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Jsomotmp = IAND(Jsomo(ii),mask) + tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Jsomotmp = IAND(Jsomo(ii),mask) + tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) case (3) ! SOMO -> VMO !print *,"type -> SOMO -> VMO" !Isomo = IEOR(Isomo,Jsomo) if(p.LT.q) then - mask = ISHFT(1_8,p) - 1 - Isomo = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomo = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + !mask = ISHFT(1_8,p) - 1 + !Isomo = IAND(Isomo,mask) + !pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + !mask = ISHFT(1_8,q) - 1 + !Jsomo = IAND(Jsomo,mask) + !qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Isomotmp = IAND(Isomo(ii),mask) + tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + end do + mask = ISHFT(1_8,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Jsomotmp = IAND(Jsomo(ii),mask) + tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1 else - mask = ISHFT(1_8,p) - 1 - Isomo = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 - mask = ISHFT(1_8,q) - 1 - Jsomo = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + !mask = ISHFT(1_8,p) - 1 + !Isomo = IAND(Isomo,mask) + !pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + !mask = ISHFT(1_8,q) - 1 + !Jsomo = IAND(Jsomo,mask) + !qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Isomotmp = IAND(Isomo(ii),mask) + tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1 + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Jsomotmp = IAND(Jsomo(ii),mask) + tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) endif case (4) ! DOMO -> SOMO @@ -1314,19 +1420,67 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod !print *,"type -> DOMO -> SOMO" !Isomo = IEOR(Isomo,Jsomo) if(p.LT.q) then - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + !mask = ISHFT(1_8,p) - 1 + !Jsomo = IAND(Jsomo,mask) + !pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + !mask = ISHFT(1_8,q) - 1 + !Isomo = IAND(Isomo,mask) + !qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Jsomotmp = IAND(Jsomo(ii),mask) + tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Isomotmp = IAND(Isomo(ii),mask) + tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1 else - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + !mask = ISHFT(1_8,p) - 1 + !Jsomo = IAND(Jsomo,mask) + !pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + !mask = ISHFT(1_8,q) - 1 + !Isomo = IAND(Isomo,mask) + !qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Jsomotmp = IAND(Jsomo(ii),mask) + tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1 + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + mask = ISHFT(1_bit_kind,-1)-1_bit_kind + Isomotmp = IAND(Isomo(ii),mask) + tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) endif case default print *,"something is wrong in convertOrbIdsToModelSpaceIds" @@ -1415,7 +1569,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze allocate(diag_energies(n_CSF)) call calculate_preconditioner_cfg(diag_energies) - !print *," diag energy =",diag_energies(1) + print *," diag energy =",diag_energies(1) MS = 0 norm_coef_cfg=0.d0 @@ -1590,6 +1744,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze call omp_set_lock(lock(jj)) do kk = 1,n_st psi_out(kk,jj) = psi_out(kk,jj) + meCC1 * psi_in(kk,ii) + print *,"jj=",jj,'psi_out(kk)=',psi_out(kk,jj) enddo call omp_unset_lock(lock(jj)) enddo @@ -1666,7 +1821,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze nconnectedI, excitationIds, excitationTypes, diagfactors) !if(i .EQ. 1) then - ! print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' kcfgDOMO=',alphas_Icfg(1,2,k),' ',POPCNT(alphas_Icfg(1,2,k)) + ! print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),alphas_Icfg(2,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' & + ! kcfgDOMO=',alphas_Icfg(1,2,k),alphas_Icfg(2,2,k),' ',POPCNT(alphas_Icfg(1,2,k)), " NconnectedI=",nconnectedI !endif @@ -1682,9 +1838,13 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze NSOMOI = getNSOMO(connectedI_alpha(:,:,j)) p = excitationIds(1,j) q = excitationIds(2,j) + !print *,"j=",j, " p=",p," q=",q extype = excitationTypes(j) call convertOrbIdsToModelSpaceIds(alphas_Icfg(1,1,k), connectedI_alpha(1,1,j), p, q, extype, pmodel, qmodel) ! for E_pp E_rs and E_ppE_rr case + !if(k.eq.722)then + ! print *,"j=",j," k=",k,"p=",p,"q=",q,"NSOMOalpha=",NSOMOalpha, "pmodel=",pmodel,"qmodel=",qmodel, "extype=",extype + !endif 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 @@ -1692,6 +1852,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze rowsTKI = rowsikpq enddo + !if(i.eq.1)then + ! print *,"n_st=",n_st,"rowsTKI=",rowsTKI, " nconnectedI=",nconnectedI, & + ! "totcolsTKI=",totcolsTKI + !endif allocate(TKI(n_st,rowsTKI,totcolsTKI)) ! coefficients of CSF ! Initialize the integral container ! dims : (totcolsTKI, nconnectedI) @@ -1721,10 +1885,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) & * psi_in(kk,idxs_connectedI_alpha(j)+m-1) enddo - !if(i.eq.1) then - ! print *,AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) - !endif enddo + !if(i.eq.1) then + ! print *,"j=",j,"psi_in=",psi_in(1,idxs_connectedI_alpha(j)+m-1) + !endif enddo diagfactors_0 = diagfactors(j)*0.5d0 @@ -1763,16 +1927,24 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze rowsTKI = rowsikpq CCmattmp = 0.d0 + !if(i.eq.1)then + ! print *,"\t n_st=",n_st," colsikpq=",colsikpq," rowsTKI=",rowsTKI,& + ! " | ",size(TKIGIJ,1),size(AIJpqContainer,1),size(CCmattmp,1) + !endif call dgemm('N','N', n_st, colsikpq, rowsTKI, 1.d0, & TKIGIJ(1,1,j), size(TKIGIJ,1), & AIJpqContainer(1,1,pmodel,qmodel,extype,NSOMOalpha), & size(AIJpqContainer,1), 0.d0, & CCmattmp, size(CCmattmp,1) ) + !print *,"j=",j,"colsikpq=",colsikpq, "sizeTIG=",size(TKIGIJ,1),"sizeaijpq=",size(AIJpqContainer,1) do m = 1,colsikpq call omp_set_lock(lock(idxs_connectedI_alpha(j)+m-1)) do kk = 1,n_st psi_out(kk,idxs_connectedI_alpha(j)+m-1) += CCmattmp(kk,m) + !if(dabs(CCmattmp(kk,m)).gt.1e-10)then + ! print *, CCmattmp(kk,m), " | ",idxs_connectedI_alpha(j)+m-1 + !end if enddo call omp_unset_lock(lock(idxs_connectedI_alpha(j)+m-1)) enddo diff --git a/src/davidson/diagonalization_hcfg.irp.f b/src/davidson/diagonalization_hcfg.irp.f index 659602a1..b88c188d 100644 --- a/src/davidson/diagonalization_hcfg.irp.f +++ b/src/davidson/diagonalization_hcfg.irp.f @@ -112,6 +112,8 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N double precision, allocatable :: U(:,:), U_csf(:,:), overlap(:,:) double precision, allocatable :: tmpU(:,:), tmpW(:,:) double precision, pointer :: W(:,:), W_csf(:,:) + double precision, pointer :: W2(:,:), W_csf2(:,:) + double precision, allocatable :: U2(:,:), U_csf2(:,:) logical :: disk_based double precision :: energy_shift(N_st_diag_in*davidson_sze_max) @@ -234,12 +236,15 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N call c_f_pointer(ptr_w, W_csf, (/sze_csf,N_st_diag*itermax/)) else allocate(W(sze,N_st_diag),W_csf(sze_csf,N_st_diag*itermax)) + allocate(W2(sze,N_st_diag),W_csf2(sze_csf,N_st_diag*itermax)) endif allocate( & ! Large U(sze,N_st_diag), & + U2(sze,N_st_diag), & U_csf(sze_csf,N_st_diag*itermax), & + U_csf2(sze_csf,N_st_diag*itermax), & ! Small h(N_st_diag*itermax,N_st_diag*itermax), & @@ -324,8 +329,8 @@ 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 =0.0d0 - !tmpU(1,2)=1.0d0 + tmpU =0.0d0 + tmpU(1,1)=1.0d0 double precision :: irp_rdtsc double precision :: ticks_0, ticks_1 integer*8 :: irp_imax @@ -340,19 +345,19 @@ 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)) - ! 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 + 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)) + !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) deallocate(tmpU) endif