mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-25 13:03:28 +01:00
Fixed OpenMP
This commit is contained in:
parent
622fc1531c
commit
acef0b913b
@ -287,7 +287,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
|
|||||||
call write_int(6,nproc_target,'Number of threads for PT2')
|
call write_int(6,nproc_target,'Number of threads for PT2')
|
||||||
call write_double(6,mem,'Memory (Gb)')
|
call write_double(6,mem,'Memory (Gb)')
|
||||||
|
|
||||||
call omp_set_nested(.false.)
|
call omp_set_max_active_levels(1)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
print '(A)', '========== ======================= ===================== ===================== ==========='
|
print '(A)', '========== ======================= ===================== ===================== ==========='
|
||||||
|
|||||||
@ -56,6 +56,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
|
|||||||
|
|
||||||
p = 0
|
p = 0
|
||||||
q = 0
|
q = 0
|
||||||
|
if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1'
|
||||||
do i=idxI,end_index
|
do i=idxI,end_index
|
||||||
Isomo = Ialpha(1,1)
|
Isomo = Ialpha(1,1)
|
||||||
Idomo = Ialpha(1,2)
|
Idomo = Ialpha(1,2)
|
||||||
@ -71,16 +72,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
|
|||||||
ndiffDOMO = POPCNT(diffDOMO)
|
ndiffDOMO = POPCNT(diffDOMO)
|
||||||
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
|
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
|
||||||
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
|
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
|
||||||
!print *,"-I--i=",i,ndiffSOMO,ndiffDOMO,nxordiffSOMODOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO
|
|
||||||
!if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle
|
|
||||||
!print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO
|
|
||||||
!if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
|
|
||||||
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
|
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
|
||||||
!call debug_spindet(Isomo,1)
|
|
||||||
!call debug_spindet(Idomo,1)
|
|
||||||
!print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration
|
|
||||||
!call debug_spindet(Jsomo,1)
|
|
||||||
!call debug_spindet(Jdomo,1)
|
|
||||||
select case(ndiffDOMO)
|
select case(ndiffDOMO)
|
||||||
case (0)
|
case (0)
|
||||||
! SOMO -> VMO
|
! SOMO -> VMO
|
||||||
|
|||||||
@ -232,8 +232,7 @@ end subroutine get_phase_qp_to_cfg
|
|||||||
print *,"Norm det=",norm_det1, size(psi_coef_config,1), " Dim csf=", countcsf
|
print *,"Norm det=",norm_det1, size(psi_coef_config,1), " Dim csf=", countcsf
|
||||||
!$OMP END MASTER
|
!$OMP END MASTER
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
call omp_set_max_active_levels(4)
|
||||||
call omp_set_nested(.True.)
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -987,12 +986,6 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
allocate(excitationIds_single(2,max(sze,100)))
|
allocate(excitationIds_single(2,max(sze,100)))
|
||||||
allocate(excitationTypes_single(max(sze,100)))
|
allocate(excitationTypes_single(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)))
|
|
||||||
|
|
||||||
!!! Single Excitations !!!
|
!!! Single Excitations !!!
|
||||||
|
|
||||||
@ -1059,7 +1052,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
endi = psi_config_data(i,2)
|
endi = psi_config_data(i,2)
|
||||||
NSOMOI = getNSOMO(Icfg)
|
NSOMOI = getNSOMO(Icfg)
|
||||||
|
|
||||||
call generate_all_singles_cfg_with_type(Icfg,singlesI,idxs_singlesI,excitationIds_single, &
|
call generate_all_singles_cfg_with_type(Icfg,singlesI,idxs_singlesI,excitationIds_single,&
|
||||||
excitationTypes_single,nsinglesI,N_int)
|
excitationTypes_single,nsinglesI,N_int)
|
||||||
|
|
||||||
do j = 1,nsinglesI
|
do j = 1,nsinglesI
|
||||||
@ -1089,17 +1082,17 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
endj = psi_config_data(idxI,2)
|
endj = psi_config_data(idxI,2)
|
||||||
|
|
||||||
!!! One-electron contribution !!!
|
!!! One-electron contribution !!!
|
||||||
do kk = 1,n_st
|
|
||||||
cnti = 0
|
|
||||||
do ii = starti, endi
|
do ii = starti, endi
|
||||||
cnti += 1
|
cnti = ii-starti+1
|
||||||
cntj = 0
|
|
||||||
do jj = startj, endj
|
do jj = startj, endj
|
||||||
cntj += 1
|
cntj = jj-startj+1
|
||||||
!meCC1 = AIJpqContainer(NSOMOI,extype,pmodel,qmodel,cnti,cntj)
|
!meCC1 = AIJpqContainer(NSOMOI,extype,pmodel,qmodel,cnti,cntj)
|
||||||
meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)
|
meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)
|
||||||
|
call omp_set_lock(lock(jj))
|
||||||
|
do kk = 1,n_st
|
||||||
psi_out(kk,jj) += meCC1 * psi_in(kk,ii) * h_core_ri(p,q)
|
psi_out(kk,jj) += meCC1 * psi_in(kk,ii) * h_core_ri(p,q)
|
||||||
enddo
|
enddo
|
||||||
|
call omp_unset_lock(lock(jj))
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
@ -1111,12 +1104,23 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
nholes -= 1
|
nholes -= 1
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
! enddo
|
enddo
|
||||||
! !$OMP END DO
|
!$OMP END DO
|
||||||
|
deallocate(singlesI)
|
||||||
|
deallocate(idxs_singlesI)
|
||||||
|
deallocate(excitationIds_single)
|
||||||
|
deallocate(excitationTypes_single)
|
||||||
|
|
||||||
! Loop over all selected configurations
|
allocate(alphas_Icfg(N_INT,2,max(sze,100)))
|
||||||
! !$OMP DO SCHEDULE(dynamic,128)
|
allocate(connectedI_alpha(N_INT,2,max(sze,100)))
|
||||||
! do i = istart_cfg,iend_cfg
|
allocate(idxs_connectedI_alpha(max(sze,100)))
|
||||||
|
allocate(excitationIds(2,max(sze,100)))
|
||||||
|
allocate(excitationTypes(max(sze,100)))
|
||||||
|
allocate(diagfactors(max(sze,100)))
|
||||||
|
|
||||||
|
! Loop over all selected configurations
|
||||||
|
!$OMP DO SCHEDULE(dynamic,16)
|
||||||
|
do i = istart_cfg,iend_cfg
|
||||||
|
|
||||||
! if Seniority_range > 8 then
|
! if Seniority_range > 8 then
|
||||||
! continue
|
! continue
|
||||||
@ -1134,7 +1138,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
! test if size(alphas_Icfg,1) < Nmo**2) then deallocate + allocate
|
! test if size(alphas_Icfg,1) < Nmo**2) then deallocate + allocate
|
||||||
!call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
|
!call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
|
||||||
Nalphas_Icfg = NalphaIcfg_list(i)
|
Nalphas_Icfg = NalphaIcfg_list(i)
|
||||||
alphas_Icfg(1:n_int,1:2,1:Nalphas_Icfg) = alphasIcfg_list(1:n_int,1:2,i,1:Nalphas_Icfg)
|
alphas_Icfg(1:N_int,1:2,1:Nalphas_Icfg) = alphasIcfg_list(1:n_int,1:2,i,1:Nalphas_Icfg)
|
||||||
|
|
||||||
! TODO : remove doubly excited for return
|
! TODO : remove doubly excited for return
|
||||||
! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
|
! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
|
||||||
@ -1142,14 +1146,11 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
! Now generate all singly excited with respect to a given alpha CFG
|
! Now generate all singly excited with respect to a given alpha CFG
|
||||||
call obtain_connected_I_foralpha(i,alphas_Icfg(1,1,k),connectedI_alpha,idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors)
|
call obtain_connected_I_foralpha(i,alphas_Icfg(1,1,k),connectedI_alpha,idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors)
|
||||||
|
|
||||||
if(nconnectedI .EQ. 0) then
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
totcolsTKI = 0
|
totcolsTKI = 0
|
||||||
rowsTKI = -1
|
rowsTKI = -1
|
||||||
do j = 1,nconnectedI
|
do j = 1,nconnectedI
|
||||||
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
|
NSOMOalpha = getNSOMO(alphas_Icfg(1,1,k))
|
||||||
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
|
NSOMOI = getNSOMO(connectedI_alpha(1,1,j))
|
||||||
p = excitationIds(1,j)
|
p = excitationIds(1,j)
|
||||||
q = excitationIds(2,j)
|
q = excitationIds(2,j)
|
||||||
extype = excitationTypes(j)
|
extype = excitationTypes(j)
|
||||||
@ -1161,12 +1162,12 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
||||||
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
||||||
totcolsTKI += colsikpq
|
totcolsTKI += colsikpq
|
||||||
if(rowsTKI .LT. rowsikpq .AND. rowsTKI .NE. -1) then
|
! if(rowsTKI .LT. rowsikpq .AND. rowsTKI .NE. -1) then
|
||||||
print *,">",j,"Something is wrong in sigma-vector", rowsTKI, rowsikpq, "(p,q)=",pmodel,qmodel,"ex=",extype,"na=",NSOMOalpha," nI=",NSOMOI
|
! print *,">",j,"Something is wrong in sigma-vector", rowsTKI, rowsikpq, "(p,q)=",pmodel,qmodel,"ex=",extype,"na=",NSOMOalpha," nI=",NSOMOI
|
||||||
!rowsTKI = rowsikpq
|
! !rowsTKI = rowsikpq
|
||||||
else
|
! else
|
||||||
rowsTKI = rowsikpq
|
rowsTKI = rowsikpq
|
||||||
endif
|
! endif
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
allocate(TKI(n_st,rowsTKI,totcolsTKI)) ! coefficients of CSF
|
allocate(TKI(n_st,rowsTKI,totcolsTKI)) ! coefficients of CSF
|
||||||
@ -1177,31 +1178,21 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
|
|
||||||
totcolsTKI = 0
|
totcolsTKI = 0
|
||||||
do j = 1,nconnectedI
|
do j = 1,nconnectedI
|
||||||
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
|
NSOMOalpha = getNSOMO(alphas_Icfg(1,1,k))
|
||||||
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
|
NSOMOI = getNSOMO(connectedI_alpha(1,1,j))
|
||||||
p = excitationIds(1,j)
|
p = excitationIds(1,j)
|
||||||
q = excitationIds(2,j)
|
q = excitationIds(2,j)
|
||||||
extype = excitationTypes(j)
|
extype = excitationTypes(j)
|
||||||
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
|
call convertOrbIdsToModelSpaceIds(alphas_Icfg(1,1,k), connectedI_alpha(1,1,j), p, q, extype, pmodel, qmodel)
|
||||||
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
||||||
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
||||||
allocate(CCmattmp(colsikpq,n_st))
|
|
||||||
!do kk = 1,n_st
|
|
||||||
!do m = 1,colsikpq
|
|
||||||
! CCmattmp(m,kk) = psi_in(idxs_connectedI_alpha(j)+m-1,kk)
|
|
||||||
!enddo
|
|
||||||
!enddo
|
|
||||||
do m = 1,colsikpq
|
do m = 1,colsikpq
|
||||||
do l = 1,rowsTKI
|
do l = 1,rowsTKI
|
||||||
do kk = 1,n_st
|
do kk = 1,n_st
|
||||||
!tmpvar = CCmattmp(m,kk)
|
|
||||||
!TKI(kk,l,totcolsTKI+m) = AIJpqContainer(NSOMOalpha,extype,pmodel,qmodel,l,m) * tmpvar
|
|
||||||
!TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) * tmpvar
|
|
||||||
TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) * psi_in(kk,idxs_connectedI_alpha(j)+m-1)
|
TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) * psi_in(kk,idxs_connectedI_alpha(j)+m-1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
deallocate(CCmattmp)
|
|
||||||
do m = 1,colsikpq
|
do m = 1,colsikpq
|
||||||
do l = 1,nconnectedI
|
do l = 1,nconnectedI
|
||||||
! <ij|kl> = (ik|jl)
|
! <ij|kl> = (ik|jl)
|
||||||
@ -1229,47 +1220,28 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
|
|
||||||
! Do big BLAS
|
! Do big BLAS
|
||||||
! TODO TKI, size(TKI,1)*size(TKI,2)
|
! TODO TKI, size(TKI,1)*size(TKI,2)
|
||||||
call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0, &
|
call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0,&
|
||||||
TKI, size(TKI,1)*size(TKI,2), GIJpqrs, size(GIJpqrs,1), 0.d0, &
|
TKI, size(TKI,1)*size(TKI,2), GIJpqrs, size(GIJpqrs,1), 0.d0,&
|
||||||
TKIGIJ , size(TKIGIJ,1)*size(TKIGIJ,2) )
|
TKIGIJ , size(TKIGIJ,1)*size(TKIGIJ,2) )
|
||||||
|
|
||||||
!print *,"DIMs = ",rowsTKI,n_st,totcolsTKI,nconnectedI
|
|
||||||
!print *,"TKI mat"
|
|
||||||
!do kk=1,n_st
|
|
||||||
! do j=1,totcolsTKI
|
|
||||||
! print *,TKI(:,kk,j)
|
|
||||||
! enddo
|
|
||||||
! print *,"--"
|
|
||||||
!enddo
|
|
||||||
|
|
||||||
!print *,"TKIGIJ mat"
|
|
||||||
!do kk=1,n_st
|
|
||||||
! do j=1,nconnectedI
|
|
||||||
! print *,TKIGIJ(:,kk,j)
|
|
||||||
! enddo
|
|
||||||
! print *,"--"
|
|
||||||
!enddo
|
|
||||||
|
|
||||||
|
|
||||||
! Collect the result
|
! Collect the result
|
||||||
totcolsTKI = 0
|
totcolsTKI = 0
|
||||||
do j = 1,nconnectedI
|
do j = 1,nconnectedI
|
||||||
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
|
NSOMOalpha = getNSOMO(alphas_Icfg(1,1,k))
|
||||||
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
|
NSOMOI = getNSOMO(connectedI_alpha(1,1,j))
|
||||||
p = excitationIds(1,j)
|
p = excitationIds(1,j)
|
||||||
q = excitationIds(2,j)
|
q = excitationIds(2,j)
|
||||||
extype = excitationTypes(j)
|
extype = excitationTypes(j)
|
||||||
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
|
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
|
||||||
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
||||||
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
||||||
!print *,">j=",j,rowsikpq,colsikpq, ">>",totcolsTKI,",",idxs_connectedI_alpha(j)
|
|
||||||
do m = 1,colsikpq
|
do m = 1,colsikpq
|
||||||
!tmpvar = psi_out(kk,idxs_connectedI_alpha(j)+m-1)
|
|
||||||
do l = 1,rowsTKI
|
do l = 1,rowsTKI
|
||||||
call omp_set_lock(lock(idxs_connectedI_alpha(j)+m-1))
|
call omp_set_lock(lock(idxs_connectedI_alpha(j)+m-1))
|
||||||
do kk = 1,n_st
|
do kk = 1,n_st
|
||||||
!tmpvar += AIJpqContainer(NSOMOalpha,extype,pmodel,qmodel,l,m) * TKIGIJ(kk,l,j)
|
psi_out(kk,idxs_connectedI_alpha(j)+m-1) = psi_out(kk,idxs_connectedI_alpha(j)+m-1) + &
|
||||||
psi_out(kk,idxs_connectedI_alpha(j)+m-1) += AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) * TKIGIJ(kk,l,j)
|
AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) * TKIGIJ(kk,l,j)
|
||||||
enddo
|
enddo
|
||||||
call omp_unset_lock(lock(idxs_connectedI_alpha(j)+m-1))
|
call omp_unset_lock(lock(idxs_connectedI_alpha(j)+m-1))
|
||||||
enddo
|
enddo
|
||||||
@ -1304,354 +1276,9 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
|
|||||||
|
|
||||||
!$OMP end parallel
|
!$OMP end parallel
|
||||||
call omp_set_max_active_levels(4)
|
call omp_set_max_active_levels(4)
|
||||||
|
! do i=1,sze
|
||||||
|
! call omp_deinit_lock(lock(i))
|
||||||
|
! enddo
|
||||||
|
|
||||||
end subroutine calculate_sigma_vector_cfg_nst_naive_store
|
end subroutine calculate_sigma_vector_cfg_nst_naive_store
|
||||||
|
|
||||||
subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep)
|
|
||||||
implicit none
|
|
||||||
use bitmasks
|
|
||||||
BEGIN_DOC
|
|
||||||
! Documentation for sigma-vector calculation
|
|
||||||
!
|
|
||||||
! Calculates the result of the
|
|
||||||
! application of the hamiltonian to the
|
|
||||||
! wavefunction in CFG basis once
|
|
||||||
! TODO : Things prepare outside this routine
|
|
||||||
! 1. Touch the providers for
|
|
||||||
! a. ApqIJ containers
|
|
||||||
! b. DET to CSF transformation matrices
|
|
||||||
! 2. DET to CSF transcormation
|
|
||||||
! 2. CSF to DET back transcormation
|
|
||||||
! returns : psi_coef_out_det :
|
|
||||||
END_DOC
|
|
||||||
integer,intent(in) :: sze, istart,iend, istep, ishift, n_st
|
|
||||||
real*8,intent(in) :: psi_in(sze,n_st)
|
|
||||||
real*8,intent(out) :: psi_out(sze,n_st)
|
|
||||||
integer(bit_kind) :: Icfg(N_INT,2)
|
|
||||||
integer :: i,j,k,l,p,q,noccp,noccq, ii, jj, m, n, idxI, kk, nocck,orbk
|
|
||||||
integer(bit_kind),dimension(:,:,:),allocatable :: alphas_Icfg
|
|
||||||
integer(bit_kind),dimension(:,:,:),allocatable :: singlesI
|
|
||||||
integer(bit_kind),dimension(:,:,:),allocatable :: connectedI_alpha
|
|
||||||
integer,dimension(:),allocatable :: idxs_singlesI
|
|
||||||
integer,dimension(:),allocatable :: idxs_connectedI_alpha
|
|
||||||
integer,dimension(:,:),allocatable :: excitationIds_single
|
|
||||||
integer,dimension(:),allocatable :: excitationTypes_single
|
|
||||||
integer,dimension(:,:),allocatable :: excitationIds
|
|
||||||
integer,dimension(:),allocatable :: excitationTypes
|
|
||||||
real*8,dimension(:),allocatable :: diagfactors
|
|
||||||
integer :: nholes
|
|
||||||
integer :: nvmos
|
|
||||||
integer :: listvmos(mo_num)
|
|
||||||
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
|
|
||||||
integer :: listholes(mo_num)
|
|
||||||
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
|
|
||||||
integer :: Nalphas_Icfg, nconnectedI, rowsikpq, colsikpq, nsinglesI
|
|
||||||
integer :: extype,NSOMOalpha,NSOMOI,NSOMOJ,pmodel,qmodel
|
|
||||||
integer :: getNSOMO
|
|
||||||
integer :: totcolsTKI
|
|
||||||
integer :: rowsTKI
|
|
||||||
integer :: noccpp
|
|
||||||
integer :: istart_cfg, iend_cfg
|
|
||||||
integer*8 :: MS, Isomo, Idomo, Jsomo, Jdomo, Ialpha, Ibeta
|
|
||||||
integer :: moi, moj, mok, mol, starti, endi, startj, endj, cnti, cntj, cntk
|
|
||||||
real*8 :: norm_coef_cfg, fac2eints
|
|
||||||
real*8 :: norm_coef_det
|
|
||||||
real*8 :: meCC1, meCC2, diagfac
|
|
||||||
real*8,dimension(:,:,:),allocatable :: TKI
|
|
||||||
real*8,dimension(:,:),allocatable :: GIJpqrs
|
|
||||||
real*8,dimension(:,:,:),allocatable :: TKIGIJ
|
|
||||||
real*8, external :: mo_two_e_integral
|
|
||||||
real*8, external :: get_two_e_integral
|
|
||||||
real*8 :: diag_energies(n_CSF)
|
|
||||||
|
|
||||||
! allocate
|
|
||||||
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
|
|
||||||
call calculate_preconditioner_cfg(diag_energies)
|
|
||||||
|
|
||||||
MS = 0
|
|
||||||
norm_coef_cfg=0.d0
|
|
||||||
|
|
||||||
psi_out=0.d0
|
|
||||||
|
|
||||||
istart_cfg = psi_csf_to_config_data(istart)
|
|
||||||
iend_cfg = psi_csf_to_config_data(iend)
|
|
||||||
|
|
||||||
|
|
||||||
!!! Single Excitations !!!
|
|
||||||
do i=istart_cfg,iend_cfg
|
|
||||||
print *,"I=",i
|
|
||||||
|
|
||||||
! if Seniority_range > 8 then
|
|
||||||
! continue
|
|
||||||
! else
|
|
||||||
! cycle
|
|
||||||
|
|
||||||
Icfg(1,1) = psi_configuration(1,1,i)
|
|
||||||
Icfg(1,2) = psi_configuration(1,2,i)
|
|
||||||
Isomo = Icfg(1,1)
|
|
||||||
Idomo = Icfg(1,2)
|
|
||||||
NSOMOI = getNSOMO(Icfg)
|
|
||||||
|
|
||||||
! find out all pq holes possible
|
|
||||||
nholes = 0
|
|
||||||
! holes in SOMO
|
|
||||||
! list_act
|
|
||||||
! list_core
|
|
||||||
! list_core_inact
|
|
||||||
! bitmasks
|
|
||||||
!do k = n_core_orb+1,n_core_orb + n_act_orb
|
|
||||||
do k = 1,mo_num
|
|
||||||
if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
|
|
||||||
nholes += 1
|
|
||||||
listholes(nholes) = k
|
|
||||||
holetype(nholes) = 1
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
! holes in DOMO
|
|
||||||
!do k = n_core_orb+1,n_core_orb + n_act_orb
|
|
||||||
do k = 1,mo_num
|
|
||||||
if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
|
|
||||||
nholes += 1
|
|
||||||
listholes(nholes) = k
|
|
||||||
holetype(nholes) = 2
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! find vmos
|
|
||||||
listvmos = -1
|
|
||||||
vmotype = -1
|
|
||||||
nvmos = 0
|
|
||||||
!do k = n_core_orb+1,n_core_orb + n_act_orb
|
|
||||||
do k = 1,mo_num
|
|
||||||
!print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
|
|
||||||
if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0) then
|
|
||||||
nvmos += 1
|
|
||||||
listvmos(nvmos) = k
|
|
||||||
vmotype(nvmos) = 0
|
|
||||||
else if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0 ) then
|
|
||||||
nvmos += 1
|
|
||||||
listvmos(nvmos) = k
|
|
||||||
vmotype(nvmos) = 1
|
|
||||||
end if
|
|
||||||
enddo
|
|
||||||
|
|
||||||
|
|
||||||
! Icsf ids
|
|
||||||
starti = psi_config_data(i,1)
|
|
||||||
endi = psi_config_data(i,2)
|
|
||||||
NSOMOI = getNSOMO(Icfg)
|
|
||||||
|
|
||||||
call generate_all_singles_cfg_with_type(Icfg,singlesI,idxs_singlesI,excitationIds_single, &
|
|
||||||
excitationTypes_single,nsinglesI,N_int)
|
|
||||||
|
|
||||||
do j = 1,nsinglesI
|
|
||||||
idxI = idxs_singlesI(j)
|
|
||||||
NSOMOJ = getNSOMO(singlesI(:,:,j))
|
|
||||||
p = excitationIds_single(1,j)
|
|
||||||
q = excitationIds_single(2,j)
|
|
||||||
extype = excitationTypes_single(j)
|
|
||||||
! Off diagonal terms
|
|
||||||
call convertOrbIdsToModelSpaceIds(Icfg, singlesI(:,:,j), p, q, extype, pmodel, qmodel)
|
|
||||||
Jsomo = singlesI(1,1,j)
|
|
||||||
Jdomo = singlesI(1,2,j)
|
|
||||||
|
|
||||||
! Add the hole on J
|
|
||||||
if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
|
|
||||||
nholes += 1
|
|
||||||
listholes(nholes) = q
|
|
||||||
holetype(nholes) = 1
|
|
||||||
endif
|
|
||||||
if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
|
|
||||||
nholes += 1
|
|
||||||
listholes(nholes) = q
|
|
||||||
holetype(nholes) = 2
|
|
||||||
endif
|
|
||||||
|
|
||||||
startj = psi_config_data(idxI,1)
|
|
||||||
endj = psi_config_data(idxI,2)
|
|
||||||
|
|
||||||
!!! One-electron contribution !!!
|
|
||||||
do kk = 1,n_st
|
|
||||||
cnti = 0
|
|
||||||
do ii = starti, endi
|
|
||||||
cnti += 1
|
|
||||||
cntj = 0
|
|
||||||
do jj = startj, endj
|
|
||||||
cntj += 1
|
|
||||||
meCC1 = AIJpqContainer(NSOMOI,extype,pmodel,qmodel,cnti,cntj)
|
|
||||||
psi_out(jj,kk) += meCC1 * psi_in(ii,kk) * h_core_ri(p,q)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! Undo setting in listholes
|
|
||||||
if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
|
|
||||||
nholes -= 1
|
|
||||||
endif
|
|
||||||
if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
|
|
||||||
nholes -= 1
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!!! Double Excitations !!!
|
|
||||||
|
|
||||||
! Loop over all selected configurations
|
|
||||||
do i = istart_cfg,iend_cfg
|
|
||||||
|
|
||||||
! if Seniority_range > 8 then
|
|
||||||
! continue
|
|
||||||
! else
|
|
||||||
! cycle
|
|
||||||
|
|
||||||
Icfg(1,1) = psi_configuration(1,1,i)
|
|
||||||
Icfg(1,2) = psi_configuration(1,2,i)
|
|
||||||
starti = psi_config_data(i,1)
|
|
||||||
endi = psi_config_data(i,2)
|
|
||||||
|
|
||||||
! Returns all unique (checking the past) singly excited cfgs connected to I
|
|
||||||
Nalphas_Icfg = 0
|
|
||||||
! TODO:
|
|
||||||
! test if size(alphas_Icfg,1) < Nmo**2) then deallocate + allocate
|
|
||||||
call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
|
|
||||||
! TODO : remove doubly excited for return
|
|
||||||
! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
|
|
||||||
do k = 1,Nalphas_Icfg
|
|
||||||
! Now generate all singly excited with respect to a given alpha CFG
|
|
||||||
call obtain_connected_I_foralpha(i,alphas_Icfg(1,1,k),connectedI_alpha,idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors)
|
|
||||||
|
|
||||||
if(nconnectedI .EQ. 0) then
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
totcolsTKI = 0
|
|
||||||
rowsTKI = -1
|
|
||||||
do j = 1,nconnectedI
|
|
||||||
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
|
|
||||||
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
|
|
||||||
p = excitationIds(1,j)
|
|
||||||
q = excitationIds(2,j)
|
|
||||||
extype = excitationTypes(j)
|
|
||||||
call convertOrbIdsToModelSpaceIds(alphas_Icfg(1,1,k), connectedI_alpha(1,1,j), p, q, extype, pmodel, qmodel)
|
|
||||||
! for E_pp E_rs and E_ppE_rr case
|
|
||||||
if(p.EQ.q) then
|
|
||||||
NSOMOalpha = NSOMOI
|
|
||||||
endif
|
|
||||||
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
|
||||||
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
|
||||||
totcolsTKI += colsikpq
|
|
||||||
if(rowsTKI .LT. rowsikpq .AND. rowsTKI .NE. -1) then
|
|
||||||
print *,">",j,"Something is wrong in sigma-vector", rowsTKI, rowsikpq, "(p,q)=",pmodel,qmodel,"ex=",extype,"na=",NSOMOalpha," nI=",NSOMOI
|
|
||||||
!rowsTKI = rowsikpq
|
|
||||||
else
|
|
||||||
rowsTKI = rowsikpq
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
allocate(TKI(rowsTKI,n_st,totcolsTKI)) ! coefficients of CSF
|
|
||||||
! Initialize the inegral container
|
|
||||||
! dims : (totcolsTKI, nconnectedI)
|
|
||||||
allocate(GIJpqrs(totcolsTKI,nconnectedI)) ! gpqrs
|
|
||||||
allocate(TKIGIJ(rowsTKI,n_st,nconnectedI)) ! gpqrs
|
|
||||||
|
|
||||||
totcolsTKI = 0
|
|
||||||
do j = 1,nconnectedI
|
|
||||||
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
|
|
||||||
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
|
|
||||||
p = excitationIds(1,j)
|
|
||||||
q = excitationIds(2,j)
|
|
||||||
extype = excitationTypes(j)
|
|
||||||
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
|
|
||||||
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
|
||||||
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
|
|
||||||
do kk = 1,n_st
|
|
||||||
do l = 1,rowsTKI
|
|
||||||
do m = 1,colsikpq
|
|
||||||
TKI(l,kk,totcolsTKI+m) = AIJpqContainer(NSOMOalpha,extype,pmodel,qmodel,l,m) * psi_in(idxs_connectedI_alpha(j)+m-1,kk)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
do m = 1,colsikpq
|
|
||||||
do l = 1,nconnectedI
|
|
||||||
! <ij|kl> = (ik|jl)
|
|
||||||
moi = excitationIds(1,j) ! p
|
|
||||||
mok = excitationIds(2,j) ! q
|
|
||||||
moj = excitationIds(2,l) ! s
|
|
||||||
mol = excitationIds(1,l) ! r
|
|
||||||
if(moi.EQ.mok .AND. moj.EQ.mol)then
|
|
||||||
diagfac = diagfactors(j)
|
|
||||||
diagfac *= diagfactors(l)
|
|
||||||
!print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
|
|
||||||
GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
|
|
||||||
else
|
|
||||||
diagfac = diagfactors(j)*diagfactors(l)
|
|
||||||
!print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
|
|
||||||
GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
|
|
||||||
!endif
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
totcolsTKI += colsikpq
|
|
||||||
enddo
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
! Do big BLAS
|
|
||||||
! TODO TKI, size(TKI,1)*size(TKI,2)
|
|
||||||
call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0, &
|
|
||||||
TKI, size(TKI,1)*n_st, GIJpqrs, size(GIJpqrs,1), 0.d0, &
|
|
||||||
TKIGIJ , size(TKIGIJ,1)*n_st )
|
|
||||||
|
|
||||||
|
|
||||||
! Collect the result
|
|
||||||
totcolsTKI = 0
|
|
||||||
do j = 1,nconnectedI
|
|
||||||
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
|
|
||||||
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
|
|
||||||
p = excitationIds(1,j)
|
|
||||||
q = excitationIds(2,j)
|
|
||||||
extype = excitationTypes(j)
|
|
||||||
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
|
|
||||||
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
|
|
||||||
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
|
|
||||||
enddo
|
|
||||||
totcolsTKI += colsikpq
|
|
||||||
enddo
|
|
||||||
|
|
||||||
deallocate(TKI) ! coefficients of CSF
|
|
||||||
! Initialize the inegral container
|
|
||||||
! dims : (totcolsTKI, nconnectedI)
|
|
||||||
deallocate(GIJpqrs) ! gpqrs
|
|
||||||
deallocate(TKIGIJ) ! gpqrs
|
|
||||||
|
|
||||||
enddo ! loop over alphas
|
|
||||||
enddo ! loop over I
|
|
||||||
|
|
||||||
|
|
||||||
! Add the diagonal contribution
|
|
||||||
do kk=1,n_st
|
|
||||||
do i = 1,n_CSF
|
|
||||||
psi_out(i,kk) += 1.0d0*diag_energies(i)*psi_in(i,kk)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
|
|
||||||
end subroutine calculate_sigma_vector_cfg_nst
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user