diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c index be4653a3..d55da188 100644 --- a/src/csf/cfgCI_utils.c +++ b/src/csf/cfgCI_utils.c @@ -342,7 +342,11 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl ************************************/ + //printf(" --- In convet ----\n"); convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI); + //printf(" --- done bf det basis ---- row=%d col=%d\n",rowsbftodetI,colsbftodetI); + + //printRealMatrix(bftodetmatrixI,rowsbftodetI,colsbftodetI); int rowsI = 0; int colsI = 0; @@ -350,6 +354,8 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO); getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO); + //printf("Overlap matrix\n"); + //printRealMatrix(overlapMatrixI,rowsI,colsI); /*********************************** Get Orthonormalization Matrix @@ -359,6 +365,9 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl gramSchmidt(overlapMatrixI, rowsI, colsI, orthoMatrixI); + //printf("Ortho matrix\n"); + //printRealMatrix(orthoMatrixI,rowsI,colsI); + /*********************************** Get Final CSF to Det Matrix ************************************/ @@ -1367,7 +1376,9 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int * int ndets = 0; int NBF = 0; double dNSOMO = NSOMO*1.0; + // MS = alpha_num - beta_num double nalpha = (NSOMO + MS)/2.0; + //printf(" in convertbftodet : MS=%d nalpha=%3.2f\n",MS,nalpha); ndets = (int)binom(dNSOMO, nalpha); Tree dettree = (Tree){ .rootNode = NULL, .NBF = -1 }; diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index 75f6e539..9d6e5219 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -24,7 +24,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) double precision, intent(out) :: psi_coef_cfg_out(n_CSF,N_st) integer*8 :: Isomo, Idomo, mask integer(bit_kind) :: Ialpha(N_int) ,Ibeta(N_int) - integer :: rows, cols, i, j, k + integer :: rows, cols, i, j, k, salpha integer :: startdet, enddet integer :: ndetI integer :: getNSOMO @@ -65,13 +65,15 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) enddo if(iand(s,1) .EQ. 0) then - bfIcfg = max(1,nint((binom(s,s/2)-binom(s,(s/2)+1)))) + salpha = (s + MS)/2 + bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1)))) else - bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) + bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1)))) endif ! 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) @@ -99,7 +101,7 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) double precision,intent(in) :: psi_coef_cfg_in(n_CSF,N_st) double precision,intent(out) :: psi_coef_det(N_det,N_st) double precision :: tmp_psi_coef_det(maxDetDimPerBF,N_st) - integer :: s, bfIcfg + integer :: s, bfIcfg, salpha integer :: countcsf integer(bit_kind) :: Ialpha(N_int), Ibeta(N_int) integer :: rows, cols, i, j, k @@ -110,6 +112,8 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) double precision,allocatable :: tempCoeff (:,:) double precision :: phasedet integer :: idx + integer MS + MS = elec_alpha_num-elec_beta_num countcsf = 0 @@ -123,7 +127,8 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) 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)))) allocate(tempCoeff(bfIcfg,N_st)) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index f631d1a8..60c4631e 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -20,12 +20,15 @@ ! The maximum number of SOMOs for the current calculation. ! required for the calculation of prototype arrays. END_DOC + integer MS, ialpha + MS = elec_alpha_num-elec_beta_num NSOMOMax = min(elec_num, cfg_nsomo_max + 2) NSOMOMin = max(0,cfg_nsomo_min-2) ! Note that here we need NSOMOMax + 2 sizes - NCSFMax = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2)-binom(NSOMOMax,((NSOMOMax+1)/2)+1)))) ! TODO: NCSFs for MS=0 + ialpha = (NSOMOMax + MS)/2 + NCSFMax = max(1,nint((binom(NSOMOMax,ialpha)-binom(NSOMOMax,ialpha+1)))) ! TODO: NCSFs for MS=0 (CHECK) NBFMax = NCSFMax - maxDetDimPerBF = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2)))) + maxDetDimPerBF = max(1,nint((binom(NSOMOMax,ialpha)))) NMO = n_act_orb integer i,j,k,l integer startdet,enddet @@ -34,11 +37,9 @@ integer dimcsfpercfg integer detDimperBF real*8 :: coeff, binom1, binom2 - integer MS integer ncfgpersomo real*8, external :: logabsgamma detDimperBF = 0 - MS = elec_alpha_num-elec_beta_num ! number of cfgs = number of dets for 0 somos n_CSF = cfg_seniority_index(NSOMOMin)-1 ncfgprev = cfg_seniority_index(NSOMOMin) @@ -105,7 +106,11 @@ ! endif !enddo n_CSF = 0 - ncfgprev = cfg_seniority_index(0) ! should be 1 + print *," -9(((((((((((((( NSOMOMin=",NSOMOMin + ncfgprev = cfg_seniority_index(NSOMOMin) ! can be -1 + if(ncfgprev.eq.-1)then + ncfgprev=1 + endif do i=NSOMOMin,NSOMOMax+2,2 !k=0 !do while((cfg_seniority_index(i+2+k) .eq. -1) .and. (k.le.NSOMOMax)) @@ -126,9 +131,11 @@ dimcsfpercfg = 2 else if(iand(MS,1) .EQ. 0) then - dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) + ialpha = (i + MS)/2 + dimcsfpercfg = max(1,nint((binom(i,ialpha)-binom(i,ialpha+1)))) else - dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) + ialpha = (i + MS)/2 + dimcsfpercfg = max(1,nint((binom(i,ialpha)-binom(i,ialpha+1)))) endif endif n_CSF += ncfg*dimcsfpercfg @@ -212,7 +219,7 @@ end subroutine get_phase_qp_to_cfg integer(bit_kind) :: Ialpha(N_int),Ibeta(N_int) integer :: rows, cols, i, j, k integer :: startdet, enddet, idx - integer*8 MS + integer*8 MS, salpha integer ndetI integer :: getNSOMO real*8,dimension(:,:),allocatable :: tempBuffer @@ -231,8 +238,11 @@ end subroutine get_phase_qp_to_cfg 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)))) + salpha = (i+MS)/2 + bfIcfg = max(1,nint((binom(i,salpha)-binom(i,salpha+1)))) + ndetI = max(1,nint((binom(i,salpha)))) + !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) @@ -1364,6 +1374,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 allocate(lock(sze)) do i=1,sze call omp_init_lock(lock(i)) diff --git a/src/csf/tree_utils.c b/src/csf/tree_utils.c index 1266b890..26cb15e0 100644 --- a/src/csf/tree_utils.c +++ b/src/csf/tree_utils.c @@ -311,3 +311,13 @@ void callBlasMatxMat(double *A, int rowA, int colA, double *B, int rowB, int col break; } } + +void printRealMatrix(double *orthoMatrix, int rows, int cols){ + int i,j; + for(i=0;i