9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-13 14:34:30 +02:00
qp2/src/csf/sigma_vector.irp.f
2022-03-24 16:18:51 +01:00

707 lines
21 KiB
Fortran

real*8 function logabsgamma(x)
implicit none
real*8, intent(in) :: x
logabsgamma = 1.d32 ! Avoid floating point exception
if (x>0.d0) then
logabsgamma = log(abs(gamma(x)))
endif
end function logabsgamma
BEGIN_PROVIDER [ integer, NSOMOMax]
&BEGIN_PROVIDER [ integer, NCSFMax]
&BEGIN_PROVIDER [ integer*8, NMO]
&BEGIN_PROVIDER [ integer, NBFMax]
&BEGIN_PROVIDER [ integer, n_CSF]
&BEGIN_PROVIDER [ integer, maxDetDimPerBF]
implicit none
BEGIN_DOC
! Documentation for NSOMOMax
! The maximum number of SOMOs for the current calculation.
! required for the calculation of prototype arrays.
END_DOC
NSOMOMax = min(elec_num, cfg_nsomo_max + 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
NBFMax = NCSFMax
maxDetDimPerBF = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2))))
NMO = n_act_orb
integer i,j,k,l
integer startdet,enddet
integer ncfg,ncfgprev
integer NSOMO
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 = 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
! l = i+k+2
! Loop over l with a constraint to ensure that l <= size(cfg_seniority_index,1)-1
! Old version commented just below
do l = min(size(cfg_seniority_index,1)-1, i+2), size(cfg_seniority_index,1)-1, 2
if (cfg_seniority_index(l) >= ncfgpersomo) then
ncfgpersomo = cfg_seniority_index(l)
endif
enddo
!k = 0
!if ((i+2+k) < size(cfg_seniority_index,1)) then
! do while(cfg_seniority_index(i+2+k) < ncfgpersomo)
! k = k + 2
! if ((i+2+k) >= size(cfg_seniority_index,1)) then
! exit
! endif
! ncfgpersomo = cfg_seniority_index(i+2+k)
! enddo
!endif
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
! l = i+k+2
! Loop over l with a constraint to ensure that l <= size(cfg_seniority_index,1)-1
! Old version commented just below
do l = min(size(cfg_seniority_index,1)-1, i+2), size(cfg_seniority_index,1)-1, 2
if (cfg_seniority_index(l) >= ncfgprev) then
ncfgprev = cfg_seniority_index(l)
endif
enddo
!k = 0
!if ((i+2+k) < size(cfg_seniority_index,1)) then
! do while(cfg_seniority_index(i+2+k) < ncfgprev)
! k = k + 2
! if ((i+2+k) >= size(cfg_seniority_index,1)) then
! exit
! endif
! ncfgprev = cfg_seniority_index(i+2+k)
! enddo
!endif
endif
enddo
END_PROVIDER
subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout)
use bitmasks
implicit none
BEGIN_DOC
! Documentation for get_phase_qp_to_cfg
!
! This function converts from (aaaa)(bbbb)
! notation to (ab)(ab)(ab)(ab)
! notation.
! The cfgCI code works in (ab)(ab)(ab)(ab)
! notation throughout.
END_DOC
integer(bit_kind),intent(in) :: Ialpha(N_int)
integer(bit_kind),intent(in) :: Ibeta(N_int)
real*8,intent(out) :: phaseout
integer(bit_kind) :: mask, deta(N_int), detb(N_int)
integer :: nbetas
integer :: k
! Initialize deta and detb
deta = Ialpha
detb = Ibeta
! Find how many alpha electrons there are in all the N_ints
integer :: Na(N_int)
do k=1,N_int
Na(k) = popcnt(deta(k))
enddo
integer :: shift, ipos, nperm
phaseout = 1.d0
do k=1,N_int
do while(detb(k) /= 0_bit_kind)
! Find the lowest beta electron and clear it
ipos = trailz(detb(k))
detb(k) = ibclr(detb(k),ipos)
! Create a mask will all MOs higher than the beta electron
mask = not(shiftl(1_bit_kind,ipos + 1) - 1_bit_kind)
! Apply the mask to the alpha string to count how many electrons to cross
nperm = popcnt( iand(mask, deta(k)) )
! Count how many alpha electrons are above the beta electron in the other integers
nperm = nperm + sum(Na(k+1:N_int))
if (iand(nperm,1) == 1) then
phaseout = -phaseout
endif
enddo
enddo
end subroutine get_phase_qp_to_cfg
BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,2)]
&BEGIN_PROVIDER [ integer, rowsmax]
&BEGIN_PROVIDER [ integer, colsmax]
use cfunctions
implicit none
BEGIN_DOC
! Documentation for AIJpqMatrixList
! The prototype matrix containing the <I|E_{pq}|J>
! matrices for each I,J somo pair and orb ids.
END_DOC
integer i,j,k,l
integer*8 Isomo, Jsomo, tmpsomo
Isomo = 0
Jsomo = 0
integer rows, cols, nsomoi, nsomoj
rows = -1
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
!allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
! Type
! 1. SOMO -> SOMO
do i = 2-iand(nsomomin,1), NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1
do j = i-2,i-2, 2
Jsomo = ISHFT(1_8,j)-1
if(j .GT. NSOMOMax .OR. j .LT. 0) then
cycle
end if
do k = 1,i
do l = 1,i
! Define Jsomo
if(k.NE.l)then
Jsomo = IBCLR(Isomo, k-1)
Jsomo = IBCLR(Jsomo, l-1)
nsomoi = i
nsomoj = j
else
Isomo = ISHFT(1_8,i)-1
Jsomo = ISHFT(1_8,i)-1
nsomoi = i
nsomoj = i
endif
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
if(rowsmax .LT. rows) then
rowsmax = rows
end if
if(colsmax .LT. cols) then
colsmax = cols
end if
! i -> j
AIJpqMatrixDimsList(nsomoi,nsomoj,1,k,l,1) = rows
AIJpqMatrixDimsList(nsomoi,nsomoj,1,k,l,2) = cols
end do
end do
end do
end do
! Type
! 2. DOMO -> VMO
do i = 0+iand(nsomomin,1), NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1
tmpsomo = ISHFT(1_8,i+2)-1
do j = i+2,i+2, 2
Jsomo = ISHFT(1_8,j)-1
if(j .GT. NSOMOMax .OR. j .LT. 0) then
cycle
end if
do k = 1,j
do l = 1,j
if(k .NE. l) then
Isomo = IBCLR(tmpsomo,k-1)
Isomo = IBCLR(Isomo,l-1)
! Define Jsomo
Jsomo = ISHFT(1_8,j)-1;
nsomoi = i
nsomoj = j
else
Isomo = ISHFT(1_8,j)-1
Isomo = IBCLR(Isomo,1-1)
Isomo = IBCLR(Isomo,j-1)
Jsomo = ISHFT(1_8,j)-1
Isomo = ISHFT(1_8,j)-1
nsomoi = j
nsomoj = j
endif
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
if(rowsmax .LT. rows) then
rowsmax = rows
end if
if(colsmax .LT. cols) then
colsmax = cols
end if
! i -> j
AIJpqMatrixDimsList(nsomoi,nsomoj,2,k,l,1) = rows
AIJpqMatrixDimsList(nsomoi,nsomoj,2,k,l,2) = cols
end do
end do
end do
end do
! Type
! 3. SOMO -> VMO
!print *,"Doing SOMO->VMO"
do i = 2-iand(nsomomin,1), NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1
do j = i,i, 2
Jsomo = ISHFT(1_8,j)-1
if(j .GT. NSOMOMax .OR. j .LE. 0) then
cycle
end if
do k = 1,i
do l = 1,i
if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,l-1)
Jsomo = ISHFT(1_8,j+1)-1
Jsomo = IBCLR(Jsomo,k-1)
else
Isomo = ISHFT(1_8,i)-1
Jsomo = ISHFT(1_8,j)-1
endif
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
if(rowsmax .LT. rows) then
rowsmax = rows
end if
if(colsmax .LT. cols) then
colsmax = cols
end if
! i -> j
AIJpqMatrixDimsList(i,j,3,k,l,1) = rows
AIJpqMatrixDimsList(i,j,3,k,l,2) = cols
end do
end do
end do
end do
! Type
! 4. DOMO -> VMO
!print *,"Doing DOMO->SOMO"
do i = 2-iand(nsomomin,1), NSOMOMax, 2
do j = i,i, 2
if(j .GT. NSOMOMax .OR. j .LE. 0) then
cycle
end if
do k = 1,i
do l = 1,i
if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,k+1-1)
Jsomo = ISHFT(1_8,j+1)-1
Jsomo = IBCLR(Jsomo,l-1)
else
Isomo = ISHFT(1_8,i)-1
Jsomo = ISHFT(1_8,j)-1
endif
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
if(rowsmax .LT. rows) then
rowsmax = rows
end if
if(colsmax .LT. cols) then
colsmax = cols
end if
! i -> j
AIJpqMatrixDimsList(i,j,4,k,l,1) = rows
AIJpqMatrixDimsList(i,j,4,k,l,2) = cols
end do
end do
end do
end do
END_PROVIDER
BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,NBFMax,NBFMax)]
use cfunctions
implicit none
BEGIN_DOC
! Documentation for AIJpqMatrixList
! The prototype matrix containing the <I|E_{pq}|J>
! matrices for each I,J somo pair and orb ids.
!
! Due to the different types of excitations which
! include DOMOs and VMOs two prototype DOMOs and two
! prototype VMOs are needed. Therefore
! the 4th and 5th dimensions are NSOMOMax+4 and NSOMOMax+4
! respectively.
!
! The type of excitations are ordered as follows:
! Type 1 - SOMO -> SOMO
! Type 2 - DOMO -> VMO
! Type 3 - SOMO -> VMO
! Type 4 - DOMO -> SOMO
END_DOC
integer i,j,k,l, orbp, orbq, ri, ci
orbp = 0
orbq = 0
integer*8 Isomo, Jsomo, tmpsomo
Isomo = 0
Jsomo = 0
integer rows, cols, nsomoi, nsomoj
rows = -1
cols = -1
integer*8 MS
MS = 0
touch AIJpqMatrixDimsList
real*8,dimension(:,:),allocatable :: meMatrix
integer maxdim
!maxdim = max(rowsmax,colsmax)
! allocate matrix
!allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
! Type
! 1. SOMO -> SOMO
do i = 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
do k = 1,i
do l = 1,i
! Define Jsomo
if(k .NE. l) then
Jsomo = IBCLR(Isomo, k-1)
Jsomo = IBCLR(Jsomo, l-1)
nsomoi = i
nsomoj = j
else
Isomo = ISHFT(1_8,i)-1
Jsomo = ISHFT(1_8,i)-1
nsomoi = i
nsomoj = i
endif
AIJpqContainer(nsomoi,nsomoj,1,k,l,:,:) = 0.0d0
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
orbp = k
orbq = l
allocate(meMatrix(rows,cols))
meMatrix = 0.0d0
! fill matrix
call getApqIJMatrixDriver(Isomo, &
Jsomo, &
orbp, &
orbq, &
MS, &
NMO, &
meMatrix, &
rows, &
cols)
! i -> j
do ri = 1,rows
do ci = 1,cols
AIJpqContainer(nsomoi,nsomoj,1,k,l,ri,ci) = meMatrix(ri, ci)
end do
end do
deallocate(meMatrix)
end do
end do
end do
end do
! Type
! 2. DOMO -> VMO
do i = 0, NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1
tmpsomo = ISHFT(1_8,i+2)-1
do j = i+2,i+2, 2
if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
Jsomo = ISHFT(1_8,j)-1
do k = 1,j
do l = 1,j
if(k .NE. l) then
Isomo = IBCLR(tmpsomo,k-1)
Isomo = IBCLR(Isomo,l-1)
! Define Jsomo
Jsomo = ISHFT(1_8,j)-1;
nsomoi = i
nsomoj = j
else
Isomo = ISHFT(1_8,j)-1
Isomo = IBCLR(Isomo,1-1)
Isomo = IBCLR(Isomo,j-1)
Jsomo = ISHFT(1_8,j)-1
Isomo = ISHFT(1_8,j)-1
nsomoi = j
nsomoj = j
endif
AIJpqContainer(nsomoi,nsomoj,2,k,l,:,:) = 0.0d0
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
orbp = k
orbq = l
allocate(meMatrix(rows,cols))
meMatrix = 0.0d0
! fill matrix
call getApqIJMatrixDriver(Isomo, &
Jsomo, &
orbp, &
orbq, &
MS, &
NMO, &
meMatrix, &
rows, &
cols)
! i -> j
do ri = 1,rows
do ci = 1,cols
AIJpqContainer(nsomoi,nsomoj,2,k,l,ri,ci) = meMatrix(ri, ci)
end do
end do
deallocate(meMatrix)
end do
end do
end do
end do
! Type
! 3. SOMO -> VMO
do i = 2, NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1
do j = i,i, 2
Jsomo = ISHFT(1_8,j)-1
if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
do k = 1,i
do l = 1,i
if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,l-1)
Jsomo = ISHFT(1_8,j+1)-1
Jsomo = IBCLR(Jsomo,k-1)
else
Isomo = ISHFT(1_8,i)-1
Jsomo = ISHFT(1_8,j)-1
endif
AIJpqContainer(i,j,3,k,l,:,:) = 0.0d0
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
orbp = k
orbq = l
allocate(meMatrix(rows,cols))
meMatrix = 0.0d0
! fill matrix
call getApqIJMatrixDriver(Isomo, &
Jsomo, &
orbp, &
orbq, &
MS, &
NMO, &
meMatrix, &
rows, &
cols)
! i -> j
do ri = 1,rows
do ci = 1,cols
AIJpqContainer(i,j,3,k,l,ri,ci) = meMatrix(ri, ci)
end do
end do
deallocate(meMatrix)
end do
end do
end do
end do
! Type
! 4. DOMO -> SOMO
do i = 2, NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1
do j = i,i, 2
Jsomo = ISHFT(1_8,i)-1
if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
do k = 1,i
do l = 1,i
if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,k-1)
Jsomo = ISHFT(1_8,j+1)-1
Jsomo = IBCLR(Jsomo,l+1-1)
else
Isomo = ISHFT(1_8,i)-1
Jsomo = ISHFT(1_8,j)-1
endif
AIJpqContainer(i,j,4,k,l,:,:) = 0.0d0
call getApqIJMatrixDims(Isomo, &
Jsomo, &
MS, &
rows, &
cols)
orbp = k
orbq = l
allocate(meMatrix(rows,cols))
meMatrix = 0.0d0
! fill matrix
call getApqIJMatrixDriver(Isomo, &
Jsomo, &
orbp, &
orbq, &
MS, &
NMO, &
meMatrix, &
rows, &
cols)
! i -> j
do ri = 1,rows
do ci = 1,cols
AIJpqContainer(i,j,4,k,l,ri,ci) = meMatrix(ri, ci)
end do
end do
deallocate(meMatrix)
end do
end do
end do
end do
END_PROVIDER
!!!!!!
BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)]
&BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF)]
&BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)]
use cfunctions
use bitmasks
implicit none
BEGIN_DOC
! Documentation for DetToCSFTransformationMatrix
! Provides the matrix of transformatons for the
! conversion between determinant to CSF basis (in BFs)
END_DOC
integer(bit_kind) :: mask(N_int), Ialpha(N_int),Ibeta(N_int)
integer :: rows, cols, i, j, k
integer :: startdet, enddet
integer*8 MS, Isomo, Idomo
integer ndetI
integer :: getNSOMO
real*8,dimension(:,:),allocatable :: tempBuffer
real*8,dimension(:),allocatable :: tempCoeff
real*8 :: norm_det1, phasedet
norm_det1 = 0.d0
MS = elec_alpha_num - elec_beta_num
! initialization
psi_coef_config = 0.d0
DetToCSFTransformationMatrix(0,:,:) = 1.d0
do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2
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))))
allocate(tempBuffer(bfIcfg,ndetI))
call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer)
DetToCSFTransformationMatrix(i,1:bfIcfg,1:ndetI) = tempBuffer(1:bfIcfg,1:ndetI)
deallocate(tempBuffer)
enddo
integer s, bfIcfg
integer countcsf
countcsf = 0
integer countdet
countdet = 0
integer idx
integer istate
istate = 1
phasedet = 1.0d0
do i = 1,N_configuration
startdet = psi_configuration_to_psi_det(1,i)
enddet = psi_configuration_to_psi_det(2,i)
ndetI = enddet-startdet+1
allocate(tempCoeff(ndetI))
countdet = 1
do j = startdet, enddet
idx = psi_configuration_to_psi_det_data(j)
Ialpha(:) = psi_det(:,1,idx)
Ibeta(:) = psi_det(:,2,idx)
call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet)
tempCoeff(countdet) = psi_coef(idx, istate)*phasedet
norm_det1 += tempCoeff(countdet)*tempCoeff(countdet)
countdet += 1
enddo
s = 0
do k=1,N_int
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))))
! perhaps blocking with CFGs of same seniority
! can be more efficient
allocate(tempBuffer(bfIcfg,ndetI))
tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI)
call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1), size(psi_coef_config,1))
!call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1)
deallocate(tempCoeff)
deallocate(tempBuffer)
psi_config_data(i,1) = countcsf + 1
countcsf += bfIcfg
psi_config_data(i,2) = countcsf
enddo
END_PROVIDER