From 77949a12a05e76f636e3a94314b5ab04b31eb668 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Fri, 12 Mar 2021 19:24:02 +0100 Subject: [PATCH] Added dgemv in sigma vector. --- src/csf/configuration_CI_sigma_helpers.irp.f | 12 +-- src/csf/configurations.irp.f | 9 ++- src/csf/sigma_vector.irp.f | 84 +++++++++++++------- 3 files changed, 69 insertions(+), 36 deletions(-) diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index 0e06817a..40490a33 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -114,14 +114,14 @@ use bitmasks Jsomo = IBCLR(Isomo,p-1) Jsomo = IBSET(Jsomo,q-1) Jdomo = Idomo - kstart = max(0,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) kend = idxI-1 else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then ! SOMO -> SOMO Jsomo = IBCLR(Isomo,p-1) Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBSET(Idomo,q-1) - kstart = max(0,cfg_seniority_index(max(0,Nsomo_I-4))) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) kend = idxI-1 else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then ! DOMO -> VMO @@ -136,7 +136,7 @@ use bitmasks Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBCLR(Idomo,p-1) Jdomo = IBSET(Jdomo,q-1) - kstart = max(0,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) kend = idxI-1 else print*,"Something went wrong in obtain_associated_alphaI" @@ -386,14 +386,14 @@ END_PROVIDER Jsomo = IBCLR(Isomo,p-1) Jsomo = IBSET(Jsomo,q-1) Jdomo = Idomo - kstart = max(0,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) kend = idxI-1 else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then ! SOMO -> SOMO Jsomo = IBCLR(Isomo,p-1) Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBSET(Idomo,q-1) - kstart = max(0,cfg_seniority_index(max(0,Nsomo_I-4))) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) kend = idxI-1 else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then ! DOMO -> VMO @@ -408,7 +408,7 @@ END_PROVIDER Jsomo = IBCLR(Jsomo,q-1) Jdomo = IBCLR(Idomo,p-1) Jdomo = IBSET(Jdomo,q-1) - kstart = max(0,cfg_seniority_index(max(0,Nsomo_I-2))) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) kend = idxI-1 else print*,"Something went wrong in obtain_associated_alphaI" diff --git a/src/csf/configurations.irp.f b/src/csf/configurations.irp.f index 336691c2..8a398a66 100644 --- a/src/csf/configurations.irp.f +++ b/src/csf/configurations.irp.f @@ -460,16 +460,19 @@ END_PROVIDER BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num) ] &BEGIN_PROVIDER [ integer, cfg_nsomo_max ] +&BEGIN_PROVIDER [ integer, cfg_nsomo_min ] implicit none BEGIN_DOC ! Returns the index in psi_configuration of the first cfg with ! the requested seniority ! ! cfg_nsomo_max : Max number of SOMO in the current wave function + ! cfg_nsomo_min : Min number of SOMO in the current wave function END_DOC - integer :: i, k, s, sold + integer :: i, k, s, sold, soldmin cfg_seniority_index(:) = -1 sold = -1 + soldmin = 2000 cfg_nsomo_max = 0 do i=1,N_configuration s = 0 @@ -482,6 +485,10 @@ END_PROVIDER cfg_seniority_index(s) = i cfg_nsomo_max = s endif + if (soldmin .GT. s ) then + soldmin = s + cfg_nsomo_min = s + endif enddo END_PROVIDER diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 377740b2..8dc94ce2 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1,4 +1,5 @@ BEGIN_PROVIDER [ integer, NSOMOMax] + &BEGIN_PROVIDER [ integer, NSOMOMin] &BEGIN_PROVIDER [ integer, NCSFMax] &BEGIN_PROVIDER [ integer*8, NMO] &BEGIN_PROVIDER [ integer, NBFMax] @@ -11,6 +12,7 @@ ! required for the calculation of prototype arrays. END_DOC NSOMOMax = min(elec_num, cfg_nsomo_max + 2) + NSOMOMin = cfg_nsomo_min ! 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 NBFMax = NCSFMax @@ -27,11 +29,13 @@ integer ncfgpersomo detDimperBF = 0 MS = elec_alpha_num-elec_beta_num - !print *,"NSOMOMax=",NSOMOMax, cfg_seniority_index(0) + !print *,"NSOMOMax=",NSOMOMax, cfg_seniority_index(NSOMOMin) + !print *,"NSOMOMin=",NSOMOMin ! number of cfgs = number of dets for 0 somos - n_CSF = cfg_seniority_index(0)-1 - ncfgprev = cfg_seniority_index(0) - do i = 0-iand(MS,1)+2, NSOMOMax,2 + n_CSF = cfg_seniority_index(NSOMOMin)-1 + ncfgprev = cfg_seniority_index(NSOMOMin) + !do i = 0-iand(MS,1)+2, NSOMOMax,2 + do i = NSOMOMin, NSOMOMax,2 if(cfg_seniority_index(i) .EQ. -1)then ncfgpersomo = N_configuration + 1 else @@ -163,6 +167,7 @@ end subroutine get_phase_qp_to_cfg countdet = 0 integer istate istate = 1 + psi_csf_to_config_data(1) = 1 phasedet = 1.0d0 do i = 1,N_configuration startdet = psi_configuration_to_psi_det(1,i) @@ -208,15 +213,17 @@ end subroutine get_phase_qp_to_cfg deallocate(tempCoeff) deallocate(tempBuffer) psi_config_data(i,1) = countcsf + 1 + do k=1,bfIcfg + psi_csf_to_config_data(countcsf+k) = i + enddo 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 - BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)] + BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (NSOMOMin:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)] &BEGIN_PROVIDER [ integer, rowsmax] &BEGIN_PROVIDER [ integer, colsmax] use cfunctions @@ -235,7 +242,6 @@ end subroutine get_phase_qp_to_cfg cols = -1 integer*8 MS MS = elec_alpha_num-elec_beta_num - integer nsomomin nsomomin = elec_alpha_num-elec_beta_num rowsmax = 0 colsmax = 0 @@ -246,7 +252,7 @@ end subroutine get_phase_qp_to_cfg !print *,"Doing SOMO->SOMO" AIJpqMatrixDimsList(0,1,1,1,1) = 1 AIJpqMatrixDimsList(0,1,1,1,2) = 1 - do i = 2-iand(nsomomin,1), NSOMOMax, 2 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i-2,i-2, 2 Jsomo = ISHFT(1_8,j)-1 @@ -292,7 +298,7 @@ end subroutine get_phase_qp_to_cfg !print *,"Doing DOMO->VMO" AIJpqMatrixDimsList(0,2,1,1,1) = 1 AIJpqMatrixDimsList(0,2,1,1,2) = 1 - do i = 0+iand(nsomomin,1), NSOMOMax, 2 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 do j = i+2,i+2, 2 @@ -344,7 +350,7 @@ end subroutine get_phase_qp_to_cfg !print *,"Doing SOMO->VMO" AIJpqMatrixDimsList(0,3,1,1,1) = 1 AIJpqMatrixDimsList(0,3,1,1,2) = 1 - do i = 2-iand(nsomomin,1), NSOMOMax, 2 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 Jsomo = ISHFT(1_8,j)-1 @@ -386,7 +392,7 @@ end subroutine get_phase_qp_to_cfg !print *,"Doing DOMO->SOMO" AIJpqMatrixDimsList(0,4,1,1,1) = 1 AIJpqMatrixDimsList(0,4,1,1,2) = 1 - do i = 2-iand(nsomomin,1), NSOMOMax, 2 + do i = NSOMOMin, NSOMOMax, 2 do j = i,i, 2 if(j .GT. NSOMOMax .OR. j .LE. 0) then cycle @@ -425,7 +431,6 @@ end subroutine get_phase_qp_to_cfg END_PROVIDER BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,NBFMax,NBFMax)] - !BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,-1:1,4,NSOMOMax+1,-1:1,NBFMax,NBFMax)] use cfunctions implicit none BEGIN_DOC @@ -467,7 +472,7 @@ end subroutine get_phase_qp_to_cfg ! 1. SOMO -> SOMO !print *,"Doing SOMO -> SOMO" AIJpqContainer(0,1,1,1,1,1) = 1.0d0 - do i = 2, NSOMOMax, 2 + do i = NSOMOMin+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 @@ -530,7 +535,7 @@ end subroutine get_phase_qp_to_cfg ! 2. DOMO -> VMO !print *,"Doing DOMO -> VMO" AIJpqContainer(0,2,1,1,1,1) = 1.0d0 - do i = 0, NSOMOMax, 2 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 do j = i+2,i+2, 2 @@ -596,7 +601,7 @@ end subroutine get_phase_qp_to_cfg ! Type ! 3. SOMO -> VMO !print *,"Doing SOMO -> VMO" - do i = 2, NSOMOMax, 2 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 Jsomo = ISHFT(1_8,j)-1 @@ -655,7 +660,7 @@ end subroutine get_phase_qp_to_cfg ! 4. DOMO -> SOMO !print *,"Doing DOMO -> SOMO" AIJpqContainer(0,4,1,1,1,1) = 1.0d0 - do i = 2, NSOMOMax, 2 + do i = NSOMOMin+2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 Jsomo = ISHFT(1_8,i)-1 @@ -902,6 +907,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze real*8,dimension(:,:,:),allocatable :: TKI real*8,dimension(:,:),allocatable :: GIJpqrs real*8,dimension(:,:,:),allocatable :: TKIGIJ + real*8,dimension(:),allocatable :: psi_out_tmp + real*8,dimension(:,:),allocatable :: CCmattmp real*8, external :: mo_two_e_integral real*8, external :: get_two_e_integral real*8 :: diag_energies(n_CSF) @@ -1187,11 +1194,29 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) !print *,">j=",j,rowsikpq,colsikpq, ">>",totcolsTKI,",",idxs_connectedI_alpha(j) do kk = 1,n_st + !do m = 1,colsikpq + ! do l = 1,rowsTKI + ! psi_out(idxs_connectedI_alpha(j)+m-1,kk) += AIJpqContainer(NSOMOalpha,extype,pmodel,qmodel,l,m) * TKIGIJ(l,kk,j) + ! enddo + !enddo + allocate(psi_out_tmp(colsikpq)) + !allocate(CCmattmp(rowsTKI,colsikpq)) + !do m=1,colsikpq + ! do l=1,rowsTKI + ! CCmattmp(l,m) = AIJpqContainer(NSOMOalpha,extype,pmodel,qmodel,l,m) + ! enddo + !enddo + !call dgemv('T',rowsTKI, colsikpq, 1.d0, & + ! CCmattmp, size(CCmattmp,1), TKIGIJ(:,kk,j), 1, 0.d0, & + ! psi_out_tmp, 1) + call dgemv('T',rowsTKI, colsikpq, 1.d0, & + AIJpqContainer(NSOMOalpha,extype,pmodel,qmodel,1:rowsTKI,1:colsikpq), rowsTKI, TKIGIJ(:,kk,j), 1, 0.d0, & + psi_out_tmp, 1) do m = 1,colsikpq - do l = 1,rowsTKI - psi_out(idxs_connectedI_alpha(j)+m-1,kk) += AIJpqContainer(NSOMOalpha,extype,pmodel,qmodel,l,m) * TKIGIJ(l,kk,j) - enddo + psi_out(idxs_connectedI_alpha(j)+m-1,kk) += psi_out_tmp(m) enddo + deallocate(psi_out_tmp) + !deallocate(CCmattmp) enddo totcolsTKI += colsikpq enddo @@ -1274,16 +1299,16 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie real*8 :: diag_energies(n_CSF) ! allocate - allocate(alphas_Icfg(N_INT,2,max(sze,100))) - allocate(singlesI(N_INT,2,max(sze,100))) - allocate(connectedI_alpha(N_INT,2,max(sze,100))) - allocate(idxs_singlesI(max(sze,100))) - allocate(idxs_connectedI_alpha(max(sze,100))) - allocate(excitationIds_single(2,max(sze,100))) - allocate(excitationTypes_single(max(sze,100))) - allocate(excitationIds(2,max(sze,100))) - allocate(excitationTypes(max(sze,100))) - allocate(diagfactors(max(sze,100))) + allocate(alphas_Icfg(N_INT,2,max(sze/2,100))) + allocate(singlesI(N_INT,2,max(sze/2,100))) + allocate(connectedI_alpha(N_INT,2,max(sze/2,100))) + allocate(idxs_singlesI(max(sze/2,100))) + allocate(idxs_connectedI_alpha(max(sze/2,100))) + allocate(excitationIds_single(2,max(sze/2,100))) + allocate(excitationTypes_single(max(sze/2,100))) + allocate(excitationIds(2,max(sze/2,100))) + allocate(excitationTypes(max(sze/2,100))) + allocate(diagfactors(max(sze/2,100))) !print *," sze = ",sze @@ -1300,6 +1325,7 @@ subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, ie !!! Single Excitations !!! do i=istart_cfg,iend_cfg + print *,"I=",i ! if Seniority_range > 8 then ! continue