LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-02 09:58:24 +01:00
Code Releases Activity

Merge branch 'dev' into guix

Browse Source
This commit is contained in:
Anthony Scemama 2022-11-01 17:04:31 +01:00
commit 6d2f444536
16 changed files with 3562 additions and 389 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/csf/cfgCI_interface.f90
View File

@ -46,6 +46,24 @@ module cfunctions
real (kind=C_DOUBLE ),intent(out) :: csftodetmatrix(rowsmax,colsmax)
end subroutine getCSFtoDETTransformationMatrix
end interface
interface
subroutine gramSchmidt(A, m, n, B) bind(C, name='gramSchmidt')
import C_INT32_T, C_INT64_T, C_DOUBLE
integer(kind=C_INT32_T),value,intent(in) :: m
integer(kind=C_INT32_T),value,intent(in) :: n
real (kind=C_DOUBLE ),intent(in) :: A(m,n)
real (kind=C_DOUBLE ),intent(out) :: B(m,n)
end subroutine gramSchmidt
end interface
interface
subroutine gramSchmidt_qp(A, m, n, B) bind(C, name='gramSchmidt_qp')
import C_INT32_T, C_INT64_T, C_DOUBLE
integer(kind=C_INT32_T),value,intent(in) :: m
integer(kind=C_INT32_T),value,intent(in) :: n
real (kind=C_DOUBLE ),intent(in) :: A(m,n)
real (kind=C_DOUBLE ),intent(out) :: B(m,n)
end subroutine gramSchmidt_qp
end interface
end module cfunctions
subroutine f_dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) &

110
src/csf/cfgCI_utils.c
View File

@ -1,5 +1,6 @@
#include <stdint.h>
#include <stdio.h>
#include <assert.h>
#include "tree_utils.h"
void int_to_bin_digit(int64_t in, int count, int* out)
@ -28,19 +29,19 @@ void getncsfs1(int *inpnsomo, int *inpms, int *outncsfs){
int nsomo = *inpnsomo;
int ms = *inpms;
int nparcoupl = (nsomo + ms)/2;
*outncsfs = binom(nsomo, nparcoupl);
*outncsfs = binom((double)nsomo, (double)nparcoupl);
}
void getncsfs(int NSOMO, int MS, int *outncsfs){
int nparcoupl = (NSOMO + MS)/2;
int nparcouplp1 = ((NSOMO + MS)/2)+1;
int nparcoupl = (NSOMO + MS)/2; // n_alpha
int nparcouplp1 = ((NSOMO + MS)/2)+1; // n_alpha + 1
double tmpndets=0.0;
if(NSOMO == 0){
(*outncsfs) = 1;
return;
}
tmpndets = binom(NSOMO, nparcoupl);
(*outncsfs) = round(tmpndets - binom(NSOMO, nparcouplp1));
tmpndets = binom((double)NSOMO, (double)nparcoupl);
(*outncsfs) = round(tmpndets - binom((double)NSOMO, (double)nparcouplp1));
}
#include <stdint.h>
@ -252,6 +253,26 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
buildTreeDriver(bftree, *NSOMO, MS, NBF);
}
void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols);
void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
int i,j;
//for(j=0;j<cols;++j){
// for(i=0;i<rows;++i){
// printf(" %3.2f ",overlapMatrix[j*rows + i]);
// }
// printf("\n");
//}
// Call the function ortho_qr from qp
ortho_qr_csf(overlapMatrix, rows, orthoMatrix, rows, cols);
//for(j=0;j<cols;++j){
// for(i=0;i<rows;++i){
// printf(" %3.2f ",orthoMatrix[j*rows + i]);
// }
// printf("\n");
//}
}
void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
// vector
@ -341,8 +362,12 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
Get BFtoDeterminant Matrix
************************************/
printf("In convertcsftodet\n");
//printf(" --- In convet ----\n");
convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
//printf(" --- done bf det basis ---- row=%d col=%d\n",rowsbftodetI,colsbftodetI);
//printRealMatrix(bftodetmatrixI,rowsbftodetI,colsbftodetI);
int rowsI = 0;
int colsI = 0;
@ -350,6 +375,8 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
//getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
//printf("Overlap matrix\n");
//printRealMatrix(overlapMatrixI,rowsI,colsI);
/***********************************
Get Orthonormalization Matrix
@ -359,6 +386,9 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
gramSchmidt(overlapMatrixI, rowsI, colsI, orthoMatrixI);
//printf("Ortho matrix\n");
//printRealMatrix(orthoMatrixI,rowsI,colsI);
/***********************************
Get Final CSF to Det Matrix
************************************/
@ -1340,11 +1370,11 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
for(int i = 0; i < npairs; i++){
for(int j = 0; j < NSOMO; j++) {
inpdet[j] = detslist[i*NSOMO + j];
printf(" %d ",inpdet[j]);
//printf(" %d ",inpdet[j]);
}
printf("\n");
//printf("\n");
findAddofDetDriver(dettree, NSOMO, inpdet, &addr);
printf("(%d) - addr = %d\n",i,addr);
//printf("(%d) - addr = %d\n",i,addr);
// Calculate the phase for cfg to QP2 conversion
//get_phase_cfg_to_qp_inpList(inpdet, NSOMO, &phase_cfg_to_qp);
//rowvec[addr] = 1.0 * phaselist[i]*phase_cfg_to_qp/sqrt(fac);
@ -1363,12 +1393,23 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *rows, int *cols){
int NSOMO=0;
//printf("before getSetBits Isomo=%ld, NSOMO=%ld\n",Isomo,NSOMO);
getSetBits(Isomo, &NSOMO);
//printf("Isomo=%ld, NSOMO=%ld\n",Isomo,NSOMO);
int ndets = 0;
int NBF = 0;
double dNSOMO = NSOMO*1.0;
double nalpha = (NSOMO + MS)/2.0;
ndets = (int)binom(dNSOMO, nalpha);
//double dNSOMO = NSOMO*1.0;
// MS = alpha_num - beta_num
int nalpha = (NSOMO + MS)/2;
//printf(" in convertbftodet : MS=%d nalpha=%3.2f\n",MS,nalpha);
//ndets = (int)binom(dNSOMO, nalpha);
if(NSOMO > 0){
ndets = (int)binom((double)NSOMO, (double)nalpha);
}
else if(NSOMO == 0){
ndets = 1;
}
else printf("Something is wrong in calcMEdetpair\n");
Tree dettree = (Tree){ .rootNode = NULL, .NBF = -1 };
dettree.rootNode = malloc(sizeof(Node));
@ -1389,16 +1430,6 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
}
else{
//int addr = -1;
//int inpdet[NSOMO];
//inpdet[0] = 1;
//inpdet[1] = 1;
//inpdet[2] = 1;
//inpdet[3] = 0;
//inpdet[4] = 0;
//inpdet[5] = 0;
//findAddofDetDriver(&dettree, NSOMO, inpdet, &addr);
int detlist[ndets];
getDetlistDriver(&dettree, NSOMO, detlist);
@ -1411,6 +1442,9 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
generateAllBFs(Isomo, MS, &bftree, &NBF, &NSOMO);
// Initialize transformation matrix
//printf("MS=%d NBF=%d ndets=%d NSOMO=%d\n",MS,NBF,ndets,NSOMO);
assert( NBF > 0);
assert( ndets > 0);
(*bftodetmatrixptr) = malloc(NBF*ndets*sizeof(double));
(*rows) = NBF;
(*cols) = ndets;
@ -1465,9 +1499,10 @@ void convertBFtoDetBasisWithArrayDims(int64_t Isomo, int MS, int rowsmax, int co
getSetBits(Isomo, &NSOMO);
int ndets = 0;
int NBF = 0;
double dNSOMO = NSOMO*1.0;
double nalpha = (NSOMO + MS)/2.0;
ndets = (int)binom(dNSOMO, nalpha);
//double dNSOMO = NSOMO*1.0;
//double nalpha = (NSOMO + MS)/2.0;
int nalpha = (NSOMO + MS)/2;
ndets = (int)binom((double)NSOMO, (double)nalpha);
Tree dettree = (Tree){ .rootNode = NULL, .NBF = -1 };
dettree.rootNode = malloc(sizeof(Node));
@ -1551,6 +1586,7 @@ void getApqIJMatrixDims(int64_t Isomo, int64_t Jsomo, int64_t MS, int32_t *rowso
getncsfs(NSOMOJ, MS, &NBFJ);
(*rowsout) = NBFI;
(*colsout) = NBFJ;
//exit(0);
}
void getApqIJMatrixDriver(int64_t Isomo, int64_t Jsomo, int orbp, int orbq, int64_t MS, int64_t NMO, double **CSFICSFJApqIJptr, int *rowsout, int *colsout){
@ -1669,6 +1705,7 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
int rowsbftodetI, colsbftodetI;
//printf(" 1Calling convertBFtoDetBasis Isomo=%ld MS=%ld\n",Isomo,MS);
convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
// Fill matrix
@ -1676,8 +1713,14 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
int colsI = 0;
//getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
//getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
//printf("Isomo=%ld\n",Isomo);
getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
if(Isomo == 0){
rowsI = 1;
colsI = 1;
}
//printf("Isomo=%ld\n",Isomo);
orthoMatrixI = malloc(rowsI*colsI*sizeof(double));
@ -1689,6 +1732,7 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
int rowsbftodetJ, colsbftodetJ;
//printf(" 2Calling convertBFtoDetBasis Jsomo=%ld MS=%ld\n",Jsomo,MS);
convertBFtoDetBasis(Jsomo, MS, &bftodetmatrixJ, &rowsbftodetJ, &colsbftodetJ);
int rowsJ = 0;
@ -1696,6 +1740,10 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
// Fill matrix
//getOverlapMatrix(Jsomo, MS, &overlapMatrixJ, &rowsJ, &colsJ, &NSOMO);
getOverlapMatrix_withDet(bftodetmatrixJ, rowsbftodetJ, colsbftodetJ, Jsomo, MS, &overlapMatrixJ, &rowsJ, &colsJ, &NSOMO);
if(Jsomo == 0){
rowsJ = 1;
colsJ = 1;
}
orthoMatrixJ = malloc(rowsJ*colsJ*sizeof(double));
@ -1713,18 +1761,25 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
int transA=false;
int transB=false;
//printf("1Calling blas\n");
//printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsbftodetI,colsbftodetI,rowsA,colsA);
callBlasMatxMat(bftodetmatrixI, rowsbftodetI, colsbftodetI, ApqIJ, rowsA, colsA, bfIApqIJ, transA, transB);
//printf("done\n");
// now transform I in csf basis
double *CSFIApqIJ = malloc(rowsI*colsA*sizeof(double));
transA = false;
transB = false;
//printf("2Calling blas\n");
//printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,colsI,colsI,colsA);
callBlasMatxMat(orthoMatrixI, rowsI, colsI, bfIApqIJ, colsI, colsA, CSFIApqIJ, transA, transB);
// now transform J in BF basis
double *CSFIbfJApqIJ = malloc(rowsI*rowsbftodetJ*sizeof(double));
transA = false;
transB = true;
//printf("3Calling blas\n");
//printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,colsA,rowsbftodetJ,colsbftodetJ);
callBlasMatxMat(CSFIApqIJ, rowsI, colsA, bftodetmatrixJ, rowsbftodetJ, colsbftodetJ, CSFIbfJApqIJ, transA, transB);
// now transform J in CSF basis
@ -1735,13 +1790,14 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
double *tmpCSFICSFJApqIJ = malloc(rowsI*rowsJ*sizeof(double));
transA = false;
transB = true;
//printf("4Calling blas\n");
//printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,rowsbftodetJ,rowsJ,colsJ);
callBlasMatxMat(CSFIbfJApqIJ, rowsI, rowsbftodetJ, orthoMatrixJ, rowsJ, colsJ, tmpCSFICSFJApqIJ, transA, transB);
// Transfer to actual buffer in Fortran order
for(int i = 0; i < rowsI; i++)
for(int j = 0; j < rowsJ; j++)
CSFICSFJApqIJ[j*rowsI + i] = tmpCSFICSFJApqIJ[i*rowsJ + j];
// Garbage collection
free(overlapMatrixI);
free(overlapMatrixJ);

684
src/csf/configuration_CI_sigma_helpers.irp.f
View File

@ -1,3 +1,592 @@
use bitmasks
BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*(mo_num))]
&BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ]
implicit none
!use bitmasks
BEGIN_DOC
! Documentation for alphasI
! Returns the associated alpha's for
! the input configuration Icfg.
END_DOC
integer :: idxI ! The id of the Ith CFG
integer(bit_kind) :: Icfg(N_int,2)
integer :: NalphaIcfg
logical,dimension(:,:),allocatable :: tableUniqueAlphas
integer :: listholes(mo_num)
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
integer :: nholes
integer :: nvmos
integer :: listvmos(mo_num)
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
integer*8 :: Idomo
integer*8 :: Isomo
integer*8 :: Jdomo
integer*8 :: Jsomo
integer*8 :: diffSOMO
integer*8 :: diffDOMO
integer*8 :: xordiffSOMODOMO
integer :: ndiffSOMO
integer :: ndiffDOMO
integer :: nxordiffSOMODOMO
integer :: ndiffAll
integer :: i,ii
integer :: j,jj
integer :: k,kk
integer :: kstart
integer :: kend
integer :: Nsomo_I
integer :: hole
integer :: p
integer :: q
integer :: countalphas
logical :: pqAlreadyGenQ
logical :: pqExistsQ
logical :: ppExistsQ
integer*8 :: MS
double precision :: t0, t1
call wall_time(t0)
MS = elec_alpha_num-elec_beta_num
allocate(tableUniqueAlphas(mo_num,mo_num))
NalphaIcfg_list = 0
do idxI = 1, N_configuration
Icfg = psi_configuration(:,:,idxI)
Isomo = iand(act_bitmask(1,1),Icfg(1,1))
Idomo = iand(act_bitmask(1,1),Icfg(1,2))
! find out all pq holes possible
nholes = 0
! holes in SOMO
do ii = 1,n_act_orb
i = list_act(ii)
if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = i
holetype(nholes) = 1
endif
end do
! holes in DOMO
do ii = 1,n_act_orb
i = list_act(ii)
if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = i
holetype(nholes) = 2
endif
end do
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
do ii = 1,n_act_orb
i = list_act(ii)
if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then
if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then
nvmos += 1
listvmos(nvmos) = i
vmotype(nvmos) = 1
else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then
nvmos += 1
listvmos(nvmos) = i
vmotype(nvmos) = 2
end if
end if
end do
tableUniqueAlphas = .FALSE.
! Now find the allowed (p,q) excitations
Isomo = iand(act_bitmask(1,1),Icfg(1,1))
Idomo = iand(act_bitmask(1,1),Icfg(1,2))
Nsomo_I = POPCNT(Isomo)
if(Nsomo_I .EQ. 0) then
kstart = 1
else
kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))
endif
kend = idxI-1
do i = 1,nholes
p = listholes(i)
do j = 1,nvmos
q = listvmos(j)
if(p .EQ. q) cycle
if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
! SOMO -> VMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = Idomo
kstart = max(1,cfg_seniority_index(max(NSOMOMin,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)
! Check for Minimal alpha electrons (MS)
if(POPCNT(Jsomo).ge.MS)then
kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
kend = idxI-1
else
cycle
endif
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
! DOMO -> VMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
kstart = cfg_seniority_index(Nsomo_I)
kend = idxI-1
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
! DOMO -> SOMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
Jdomo = IBSET(Jdomo,q-1)
kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
kend = idxI-1
else
print*,"Something went wrong in obtain_associated_alphaI"
endif
! Check for Minimal alpha electrons (MS)
if(POPCNT(Jsomo).lt.MS)then
cycle
endif
! Again, we don't have to search from 1
! we just use seniority to find the
! first index with NSOMO - 2 to NSOMO + 2
! this is what is done in kstart, kend
pqAlreadyGenQ = .FALSE.
! First check if it can be generated before
do k = kstart, kend
diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)))
ndiffSOMO = POPCNT(diffSOMO)
if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle
diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
!if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
pqAlreadyGenQ = .TRUE.
ppExistsQ = .TRUE.
EXIT
endif
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
pqAlreadyGenQ = .TRUE.
EXIT
endif
end do
if(pqAlreadyGenQ) cycle
pqExistsQ = .FALSE.
if(.NOT. pqExistsQ) then
tableUniqueAlphas(p,q) = .TRUE.
endif
end do
end do
!print *,tableUniqueAlphas(:,:)
! prune list of alphas
Isomo = Icfg(1,1)
Idomo = Icfg(1,2)
Jsomo = Icfg(1,1)
Jdomo = Icfg(1,2)
NalphaIcfg = 0
do i = 1, nholes
p = listholes(i)
do j = 1, nvmos
q = listvmos(j)
if(p .EQ. q) cycle
if(tableUniqueAlphas(p,q)) then
if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
! SOMO -> VMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = Idomo
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)
if(POPCNT(Jsomo).ge.MS)then
kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
kend = idxI-1
else
cycle
endif
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
! DOMO -> VMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
! DOMO -> SOMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
Jdomo = IBSET(Jdomo,q-1)
else
print*,"Something went wrong in obtain_associated_alphaI"
endif
! SOMO
!print *,i,j,"|",NalphaIcfg, Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
if(POPCNT(Jsomo) .ge. NSOMOMin) then
NalphaIcfg += 1
alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo
alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
NalphaIcfg_list(idxI) = NalphaIcfg
endif
endif
end do
end do
! Check if this Icfg has been previously generated as a mono
ppExistsQ = .False.
Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
do k = kstart, idxI-1
diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k)))
ndiffSOMO = POPCNT(diffSOMO)
if (ndiffSOMO /= 2) cycle
diffDOMO = IEOR(Idomo,iand(act_bitmask(1,1),psi_configuration(1,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then
ppExistsQ = .TRUE.
EXIT
endif
end do
! Diagonal part (pp,qq)
if(nholes > 0 .AND. (.NOT. ppExistsQ))then
! SOMO
if(POPCNT(Jsomo) .ge. NSOMOMin) then
NalphaIcfg += 1
alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1)
alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2)
NalphaIcfg_list(idxI) = NalphaIcfg
endif
endif
NalphaIcfg = 0
enddo ! end loop idxI
call wall_time(t1)
print *, 'Preparation : ', t1 - t0
END_PROVIDER
subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg)
implicit none
use bitmasks
BEGIN_DOC
! Documentation for alphasI
! Returns the associated alpha's for
! the input configuration Icfg.
END_DOC
integer,intent(in) :: idxI ! The id of the Ith CFG
integer(bit_kind),intent(in) :: Icfg(N_int,2)
integer,intent(out) :: NalphaIcfg
integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*)
logical,dimension(:,:),allocatable :: tableUniqueAlphas
integer :: listholes(mo_num)
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
integer :: nholes
integer :: nvmos
integer :: listvmos(mo_num)
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
integer*8 :: Idomo
integer*8 :: Isomo
integer*8 :: Jdomo
integer*8 :: Jsomo
integer*8 :: diffSOMO
integer*8 :: diffDOMO
integer*8 :: xordiffSOMODOMO
integer :: ndiffSOMO
integer :: ndiffDOMO
integer :: nxordiffSOMODOMO
integer :: ndiffAll
integer :: i, ii
integer :: j, jj
integer :: k, kk
integer :: kstart
integer :: kend
integer :: Nsomo_I
integer :: hole
integer :: p
integer :: q
integer :: countalphas
logical :: pqAlreadyGenQ
logical :: pqExistsQ
logical :: ppExistsQ
Isomo = iand(act_bitmask(1,1),Icfg(1,1))
Idomo = iand(act_bitmask(1,1),Icfg(1,2))
!print*,"Input cfg"
!call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1)
! find out all pq holes possible
nholes = 0
! holes in SOMO
do ii = 1,n_act_orb
i = list_act(ii)
if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = i
holetype(nholes) = 1
endif
end do
! holes in DOMO
do ii = 1,n_act_orb
i = list_act(ii)
if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = i
holetype(nholes) = 2
endif
end do
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
do ii = 1,n_act_orb
i = list_act(ii)
!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,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then
nvmos += 1
listvmos(nvmos) = i
vmotype(nvmos) = 1
else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then
nvmos += 1
listvmos(nvmos) = i
vmotype(nvmos) = 2
end if
end do
!print *,"Nvmo=",nvmos
!print *,listvmos
!print *,vmotype
allocate(tableUniqueAlphas(mo_num,mo_num))
tableUniqueAlphas = .FALSE.
! Now find the allowed (p,q) excitations
Isomo = iand(act_bitmask(1,1),Icfg(1,1))
Idomo = iand(act_bitmask(1,1),Icfg(1,2))
Nsomo_I = POPCNT(Isomo)
if(Nsomo_I .EQ. 0) then
kstart = 1
else
kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))
endif
kend = idxI-1
!print *,"Isomo"
!call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1)
!print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI
!do i = 1,nholes
! print *,i,"->",listholes(i)
!enddo
!do i = 1,nvmos
! print *,i,"->",listvmos(i)
!enddo
do i = 1,nholes
p = listholes(i)
do j = 1,nvmos
q = listvmos(j)
if(p .EQ. q) cycle
if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
! SOMO -> VMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = Idomo
kstart = max(1,cfg_seniority_index(max(NSOMOMin,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(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
kend = idxI-1
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
! DOMO -> VMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
kstart = cfg_seniority_index(Nsomo_I)
kend = idxI-1
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
! DOMO -> SOMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
Jdomo = IBSET(Jdomo,q-1)
kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
kend = idxI-1
else
print*,"Something went wrong in obtain_associated_alphaI"
endif
! Again, we don't have to search from 1
! we just use seniortiy to find the
! first index with NSOMO - 2 to NSOMO + 2
! this is what is done in kstart, kend
pqAlreadyGenQ = .FALSE.
! First check if it can be generated before
do k = kstart, kend
diffSOMO = IEOR(Jsomo,iand(act_bitmask(1,1),psi_configuration(1,1,k)))
ndiffSOMO = POPCNT(diffSOMO)
if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle
diffDOMO = IEOR(Jdomo,iand(act_bitmask(1,1),psi_configuration(1,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
!if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
pqAlreadyGenQ = .TRUE.
ppExistsQ = .TRUE.
EXIT
endif
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
pqAlreadyGenQ = .TRUE.
!EXIT
!ppExistsQ = .TRUE.
!print *,i,k,ndiffSOMO,ndiffDOMO
!call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1)
!call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1)
!call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1)
EXIT
endif
end do
!print *,"(,",p,",",q,")",pqAlreadyGenQ
if(pqAlreadyGenQ) cycle
pqExistsQ = .FALSE.
! now check if this exists in the selected list
!do k = idxI+1, N_configuration
! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k))
! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k))
! ndiffSOMO = POPCNT(diffSOMO)
! ndiffDOMO = POPCNT(diffDOMO)
! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
! pqExistsQ = .TRUE.
! EXIT
! endif
!end do
if(.NOT. pqExistsQ) then
tableUniqueAlphas(p,q) = .TRUE.
!print *,p,q
!call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1)
endif
end do
end do
!print *,tableUniqueAlphas(:,:)
! prune list of alphas
Isomo = Icfg(1,1)
Idomo = Icfg(1,2)
Jsomo = Icfg(1,1)
Jdomo = Icfg(1,2)
NalphaIcfg = 0
do i = 1, nholes
p = listholes(i)
do j = 1, nvmos
q = listvmos(j)
if(p .EQ. q) cycle
if(tableUniqueAlphas(p,q)) then
if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
! SOMO -> VMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = Idomo
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)
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
! DOMO -> VMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
! DOMO -> SOMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
Jdomo = IBSET(Jdomo,q-1)
else
print*,"Something went wrong in obtain_associated_alphaI"
endif
! SOMO
NalphaIcfg += 1
!print *,i,j,"|",NalphaIcfg
alphasIcfg(1,1,NalphaIcfg) = Jsomo
alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
!print *,"I = ",idxI, " Na=",NalphaIcfg," - ",Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
endif
end do
end do
! Check if this Icfg has been previously generated as a mono
ppExistsQ = .False.
Isomo = iand(act_bitmask(1,1),Icfg(1,1))
Idomo = iand(act_bitmask(1,1),Icfg(1,2))
do k = 1, idxI-1
diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k)))
diffDOMO = IEOR(Idomo,iand(act_bitmask(1,1),psi_configuration(1,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffSOMO = POPCNT(diffSOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
ppExistsQ = .TRUE.
EXIT
endif
end do
! Diagonal part (pp,qq)
if(nholes > 0 .AND. (.NOT. ppExistsQ))then
! SOMO
NalphaIcfg += 1
!print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg
!call debug_spindet(Idomo,1)
!call debug_spindet(Jdomo,1)
alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1)
alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2)
endif
end subroutine
function getNSOMO(Icfg) result(NSOMO)
implicit none
integer(bit_kind),intent(in) :: Icfg(N_int,2)
@ -8,98 +597,3 @@
NSOMO += POPCNT(Icfg(i,1))
enddo
end function getNSOMO
subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmodel)
implicit none
BEGIN_DOC
! This function converts the orbital ids
! in real space to those used in model space
! in order to identify the matrices required
! for the calculation of MEs.
!
! 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(bit_kind),intent(in) :: Ialpha(N_int,2)
integer(bit_kind),intent(in) :: Jcfg(N_int,2)
integer,intent(in) :: p,q
integer,intent(in) :: extype
integer,intent(out) :: pmodel,qmodel
integer*8 :: Isomo
integer*8 :: Idomo
integer*8 :: Jsomo
integer*8 :: Jdomo
integer*8 :: mask
integer*8 :: Isomotmp
integer*8 :: Jsomotmp
integer :: pos0,pos0prev
! TODO Flag (print) when model space indices is > 64
Isomo = Ialpha(1,1)
Idomo = Ialpha(1,2)
Jsomo = Jcfg(1,1)
Jdomo = Jcfg(1,2)
pos0prev = 0
pmodel = p
qmodel = q
if(p .EQ. q) then
pmodel = 1
qmodel = 1
else
!print *,"input pq=",p,q,"extype=",extype
!call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1)
!call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1)
select case(extype)
case (1)
! SOMO -> SOMO
! remove all domos
!print *,"type -> SOMO -> SOMO"
mask = ISHFT(1_8,p) - 1
Isomotmp = IAND(Isomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
mask = ISHFT(1_8,q) - 1
Isomotmp = IAND(Isomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
case (2)
! DOMO -> VMO
! remove all domos except one at p
!print *,"type -> DOMO -> VMO"
mask = ISHFT(1_8,p) - 1
Jsomotmp = IAND(Jsomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
mask = ISHFT(1_8,q) - 1
Jsomotmp = IAND(Jsomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
case (3)
! SOMO -> VMO
!print *,"type -> SOMO -> VMO"
!Isomo = IEOR(Isomo,Jsomo)
mask = ISHFT(1_8,p) - 1
Isomo = IAND(Isomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
mask = ISHFT(1_8,q) - 1
Jsomo = IAND(Jsomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
case (4)
! DOMO -> SOMO
! remove all domos except one at p
!print *,"type -> DOMO -> SOMO"
!Isomo = IEOR(Isomo,Jsomo)
mask = ISHFT(1_8,p) - 1
Jsomo = IAND(Jsomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
mask = ISHFT(1_8,q) - 1
Isomo = IAND(Isomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
case default
print *,"something is wrong in convertOrbIdsToModelSpaceIds"
end select
endif
!print *,p,q,"model ids=",pmodel,qmodel
end subroutine convertOrbIdsToModelSpaceIds

121
src/csf/configurations.irp.f
View File

@ -458,8 +458,9 @@ end
END_PROVIDER
BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num) ]
BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num+2) ]
&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
@ -467,9 +468,10 @@ END_PROVIDER
!
! cfg_nsomo_max : Max 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 +484,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
@ -743,41 +749,112 @@ BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_d
enddo
END_PROVIDER
subroutine binary_search_cfg(cfgInp,addcfg)
subroutine binary_search_cfg(cfgInp,addcfg,bit_tmp)
use bitmasks
implicit none
BEGIN_DOC
! Documentation for binary_search
!
! Does a binary search to find
!
! Does a binary search to find
! the address of a configuration in a list of
! configurations.
END_DOC
integer(bit_kind), intent(in) :: cfgInp(N_int,2)
integer , intent(out) :: addcfg
integer :: i,j,k,r,l
integer*8 :: key, key2
logical :: found
!integer*8, allocatable :: bit_tmp(:)
!integer*8, external :: configuration_search_key
integer*8, intent(in) :: bit_tmp(0:N_configuration+1)
!allocate(bit_tmp(0:N_configuration))
!bit_tmp(0) = 0
do i=1,N_configuration
!bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int)
found = .True.
do k=1,N_int
found = found .and. (psi_configuration(k,1,i) == cfgInp(k,1)) &
.and. (psi_configuration(k,2,i) == cfgInp(k,2))
enddo
if (found) then
addcfg = i
exit
logical :: found
integer :: l, r, j, k
integer*8 :: key
integer*8, external :: configuration_search_key
key = configuration_search_key(cfgInp,N_int)
! Binary search
l = 0
r = N_configuration+1
IRP_IF WITHOUT_SHIFTRL
j = ishft(r-l,-1)
IRP_ELSE
j = shiftr(r-l,1)
IRP_ENDIF
do while (j>=1)
j = j+l
if (bit_tmp(j) == key) then
! Find 1st element which matches the key
if (j > 1) then
do while (j>1 .and. bit_tmp(j-1) == key)
j = j-1
enddo
endif
! Find correct element matching the key
do while (bit_tmp(j) == key)
found = .True.
do k=1,N_int
found = found .and. (psi_configuration(k,1,j) == cfgInp(k,1))&
.and. (psi_configuration(k,2,j) == cfgInp(k,2))
enddo
if (found) then
addcfg = j
return
endif
j = j+1
enddo
addcfg = -1
return
else if (bit_tmp(j) > key) then
r = j
else
l = j
endif
IRP_IF WITHOUT_SHIFTRL
j = ishft(r-l,-1)
IRP_ELSE
j = shiftr(r-l,1)
IRP_ENDIF
enddo
addcfg = -1
return
end subroutine
!subroutine binary_search_cfg(cfgInp,addcfg)
! use bitmasks
! implicit none
! BEGIN_DOC
! ! Documentation for binary_search
! !
! ! Does a binary search to find
! ! the address of a configuration in a list of
! ! configurations.
! END_DOC
! integer(bit_kind), intent(in) :: cfgInp(N_int,2)
! integer , intent(out) :: addcfg
! integer :: i,j,k,r,l
! integer*8 :: key, key2
! logical :: found
! !integer*8, allocatable :: bit_tmp(:)
! !integer*8, external :: configuration_search_key
!
! !allocate(bit_tmp(0:N_configuration))
! !bit_tmp(0) = 0
! do i=1,N_configuration
! !bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int)
! found = .True.
! do k=1,N_int
! found = found .and. (psi_configuration(k,1,i) == cfgInp(k,1)) &
! .and. (psi_configuration(k,2,i) == cfgInp(k,2))
! enddo
! if (found) then
! addcfg = i
! exit
! endif
! enddo
!
!end subroutine
!
BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ]
&BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]

15
src/csf/conversion.irp.f
View File

@ -24,7 +24,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
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 :: rows, cols, i, j, k, salpha
integer :: startdet, enddet
integer :: ndetI
integer :: getNSOMO
@ -65,9 +65,11 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
enddo
if(iand(s,1) .EQ. 0) then
bfIcfg = max(1,nint((binom(s,s/2)-binom(s,(s/2)+1))))
salpha = (s + MS)/2
bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
else
bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))
salpha = (s + MS)/2
bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
endif
! perhaps blocking with CFGs of same seniority
@ -99,7 +101,7 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
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 :: s, bfIcfg, salpha
integer :: countcsf
integer(bit_kind) :: Ialpha(N_int), Ibeta(N_int)
integer :: rows, cols, i, j, k
@ -110,6 +112,8 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
double precision,allocatable :: tempCoeff (:,:)
double precision :: phasedet
integer :: idx
integer MS
MS = elec_alpha_num-elec_beta_num
countcsf = 0
@ -123,7 +127,8 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
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))))
salpha = (s + MS)/2
bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
allocate(tempCoeff(bfIcfg,N_st))

15
src/csf/create_excitations.irp.f
View File

@ -226,7 +226,7 @@ subroutine generate_all_singles_cfg(cfg,singles,n_singles,Nint)
enddo
end
subroutine generate_all_singles_cfg_with_type(cfgInp,singles,idxs_singles,pq_singles,ex_type_singles,n_singles,Nint)
subroutine generate_all_singles_cfg_with_type(bit_tmp,cfgInp,singles,idxs_singles,pq_singles,ex_type_singles,n_singles,Nint)
implicit none
use bitmasks
BEGIN_DOC
@ -238,6 +238,7 @@ subroutine generate_all_singles_cfg_with_type(cfgInp,singles,idxs_singles,pq_sin
! ex_type_singles : on output contains type of excitations :
!
END_DOC
integer*8, intent(in) :: bit_tmp(0:N_configuration+1)
integer, intent(in) :: Nint
integer, intent(inout) :: n_singles
integer, intent(out) :: idxs_singles(*)
@ -248,20 +249,26 @@ subroutine generate_all_singles_cfg_with_type(cfgInp,singles,idxs_singles,pq_sin
integer(bit_kind) :: Jdet(Nint,2)
integer :: i,k, n_singles_ma, i_hole, i_particle, ex_type, addcfg
integer :: ii,kk
integer(bit_kind) :: single(Nint,2)
logical :: i_ok
n_singles = 0
!TODO
!Make list of Somo and Domo for holes
!Make list of Unocc and Somo for particles
do i_hole = 1+n_core_orb, n_core_orb + n_act_orb
do i_particle = 1+n_core_orb, n_core_orb + n_act_orb
!do i_hole = 1+n_core_orb, n_core_orb + n_act_orb
do ii = 1, n_act_orb
i_hole = list_act(ii)