diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c index c11b38dc..a90fbc1f 100644 --- a/src/csf/cfgCI_utils.c +++ b/src/csf/cfgCI_utils.c @@ -341,7 +341,6 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl Get BFtoDeterminant Matrix ************************************/ - printf("In convertcsftodet\n"); convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI); int rowsI = 0; diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 7a84cfc2..80e11cfc 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -40,49 +40,68 @@ detDimperBF = 0 MS = elec_alpha_num-elec_beta_num ! number of cfgs = number of dets for 0 somos - n_CSF = 0 - ncfgprev = cfg_seniority_index(0) - ncfgpersomo = ncfgprev - do i = iand(MS,1), NSOMOMax-2,2 - if(cfg_seniority_index(i) .EQ. -1) then - cycle - endif - if(cfg_seniority_index(i+2) .EQ. -1) then - ncfgpersomo = N_configuration + 1 - else - if(cfg_seniority_index(i+2) > ncfgpersomo) then - ncfgpersomo = cfg_seniority_index(i+2) - else - k = 0 - do while(cfg_seniority_index(i+2+k) < ncfgpersomo) - k = k + 2 - ncfgpersomo = cfg_seniority_index(i+2+k) - enddo - endif - endif - ncfg = ncfgpersomo - ncfgprev - if(i .EQ. 0 .OR. i .EQ. 1) then - dimcsfpercfg = 1 - elseif( i .EQ. 3) then - dimcsfpercfg = 2 - else - if(iand(MS,1) .EQ. 0) then - dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) - else - dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) - endif - endif - n_CSF += ncfg * dimcsfpercfg - if(cfg_seniority_index(i+2) > ncfgprev) then - ncfgprev = cfg_seniority_index(i+2) - else - k = 0 - do while(cfg_seniority_index(i+2+k) < ncfgprev) - k = k + 2 - ncfgprev = cfg_seniority_index(i+2+k) - enddo - endif + n_CSF = cfg_seniority_index(NSOMOMin)-1 + print *,"start=",n_CSF + ncfgprev = cfg_seniority_index(NSOMOMin) + !do i = 0-iand(MS,1)+2, NSOMOMax,2 + do i = NSOMOMin+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 + print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", n_CSF + !if(cfg_seniority_index(i+2) == -1) EXIT + !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF + ncfgprev = cfg_seniority_index(i) enddo + !n_CSF = 0 + !ncfgprev = cfg_seniority_index(0) + !ncfgpersomo = ncfgprev + !do i = iand(MS,1), NSOMOMax-2,2 + ! if(cfg_seniority_index(i) .EQ. -1) then + ! cycle + ! endif + ! if(cfg_seniority_index(i+2) .EQ. -1) then + ! ncfgpersomo = N_configuration + 1 + ! else + ! if(cfg_seniority_index(i+2) > ncfgpersomo) then + ! ncfgpersomo = cfg_seniority_index(i+2) + ! else + ! k = 0 + ! do while(cfg_seniority_index(i+2+k) < ncfgpersomo) + ! k = k + 2 + ! ncfgpersomo = cfg_seniority_index(i+2+k) + ! enddo + ! endif + ! endif + ! ncfg = ncfgpersomo - ncfgprev + ! if(i .EQ. 0 .OR. i .EQ. 1) then + ! dimcsfpercfg = 1 + ! elseif( i .EQ. 3) then + ! dimcsfpercfg = 2 + ! else + ! if(iand(MS,1) .EQ. 0) then + ! dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) + ! else + ! dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) + ! endif + ! endif + ! n_CSF += ncfg * dimcsfpercfg + ! if(cfg_seniority_index(i+2) > ncfgprev) then + ! ncfgprev = cfg_seniority_index(i+2) + ! else + ! k = 0 + ! do while(cfg_seniority_index(i+2+k) < ncfgprev) + ! k = k + 2 + ! ncfgprev = cfg_seniority_index(i+2+k) + ! enddo + ! endif + !enddo END_PROVIDER @@ -1477,14 +1496,14 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze deallocate(excitationIds_single) deallocate(excitationTypes_single) - allocate(listconnectedJ(N_INT,2,max(sze,100))) - allocate(alphas_Icfg(N_INT,2,max(sze,100))) - allocate(connectedI_alpha(N_INT,2,max(sze,100))) - allocate(idxs_connectedI_alpha(max(sze,100))) - allocate(excitationIds(2,max(sze,100))) - allocate(excitationTypes(max(sze,100))) - allocate(diagfactors(max(sze,100))) - allocate(idslistconnectedJ(max(sze,100))) + allocate(listconnectedJ(N_INT,2,max(sze,1000))) + allocate(alphas_Icfg(N_INT,2,max(sze,1000))) + allocate(connectedI_alpha(N_INT,2,max(sze,1000))) + allocate(idxs_connectedI_alpha(max(sze,1000))) + allocate(excitationIds(2,max(sze,1000))) + allocate(excitationTypes(max(sze,1000))) + allocate(diagfactors(max(sze,1000))) + allocate(idslistconnectedJ(max(sze,1000))) allocate(CCmattmp(n_st,NBFmax)) ! Loop over all selected configurations @@ -1643,11 +1662,13 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze do i = 1,n_CSF do kk=1,n_st psi_out(kk,i) += diag_energies(i)*psi_in(kk,i) + print *,psi_in(kk,i)," | ",psi_out(kk,i) enddo enddo !$OMP END DO !$OMP END PARALLEL + stop call omp_set_max_active_levels(4) deallocate(diag_energies) diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index b6d5a3f8..2d3ed8b1 100644 --- a/src/davidson/diagonalization_hcsf_dressed.irp.f +++ b/src/davidson/diagonalization_hcsf_dressed.irp.f @@ -267,8 +267,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N do k=N_st+1,N_st_diag do i=1,sze - call random_number(r1) - call random_number(r2) + !call random_number(r1) + !call random_number(r2) + r1 = 0.5 + r2 = 0.5 r1 = dsqrt(-2.d0*dlog(r1)) r2 = dtwo_pi*r2 u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) @@ -286,8 +288,8 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N enddo ! Make random verctors eigenstates of S2 - call convertWFfromDETtoCSF(N_st_diag,U(1,1),U_csf(1,1)) - call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1)) + call convertWFfromDETtoCSF(N_st_diag,U,U_csf) + !call convertWFfromCSFtoDET(N_st_diag,U_csf,U) do while (.not.converged) itertot = itertot+1 @@ -330,8 +332,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N do kk=1,N_st_diag do ii=1,sze_csf tmpU(kk,ii) = U_csf(ii,shift+kk) + tmpU(kk,ii) = 0.0d0 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 @@ -479,24 +483,24 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N ! Compute residual vector and davidson step ! ----------------------------------------- - if (without_diagonal) then - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) - do k=1,N_st_diag - do i=1,sze - U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) & - /max(H_jj(i) - lambda (k),1.d-2) - enddo + !if (without_diagonal) then + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag + do i=1,sze + U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) & + /max(H_jj(i) - lambda (k),1.d-2) enddo - !$OMP END PARALLEL DO - else - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) - do k=1,N_st_diag - do i=1,sze - U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) - enddo - enddo - !$OMP END PARALLEL DO - endif + enddo + !$OMP END PARALLEL DO + !else + ! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + ! do k=1,N_st_diag + ! do i=1,sze + ! U(i,k) = (lambda(k) * U(i,k) - W(i,k) ) + ! enddo + ! enddo + ! !$OMP END PARALLEL DO + !endif do k=1,N_st residual_norm(k) = u_dot_u(U(1,k),sze) @@ -543,6 +547,15 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N ! Re-contract U ! ------------- + call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0, & + W_csf, size(W_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze_csf + W_csf(i,k) = u_in(i,k) + enddo + enddo + call convertWFfromCSFtoDET(N_st_diag,W_csf,W) + call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0, & U_csf, size(U_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) do k=1,N_st_diag