PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-17 00:20:41 +02:00
Code Releases Activity

Merge pull request #219 from v1j4y/fix_csf_merge

Browse Source 
Fix csf merge
This commit is contained in:
Anthony Scemama 2022-11-03 15:51:13 +01:00 committed by GitHub
commit f82e9d358f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 121 additions and 119 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/csf/cfgCI_utils.c
View File

@ -255,23 +255,24 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols); 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){ // QR to orthogonalize CSFs does not work
int i,j; //void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
//for(j=0;j<cols;++j){ // int i,j;
// for(i=0;i<rows;++i){ // //for(j=0;j<cols;++j){
// printf(" %3.2f ",overlapMatrix[j*rows + i]); // // 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");
// //}
//} //}
// 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){ void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix){

14
src/csf/sigma_vector.irp.f
View File

@ -112,7 +112,7 @@
! endif ! endif
!enddo !enddo
n_CSF = 0 n_CSF = 0
print *," -9(((((((((((((( NSOMOMin=",NSOMOMin !print *," -9(((((((((((((( NSOMOMin=",NSOMOMin
ncfgprev = cfg_seniority_index(NSOMOMin) ! can be -1 ncfgprev = cfg_seniority_index(NSOMOMin) ! can be -1
if(ncfgprev.eq.-1)then if(ncfgprev.eq.-1)then
ncfgprev=1 ncfgprev=1
@ -145,10 +145,10 @@
endif endif
endif endif
n_CSF += ncfg*dimcsfpercfg n_CSF += ncfg*dimcsfpercfg
print *," i=",i," dimcsf=",dimcsfpercfg," ncfg=",ncfg, " ncfgprev=",ncfgprev, " senor=",cfg_seniority_index(i) !print *," i=",i," dimcsf=",dimcsfpercfg," ncfg=",ncfg, " ncfgprev=",ncfgprev, " senor=",cfg_seniority_index(i)
ncfgprev = cfg_seniority_index(i+2) ncfgprev = cfg_seniority_index(i+2)
end do end do
print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration !print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration
END_PROVIDER END_PROVIDER
@ -240,7 +240,7 @@ end subroutine get_phase_qp_to_cfg
! initialization ! initialization
psi_coef_config = 0.d0 psi_coef_config = 0.d0
DetToCSFTransformationMatrix(0,:,:) = 1.d0 DetToCSFTransformationMatrix(0,:,:) = 1.d0
do i = 2-iand(MS,1), NSOMOMax,2 do i = 2-iand(MS,1_8), NSOMOMax,2
Isomo = IBSET(0_8, i) - 1_8 Isomo = IBSET(0_8, i) - 1_8
! rows = Ncsfs ! rows = Ncsfs
! cols = Ndets ! cols = Ndets
@ -342,8 +342,8 @@ end subroutine get_phase_qp_to_cfg
nsomomin = elec_alpha_num-elec_beta_num nsomomin = elec_alpha_num-elec_beta_num
rowsmax = 0 rowsmax = 0
colsmax = 0 colsmax = 0
print *,"NSOMOMax = ",NSOMOMax !print *,"NSOMOMax = ",NSOMOMax
print *,"NSOMOMin = ",NSOMOMin !print *,"NSOMOMin = ",NSOMOMin
!allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
! Type ! Type
! 1. SOMO -> SOMO ! 1. SOMO -> SOMO
@ -525,7 +525,7 @@ end subroutine get_phase_qp_to_cfg
end do end do
end do end do
end do end do
print *,"Rowsmax=",rowsmax," Colsmax=",colsmax !print *,"Rowsmax=",rowsmax," Colsmax=",colsmax
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ real*8, AIJpqContainer, (NBFMax,NBFmax,NSOMOMax+1,NSOMOMax+1,4,NSOMOMin:NSOMOMax)] BEGIN_PROVIDER [ real*8, AIJpqContainer, (NBFMax,NBFmax,NSOMOMax+1,NSOMOMax+1,4,NSOMOMin:NSOMOMax)]

173
src/utils/linear_algebra.irp.f
View File

@ -1133,102 +1133,103 @@ subroutine ortho_svd(A,LDA,m,n)
end end
subroutine ortho_qr_withB(A,LDA,B,m,n) ! QR to orthonormalize CSFs does not work :-(
implicit none !subroutine ortho_qr_withB(A,LDA,B,m,n)
BEGIN_DOC ! implicit none
! Orthogonalization using Q.R factorization ! BEGIN_DOC
! ! Orthogonalization using Q.R factorization
! !
! ! A : Overlap Matrix
! !
! ! LDA : leftmost dimension of A
! !
! ! m : Number of rows of A
! !
! ! n : Number of columns of A
! !
! ! B : Output orthogonal basis
! !
! END_DOC
! integer, intent(in) :: m,n, LDA
! double precision, intent(inout) :: A(LDA,n)
! double precision, intent(inout) :: B(LDA,n)
! !
! A : Overlap Matrix ! integer :: LWORK, INFO
! integer, allocatable :: jpvt(:)
! double precision, allocatable :: TAU(:), WORK(:)
! double precision, allocatable :: C(:,:)
! double precision :: norm
! integer :: i,j
! !
! LDA : leftmost dimension of A ! allocate (TAU(min(m,n)), WORK(1))
! allocate (jpvt(n))
! !print *," In function ortho"
! B = A
! !
! m : Number of rows of A ! jpvt(1:n)=1
! !
! n : Number of columns of A
!
! B : Output orthogonal basis
!
END_DOC
integer, intent(in) :: m,n, LDA
double precision, intent(inout) :: A(LDA,n)
double precision, intent(inout) :: B(LDA,n)
integer :: LWORK, INFO
integer, allocatable :: jpvt(:)
double precision, allocatable :: TAU(:), WORK(:)
double precision, allocatable :: C(:,:)
double precision :: norm
integer :: i,j
allocate (TAU(min(m,n)), WORK(1))
allocate (jpvt(n))
!print *," In function ortho"
B = A
jpvt(1:n)=1
LWORK=-1
call dgeqp3( m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
LWORK=max(n,int(WORK(1)))
deallocate(WORK)
allocate(WORK(LWORK))
call dgeqp3(m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
print *,A
print *,jpvt
deallocate(WORK,TAU)
!stop
! LWORK=-1 ! LWORK=-1
!call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO ) ! call dgeqp3( m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
!
! ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 ! ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
! LWORK=max(n,int(WORK(1))) ! LWORK=max(n,int(WORK(1)))
!deallocate(WORK)
!allocate(WORK(LWORK))
!call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
!LWORK=-1
!call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
!! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
!LWORK=max(n,int(WORK(1)))
!deallocate(WORK)
!allocate(WORK(LWORK))
!call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
! !
!allocate(C(LDA,n)) ! deallocate(WORK)
!call dgemm('N','N',m,n,n,1.0d0,B,LDA,A,LDA,0.0d0,C,LDA) ! allocate(WORK(LWORK))
!norm = 0.0d0 ! call dgeqp3(m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
!B = 0.0d0 ! print *,A
!!print *,C ! print *,jpvt
!do i=1,m
! norm = 0.0d0
! do j=1,n
! norm = norm + C(j,i)*C(j,i)
! end do
! norm = 1.0d0/dsqrt(norm)
! do j=1,n
! B(j,i) = C(j,i)
! end do
!end do
!print *,B
! deallocate(WORK,TAU) ! deallocate(WORK,TAU)
end ! !stop
!
! !LWORK=-1
! !call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
! !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
! !LWORK=max(n,int(WORK(1)))
!
! !deallocate(WORK)
! !allocate(WORK(LWORK))
! !call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
!
! !LWORK=-1
! !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
! !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
! !LWORK=max(n,int(WORK(1)))
!
! !deallocate(WORK)
! !allocate(WORK(LWORK))
! !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
! !
! !allocate(C(LDA,n))
! !call dgemm('N','N',m,n,n,1.0d0,B,LDA,A,LDA,0.0d0,C,LDA)
! !norm = 0.0d0
! !B = 0.0d0
! !!print *,C
! !do i=1,m
! ! norm = 0.0d0
! ! do j=1,n
! ! norm = norm + C(j,i)*C(j,i)
! ! end do
! ! norm = 1.0d0/dsqrt(norm)
! ! do j=1,n
! ! B(j,i) = C(j,i)
! ! end do
! !end do
! !print *,B
!
!
! !deallocate(WORK,TAU)
!end
subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf") !subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf")
use iso_c_binding ! use iso_c_binding
integer(c_int32_t), value :: LDA ! integer(c_int32_t), value :: LDA
integer(c_int32_t), value :: m ! integer(c_int32_t), value :: m
integer(c_int32_t), value :: n ! integer(c_int32_t), value :: n
integer(c_int16_t) :: A(LDA,n) ! integer(c_int16_t) :: A(LDA,n)
integer(c_int16_t) :: B(LDA,n) ! integer(c_int16_t) :: B(LDA,n)
call ortho_qr_withB(A,LDA,B,m,n) ! call ortho_qr_withB(A,LDA,B,m,n)
end subroutine ortho_qr_csf !end subroutine ortho_qr_csf
subroutine ortho_qr(A,LDA,m,n) subroutine ortho_qr(A,LDA,m,n)
implicit none implicit none