mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-10-13 03:11:31 +02:00
First CSF davidson working
This commit is contained in:
parent
554579492b
commit
b87e87b740
49
src/csf/cfgCI_interface.f90
Normal file
49
src/csf/cfgCI_interface.f90
Normal file
@ -0,0 +1,49 @@
|
||||
module cfunctions
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
interface
|
||||
subroutine printcfglist(nint, ncfgs, cfglist) bind(C, name='printCFGList')
|
||||
import C_INT32_T, C_INT64_T
|
||||
integer(kind=C_INT32_T) :: nint
|
||||
integer(kind=C_INT32_T) :: ncfgs
|
||||
integer(kind=C_INT64_T) :: cfglist(nint,2,ncfgs)
|
||||
end subroutine printcfglist
|
||||
end interface
|
||||
interface
|
||||
subroutine getApqIJMatrixDims(Isomo, Jsomo, MS, rowsout, colsout) &
|
||||
bind(C, name='getApqIJMatrixDims')
|
||||
import C_INT32_T, C_INT64_T
|
||||
integer(kind=C_INT64_T),value,intent(in) :: Isomo ! CSFI
|
||||
integer(kind=C_INT64_T),value,intent(in) :: Jsomo ! CSFJ
|
||||
integer(kind=C_INT64_T),value,intent(in) :: MS ! Ms = 2*Spin
|
||||
integer(kind=C_INT32_T),intent(out):: rowsout
|
||||
integer(kind=C_INT32_T),intent(out):: colsout
|
||||
end subroutine getApqIJMatrixDims
|
||||
end interface
|
||||
interface
|
||||
subroutine getApqIJMatrixDriver(Isomo, Jsomo, orbp, orbq, &
|
||||
MS, NMO, CSFICSFJApqIJ, rowsmax, colsmax) bind(C, name='getApqIJMatrixDriverArrayInp')
|
||||
import C_INT32_T, C_INT64_T, C_DOUBLE
|
||||
integer(kind=C_INT64_T),value,intent(in) :: Isomo
|
||||
integer(kind=C_INT64_T),value,intent(in) :: Jsomo
|
||||
integer(kind=C_INT32_T),value,intent(in) :: orbp
|
||||
integer(kind=C_INT32_T),value,intent(in) :: orbq
|
||||
integer(kind=C_INT64_T),value,intent(in) :: MS
|
||||
integer(kind=C_INT64_T),value,intent(in) :: NMO
|
||||
integer(kind=C_INT32_T),intent(in) :: rowsmax
|
||||
integer(kind=C_INT32_T),intent(in) :: colsmax
|
||||
real (kind=C_DOUBLE ),intent(out) :: CSFICSFJApqIJ(rowsmax,colsmax)
|
||||
!integer(kind=C_INT32_T),dimension(rowApqIJ,colApqIJ) :: ApqIJ
|
||||
end subroutine getApqIJMatrixDriver
|
||||
end interface
|
||||
interface
|
||||
subroutine getCSFtoDETTransformationMatrix(Isomo,&
|
||||
MS, rowsmax, colsmax, csftodetmatrix) bind(C, name='convertCSFtoDetBasis')
|
||||
import C_INT32_T, C_INT64_T, C_DOUBLE
|
||||
integer(kind=C_INT64_T),value,intent(in) :: Isomo
|
||||
integer(kind=C_INT64_T),value,intent(in) :: MS
|
||||
integer(kind=C_INT32_T),intent(in) :: rowsmax
|
||||
integer(kind=C_INT32_T),intent(in) :: colsmax
|
||||
real (kind=C_DOUBLE ),intent(out) :: csftodetmatrix(rowsmax,colsmax)
|
||||
end subroutine getCSFtoDETTransformationMatrix
|
||||
end interface
|
||||
end module cfunctions
|
1763
src/csf/cfgCI_utils.c
Normal file
1763
src/csf/cfgCI_utils.c
Normal file
File diff suppressed because it is too large
Load Diff
@ -1,4 +1,13 @@
|
||||
use bitmasks
|
||||
|
||||
BEGIN_PROVIDER [ integer, spin_multiplicity ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! n_alpha - n_beta + 1
|
||||
END_DOC
|
||||
spin_multiplicity = elec_alpha_num - elec_beta_num + 1
|
||||
END_PROVIDER
|
||||
|
||||
subroutine configuration_of_det(d,o,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
29
src/csf/connected_to_ref.irp.f
Normal file
29
src/csf/connected_to_ref.irp.f
Normal file
@ -0,0 +1,29 @@
|
||||
integer*8 function configuration_search_key(cfg,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns an integer*8 corresponding to a determinant index for searching.
|
||||
! The left-most 8 bits contain the number of open shells+1. This ensures that the CSF
|
||||
! are packed with the same seniority.
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: cfg(Nint,2)
|
||||
integer :: i, n_open_shells
|
||||
integer*8 :: mask
|
||||
|
||||
i = shiftr(elec_alpha_num, bit_kind_shift)+1
|
||||
configuration_search_key = int(shiftr(ior(cfg(i,1),cfg(i,2)),1)+sum(cfg),8)
|
||||
|
||||
mask = X'00FFFFFFFFFFFFFF'
|
||||
configuration_search_key = iand(mask,configuration_search_key)
|
||||
|
||||
n_open_shells = 1
|
||||
do i=1,Nint
|
||||
if (cfg(i,1) == 0_bit_kind) cycle
|
||||
n_open_shells = n_open_shells + popcnt(cfg(i,1))
|
||||
enddo
|
||||
mask = n_open_shells
|
||||
mask = shiftl(mask,56)
|
||||
configuration_search_key = ior (mask,configuration_search_key)
|
||||
|
||||
end
|
148
src/csf/conversion.irp.f
Normal file
148
src/csf/conversion.irp.f
Normal file
@ -0,0 +1,148 @@
|
||||
subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
|
||||
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, intent(in) :: N_st
|
||||
double precision, intent(in) :: psi_coef_det_in(N_det,N_st)
|
||||
double precision, intent(out) :: psi_coef_cfg_out(n_CSF,N_st)
|
||||
integer*8 :: Isomo, Idomo, mask
|
||||
integer(bit_kind) :: Ialpha(N_int) ,Ibeta(N_int)
|
||||
integer :: rows, cols, i, j, k
|
||||
integer :: startdet, enddet
|
||||
integer :: ndetI
|
||||
integer :: getNSOMO
|
||||
double precision,allocatable :: tempBuffer(:,:)
|
||||
double precision,allocatable :: tempCoeff(:,:)
|
||||
double precision :: phasedet
|
||||
integer :: idx
|
||||
|
||||
! initialization
|
||||
psi_coef_cfg_out(:,1) = 0.d0
|
||||
|
||||
integer s, bfIcfg
|
||||
integer countcsf
|
||||
countcsf = 0
|
||||
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,N_st))
|
||||
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)
|
||||
do k=1,N_st
|
||||
tempCoeff(j-startdet+1,k) = psi_coef_det_in(idx, k)*phasedet
|
||||
enddo
|
||||
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, N_st, ndetI, 1.d0, tempBuffer, size(tempBuffer,1),&
|
||||
tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_cfg_out(countcsf+1,1),&
|
||||
size(psi_coef_cfg_out,1))
|
||||
|
||||
deallocate(tempCoeff)
|
||||
deallocate(tempBuffer)
|
||||
countcsf += bfIcfg
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Documentation for convertCSFtoDET
|
||||
! This function converts the wavefunction
|
||||
! in CFG basis to DET basis using the
|
||||
! transformation matrix provided before.
|
||||
END_DOC
|
||||
integer, intent(in) :: N_st
|
||||
double precision,intent(in) :: psi_coef_cfg_in(n_CSF,N_st)
|
||||
double precision,intent(out) :: psi_coef_det(N_det,N_st)
|
||||
double precision :: tmp_psi_coef_det(maxDetDimPerBF,N_st)
|
||||
integer :: s, bfIcfg
|
||||
integer :: countcsf
|
||||
integer(bit_kind) :: Ialpha(N_int), Ibeta(N_int)
|
||||
integer :: rows, cols, i, j, k
|
||||
integer :: startdet, enddet
|
||||
integer :: ndetI
|
||||
integer :: getNSOMO
|
||||
double precision,allocatable :: tempBuffer(:,:)
|
||||
double precision,allocatable :: tempCoeff (:,:)
|
||||
double precision :: phasedet
|
||||
integer :: idx
|
||||
|
||||
countcsf = 0
|
||||
|
||||
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
|
||||
|
||||
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))))
|
||||
|
||||
allocate(tempCoeff(bfIcfg,N_st))
|
||||
|
||||
do k=1,N_st
|
||||
do j = 1,bfIcfg
|
||||
tempCoeff(j,k) = psi_coef_cfg_in(countcsf+j,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
countcsf += bfIcfg
|
||||
! perhaps blocking with CFGs of same seniority
|
||||
! can be more efficient
|
||||
allocate(tempBuffer(bfIcfg,ndetI))
|
||||
tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI)
|
||||
|
||||
call dgemm('T','N', ndetI, N_st, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1),&
|
||||
tempCoeff, size(tempCoeff,1), 0.d0, tmp_psi_coef_det, &
|
||||
size(tmp_psi_coef_det,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)
|
||||
do k=1,N_st
|
||||
psi_coef_det(idx,k) = tmp_psi_coef_det(j-startdet+1,k) * phasedet
|
||||
enddo
|
||||
enddo
|
||||
|
||||
deallocate(tempCoeff)
|
||||
deallocate(tempBuffer)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
641
src/csf/sigma_vector.irp.f
Normal file
641
src/csf/sigma_vector.irp.f
Normal file
@ -0,0 +1,641 @@
|
||||
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
|
||||
integer MS
|
||||
integer ncfgpersomo
|
||||
detDimperBF = 0
|
||||
MS = elec_alpha_num-elec_beta_num
|
||||
! 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
|
||||
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
|
||||
!if(cfg_seniority_index(i+2) == -1) EXIT
|
||||
!if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF
|
||||
ncfgprev = cfg_seniority_index(i)
|
||||
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(N_int), deta(N_int), detb(N_int)
|
||||
integer :: nbetas
|
||||
integer :: count, k
|
||||
if (N_int >1 ) then
|
||||
stop 'TODO: get_phase_qp_to_cfg '
|
||||
endif
|
||||
|
||||
nbetas = 0
|
||||
mask = 0_bit_kind
|
||||
count = 0
|
||||
deta = Ialpha
|
||||
detb = Ibeta
|
||||
! remove the domos
|
||||
mask = IAND(deta,detb)
|
||||
deta = IEOR(deta,mask)
|
||||
detb = IEOR(detb,mask)
|
||||
mask = 0
|
||||
phaseout = 1.0
|
||||
k = 1
|
||||
do while((deta(k)).GT.0_8)
|
||||
mask(k) = ISHFT(1_8,count)
|
||||
if(POPCNT(IAND(deta(k),mask(k))).EQ.1)then
|
||||
if(IAND(nbetas,1).EQ.0) then
|
||||
phaseout *= 1.0d0
|
||||
else
|
||||
phaseout *= -1.0d0
|
||||
endif
|
||||
deta(k) = IEOR(deta(k),mask(k))
|
||||
else
|
||||
if(POPCNT(IAND(detb(k),mask(k))).EQ.1) then
|
||||
nbetas += 1
|
||||
detb(k) = IEOR(detb(k),mask(k))
|
||||
endif
|
||||
endif
|
||||
count += 1
|
||||
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
|
309
src/csf/tree_utils.c
Normal file
309
src/csf/tree_utils.c
Normal file
@ -0,0 +1,309 @@
|
||||
#include "tree_utils.h"
|
||||
|
||||
void buildTree(Tree *bftree,
|
||||
Node **inode,
|
||||
int isomo,
|
||||
int izeros,
|
||||
int icpl,
|
||||
int NSOMOMax,
|
||||
int MSmax){
|
||||
|
||||
// Find the maximum parallel couplings 0
|
||||
// the maximum anti-parallel couplings 1
|
||||
int zeromax = MSmax + (NSOMOMax-MSmax)/2;
|
||||
int onemax = NSOMOMax - zeromax;
|
||||
|
||||
// Exit condition
|
||||
if(isomo > NSOMOMax || icpl < 0 || izeros > zeromax ) return;
|
||||
|
||||
// If we find a valid BF assign its address
|
||||
if(isomo == NSOMOMax){
|
||||
(*inode)->addr = bftree->rootNode->addr;
|
||||
bftree->rootNode->addr += 1;
|
||||
return;
|
||||
}
|
||||
|
||||
// Call 0 branch
|
||||
if(((*inode)->C0) == NULL && izeros+1 <= zeromax){
|
||||
((*inode)->C0) = malloc(sizeof(Node));
|
||||
(*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
|
||||
buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
|
||||
}
|
||||
else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
|
||||
|
||||
// Call 1 branch
|
||||
if(((*inode)->C1) == NULL && icpl-1 >= 0){
|
||||
((*inode)->C1) = malloc(sizeof(Node));
|
||||
(*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
|
||||
buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
|
||||
}
|
||||
else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void buildTreeDriver(Tree *bftree, int NSOMO, int MS, int *NBF){
|
||||
int isomo = 0; // counts the total number of SOMO's
|
||||
int izeros= 0; // Counts the total number of parallel coupings (i.e. 0's)
|
||||
int icpl = 0; // keep track of the ith ms (cannot be -ve)
|
||||
int addr = 0; // Counts the total BF's
|
||||
|
||||
buildTree(bftree, &(bftree->rootNode), isomo, izeros, icpl, NSOMO, MS);
|
||||
|
||||
*NBF = bftree->rootNode->addr;
|
||||
}
|
||||
|
||||
|
||||
void getIthBF(Node *inode, int isomo, bool foundBF, int NSOMOMax, int getaddr, int *vecBF){
|
||||
// Exit condition
|
||||
if(foundBF) return;
|
||||
if(isomo > NSOMOMax) return;
|
||||
if(inode == NULL) return;
|
||||
|
||||
if(isomo == NSOMOMax){
|
||||
if(inode->addr == getaddr){
|
||||
for(int i = NSOMOMax-1; i > -1; i--){
|
||||
vecBF[i] = inode->cpl;
|
||||
inode = inode->PREV;
|
||||
}
|
||||
foundBF = true;
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
// Recurse to C0
|
||||
if(inode->C0 != NULL){
|
||||
getIthBF(inode->C0, isomo+1, foundBF, NSOMOMax, getaddr, vecBF);
|
||||
}
|
||||
// Recurse to C1
|
||||
if(inode->C1 != NULL){
|
||||
getIthBF(inode->C1, isomo+1, foundBF, NSOMOMax, getaddr, vecBF);
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void getIthBFDriver(Tree *bftree, int NSOMOMax, int getaddr, int *vecBF){
|
||||
int isomo = 0;
|
||||
bool foundBF = false;
|
||||
getIthBF((bftree->rootNode), isomo, foundBF, NSOMOMax, getaddr, vecBF);
|
||||
}
|
||||
|
||||
void genDets(Tree *dettree,
|
||||
Node **inode,
|
||||
int isomo,
|
||||
int izeros,
|
||||
int icpl,
|
||||
int NSOMOMax,
|
||||
int MSmax){
|
||||
|
||||
// Find the maximum parallel couplings 0
|
||||
// the maximum anti-parallel couplings 1
|
||||
int zeromax = MSmax + (NSOMOMax-MSmax)/2;
|
||||
int onemax = NSOMOMax - zeromax;
|
||||
|
||||
// Exit condition
|
||||
if(isomo > NSOMOMax || izeros > zeromax || abs(icpl) > onemax) return;
|
||||
|
||||
// If we find a valid BF assign its address
|
||||
if(isomo == NSOMOMax){
|
||||
(*inode)->addr = dettree->rootNode->addr;
|
||||
dettree->rootNode->addr += 1;
|
||||
return;
|
||||
}
|
||||
|
||||
// Call 0 branch
|
||||
if(((*inode)->C0) == NULL && izeros+1 <= zeromax){
|
||||
((*inode)->C0) = malloc(sizeof(Node));
|
||||
(*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
|
||||
genDets(dettree, &(*inode)->C0, isomo+1, izeros+1, icpl+0, NSOMOMax, MSmax);
|
||||
}
|
||||
else genDets(dettree, &(*inode)->C0, isomo+1, izeros+1, icpl+0, NSOMOMax, MSmax);
|
||||
|
||||
// Call 1 branch
|
||||
if(((*inode)->C1) == NULL && abs(icpl+1) <= onemax){
|
||||
((*inode)->C1) = malloc(sizeof(Node));
|
||||
(*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
|
||||
genDets(dettree, &(*inode)->C1, isomo+1, izeros+0, icpl+1, NSOMOMax, MSmax);
|
||||
}
|
||||
else genDets(dettree, &(*inode)->C1, isomo+1, izeros+0, icpl+1, NSOMOMax, MSmax);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void genDetsDriver(Tree *dettree, int NSOMO, int MS, int *Ndets){
|
||||
int isomo = 0; // counts the total number of SOMO's
|
||||
int izeros= 0; // Counts the total number of parallel coupings (i.e. 0's)
|
||||
int icpl = 0; // keep track of the ith ms (cannot be -ve)
|
||||
int addr = 0; // Counts the total BF's
|
||||
|
||||
genDets(dettree, &(dettree->rootNode), isomo, izeros, icpl, NSOMO, MS);
|
||||
|
||||
*Ndets = dettree->rootNode->addr;
|
||||
}
|
||||
|
||||
void getIthDet(Node *inode, int isomo, bool foundBF, int NSOMOMax, int getaddr, int *vecBF){
|
||||
// Exit condition
|
||||
if(foundBF) return;
|
||||
if(isomo > NSOMOMax) return;
|
||||
if(inode == NULL) return;
|
||||
|
||||
if(isomo == NSOMOMax){
|
||||
if(inode->addr == getaddr){
|
||||
for(int i = NSOMOMax-1; i > -1; i--){
|
||||
vecBF[i] = inode->cpl;
|
||||
inode = inode->PREV;
|
||||
}
|
||||
foundBF = true;
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
// Recurse to C0
|
||||
if(inode->C0 != NULL){
|
||||
getIthDet(inode->C0, isomo+1, foundBF, NSOMOMax, getaddr, vecBF);
|
||||
}
|
||||
// Recurse to C1
|
||||
if(inode->C1 != NULL){
|
||||
getIthDet(inode->C1, isomo+1, foundBF, NSOMOMax, getaddr, vecBF);
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void getIthDetDriver(Tree *dettree, int NSOMOMax, int getaddr, int *vecBF){
|
||||
int isomo = 0;
|
||||
bool foundBF = false;
|
||||
getIthDet((dettree->rootNode), isomo, foundBF, NSOMOMax, getaddr, vecBF);
|
||||
}
|
||||
|
||||
void findAddofDet(Node *inode, int isomo, bool foundDet, int NSOMOMax, int *inpdet, int *addr){
|
||||
// Exit condition
|
||||
if(foundDet) return;
|
||||
if(isomo == NSOMOMax){
|
||||
foundDet = true;
|
||||
*addr = inode->addr;
|
||||
return;
|
||||
}
|
||||
|
||||
// Recurse to C0
|
||||
if(inpdet[isomo] == 0){
|
||||
if(inode->C0 != NULL){
|
||||
findAddofDet(inode->C0, isomo+1, foundDet, NSOMOMax, inpdet, addr);
|
||||
}
|
||||
else{
|
||||
*addr = -1;
|
||||
return;
|
||||
}
|
||||
}
|
||||
else{
|
||||
// Recurse to C1
|
||||
if(inode->C1 != NULL){
|
||||
findAddofDet(inode->C1, isomo+1, foundDet, NSOMOMax, inpdet, addr);
|
||||
}
|
||||
else{
|
||||
*addr = -1;
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void findAddofDetDriver(Tree *dettree, int NSOMOMax, int *inpdet, int *addr){
|
||||
*addr = -1;
|
||||
int isomo = 0;
|
||||
bool foundDet = false;
|
||||
findAddofDet((dettree->rootNode), isomo, foundDet, NSOMOMax, inpdet, addr);
|
||||
}
|
||||
|
||||
void getDetlist(Node *inode, int isomo, int NSOMOMax, int *vecBF, int *detlist){
|
||||
// Exit condition
|
||||
if(isomo > NSOMOMax) return;
|
||||
if(inode == NULL) return;
|
||||
|
||||
if(isomo == NSOMOMax){
|
||||
int idet=0;
|
||||
for(int k=0;k<NSOMOMax;k++){
|
||||
if(vecBF[k] == 1) idet = idet | (1<<(NSOMOMax-1-k));
|
||||
}
|
||||
detlist[inode->addr]=idet;
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
// Recurse to C0
|
||||
if(inode->C0 != NULL){
|
||||
vecBF[isomo] = 0;
|
||||
getDetlist(inode->C0, isomo+1, NSOMOMax, vecBF, detlist);
|
||||
}
|
||||
// Recurse to C1
|
||||
if(inode->C1 != NULL){
|
||||
vecBF[isomo] = 1;
|
||||
getDetlist(inode->C1, isomo+1, NSOMOMax, vecBF, detlist);
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
void getDetlistDriver(Tree *dettree, int NSOMOMax, int *detlist){
|
||||
int isomo = 0;
|
||||
int vecBF[NSOMOMax];
|
||||
if(dettree->rootNode->addr > 1) getDetlist((dettree->rootNode), isomo, NSOMOMax, vecBF, detlist);
|
||||
}
|
||||
|
||||
void genDetBasis(Tree *dettree, int Isomo, int MS, int *ndets){
|
||||
|
||||
int NSOMO=0;
|
||||
getSetBits(Isomo, &NSOMO);
|
||||
genDetsDriver(dettree, NSOMO, MS, ndets);
|
||||
|
||||
}
|
||||
|
||||
void callBlasMatxMat(double *A, int rowA, int colA, double *B, int rowB, int colB, double *C, bool transA, bool transB){
|
||||
int m = rowA;
|
||||
int k = colA;
|
||||
int n = colB;
|
||||
double alpha = 1.0;
|
||||
double beta = 0.0;
|
||||
int val = 0;
|
||||
if (transA) val |= 0x1;
|
||||
if (transB) val |= 0x2;
|
||||
|
||||
switch (val) {
|
||||
case 0: // notransA, notransB
|
||||
m = rowA;
|
||||
n = colB;
|
||||
k = colA;
|
||||
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
|
||||
m, n, k, alpha, A, k, B, n, beta, C, n);
|
||||
break;
|
||||
case 1: // transA, notransB
|
||||
m = colA;
|
||||
n = colB;
|
||||
k = rowA;
|
||||
cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
|
||||
m, n, k, alpha, A, colA, B, n, beta, C, n);
|
||||
break;
|
||||
case 2: // notransA, transB
|
||||
//m = rowA;
|
||||
//n = rowB;
|
||||
//k = colB;
|
||||
m = rowA;
|
||||
n = rowB;
|
||||
k = colA;
|
||||
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
|
||||
m, n, k, alpha, A, k, B, colB, beta, C, n);
|
||||
break;
|
||||
case 3: // transA, transB
|
||||
m = colA;
|
||||
n = rowB;
|
||||
k = rowA;
|
||||
cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans,
|
||||
m, n, k, alpha, A, colA, B, colB, beta, C, n);
|
||||
break;
|
||||
default:
|
||||
printf("Impossible !!!!\n");
|
||||
break;
|
||||
}
|
||||
}
|
87
src/csf/tree_utils.h
Normal file
87
src/csf/tree_utils.h
Normal file
@ -0,0 +1,87 @@
|
||||
#ifndef TREE_UTILS_H
|
||||
#define TREE_UTILS_H
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdbool.h>
|
||||
#include <stdint.h>
|
||||
#include <math.h>
|
||||
|
||||
typedef struct bin_node Node;
|
||||
typedef struct bin_tree Tree;
|
||||
struct bin_node {
|
||||
Node *C0;
|
||||
Node *C1;
|
||||
Node *PREV;
|
||||
int addr;
|
||||
int cpl;
|
||||
int iSOMO;
|
||||
};
|
||||
|
||||
struct bin_tree {
|
||||
Node *rootNode;
|
||||
int NBF;
|
||||
};
|
||||
|
||||
#include "/usr/include/x86_64-linux-gnu/cblas.h"
|
||||
|
||||
#define MAX_SOMO 32
|
||||
|
||||
void buildTreeDriver(Tree *bftree, int NSOMO, int MS, int *NBF);
|
||||
|
||||
void buildTree(Tree *bftree, Node **inode, int isomo, int izeros, int icpl, int NSOMOMax, int MSmax);
|
||||
|
||||
void printTreeDriver(Tree *bftree, int NSOMOMax);
|
||||
void printTree(Node *bftree, int isomo, int NSOMOMax, int *vecBF);
|
||||
|
||||
void getIthBF(Node *node, int isomo, bool foundBF, int NSOMOMax, int getaddr, int *vecBF);
|
||||
void getIthBFDriver(Tree *bftree, int NSOMOMax, int getaddr, int *vecBF);
|
||||
|
||||
void getBFIndexList(int NSOMO, int *BF1, int *IdxListBF1);
|
||||
void getIslands(int NSOMO, int *BF1, int *BF2, int *nislands, int *phasefactor);
|
||||
|
||||
void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOMO);
|
||||
void getSetBits(int64_t n, int *nsetbits);
|
||||
void getOverlapMatrix(int64_t Isomo, int64_t MS, double **overlapMatrixptr, int *rows, int *cols, int *NSOMOout);
|
||||
void getOverlapMatrix_withDet(double *bftodetmatrixI, int rowsbftodetI, int colsbftodetI, int64_t Isomo, int64_t MS, double **overlapMatrixI, int *rowsI, int *colsI, int *NSOMO);
|
||||
void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix);
|
||||
|
||||
|
||||
void calculateMETypeSOMOSOMO(int *BF1, int *BF2, int moi, int moj, double *factor, int *phasefactor);
|
||||
void getOneElMETypeSOMOSOMO(int64_t Isomo, int64_t Jsomos, int moi, int moj, int MS, double **oneElMatrixElementsptr, int *rows, int *cols);
|
||||
|
||||
/***********************
|
||||
|
||||
Determinant Tree utils
|
||||
***********************/
|
||||
|
||||
|
||||
void genDets(Tree *dettree,
|
||||
Node **inode,
|
||||
int isomo,
|
||||
int izeros,
|
||||
int icpl,
|
||||
int NSOMOMax,
|
||||
int MSmax);
|
||||
void genDetsDriver(Tree *dettree, int NSOMO, int MS, int *Ndets);
|
||||
|
||||
void getIthDet(Node *inode, int isomo, bool foundBF, int NSOMOMax, int getaddr, int *vecBF);
|
||||
void getIthDetDriver(Tree *dettree, int NSOMOMax, int getaddr, int *vecBF);
|
||||
void getDetlistDriver(Tree *dettree, int NSOMOMax, int *detlist);
|
||||
void findAddofDet(Node *inode, int isomo, bool foundDet, int NSOMOMax, int *inpdet, int *addr);
|
||||
void findAddofDetDriver(Tree *dettree, int NSOMOMax, int *inpdet, int *addr);
|
||||
|
||||
|
||||
/************************/
|
||||
|
||||
void genDetBasis(Tree *dettree, int Isomo, int MS, int *ndets);
|
||||
void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowvec);
|
||||
void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *rows, int *cols);
|
||||
|
||||
// Misc utils
|
||||
void int_to_bin_digit(int64_t in, int count, int* out);
|
||||
void printRealMatrix(double *orthoMatrix, int rows, int cols);
|
||||
void callBlasMatxMat(double *A, int rowA, int colA, double *B, int rowB, int colB, double *C, bool transA, bool transB);
|
||||
void get_phase_cfg_to_qp_inpList(int *inpdet, int NSOMO, int *phaseout);
|
||||
void get_phase_cfg_to_qp_inpInt(int inpdet, double *phaseout);
|
||||
|
||||
#endif
|
495
src/davidson/davidson_parallel_csf.irp.f
Normal file
495
src/davidson/davidson_parallel_csf.irp.f
Normal file
@ -0,0 +1,495 @@
|
||||
use bitmasks
|
||||
use f77_zmq
|
||||
|
||||
|
||||
subroutine davidson_csf_slave_inproc(i)
|
||||
implicit none
|
||||
integer, intent(in) :: i
|
||||
|
||||
call davidson_csf_run_slave(1,i)
|
||||
end
|
||||
|
||||
|
||||
subroutine davidson_csf_slave_tcp(i)
|
||||
implicit none
|
||||
integer, intent(in) :: i
|
||||
call davidson_csf_run_slave(0,i)
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine davidson_csf_run_slave(thread,iproc)
|
||||
use f77_zmq
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Slave routine for Davidson's diagonalization.
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: thread, iproc
|
||||
|
||||
integer :: worker_id, task_id, blockb
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_push_socket
|
||||
integer(ZMQ_PTR) :: zmq_socket_push
|
||||
|
||||
integer, external :: connect_to_taskserver
|
||||
integer :: doexit, send, receive
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
|
||||
doexit = 0
|
||||
if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
|
||||
doexit=1
|
||||
endif
|
||||
IRP_IF MPI
|
||||
include 'mpif.h'
|
||||
integer :: ierr
|
||||
send = doexit
|
||||
call MPI_AllReduce(send, receive, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)
|
||||
if (ierr /= MPI_SUCCESS) then
|
||||
doexit=1
|
||||
endif
|
||||
doexit = receive
|
||||
IRP_ENDIF
|
||||
if (doexit>0) then
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||