diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index e14f11bb..d981da79 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1,9 +1,9 @@ - BEGIN_PROVIDER [ integer, NSOMOMax] -&BEGIN_PROVIDER [ integer, NCSFMax] -&BEGIN_PROVIDER [ integer*8, NMO] -&BEGIN_PROVIDER [ integer, NBFMax] -&BEGIN_PROVIDER [ integer, n_CSF] -&BEGIN_PROVIDER [ integer, maxDetDimPerBF] + BEGIN_PROVIDER [ integer, NSOMOMax] + &BEGIN_PROVIDER [ integer, NCSFMax] + &BEGIN_PROVIDER [ integer*8, NMO] + &BEGIN_PROVIDER [ integer, NBFMax] + &BEGIN_PROVIDER [ integer, dimBasisCSF] + &BEGIN_PROVIDER [ integer, maxDetDimPerBF] implicit none BEGIN_DOC ! Documentation for NSOMOMax @@ -22,29 +22,38 @@ integer NSOMO integer dimcsfpercfg integer detDimperBF - real*8 :: coeff + real*8 :: coeff integer MS integer ncfgpersomo detDimperBF = 0 MS = elec_alpha_num-elec_beta_num + !print *,"NSOMOMax=",NSOMOMax, cfg_seniority_index(0) ! number of cfgs = number of dets for 0 somos - n_CSF = cfg_seniority_index(0)-1 + dimBasisCSF = cfg_seniority_index(0)-1 ncfgprev = cfg_seniority_index(0) do i = 0-iand(MS,1)+2, NSOMOMax,2 - if(cfg_seniority_index(i) .EQ. -1)then - ncfgpersomo = N_configuration + 1 - else - ncfgpersomo = cfg_seniority_index(i) - endif - ncfg = ncfgpersomo - ncfgprev - !detDimperBF = max(1,nint((binom(i,(i+1)/2)))) - dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) - n_CSF += ncfg * dimcsfpercfg - !if(cfg_seniority_index(i+2) == -1) EXIT - !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF - ncfgprev = cfg_seniority_index(i) + if(cfg_seniority_index(i) .EQ. -1)then + ncfgpersomo = N_configuration + 1 + else + ncfgpersomo = cfg_seniority_index(i) + endif + ncfg = ncfgpersomo - ncfgprev + !detDimperBF = max(1,nint((binom(i,(i+1)/2)))) + dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) + dimBasisCSF += ncfg * dimcsfpercfg + !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", dimBasisCSF + !if(cfg_seniority_index(i+2) == -1) EXIT + !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF + ncfgprev = cfg_seniority_index(i) enddo -END_PROVIDER + if(NSOMOMax .EQ. elec_num)then + ncfgpersomo = N_configuration + 1 + ncfg = ncfgpersomo - ncfgprev + dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) + dimBasisCSF += ncfg * dimcsfpercfg + !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", dimBasisCSF + endif + END_PROVIDER subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) use bitmasks @@ -99,7 +108,296 @@ end subroutine get_phase_qp_to_cfg - BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,2)] + BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)] + &BEGIN_PROVIDER [ real*8, psi_coef_config, (dimBasisCSF,1)] + &BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)] + &BEGIN_PROVIDER [ integer, psi_csf_to_config_data, (dimBasisCSF)] + use cfunctions + implicit none + BEGIN_DOC + ! Documentation for DetToCSFTransformationMatrix + ! Provides the matrix of transformatons for the + ! conversion between determinant to CSF basis (in BFs) + END_DOC + integer*8 :: Isomo, Idomo, mask, Ialpha,Ibeta + integer :: rows, cols, i, j, k + integer :: startdet, enddet + integer*8 MS + integer ndetI + integer :: getNSOMO + real*8,dimension(:,:),allocatable :: tempBuffer + real*8,dimension(:),allocatable :: tempCoeff + real*8 :: norm_det1, phasedet + norm_det1 = 0.d0 + MS = elec_alpha_num - elec_beta_num + print *,"Maxbfdim=",NBFMax + print *,"Maxdetdim=",maxDetDimPerBF + print *,"dimBasisCSF=",dimBasisCSF + print *,"N_configurations=",N_configuration + print *,"n_core_orb=",n_core_orb + ! initialization + psi_coef_config = 0.d0 + DetToCSFTransformationMatrix(0,:,:) = 1.d0 + do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2 + Isomo = IBSET(0_8, i) - 1_8 + ! rows = Ncsfs + ! cols = Ndets + bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1)))) + ndetI = max(1,nint((binom(i,(i+1)/2)))) + + allocate(tempBuffer(bfIcfg,ndetI)) + call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer) + DetToCSFTransformationMatrix(i,:bfIcfg,:ndetI) = tempBuffer + deallocate(tempBuffer) + enddo + + integer s, bfIcfg + integer countcsf + countcsf = 0 + integer countdet + countdet = 0 + integer istate + istate = 1 + phasedet = 1.0d0 + do i = 1,N_configuration + startdet = psi_configuration_to_psi_det(1,i) + enddet = psi_configuration_to_psi_det(2,i) + ndetI = enddet-startdet+1 + + allocate(tempCoeff(ndetI)) + countdet = 1 + do j = startdet, enddet + Ialpha = psi_det(1,1,psi_configuration_to_psi_det_data(j)) + Ibeta = psi_det(1,2,psi_configuration_to_psi_det_data(j)) + !call debug_spindet(Ialpha,1,1) + !call debug_spindet(Ibeta ,1,1) + call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet) + !print *,">>",Ialpha,Ibeta,phasedet + tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate)*phasedet + !tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate) + norm_det1 += tempCoeff(countdet)*tempCoeff(countdet) + countdet += 1 + enddo + + !print *,"dimcoef=",bfIcfg,norm_det1 + !call printMatrix(tempCoeff,ndetI,1) + + s = 0 + do k=1,N_int + 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)))) + + ! perhaps blocking with CFGs of same seniority + ! can be more efficient + allocate(tempBuffer(bfIcfg,ndetI)) + tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) + !print *,"csftodetdim=",bfIcfg,ndetI + !call printMatrix(tempBuffer,bfIcfg,ndetI) + + call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1,1), size(psi_coef_config,1)) + !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1) + + !call printMatrix(psi_coef_config(countcsf+1,1),bfIcfg,1) + deallocate(tempCoeff) + deallocate(tempBuffer) + psi_config_data(i,1) = countcsf + 1 + countcsf += bfIcfg + psi_config_data(i,2) = countcsf + psi_csf_to_config_data(countcsf) = i + enddo + print *,"Norm det=",norm_det1, size(psi_coef_config,1), " Dim csf=", countcsf + + END_PROVIDER + + subroutine convertWFfromDETtoCSF(psi_coef_det_in, psi_coef_cfg_out) + use cfunctions + implicit none + BEGIN_DOC + ! Documentation for DetToCSFTransformationMatrix + ! Provides the matrix of transformatons for the + ! conversion between determinant to CSF basis (in BFs) + END_DOC + integer*8 :: Isomo, Idomo, mask, Ialpha,Ibeta + integer :: rows, cols, i, j, k + integer :: startdet, enddet + integer*8 MS + integer ndetI + integer :: getNSOMO + real*8,intent(in) :: psi_coef_det_in(n_det,1) + real*8,intent(out) :: psi_coef_cfg_out(dimBasisCSF,1) + real*8,dimension(:,:),allocatable :: tempBuffer + real*8,dimension(:),allocatable :: tempCoeff + real*8 :: norm_det1, phasedet + norm_det1 = 0.d0 + MS = elec_alpha_num - elec_beta_num + print *,"Maxbfdim=",NBFMax + print *,"Maxdetdim=",maxDetDimPerBF + print *,"dimBasisCSF=",dimBasisCSF + print *,"N_configurations=",N_configuration + print *,"n_core_orb=",n_core_orb + ! initialization + psi_coef_cfg_out(:,1) = 0.d0 + + integer s, bfIcfg + integer countcsf + countcsf = 0 + integer countdet + countdet = 0 + integer istate + istate = 1 + phasedet = 1.0d0 + do i = 1,N_configuration + startdet = psi_configuration_to_psi_det(1,i) + enddet = psi_configuration_to_psi_det(2,i) + ndetI = enddet-startdet+1 + + allocate(tempCoeff(ndetI)) + countdet = 1 + do j = startdet, enddet + Ialpha = psi_det(1,1,psi_configuration_to_psi_det_data(j)) + Ibeta = psi_det(1,2,psi_configuration_to_psi_det_data(j)) + call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet) + !print *,">>",Ialpha,Ibeta,phasedet + tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate)*phasedet + !tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate) + norm_det1 += tempCoeff(countdet)*tempCoeff(countdet) + countdet += 1 + enddo + + !print *,"dimcoef=",bfIcfg,norm_det1 + !call printMatrix(tempCoeff,ndetI,1) + + s = 0 + do k=1,N_int + 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)))) + + ! perhaps blocking with CFGs of same seniority + ! can be more efficient + allocate(tempBuffer(bfIcfg,ndetI)) + tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) + !print *,"csftodetdim=",bfIcfg,ndetI + !call printMatrix(tempBuffer,bfIcfg,ndetI) + + call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_cfg_out(countcsf+1,1), size(psi_coef_cfg_out,1)) + + deallocate(tempCoeff) + deallocate(tempBuffer) + psi_config_data(i,1) = countcsf + 1 + countcsf += bfIcfg + psi_config_data(i,2) = countcsf + enddo + print *,"Norm det=",norm_det1, size(psi_coef_cfg_out,1), " Dim csf=", countcsf + + end subroutine convertWFfromDETtoCSF + + subroutine convertWFfromCSFtoDET(psi_coef_cfg_in, psi_coef_det) + implicit none + BEGIN_DOC + ! Documentation for convertCSFtoDET + ! This function converts the wavefunction + ! in CFG basis to DET basis using the + ! transformation matrix provided before. + END_DOC + real*8,intent(in) :: psi_coef_cfg_in(dimBasisCSF,1) + real*8,intent(out) :: psi_coef_det(N_det,1) + real*8 :: tmp_psi_coef_det(maxDetDimPerBF) + integer s, bfIcfg + integer countcsf + integer countdet + integer*8 :: Isomo, Idomo, Ialpha, Ibeta + integer :: rows, cols, i, j, k + integer :: startdet, enddet + integer*8 MS + integer ndetI + integer :: getNSOMO + real*8,dimension(:,:),allocatable :: tempBuffer + real*8,dimension(:),allocatable :: tempCoeff + real*8 :: phasedet + ! number of states + integer istate + istate = 1 + countcsf = 1 + countdet = 1 + print *,"in function convertWFfromCSFtoDET()" + + + do i = 1,N_configuration + startdet = psi_configuration_to_psi_det(1,i) + enddet = psi_configuration_to_psi_det(2,i) + ndetI = enddet-startdet+1 + + s = 0 + do k=1,N_int + 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)))) + + allocate(tempCoeff(bfIcfg)) + + do j = 1,bfIcfg + tempCoeff(j) = psi_coef_cfg_in(countcsf,1) + countcsf += 1 + enddo + !print *,"dimcoef=",bfIcfg + !call printMatrix(tempCoeff,bfIcfg,1) + + ! perhaps blocking with CFGs of same seniority + ! can be more efficient + allocate(tempBuffer(bfIcfg,ndetI)) + tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) + !print *,"csftodetdim=",bfIcfg,ndetI + !call printMatrix(tempBuffer,bfIcfg,ndetI) + + !call dgemm('T','N', ndetI, 1, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_det(countdet,1), size(psi_coef_det,1)) + call dgemm('T','N', ndetI, 1, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, tmp_psi_coef_det, size(tmp_psi_coef_det,1)) + + !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf,1), 1) + + !print *,"result" + !call printMatrix(tmp_psi_coef_det,ndetI,1) + + countdet = 1 + do j=startdet,enddet + Ialpha = psi_det(1,1,psi_configuration_to_psi_det_data(j)) + Ibeta = psi_det(1,2,psi_configuration_to_psi_det_data(j)) + !call debug_spindet(Ialpha,1,1) + !call debug_spindet(Ibeta ,1,1) + call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet) + !print *,">>",Ialpha,Ibeta,phasedet + psi_coef_det(psi_configuration_to_psi_det_data(j),1) = tmp_psi_coef_det(countdet)*phasedet + countdet += 1 + enddo + + deallocate(tempCoeff) + deallocate(tempBuffer) + !countdet += ndetI + enddo + + !countdet = 1 + !tmp_psi_coef_det = psi_coef_det(:,1) + !do i=1,N_configuration + ! startdet = psi_configuration_to_psi_det(1,i) + ! enddet = psi_configuration_to_psi_det(2,i) + ! ndetI = enddet-startdet+1 + ! print *,i,">>>",startdet,enddet + ! do k=1,ndetI + ! !psi_coef_det(startdet+k-1,1) = tmp_psi_coef_det(countdet) + ! psi_coef_det(countdet,1) = tmp_psi_coef_det(startdet+k-1) + ! countdet += 1 + ! enddo + !enddo + + print *,"End ncsfs=",countcsf + + end subroutine convertCSFtoDET + + BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)] &BEGIN_PROVIDER [ integer, rowsmax] &BEGIN_PROVIDER [ integer, colsmax] use cfunctions @@ -122,9 +420,13 @@ end subroutine get_phase_qp_to_cfg nsomomin = elec_alpha_num-elec_beta_num rowsmax = 0 colsmax = 0 + print *,"NSOMOMax = ",NSOMOMax !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) ! Type ! 1. SOMO -> SOMO + !print *,"Doing SOMO->SOMO" + AIJpqMatrixDimsList(0,0,1,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,1,1,1,2) = 1 do i = 2-iand(nsomomin,1), NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i-2,i-2, 2 @@ -152,6 +454,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, "SOMO->SOMO \t",i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -167,6 +470,9 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 2. DOMO -> VMO + !print *,"Doing DOMO->VMO" + AIJpqMatrixDimsList(0,0,2,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,2,1,1,2) = 1 do i = 0+iand(nsomomin,1), NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 @@ -200,6 +506,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -216,6 +523,8 @@ end subroutine get_phase_qp_to_cfg ! Type ! 3. SOMO -> VMO !print *,"Doing SOMO->VMO" + AIJpqMatrixDimsList(0,0,3,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,3,1,1,2) = 1 do i = 2-iand(nsomomin,1), NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 @@ -223,8 +532,8 @@ end subroutine get_phase_qp_to_cfg if(j .GT. NSOMOMax .OR. j .LE. 0) then cycle end if - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 Isomo = IBCLR(Isomo,l-1) @@ -239,6 +548,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -253,18 +563,20 @@ end subroutine get_phase_qp_to_cfg end do end do ! Type - ! 4. DOMO -> VMO + ! 4. DOMO -> SOMO !print *,"Doing DOMO->SOMO" + AIJpqMatrixDimsList(0,0,4,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,4,1,1,2) = 1 do i = 2-iand(nsomomin,1), NSOMOMax, 2 do j = i,i, 2 if(j .GT. NSOMOMax .OR. j .LE. 0) then cycle end if - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 - Isomo = IBCLR(Isomo,k+1-1) + Isomo = IBCLR(Isomo,k-1) Jsomo = ISHFT(1_8,j+1)-1 Jsomo = IBCLR(Jsomo,l-1) else @@ -276,6 +588,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -289,9 +602,10 @@ end subroutine get_phase_qp_to_cfg end do end do end do + print *,"Rowsmax=",rowsmax," Colsmax=",colsmax END_PROVIDER - BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,NBFMax,NBFMax)] + BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,NBFMax,NBFMax)] use cfunctions implicit none BEGIN_DOC @@ -327,13 +641,18 @@ end subroutine get_phase_qp_to_cfg integer maxdim !maxdim = max(rowsmax,colsmax) ! allocate matrix + !print *,"rowsmax =",rowsmax," colsmax=",colsmax + !print *,"NSOMOMax = ",NSOMOMax !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) ! Type ! 1. SOMO -> SOMO + !print *,"Doing SOMO -> SOMO" + AIJpqContainer(0,0,1,1,1,1,1) = 1.0d0 do i = 2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i-2,i-2, 2 if(j .GT. NSOMOMax .OR. j .LT. 0) cycle + !print *,"i,j=",i,j do k = 1,i do l = 1,i @@ -350,6 +669,10 @@ end subroutine get_phase_qp_to_cfg nsomoj = i endif + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + AIJpqContainer(nsomoi,nsomoj,1,k,l,:,:) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & @@ -371,6 +694,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax + !call printMatrix(meMatrix,rows,cols) ! i -> j do ri = 1,rows do ci = 1,cols @@ -384,6 +709,8 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 2. DOMO -> VMO + !print *,"Doing DOMO -> VMO" + AIJpqContainer(0,0,2,1,1,1,1) = 1.0d0 do i = 0, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 @@ -409,6 +736,10 @@ end subroutine get_phase_qp_to_cfg nsomoj = j endif + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + AIJpqContainer(nsomoi,nsomoj,2,k,l,:,:) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & @@ -430,6 +761,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax + !call printMatrix(meMatrix,rows,cols) ! i -> j do ri = 1,rows do ci = 1,cols @@ -443,13 +776,14 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 3. SOMO -> VMO + !print *,"Doing SOMO -> VMO" do i = 2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 Jsomo = ISHFT(1_8,j)-1 if(j .GT. NSOMOMax .OR. j .LE. 0) cycle - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 Isomo = IBCLR(Isomo,l-1) @@ -460,6 +794,10 @@ end subroutine get_phase_qp_to_cfg Jsomo = ISHFT(1_8,j)-1 endif + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + AIJpqContainer(i,j,3,k,l,:,:) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & @@ -481,6 +819,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !call printMatrix(meMatrix,rows,cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax ! i -> j do ri = 1,rows do ci = 1,cols @@ -494,18 +834,20 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 4. DOMO -> SOMO + !print *,"Doing DOMO -> SOMO" + AIJpqContainer(0,0,4,1,1,1,1) = 1.0d0 do i = 2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 Jsomo = ISHFT(1_8,i)-1 if(j .GT. NSOMOMax .OR. j .LE. 0) cycle - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 Isomo = IBCLR(Isomo,k-1) Jsomo = ISHFT(1_8,j+1)-1 - Jsomo = IBCLR(Jsomo,l+1-1) + Jsomo = IBCLR(Jsomo,l-1) else Isomo = ISHFT(1_8,i)-1 Jsomo = ISHFT(1_8,j)-1 @@ -533,6 +875,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !call printMatrix(meMatrix,rows,cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax ! i -> j do ri = 1,rows do ci = 1,cols @@ -545,94 +889,3 @@ end subroutine get_phase_qp_to_cfg end do end do END_PROVIDER - - -!!!!!! - - BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)] - &BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF)] - &BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)] - use cfunctions - use bitmasks - implicit none - BEGIN_DOC - ! Documentation for DetToCSFTransformationMatrix - ! Provides the matrix of transformatons for the - ! conversion between determinant to CSF basis (in BFs) - END_DOC - integer(bit_kind) :: mask(N_int), Ialpha(N_int),Ibeta(N_int) - integer :: rows, cols, i, j, k - integer :: startdet, enddet - integer*8 MS, Isomo, Idomo - integer ndetI - integer :: getNSOMO - real*8,dimension(:,:),allocatable :: tempBuffer - real*8,dimension(:),allocatable :: tempCoeff - real*8 :: norm_det1, phasedet - norm_det1 = 0.d0 - MS = elec_alpha_num - elec_beta_num - ! initialization - psi_coef_config = 0.d0 - DetToCSFTransformationMatrix(0,:,:) = 1.d0 - do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2 - Isomo = IBSET(0_8, i) - 1_8 - ! rows = Ncsfs - ! cols = Ndets - bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1)))) - ndetI = max(1,nint((binom(i,(i+1)/2)))) - - allocate(tempBuffer(bfIcfg,ndetI)) - call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer) - DetToCSFTransformationMatrix(i,1:bfIcfg,1:ndetI) = tempBuffer(1:bfIcfg,1:ndetI) - deallocate(tempBuffer) - enddo - - integer s, bfIcfg - integer countcsf - countcsf = 0 - integer countdet - countdet = 0 - integer idx - integer istate - istate = 1 - phasedet = 1.0d0 - do i = 1,N_configuration - startdet = psi_configuration_to_psi_det(1,i) - enddet = psi_configuration_to_psi_det(2,i) - ndetI = enddet-startdet+1 - - allocate(tempCoeff(ndetI)) - countdet = 1 - do j = startdet, enddet - idx = psi_configuration_to_psi_det_data(j) - Ialpha(:) = psi_det(:,1,idx) - Ibeta(:) = psi_det(:,2,idx) - call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet) - tempCoeff(countdet) = psi_coef(idx, istate)*phasedet - norm_det1 += tempCoeff(countdet)*tempCoeff(countdet) - countdet += 1 - enddo - - s = 0 - do k=1,N_int - 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)))) - - ! perhaps blocking with CFGs of same seniority - ! can be more efficient - allocate(tempBuffer(bfIcfg,ndetI)) - tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) - - call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1), size(psi_coef_config,1)) - !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1) - - deallocate(tempCoeff) - deallocate(tempBuffer) - psi_config_data(i,1) = countcsf + 1 - countcsf += bfIcfg - psi_config_data(i,2) = countcsf - enddo - - END_PROVIDER diff --git a/src/csf/tree_utils.h b/src/csf/tree_utils.h index eae19c1b..049b09db 100644 --- a/src/csf/tree_utils.h +++ b/src/csf/tree_utils.h @@ -22,7 +22,8 @@ struct bin_tree { int NBF; }; -#include "/export/apps/pgi/linux86-64/18.10/include/cblas.h" +//#include "/export/apps/pgi/linux86-64/18.10/include/cblas.h" +#include "cblas.h" #define MAX_SOMO 32