From 85db84df23cf98161deb7238facae6b160cdee8a Mon Sep 17 00:00:00 2001 From: v1j4y Date: Wed, 14 Sep 2022 15:07:34 +0200 Subject: [PATCH] Triplet seems to be working. Doublet needs work. --- src/csf/cfgCI_utils.c | 65 +++++++++++++------ src/csf/conversion.irp.f | 1 + src/csf/sigma_vector.irp.f | 27 +++++--- src/csf/tree_utils.c | 4 ++ src/davidson/diagonalization_hcfg.irp.f | 8 ++- .../diagonalization_hs2_dressed.irp.f | 10 ++- 6 files changed, 82 insertions(+), 33 deletions(-) diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c index 64593159..ad3af92e 100644 --- a/src/csf/cfgCI_utils.c +++ b/src/csf/cfgCI_utils.c @@ -1,5 +1,6 @@ #include #include +#include #include "tree_utils.h" void int_to_bin_digit(int64_t in, int count, int* out) @@ -28,7 +29,7 @@ void getncsfs1(int *inpnsomo, int *inpms, int *outncsfs){ int nsomo = *inpnsomo; int ms = *inpms; int nparcoupl = (nsomo + ms)/2; - *outncsfs = binom(nsomo, nparcoupl); + *outncsfs = binom((double)nsomo, (double)nparcoupl); } void getncsfs(int NSOMO, int MS, int *outncsfs){ @@ -39,8 +40,8 @@ void getncsfs(int NSOMO, int MS, int *outncsfs){ (*outncsfs) = 1; return; } - tmpndets = binom(NSOMO, nparcoupl); - (*outncsfs) = round(tmpndets - binom(NSOMO, nparcouplp1)); + tmpndets = binom((double)NSOMO, (double)nparcoupl); + (*outncsfs) = round(tmpndets - binom((double)NSOMO, (double)nparcouplp1)); } #include @@ -1392,14 +1393,23 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *rows, int *cols){ int NSOMO=0; + //printf("before getSetBits Isomo=%ld, NSOMO=%ld\n",Isomo,NSOMO); getSetBits(Isomo, &NSOMO); + //printf("Isomo=%ld, NSOMO=%ld\n",Isomo,NSOMO); int ndets = 0; int NBF = 0; - double dNSOMO = NSOMO*1.0; + //double dNSOMO = NSOMO*1.0; // MS = alpha_num - beta_num - double nalpha = (NSOMO + MS)/2.0; + int nalpha = (NSOMO + MS)/2; //printf(" in convertbftodet : MS=%d nalpha=%3.2f\n",MS,nalpha); - ndets = (int)binom(dNSOMO, nalpha); + //ndets = (int)binom(dNSOMO, nalpha); + if(NSOMO > 0){ + ndets = (int)binom((double)NSOMO, (double)nalpha); + } + else if(NSOMO == 0){ + ndets = 1; + } + else printf("Something is wrong in calcMEdetpair\n"); Tree dettree = (Tree){ .rootNode = NULL, .NBF = -1 }; dettree.rootNode = malloc(sizeof(Node)); @@ -1420,16 +1430,6 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int * } else{ - //int addr = -1; - //int inpdet[NSOMO]; - //inpdet[0] = 1; - //inpdet[1] = 1; - //inpdet[2] = 1; - //inpdet[3] = 0; - //inpdet[4] = 0; - //inpdet[5] = 0; - - //findAddofDetDriver(&dettree, NSOMO, inpdet, &addr); int detlist[ndets]; getDetlistDriver(&dettree, NSOMO, detlist); @@ -1442,6 +1442,9 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int * generateAllBFs(Isomo, MS, &bftree, &NBF, &NSOMO); // Initialize transformation matrix + //printf("MS=%d NBF=%d ndets=%d NSOMO=%d\n",MS,NBF,ndets,NSOMO); + assert( NBF > 0); + assert( ndets > 0); (*bftodetmatrixptr) = malloc(NBF*ndets*sizeof(double)); (*rows) = NBF; (*cols) = ndets; @@ -1496,9 +1499,10 @@ void convertBFtoDetBasisWithArrayDims(int64_t Isomo, int MS, int rowsmax, int co getSetBits(Isomo, &NSOMO); int ndets = 0; int NBF = 0; - double dNSOMO = NSOMO*1.0; - double nalpha = (NSOMO + MS)/2.0; - ndets = (int)binom(dNSOMO, nalpha); + //double dNSOMO = NSOMO*1.0; + //double nalpha = (NSOMO + MS)/2.0; + int nalpha = (NSOMO + MS)/2; + ndets = (int)binom((double)NSOMO, (double)nalpha); Tree dettree = (Tree){ .rootNode = NULL, .NBF = -1 }; dettree.rootNode = malloc(sizeof(Node)); @@ -1582,6 +1586,7 @@ void getApqIJMatrixDims(int64_t Isomo, int64_t Jsomo, int64_t MS, int32_t *rowso getncsfs(NSOMOJ, MS, &NBFJ); (*rowsout) = NBFI; (*colsout) = NBFJ; + //exit(0); } void getApqIJMatrixDriver(int64_t Isomo, int64_t Jsomo, int orbp, int orbq, int64_t MS, int64_t NMO, double **CSFICSFJApqIJptr, int *rowsout, int *colsout){ @@ -1700,6 +1705,7 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in int rowsbftodetI, colsbftodetI; + //printf(" 1Calling convertBFtoDetBasis Isomo=%ld MS=%ld\n",Isomo,MS); convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI); // Fill matrix @@ -1707,7 +1713,14 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in int colsI = 0; //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO); + //printf("Isomo=%ld\n",Isomo); getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO); + if(Isomo == 0){ + rowsI = 1; + colsI = 1; + } + + //printf("Isomo=%ld\n",Isomo); orthoMatrixI = malloc(rowsI*colsI*sizeof(double)); @@ -1719,6 +1732,7 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in int rowsbftodetJ, colsbftodetJ; + //printf(" 2Calling convertBFtoDetBasis Jsomo=%ld MS=%ld\n",Jsomo,MS); convertBFtoDetBasis(Jsomo, MS, &bftodetmatrixJ, &rowsbftodetJ, &colsbftodetJ); int rowsJ = 0; @@ -1726,6 +1740,10 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in // Fill matrix //getOverlapMatrix(Jsomo, MS, &overlapMatrixJ, &rowsJ, &colsJ, &NSOMO); getOverlapMatrix_withDet(bftodetmatrixJ, rowsbftodetJ, colsbftodetJ, Jsomo, MS, &overlapMatrixJ, &rowsJ, &colsJ, &NSOMO); + if(Jsomo == 0){ + rowsJ = 1; + colsJ = 1; + } orthoMatrixJ = malloc(rowsJ*colsJ*sizeof(double)); @@ -1743,18 +1761,25 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in int transA=false; int transB=false; + //printf("1Calling blas\n"); + //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsbftodetI,colsbftodetI,rowsA,colsA); callBlasMatxMat(bftodetmatrixI, rowsbftodetI, colsbftodetI, ApqIJ, rowsA, colsA, bfIApqIJ, transA, transB); + //printf("done\n"); // now transform I in csf basis double *CSFIApqIJ = malloc(rowsI*colsA*sizeof(double)); transA = false; transB = false; + //printf("2Calling blas\n"); + //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,colsI,colsI,colsA); callBlasMatxMat(orthoMatrixI, rowsI, colsI, bfIApqIJ, colsI, colsA, CSFIApqIJ, transA, transB); // now transform J in BF basis double *CSFIbfJApqIJ = malloc(rowsI*rowsbftodetJ*sizeof(double)); transA = false; transB = true; + //printf("3Calling blas\n"); + //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,colsA,rowsbftodetJ,colsbftodetJ); callBlasMatxMat(CSFIApqIJ, rowsI, colsA, bftodetmatrixJ, rowsbftodetJ, colsbftodetJ, CSFIbfJApqIJ, transA, transB); // now transform J in CSF basis @@ -1765,6 +1790,8 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in double *tmpCSFICSFJApqIJ = malloc(rowsI*rowsJ*sizeof(double)); transA = false; transB = true; + //printf("4Calling blas\n"); + //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,rowsbftodetJ,rowsJ,colsJ); callBlasMatxMat(CSFIbfJApqIJ, rowsI, rowsbftodetJ, orthoMatrixJ, rowsJ, colsJ, tmpCSFICSFJApqIJ, transA, transB); // Transfer to actual buffer in Fortran order for(int i = 0; i < rowsI; i++) diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index 9ddd4fee..494c3bfa 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -68,6 +68,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) salpha = (s + MS)/2 bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1)))) else + salpha = (s + MS)/2 bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1)))) endif diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 1c2357b2..8f6e27e8 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -23,7 +23,13 @@ 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) + print *,'cfg_nsomo_min=',cfg_nsomo_min + print *,'cfg_nsomo_max=',cfg_nsomo_max + if(AND(cfg_nsomo_min , 1) .eq. 0)then + NSOMOMin = max(0,cfg_nsomo_min-2) + else + NSOMOMin = max(1,cfg_nsomo_min-2) + endif ! Note that here we need NSOMOMax + 2 sizes ialpha = (NSOMOMax + MS)/2 NCSFMax = max(1,nint((binom(NSOMOMax,ialpha)-binom(NSOMOMax,ialpha+1)))) ! TODO: NCSFs for MS=0 (CHECK) @@ -552,7 +558,7 @@ end subroutine get_phase_qp_to_cfg rows = -1 cols = -1 integer*8 MS - MS = 0 + MS = elec_alpha_num-elec_beta_num real*8,dimension(:,:),allocatable :: meMatrix integer maxdim @@ -931,7 +937,7 @@ subroutine calculate_preconditioner_cfg(diag_energies) noccp = holetype(k) - ! core-active + ! core-virtual do l = 1, n_core_orb jj = list_core(l) core_act_contrib += noccp * (2.d0 * mo_two_e_integrals_jj(jj,p) - mo_two_e_integrals_jj_exchange(jj,p)) @@ -1387,6 +1393,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) MS = 0 norm_coef_cfg=0.d0 @@ -1549,12 +1556,15 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze !meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* h_core_ri(p,q) core_act_contrib = 0.0d0 if(p.ne.q)then - do pp=1,n_core_orb - n=list_core(pp) - core_act_contrib += 2.d0 * get_two_e_integral(p,n,q,n,mo_integrals_map) - get_two_e_integral(p,n,n,q,mo_integrals_map) - end do + do pp=1,n_core_orb + n=list_core(pp) + core_act_contrib += 2.d0 * get_two_e_integral(p,n,q,n,mo_integrals_map) - get_two_e_integral(p,n,n,q,mo_integrals_map) + end do endif meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* (h_act_ri(p,q) + core_act_contrib) + if(jj.eq.1.and.ii.eq.1)then + print *,"CC=",AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI), " p=",p," q=",q + endif call omp_set_lock(lock(jj)) do kk = 1,n_st psi_out(kk,jj) = psi_out(kk,jj) + meCC1 * psi_in(kk,ii) @@ -1578,7 +1588,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,5)=",psi_out(1,5) + !print *," singles part psi(1,1)=",psi_out(1,1) allocate(listconnectedJ(N_INT,2,max(sze,10000))) allocate(alphas_Icfg(N_INT,2,max(sze,10000))) @@ -1751,6 +1761,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze deallocate(diagfactors) !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 + !print *," psi(1,1)=",psi_out(1,1) ! Add the diagonal contribution !$OMP DO diff --git a/src/csf/tree_utils.c b/src/csf/tree_utils.c index 26cb15e0..91bae8b2 100644 --- a/src/csf/tree_utils.c +++ b/src/csf/tree_utils.c @@ -1,3 +1,4 @@ +#include #include "tree_utils.h" void buildTree(Tree *bftree, @@ -52,6 +53,7 @@ void buildTreeDriver(Tree *bftree, int NSOMO, int MS, int *NBF){ int icpl = 0; // keep track of the ith ms (cannot be -ve) int addr = 0; // Counts the total BF's + assert(bftree->rootNode->addr == 0); buildTree(bftree, &(bftree->rootNode), isomo, izeros, icpl, NSOMO, MS); *NBF = bftree->rootNode->addr; @@ -264,6 +266,8 @@ void genDetBasis(Tree *dettree, int Isomo, int MS, int *ndets){ int NSOMO=0; getSetBits(Isomo, &NSOMO); genDetsDriver(dettree, NSOMO, MS, ndets); + // Closed shell case + if(NSOMO==0) (*ndets) = 1; } diff --git a/src/davidson/diagonalization_hcfg.irp.f b/src/davidson/diagonalization_hcfg.irp.f index 08645b53..7da2e1d2 100644 --- a/src/davidson/diagonalization_hcfg.irp.f +++ b/src/davidson/diagonalization_hcfg.irp.f @@ -329,10 +329,12 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N double precision :: ticks_0, ticks_1 integer*8 :: irp_imax irp_imax = 1 - ticks_0 = irp_rdtsc() + !ticks_0 = irp_rdtsc() call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1) - ticks_1 = irp_rdtsc() - print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----" + !ticks_1 = irp_rdtsc() + !print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----" + !print *,' tmpW(1,1)=',tmpW(1,1) + !stop do kk=1,N_st_diag do ii=1,sze_csf W_csf(ii,shift+kk)=tmpW(kk,ii) diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 2a4f53d4..b554754d 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -360,10 +360,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ double precision :: ticks_0, ticks_1 integer*8 :: irp_imax irp_imax = 1 - ticks_0 = irp_rdtsc() + !ticks_0 = irp_rdtsc() + !U = 0d0 + !U(1,1)=1.0d0 call H_S2_u_0_nstates_openmp(W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze) - ticks_1 = irp_rdtsc() - print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----" + !ticks_1 = irp_rdtsc() + !print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----" + !print *,W(1,1) + !stop endif S(1:sze,shift+1:shift+N_st_diag) = real(S_d(1:sze,1:N_st_diag)) ! else