mirror of
https://github.com/pfloos/quack
synced 2024-12-22 04:14:26 +01:00
fixed conf in pCCD
This commit is contained in:
Abdallah Ammar
commit
20875a51a4
3
.gitignore
vendored
3
.gitignore
vendored
@ -1,2 +1,5 @@
|
||||
*.o
|
||||
*.
|
||||
__pycache__
|
||||
|
||||
.ninja_deps
|
||||
|
||||
@ -2,7 +2,7 @@
|
||||
! ---
|
||||
|
||||
subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, &
|
||||
maxSCF, thresh, max_diis, nBas, nOrb, nC, nO, nV, nR, Hc, ERI, ENuc, ERHF, eHF, cHF)
|
||||
maxSCF, thresh, max_diis, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF)
|
||||
|
||||
! Coupled-cluster module
|
||||
|
||||
@ -37,7 +37,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
double precision,intent(in) :: eHF(nOrb)
|
||||
double precision,intent(in) :: cHF(nBas,nOrb)
|
||||
double precision,intent(in) :: Hc(nBas,nBas)
|
||||
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
|
||||
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
|
||||
double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
|
||||
|
||||
! Local variables
|
||||
|
||||
@ -50,7 +51,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
if(doCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call CCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call CCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -66,7 +67,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
if(doDCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call DCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call DCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR, &
|
||||
ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -84,7 +86,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
if(doCCSD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -100,7 +102,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
if(dodrCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call drCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call drCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -116,7 +118,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
if(dorCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call rCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call rCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -132,7 +134,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
if(docrCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call crCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call crCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -148,7 +150,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
if(dolCCD) then
|
||||
|
||||
call wall_time(start_CC)
|
||||
call lCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
call lCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
@ -165,7 +167,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
|
||||
|
||||
call wall_time(start_CC)
|
||||
call pCCD(dotest, maxSCF, thresh, max_diis, nBas, nOrb, &
|
||||
nC, nO, nV, nR, Hc, ERI, ENuc, ERHF, eHF, cHF)
|
||||
nC, nO, nV, nR, Hc, ERI_AO, ENuc, ERHF, eHF, cHF)
|
||||
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
|
||||
1084
src/CC/pCCD.f90
1084
src/CC/pCCD.f90
File diff suppressed because it is too large
Load Diff
@ -135,131 +135,132 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
|
||||
dim_1 = (nBas - nO) * (nBas - nO - 1) / 2
|
||||
dim_2 = nO * (nO - 1) / 2
|
||||
|
||||
allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
|
||||
ERI_1 = 0.d0
|
||||
ERI_2 = 0.d0
|
||||
if((dim_1 .eq. 0) .or. (dim_2 .eq. 0)) then
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, c, d, cd, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = c+1, nBas-nR
|
||||
cd = cd + 1
|
||||
ERI_1(p,q,cd) = ERI(p,q,c,d) - ERI(p,q,d,c)
|
||||
enddo
|
||||
enddo
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = k+1, nO
|
||||
kl = kl + 1
|
||||
ERI_2(p,q,kl) = ERI(p,q,k,l) - ERI(p,q,l,k)
|
||||
end do
|
||||
end do
|
||||
enddo
|
||||
enddo
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
|
||||
ab = 0
|
||||
|
||||
do a = nO+1, nBas-nR
|
||||
do b = a+1, nBas-nR
|
||||
|
||||
ab = ab + 1
|
||||
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = c+1, nBas-nR
|
||||
|
||||
cd = cd + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab)
|
||||
end do ! d
|
||||
end do ! c
|
||||
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = k+1, nO
|
||||
|
||||
kl = kl + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
+ (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab)
|
||||
end do ! l
|
||||
end do ! k
|
||||
end do ! b
|
||||
end do ! a
|
||||
|
||||
ij = 0
|
||||
do i = nC+1, nO
|
||||
do j = i+1, nO
|
||||
|
||||
ij = ij + 1
|
||||
|
||||
cd = 0
|
||||
|
||||
do c = nO+1, nBas-nR
|
||||
do d = c+1, nBas-nR
|
||||
|
||||
cd = cd + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij)
|
||||
end do ! d
|
||||
end do ! c
|
||||
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = k+1, nO
|
||||
|
||||
kl = kl + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
+ (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij)
|
||||
end do ! l
|
||||
end do ! k
|
||||
end do ! j
|
||||
end do ! i
|
||||
end do ! p
|
||||
end do ! q
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
|
||||
0.d0, rho1(1,1,1), nBas*nBas)
|
||||
else
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
|
||||
1.d0, rho1(1,1,1), nBas*nBas)
|
||||
allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
|
||||
ERI_1 = 0.d0
|
||||
ERI_2 = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, c, d, cd, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = c+1, nBas-nR
|
||||
cd = cd + 1
|
||||
ERI_1(p,q,cd) = ERI(p,q,c,d) - ERI(p,q,d,c)
|
||||
enddo
|
||||
enddo
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = k+1, nO
|
||||
kl = kl + 1
|
||||
ERI_2(p,q,kl) = ERI(p,q,k,l) - ERI(p,q,l,k)
|
||||
end do
|
||||
end do
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
|
||||
0.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
|
||||
1.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
|
||||
0.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
|
||||
1.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
deallocate(ERI_1, ERI_2)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
|
||||
0.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
|
||||
1.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
deallocate(ERI_1, ERI_2)
|
||||
|
||||
|
||||
! !$OMP PARALLEL DEFAULT(NONE) &
|
||||
! !$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
|
||||
! !$OMP SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
|
||||
! !$OMP DO COLLAPSE(2)
|
||||
! do q = nC+1, nBas-nR
|
||||
! do p = nC+1, nBas-nR
|
||||
!
|
||||
! ab = 0
|
||||
!
|
||||
! do a = nO+1, nBas-nR
|
||||
! do b = a+1, nBas-nR
|
||||
!
|
||||
! ab = ab + 1
|
||||
!
|
||||
! cd = 0
|
||||
! do c = nO+1, nBas-nR
|
||||
! do d = c+1, nBas-nR
|
||||
!
|
||||
! cd = cd + 1
|
||||
!
|
||||
! rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
! + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab)
|
||||
! end do ! d
|
||||
! end do ! c
|
||||
!
|
||||
! kl = 0
|
||||
! do k = nC+1, nO
|
||||
! do l = k+1, nO
|
||||
!
|
||||
! kl = kl + 1
|
||||
!
|
||||
! rho1(p,q,ab) = rho1(p,q,ab) &
|
||||
! + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab)
|
||||
! end do ! l
|
||||
! end do ! k
|
||||
!
|
||||
! end do ! b
|
||||
! end do ! a
|
||||
!
|
||||
! ij = 0
|
||||
! do i = nC+1, nO
|
||||
! do j = i+1, nO
|
||||
!
|
||||
! ij = ij + 1
|
||||
!
|
||||
! cd = 0
|
||||
!
|
||||
! do c = nO+1, nBas-nR
|
||||
! do d = c+1, nBas-nR
|
||||
!
|
||||
! cd = cd + 1
|
||||
!
|
||||
! rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
! + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij)
|
||||
! end do ! d
|
||||
! end do ! c
|
||||
!
|
||||
! kl = 0
|
||||
! do k = nC+1, nO
|
||||
! do l = k+1, nO
|
||||
!
|
||||
! kl = kl + 1
|
||||
!
|
||||
! rho2(p,q,ij) = rho2(p,q,ij) &
|
||||
! + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij)
|
||||
! end do ! l
|
||||
! end do ! k
|
||||
!
|
||||
! end do ! j
|
||||
! end do ! i
|
||||
!
|
||||
! end do ! p
|
||||
! end do ! q
|
||||
! !$OMP END DO
|
||||
! !$OMP END PARALLEL
|
||||
|
||||
end if
|
||||
endif
|
||||
endif
|
||||
|
||||
!----------------------------------------------
|
||||
! alpha-beta block
|
||||
@ -270,125 +271,126 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1
|
||||
dim_1 = (nBas - nO) * (nBas - nO)
|
||||
dim_2 = nO * nO
|
||||
|
||||
allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
|
||||
ERI_1 = 0.d0
|
||||
ERI_2 = 0.d0
|
||||
if((dim_1 .eq. 0) .or. (dim_2 .eq. 0)) then
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, c, d, cd, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = nO+1, nBas-nR
|
||||
cd = cd + 1
|
||||
ERI_1(p,q,cd) = ERI(p,q,c,d)
|
||||
enddo
|
||||
enddo
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = nC+1, nO
|
||||
kl = kl + 1
|
||||
ERI_2(p,q,kl) = ERI(p,q,k,l)
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
|
||||
ab = 0
|
||||
do a = nO+1, nBas-nR
|
||||
do b = nO+1, nBas-nR
|
||||
|
||||
ab = ab + 1
|
||||
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = nO+1, nBas-nR
|
||||
|
||||
cd = cd + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab)
|
||||
end do
|
||||
end do
|
||||
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = nC+1, nO
|
||||
|
||||
kl = kl + 1
|
||||
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
ij = 0
|
||||
do i = nC+1, nO
|
||||
do j = nC+1, nO
|
||||
|
||||
ij = ij + 1
|
||||
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = nO+1, nBas-nR
|
||||
|
||||
cd = cd + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = nC+1, nO
|
||||
|
||||
kl = kl + 1
|
||||
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
else
|
||||
|
||||
allocate(ERI_1(nBas,nBas,dim_1), ERI_2(nBas,nBas,dim_2))
|
||||
ERI_1 = 0.d0
|
||||
ERI_2 = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(p, q, c, d, cd, k, l, kl) &
|
||||
!$OMP SHARED(nC, nBas, nR, nO, ERI_1, ERI_2, ERI)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do q = nC+1, nBas-nR
|
||||
do p = nC+1, nBas-nR
|
||||
cd = 0
|
||||
do c = nO+1, nBas-nR
|
||||
do d = nO+1, nBas-nR
|
||||
cd = cd + 1
|
||||
ERI_1(p,q,cd) = ERI(p,q,c,d)
|
||||
enddo
|
||||
enddo
|
||||
kl = 0
|
||||
do k = nC+1, nO
|
||||
do l = nC+1, nO
|
||||
kl = kl + 1
|
||||
ERI_2(p,q,kl) = ERI(p,q,k,l)
|
||||
end do
|
||||
end do
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
|
||||
0.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
|
||||
1.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
|
||||
0.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
|
||||
1.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
deallocate(ERI_1, ERI_2)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, &
|
||||
0.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_1, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y1(1,1), dim_2, &
|
||||
1.d0, rho1(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_1, 1.d0, &
|
||||
ERI_1(1,1,1), nBas*nBas, X2(1,1), dim_1, &
|
||||
0.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
call dgemm("N", "N", nBas*nBas, dim_2, dim_2, 1.d0, &
|
||||
ERI_2(1,1,1), nBas*nBas, Y2(1,1), dim_2, &
|
||||
1.d0, rho2(1,1,1), nBas*nBas)
|
||||
|
||||
deallocate(ERI_1, ERI_2)
|
||||
|
||||
|
||||
! !$OMP PARALLEL DEFAULT(NONE) &
|
||||
! !$OMP PRIVATE(p, q, a, b, ab, c, d, cd, i, j, ij, k, l, kl) &
|
||||
! !$OMP SHARED(nC, nBas, nR, nO, rho1, rho2, ERI, X1, Y1, X2, Y2)
|
||||
! !$OMP DO COLLAPSE(2)
|
||||
!
|
||||
! do q = nC+1, nBas-nR
|
||||
! do p = nC+1, nBas-nR
|
||||
!
|
||||
! ab = 0
|
||||
! do a = nO+1, nBas-nR
|
||||
! do b = nO+1, nBas-nR
|
||||
!
|
||||
! ab = ab + 1
|
||||
!
|
||||
! cd = 0
|
||||
! do c = nO+1, nBas-nR
|
||||
! do d = nO+1, nBas-nR
|
||||
!
|
||||
! cd = cd + 1
|
||||
!
|
||||
! rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab)
|
||||
! end do
|
||||
! end do
|
||||
!
|
||||
! kl = 0
|
||||
! do k = nC+1, nO
|
||||
! do l = nC+1, nO
|
||||
!
|
||||
! kl = kl + 1
|
||||
!
|
||||
! rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab)
|
||||
! end do
|
||||
! end do
|
||||
!
|
||||
! end do
|
||||
! end do
|
||||
!
|
||||
! ij = 0
|
||||
! do i = nC+1, nO
|
||||
! do j = nC+1, nO
|
||||
!
|
||||
! ij = ij + 1
|
||||
!
|
||||
! cd = 0
|
||||
! do c = nO+1, nBas-nR
|
||||
! do d = nO+1, nBas-nR
|
||||
!
|
||||
! cd = cd + 1
|
||||
!
|
||||
! rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij)
|
||||
! end do
|
||||
! end do
|
||||
!
|
||||
! kl = 0
|
||||
! do k = nC+1, nO
|
||||
! do l = nC+1, nO
|
||||
!
|
||||
! kl = kl + 1
|
||||
!
|
||||
! rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij)
|
||||
! end do
|
||||
! end do
|
||||
!
|
||||
! end do
|
||||
! end do
|
||||
!
|
||||
! end do
|
||||
! end do
|
||||
! !$OMP END DO
|
||||
! !$OMP END PARALLEL
|
||||
|
||||
end if
|
||||
endif
|
||||
endif
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -209,7 +209,9 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
|
||||
! Cumulant expansion !
|
||||
!--------------------!
|
||||
|
||||
call RGWC(dotest,nBas,nC,nO,nR,nS,Om,rho,eGW,Z)
|
||||
! TODO
|
||||
!call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eGW, eGW, Z)
|
||||
call RGWC(dotest, eta, nBas, nC, nO, nV, nR, nS, Om, rho, eHF, eHF, eGW, Z)
|
||||
|
||||
! Deallocate memory
|
||||
|
||||
|
||||
143
src/LR/ppLR.f90
143
src/LR/ppLR.f90
@ -1,93 +1,128 @@
|
||||
subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
|
||||
|
||||
! Solve the pp-RPA linear eigenvalue problem
|
||||
! ---
|
||||
|
||||
subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA)
|
||||
|
||||
!
|
||||
! Solve the pp-RPA linear eigenvalue problem
|
||||
!
|
||||
! right eigen-problem: H R = R w
|
||||
! left eigen-problem: H.T L = L w
|
||||
!
|
||||
! where L.T R = 1
|
||||
!
|
||||
!
|
||||
! (+C +B)
|
||||
! H = ( ) where C = C.T and D = D.T
|
||||
! (-B.T -D)
|
||||
!
|
||||
! (w1 0) (X1 X2) (+X1 +X2)
|
||||
! w = ( ), R = ( ) and L = ( )
|
||||
! (0 w2) (Y1 Y2) (-Y1 -Y2)
|
||||
!
|
||||
!
|
||||
! the normalisation condition reduces to
|
||||
!
|
||||
! X1.T X2 - Y1.T Y2 = 0
|
||||
! X1.T X1 - Y1.T Y1 = 1
|
||||
! X2.T X2 - Y2.T Y2 = 1
|
||||
!
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(in) :: TDA
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
double precision,intent(in) :: Bpp(nVV,nOO)
|
||||
double precision,intent(in) :: Cpp(nVV,nVV)
|
||||
double precision,intent(in) :: Dpp(nOO,nOO)
|
||||
logical, intent(in) :: TDA
|
||||
integer, intent(in) :: nOO, nVV
|
||||
double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO)
|
||||
double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV)
|
||||
double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO)
|
||||
double precision, intent(out) :: EcRPA
|
||||
|
||||
! Local variables
|
||||
logical :: imp_bio, verbose
|
||||
integer :: i, j, N
|
||||
double precision :: EcRPA1, EcRPA2
|
||||
double precision :: thr_d, thr_nd, thr_deg
|
||||
double precision,allocatable :: M(:,:), Z(:,:), Om(:)
|
||||
|
||||
double precision :: trace_matrix
|
||||
double precision :: EcRPA1
|
||||
double precision :: EcRPA2
|
||||
double precision,allocatable :: M(:,:)
|
||||
double precision,allocatable :: Z(:,:)
|
||||
double precision,allocatable :: Om(:)
|
||||
double precision, external :: trace_matrix
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Om1(nVV)
|
||||
double precision,intent(out) :: X1(nVV,nVV)
|
||||
double precision,intent(out) :: Y1(nOO,nVV)
|
||||
double precision,intent(out) :: Om2(nOO)
|
||||
double precision,intent(out) :: X2(nVV,nOO)
|
||||
double precision,intent(out) :: Y2(nOO,nOO)
|
||||
double precision,intent(out) :: EcRPA
|
||||
|
||||
! Memory allocation
|
||||
N = nOO + nVV
|
||||
|
||||
allocate(M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Om(nOO+nVV))
|
||||
|
||||
!-------------------------------------------------!
|
||||
! Solve the p-p eigenproblem !
|
||||
!-------------------------------------------------!
|
||||
! !
|
||||
! | C B | | X1 X2 | | w1 0 | | X1 X2 | !
|
||||
! | | | | = | | | | !
|
||||
! | -Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 | !
|
||||
! !
|
||||
!-------------------------------------------------!
|
||||
allocate(M(N,N), Z(N,N), Om(N))
|
||||
|
||||
if(TDA) then
|
||||
|
||||
X1(:,:) = +Cpp(:,:)
|
||||
Y1(:,:) = 0d0
|
||||
if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1)
|
||||
if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1)
|
||||
|
||||
X2(:,:) = 0d0
|
||||
Y2(:,:) = -Dpp(:,:)
|
||||
if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2)
|
||||
if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2)
|
||||
|
||||
else
|
||||
|
||||
! Diagonal blocks
|
||||
|
||||
! Diagonal blocks
|
||||
M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV)
|
||||
M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO)
|
||||
|
||||
! Off-diagonal blocks
|
||||
|
||||
! Off-diagonal blocks
|
||||
M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO)
|
||||
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO))
|
||||
|
||||
! call matout(nOO,nOO,Dpp)
|
||||
if((nOO .eq. 0) .or. (nVV .eq. 0)) then
|
||||
|
||||
! Diagonalize the p-p matrix
|
||||
! Diagonalize the p-p matrix
|
||||
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z)
|
||||
! Split the various quantities in p-p and h-h parts
|
||||
call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2)
|
||||
|
||||
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z)
|
||||
else
|
||||
|
||||
! Split the various quantities in p-p and h-h parts
|
||||
thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1
|
||||
thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1
|
||||
thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not
|
||||
imp_bio = .True. ! impose bi-orthogonality
|
||||
verbose = .False.
|
||||
call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose)
|
||||
|
||||
do i = 1, nOO
|
||||
Om2(i) = Om(i)
|
||||
do j = 1, nVV
|
||||
X2(j,i) = Z(j,i)
|
||||
enddo
|
||||
do j = 1, nOO
|
||||
Y2(j,i) = Z(nVV+j,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i = 1, nVV
|
||||
Om1(i) = Om(nOO+i)
|
||||
do j = 1, nVV
|
||||
X1(j,i) = M(j,nOO+i)
|
||||
enddo
|
||||
do j = 1, nOO
|
||||
Y1(j,i) = M(nVV+j,nOO+i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
|
||||
endif
|
||||
|
||||
end if
|
||||
|
||||
! Compute the RPA correlation energy
|
||||
! Compute the RPA correlation energy
|
||||
EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp))
|
||||
EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp)
|
||||
EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp)
|
||||
|
||||
EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp) )
|
||||
EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp)
|
||||
EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp)
|
||||
|
||||
if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) &
|
||||
if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then
|
||||
print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
|
||||
endif
|
||||
|
||||
deallocate(M, Z, Om)
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
@ -232,7 +232,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
|
||||
|
||||
call wall_time(start_CC)
|
||||
call RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, &
|
||||
maxSCF_CC, thresh_CC, max_diis_CC, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_MO, ENuc, ERHF, eHF, cHF)
|
||||
maxSCF_CC, thresh_CC, max_diis_CC, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
|
||||
@ -89,8 +89,8 @@ FIX_ORDER_OF_LIBS=-Wl,--start-group
|
||||
|
||||
if sys.platform in ["linux", "linux2"]:
|
||||
# compiler = compile_gfortran_linux
|
||||
# compiler = compile_ifort_linux
|
||||
compiler = compile_olympe
|
||||
compiler = compile_ifort_linux
|
||||
# compiler = compile_olympe
|
||||
elif sys.platform == "darwin":
|
||||
compiler = compile_gfortran_mac
|
||||
else:
|
||||
|
||||
@ -9,12 +9,12 @@ subroutine check_test_value(branch)
|
||||
! Local variables
|
||||
|
||||
character(len=30) :: description
|
||||
double precision :: value
|
||||
double precision :: val
|
||||
double precision :: reference
|
||||
character(len=15) :: answer
|
||||
|
||||
logical :: failed
|
||||
double precision,parameter :: cutoff = 1d-10
|
||||
double precision,parameter :: thresh = 1d-10
|
||||
|
||||
! Output variables
|
||||
|
||||
@ -45,19 +45,19 @@ subroutine check_test_value(branch)
|
||||
do
|
||||
|
||||
read(11,'(A30)',end=11) description
|
||||
read(11,'(F20.15)',end=11) value
|
||||
read(11,'(F20.15)',end=11) val
|
||||
|
||||
read(12,*,end=12)
|
||||
read(12,'(F20.15)',end=12) reference
|
||||
|
||||
if(abs(value-reference) < cutoff) then
|
||||
if(dabs(val-reference)/(1d-15+dabs(reference)) < thresh) then
|
||||
answer = '.......... :-)'
|
||||
else
|
||||
answer = '.......... :-( '
|
||||
failed = .true.
|
||||
end if
|
||||
write(*,'(1X,A1,1X,A30,1X,A1,1X,3F15.10,1X,A1,1X,A15,1X,A1)') &
|
||||
'|',description,'|',value,reference,abs(value-reference),'|',answer,'|'
|
||||
'|',description,'|',val,reference,abs(val-reference),'|',answer,'|'
|
||||
|
||||
end do
|
||||
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine dump_test_value(branch,description,value)
|
||||
subroutine dump_test_value(branch, description, val)
|
||||
|
||||
implicit none
|
||||
|
||||
@ -7,7 +7,7 @@ subroutine dump_test_value(branch,description,value)
|
||||
character(len=1),intent(in) :: branch
|
||||
character(len=*),intent(in) :: description
|
||||
|
||||
double precision,intent(in) :: value
|
||||
double precision,intent(in) :: val
|
||||
|
||||
! Local variables
|
||||
|
||||
@ -15,18 +15,19 @@ subroutine dump_test_value(branch,description,value)
|
||||
|
||||
if(branch == 'R') then
|
||||
|
||||
write(11,*) trim(description)
|
||||
write(11,'(F20.15)') value
|
||||
!write(1231597, '(A, ": ", F20.15)') '"' // trim(description) // '"', val
|
||||
write(1231597, *) trim(description)
|
||||
write(1231597, '(F20.15)') val
|
||||
|
||||
elseif(branch == 'U') then
|
||||
|
||||
write(12,*) trim(description)
|
||||
write(12,'(F20.15)') value
|
||||
write(1232584,*) trim(description)
|
||||
write(1232584,'(F20.15)') val
|
||||
|
||||
elseif(branch == 'G') then
|
||||
|
||||
write(13,*) trim(description)
|
||||
write(13,'(F20.15)') value
|
||||
write(1234181,*) trim(description)
|
||||
write(1234181,'(F20.15)') val
|
||||
|
||||
else
|
||||
|
||||
|
||||
@ -12,10 +12,10 @@ subroutine init_test(doRtest,doUtest,doGtest)
|
||||
|
||||
! Output variables
|
||||
|
||||
if(doRtest) open(unit=11,file='test/Rtest.dat')
|
||||
if(doRtest) open(unit=1231597, file='test/Rtest.dat')
|
||||
|
||||
if(doUtest) open(unit=12,file='test/Utest.dat')
|
||||
if(doUtest) open(unit=1232584, file='test/Utest.dat')
|
||||
|
||||
if(doGtest) open(unit=13,file='test/Gtest.dat')
|
||||
if(doGtest) open(unit=1234181, file='test/Gtest.dat')
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -12,10 +12,10 @@ subroutine stop_test(doRtest,doUtest,doGtest)
|
||||
|
||||
! Output variables
|
||||
|
||||
if(doRtest) close(unit=11)
|
||||
if(doRtest) close(unit=1231597)
|
||||
|
||||
if(doUtest) close(unit=12)
|
||||
if(doUtest) close(unit=1231597)
|
||||
|
||||
if(doGtest) close(unit=13)
|
||||
if(doGtest) close(unit=1234181)
|
||||
|
||||
end subroutine
|
||||
|
||||
571
src/utils/non_sym_diag.f90
Normal file
571
src/utils/non_sym_diag.f90
Normal file
@ -0,0 +1,571 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine diagonalize_nonsym_matrix(N, A, L, e_re, thr_d, thr_nd, thr_deg, imp_bio, verbose)
|
||||
|
||||
! Diagonalize a non-symmetric matrix A
|
||||
!
|
||||
! Output
|
||||
! right-eigenvectors are saved in A
|
||||
! left-eigenvectors are saved in L
|
||||
! eigenvalues are saved in e = e_re + i e_im
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: N
|
||||
logical, intent(in) :: imp_bio, verbose
|
||||
double precision, intent(in) :: thr_d, thr_nd, thr_deg
|
||||
double precision, intent(inout) :: A(N,N)
|
||||
double precision, intent(out) :: e_re(N), L(N,N)
|
||||
|
||||
integer :: i, j, ii
|
||||
integer :: lwork, info
|
||||
double precision :: accu_d, accu_nd
|
||||
integer, allocatable :: iorder(:), deg_num(:)
|
||||
double precision, allocatable :: Atmp(:,:), Ltmp(:,:), work(:), e_im(:)
|
||||
double precision, allocatable :: S(:,:)
|
||||
|
||||
if(verbose) then
|
||||
print*, ' Starting a non-Hermitian diagonalization ...'
|
||||
print*, ' Good Luck ;)'
|
||||
print*, ' imp_bio = ', imp_bio
|
||||
endif
|
||||
|
||||
! ---
|
||||
! diagonalize
|
||||
|
||||
allocate(Atmp(N,N), e_im(N))
|
||||
Atmp(1:N,1:N) = A(1:N,1:N)
|
||||
|
||||
allocate(work(1))
|
||||
lwork = -1
|
||||
call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info)
|
||||
if(info .gt. 0) then
|
||||
print*,'dgeev failed !!', info
|
||||
stop
|
||||
endif
|
||||
|
||||
lwork = max(int(work(1)), 1)
|
||||
deallocate(work)
|
||||
allocate(work(lwork))
|
||||
|
||||
call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info)
|
||||
if(info .ne. 0) then
|
||||
print*,'dgeev failed !!', info
|
||||
stop
|
||||
endif
|
||||
|
||||
deallocate(Atmp, WORK)
|
||||
|
||||
|
||||
! ---
|
||||
! check if eigenvalues are real
|
||||
|
||||
i = 1
|
||||
ii = 0
|
||||
do while(i .le. N)
|
||||
if(dabs(e_im(i)) .gt. 1.d-12) then
|
||||
ii = ii + 1
|
||||
if(verbose) then
|
||||
print*, ' Warning: complex eigenvalue !'
|
||||
print*, i, e_re(i), e_im(i)
|
||||
if(dabs(e_im(i)/e_re(i)) .lt. 1.d-6) then
|
||||
print*, ' small enouph to be igored'
|
||||
else
|
||||
print*, ' IMAGINARY PART IS SIGNIFANT !!!'
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
i = i + 1
|
||||
enddo
|
||||
|
||||
if(verbose) then
|
||||
if(ii .eq. 0) print*, ' congratulations :) eigenvalues are real-valued !!'
|
||||
endif
|
||||
|
||||
|
||||
! ---
|
||||
! track & sort the real eigenvalues
|
||||
|
||||
allocate(Atmp(N,N), Ltmp(N,N), iorder(N))
|
||||
|
||||
do i = 1, N
|
||||
iorder(i) = i
|
||||
enddo
|
||||
call quick_sort(e_re, iorder, N)
|
||||
|
||||
Atmp(:,:) = A(:,:)
|
||||
Ltmp(:,:) = L(:,:)
|
||||
do i = 1, N
|
||||
do j = 1, N
|
||||
A(j,i) = Atmp(j,iorder(i))
|
||||
L(j,i) = Ltmp(j,iorder(i))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
deallocate(Atmp, Ltmp, iorder)
|
||||
|
||||
|
||||
|
||||
|
||||
! ---
|
||||
! check bi-orthog
|
||||
|
||||
allocate(S(N,N))
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose)
|
||||
|
||||
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then
|
||||
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are normalized and bi-orthogonalized'
|
||||
endif
|
||||
|
||||
elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then
|
||||
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are not normalized but bi-orthogonalized'
|
||||
endif
|
||||
|
||||
call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.)
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose)
|
||||
|
||||
else
|
||||
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are not normalized neither bi-orthogonalized'
|
||||
endif
|
||||
|
||||
allocate(deg_num(N))
|
||||
call reorder_degen_eigvec(N, thr_deg, deg_num, e_re, L, A)
|
||||
call impose_biorthog_degen_eigvec(N, deg_num, e_re, L, A)
|
||||
deallocate(deg_num)
|
||||
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose)
|
||||
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are now normalized and bi-orthogonalized'
|
||||
endif
|
||||
elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then
|
||||
if(verbose) then
|
||||
print *, ' lapack vectors are now not normalized but bi-orthogonalized'
|
||||
endif
|
||||
call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.)
|
||||
call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose)
|
||||
else
|
||||
if(verbose) then
|
||||
print*, ' bi-orthogonalization failed !'
|
||||
endif
|
||||
if(imp_bio) then
|
||||
print*, ' bi-orthogonalization failed !'
|
||||
deallocate(S)
|
||||
stop
|
||||
endif
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
deallocate(S)
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot, verbose)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, m
|
||||
logical, intent(in) :: stop_ifnot, verbose
|
||||
double precision, intent(in) :: Vl(n,m), Vr(n,m)
|
||||
double precision, intent(in) :: thr_d, thr_nd
|
||||
double precision, intent(out) :: accu_d, accu_nd, S(m,m)
|
||||
|
||||
integer :: i, j
|
||||
double precision, allocatable :: SS(:,:)
|
||||
|
||||
if(verbose) then
|
||||
print *, ' check bi-orthogonality'
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, Vl, size(Vl, 1), Vr, size(Vr, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
accu_d = 0.d0
|
||||
accu_nd = 0.d0
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
if(i==j) then
|
||||
accu_d = accu_d + dabs(S(i,i))
|
||||
else
|
||||
accu_nd = accu_nd + S(j,i) * S(j,i)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
accu_nd = dsqrt(accu_nd) / dble(m)
|
||||
|
||||
if(verbose) then
|
||||
if((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) then
|
||||
print *, ' non bi-orthogonal vectors !'
|
||||
print *, ' accu_nd = ', accu_nd
|
||||
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
|
||||
else
|
||||
print *, ' vectors are bi-orthogonals'
|
||||
endif
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
if(stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) then
|
||||
print *, ' non bi-orthogonal vectors !'
|
||||
print *, ' accu_nd = ', accu_nd
|
||||
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
|
||||
stop
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, m
|
||||
logical, intent(in) :: stop_ifnot
|
||||
double precision, intent(in) :: thr_d, thr_nd
|
||||
double precision, intent(inout) :: Vl(n,m), Vr(n,m)
|
||||
|
||||
integer :: i, j
|
||||
double precision :: accu_d, accu_nd, s_tmp
|
||||
double precision, allocatable :: S(:,:)
|
||||
|
||||
! ---
|
||||
|
||||
allocate(S(m,m))
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, Vl, size(Vl, 1), Vr, size(Vr, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
do i = 1, m
|
||||
if(S(i,i) .lt. 0.d0) then
|
||||
do j = 1, n
|
||||
Vl(j,i) = -1.d0 * Vl(j,i)
|
||||
enddo
|
||||
S(i,i) = -S(i,i)
|
||||
endif
|
||||
enddo
|
||||
|
||||
accu_d = 0.d0
|
||||
accu_nd = 0.d0
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
if(i==j) then
|
||||
accu_d = accu_d + S(i,i)
|
||||
else
|
||||
accu_nd = accu_nd + S(j,i) * S(j,i)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
accu_nd = dsqrt(accu_nd) / dble(m)
|
||||
|
||||
! ---
|
||||
|
||||
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then
|
||||
|
||||
do i = 1, m
|
||||
if(S(i,i) <= 0.d0) then
|
||||
print *, ' negative overlap !'
|
||||
print *, i, S(i,i)
|
||||
exit
|
||||
endif
|
||||
if(dabs(S(i,i) - 1.d0) .gt. thr_d) then
|
||||
s_tmp = 1.d0 / dsqrt(S(i,i))
|
||||
do j = 1, n
|
||||
Vl(j,i) = Vl(j,i) * s_tmp
|
||||
Vr(j,i) = Vr(j,i) * s_tmp
|
||||
enddo
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, Vl, size(Vl, 1), Vr, size(Vr, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
accu_d = 0.d0
|
||||
accu_nd = 0.d0
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
if(i==j) then
|
||||
accu_d = accu_d + S(i,i)
|
||||
else
|
||||
accu_nd = accu_nd + S(j,i) * S(j,i)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
accu_nd = dsqrt(accu_nd) / dble(m)
|
||||
|
||||
deallocate(S)
|
||||
|
||||
! ---
|
||||
|
||||
if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then
|
||||
print *, accu_nd, thr_nd
|
||||
print *, dabs(accu_d-dble(m))/dble(m), thr_d
|
||||
print *, ' biorthog_binormalize failed !'
|
||||
stop
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine reorder_degen_eigvec(n, thr_deg, deg_num, e0, L0, R0)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n
|
||||
double precision, intent(in) :: thr_deg
|
||||
double precision, intent(inout) :: e0(n), L0(n,n), R0(n,n)
|
||||
integer, intent(out) :: deg_num(n)
|
||||
|
||||
logical :: complex_root
|
||||
integer :: i, j, k, m, ii, j_tmp
|
||||
double precision :: ei, ej, de
|
||||
double precision :: accu_d, accu_nd
|
||||
double precision :: e0_tmp, L0_tmp(n), R0_tmp(n)
|
||||
double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
|
||||
|
||||
do i = 1, n
|
||||
deg_num(i) = 1
|
||||
enddo
|
||||
|
||||
do i = 1, n-1
|
||||
ei = e0(i)
|
||||
|
||||
! already considered in degen vectors
|
||||
if(deg_num(i) .eq. 0) cycle
|
||||
|
||||
ii = 0
|
||||
do j = i+1, n
|
||||
ej = e0(j)
|
||||
de = dabs(ei - ej)
|
||||
|
||||
if(de .lt. thr_deg) then
|
||||
ii = ii + 1
|
||||
|
||||
j_tmp = i + ii
|
||||
|
||||
deg_num(j_tmp) = 0
|
||||
|
||||
e0_tmp = e0(j_tmp)
|
||||
e0(j_tmp) = e0(j)
|
||||
e0(j) = e0_tmp
|
||||
|
||||
L0_tmp(1:n) = L0(1:n,j_tmp)
|
||||
L0(1:n,j_tmp) = L0(1:n,j)
|
||||
L0(1:n,j) = L0_tmp(1:n)
|
||||
|
||||
R0_tmp(1:n) = R0(1:n,j_tmp)
|
||||
R0(1:n,j_tmp) = R0(1:n,j)
|
||||
R0(1:n,j) = R0_tmp(1:n)
|
||||
endif
|
||||
enddo
|
||||
|
||||
deg_num(i) = ii + 1
|
||||
enddo
|
||||
|
||||
ii = 0
|
||||
do i = 1, n
|
||||
if(deg_num(i) .gt. 1) then
|
||||
ii = ii + 1
|
||||
endif
|
||||
enddo
|
||||
|
||||
if(ii .eq. 0) then
|
||||
print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies'
|
||||
print*, ' rotations may change energy'
|
||||
stop
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, deg_num(n)
|
||||
double precision, intent(in) :: e0(n)
|
||||
double precision, intent(inout) :: L0(n,n), R0(n,n)
|
||||
|
||||
logical :: complex_root
|
||||
integer :: i, j, k, m
|
||||
double precision :: ei, ej, de
|
||||
double precision :: accu_d, accu_nd
|
||||
double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
|
||||
|
||||
!do i = 1, n
|
||||
! if(deg_num(i) .gt. 1) then
|
||||
! print *, ' degen on', i, deg_num(i), e0(i)
|
||||
! endif
|
||||
!enddo
|
||||
|
||||
! ---
|
||||
|
||||
do i = 1, n
|
||||
m = deg_num(i)
|
||||
|
||||
if(m .gt. 1) then
|
||||
|
||||
allocate(L(n,m), R(n,m), S(m,m))
|
||||
|
||||
do j = 1, m
|
||||
L(1:n,j) = L0(1:n,i+j-1)
|
||||
R(1:n,j) = R0(1:n,i+j-1)
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, L, size(L, 1), R, size(R, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
accu_nd = 0.d0
|
||||
do j = 1, m
|
||||
do k = 1, m
|
||||
if(j==k) cycle
|
||||
accu_nd = accu_nd + dabs(S(j,k))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if(accu_nd .lt. 1d-12) then
|
||||
deallocate(S, L, R)
|
||||
cycle
|
||||
endif
|
||||
|
||||
call impose_biorthog_svd(n, m, L, R)
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, L, size(L, 1), R, size(R, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
accu_nd = 0.d0
|
||||
do j = 1, m
|
||||
do k = 1, m
|
||||
if(j==k) cycle
|
||||
accu_nd = accu_nd + dabs(S(j,k))
|
||||
enddo
|
||||
enddo
|
||||
if(accu_nd .gt. 1d-12) then
|
||||
print*, ' accu_nd =', accu_nd
|
||||
print*, ' your strategy for degenerates orbitals failed !'
|
||||
print*, m, 'deg on', i
|
||||
stop
|
||||
endif
|
||||
|
||||
deallocate(S)
|
||||
|
||||
! ---
|
||||
|
||||
do j = 1, m
|
||||
L0(1:n,i+j-1) = L(1:n,j)
|
||||
R0(1:n,i+j-1) = R(1:n,j)
|
||||
enddo
|
||||
|
||||
deallocate(L, R)
|
||||
|
||||
endif
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine impose_biorthog_svd(n, m, L, R)
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, m
|
||||
double precision, intent(inout) :: L(n,m), R(n,m)
|
||||
|
||||
integer :: i, j, num_linear_dependencies
|
||||
double precision :: threshold
|
||||
double precision, allocatable :: S(:,:), tmp(:,:)
|
||||
double precision, allocatable :: U(:,:), V(:,:), Vt(:,:), D(:)
|
||||
|
||||
allocate(S(m,m))
|
||||
|
||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||
, L, size(L, 1), R, size(R, 1) &
|
||||
, 0.d0, S, size(S, 1) )
|
||||
|
||||
! ---
|
||||
|
||||
allocate(U(m,m), Vt(m,m), D(m))
|
||||
|
||||
call svd(S, m, U, m, D, Vt, m, m, m)
|
||||
|
||||
deallocate(S)
|
||||
|
||||
threshold = 1.d-6
|
||||
num_linear_dependencies = 0
|
||||
do i = 1, m
|
||||
if(abs(D(i)) <= threshold) then
|
||||
D(i) = 0.d0
|
||||
num_linear_dependencies = num_linear_dependencies + 1
|
||||
else
|
||||
D(i) = 1.d0 / dsqrt(D(i))
|
||||
endif
|
||||
enddo
|
||||
if(num_linear_dependencies > 0) then
|
||||
write(*,*) ' linear dependencies = ', num_linear_dependencies
|
||||
write(*,*) ' m = ', m
|
||||
stop
|
||||
endif
|
||||
|
||||
allocate(V(m,m))
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
V(j,i) = Vt(i,j)
|
||||
enddo
|
||||
enddo
|
||||
deallocate(Vt)
|
||||
|
||||
! ---
|
||||
|
||||
! R <-- R x V x D^{-0.5}
|
||||
! L <-- L x U x D^{-0.5}
|
||||
|
||||
do i = 1, m
|
||||
do j = 1, m
|
||||
V(j,i) = V(j,i) * D(i)
|
||||
U(j,i) = U(j,i) * D(i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
allocate(tmp(n,m))
|
||||
tmp(:,:) = R(:,:)
|
||||
call dgemm( 'N', 'N', n, m, m, 1.d0 &
|
||||
, tmp, size(tmp, 1), V, size(V, 1) &
|
||||
, 0.d0, R, size(R, 1))
|
||||
|
||||
tmp(:,:) = L(:,:)
|
||||
call dgemm( 'N', 'N', n, m, m, 1.d0 &
|
||||
, tmp, size(tmp, 1), U, size(U, 1) &
|
||||
, 0.d0, L, size(L, 1))
|
||||
|
||||
deallocate(tmp, U, V, D)
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
26
test/export_tobench.py
Normal file
26
test/export_tobench.py
Normal file
@ -0,0 +1,26 @@
|
||||
|
||||
import sys
|
||||
|
||||
|
||||
def read_quantities_from_file(filename):
|
||||
quantities = {}
|
||||
|
||||
with open(filename, 'r') as file:
|
||||
lines = file.readlines()
|
||||
for i in range(0, len(lines), 2):
|
||||
# Remove any leading or trailing whitespace/newline characters
|
||||
quantity_name = lines[i].strip()
|
||||
quantity_value = float(lines[i+1].strip())
|
||||
quantities[quantity_name] = quantity_value
|
||||
|
||||
return quantities
|
||||
|
||||
def print_quantities(quantities):
|
||||
for key, value in quantities.items():
|
||||
print(f'"{key}": {value},')
|
||||
|
||||
filename = sys.argv[1]
|
||||
|
||||
quantities = read_quantities_from_file(filename)
|
||||
print_quantities(quantities)
|
||||
|
||||
7
tests/.gitignore
vendored
Normal file
7
tests/.gitignore
vendored
Normal file
@ -0,0 +1,7 @@
|
||||
|
||||
FeatherBench.db
|
||||
FeatherBench.json
|
||||
|
||||
*.xyz
|
||||
work
|
||||
|
||||
0
tests/balance_bench.py
Normal file
0
tests/balance_bench.py
Normal file
52
tests/create_database.py
Normal file
52
tests/create_database.py
Normal file
@ -0,0 +1,52 @@
|
||||
|
||||
import argparse
|
||||
|
||||
from molecule import save_molecules_to_json, load_molecules_from_json
|
||||
from molecule import create_database, add_molecule_to_db, remove_database
|
||||
|
||||
from feather_bench import FeatherBench
|
||||
|
||||
|
||||
parser = argparse.ArgumentParser(description="Benchmark Data Sets")
|
||||
|
||||
parser.add_argument(
|
||||
'-s', '--set_type',
|
||||
choices=['light', 'medium', 'heavy'],
|
||||
default='light',
|
||||
help="Specify the type of data set: light (default), medium, or heavy."
|
||||
)
|
||||
|
||||
|
||||
args = parser.parse_args()
|
||||
|
||||
if args.set_type == 'light':
|
||||
bench = 'FeatherBench'
|
||||
bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
|
||||
elif args.set_type == 'medium':
|
||||
bench = 'BalanceBench'
|
||||
bench_title = "\n\nSelected Medium Benchmark: {}\n\n".format(bench)
|
||||
elif args.set_type == 'heavy':
|
||||
bench = 'TitanBench'
|
||||
bench_title = "\n\nSelected Heavy Benchmark: {}\n\n".format(bench)
|
||||
else:
|
||||
bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
|
||||
|
||||
|
||||
db_name = '{}.db'.format(bench)
|
||||
|
||||
|
||||
# Save molecules to JSON
|
||||
#save_molecules_to_json(FeatherBench, 'FeatherBench.json')
|
||||
|
||||
# Load molecules from JSON
|
||||
#loaded_molecules = load_molecules_from_json('FeatherBench.json')
|
||||
#print(loaded_molecules)
|
||||
|
||||
#remove_database(db_name)
|
||||
|
||||
create_database(db_name)
|
||||
for molecule in FeatherBench:
|
||||
add_molecule_to_db(db_name, molecule)
|
||||
|
||||
|
||||
|
||||
113
tests/feather_bench.py
Normal file
113
tests/feather_bench.py
Normal file
@ -0,0 +1,113 @@
|
||||
|
||||
from molecule import Molecule
|
||||
|
||||
|
||||
He = Molecule(
|
||||
name="He",
|
||||
multiplicity=1,
|
||||
geometry=[
|
||||
{"element": "He", "x": 0.0, "y": 0.0, "z": 0.0}
|
||||
],
|
||||
properties={
|
||||
"properties_rhf":{
|
||||
"6-31g": {
|
||||
"RHF energy": -2.855160426154444,
|
||||
"RHF HOMO energy": -0.914126628640145,
|
||||
"RHF LUMO energy": 1.399859335255765,
|
||||
"RHF dipole moment": 0.0,
|
||||
"MP2 correlation energy": -0.011200122909934,
|
||||
"CCD correlation energy": -0.014985063116,
|
||||
"CCSD correlation energy": -0.015001711549092,
|
||||
"drCCD correlation energy": -0.01884537385338,
|
||||
"rCCD correlation energy": -0.016836322809386,
|
||||
"crCCD correlation energy": 0.008524676641474,
|
||||
"lCCD correlation energy": -0.00808242082105,
|
||||
"CIS singlet excitation energy": 1.911193619991987,
|
||||
"CIS triplet excitation energy": 1.455852629458543,
|
||||
"phRPA correlation energy": -0.018845374128748,
|
||||
"phRPAx correlation energy": -0.015760565120758,
|
||||
"crRPA correlation energy": -0.008868581132249,
|
||||
"ppRPA correlation energy": -0.008082420814972,
|
||||
"G0F2 correlation energy": -0.011438430540104,
|
||||
"G0F2 HOMO energy": -0.882696116274599,
|
||||
"G0F2 LUMO energy": 1.383080391842522,
|
||||
"G0W0 correlation energy": -0.019314094399372,
|
||||
"G0W0 HOMO energy": -0.87053388021722,
|
||||
"G0W0 LUMO energy": 1.377171287041735,
|
||||
"evGW correlation energy": -0.019335511771337,
|
||||
"evGW HOMO energy": -0.868460640984803,
|
||||
"evGW LUMO energy": 1.376287581502582,
|
||||
"G0T0pp correlation energy": -0.008161908540634,
|
||||
"G0T0pp HOMO energy": -0.898869172597701,
|
||||
"G0T0pp LUMO energy": 1.383928087417952,
|
||||
}
|
||||
},
|
||||
"properties_uhf":{
|
||||
"6-31g": {
|
||||
}
|
||||
},
|
||||
"properties_ghf":{
|
||||
"6-31g": {
|
||||
}
|
||||
},
|
||||
}
|
||||
)
|
||||
|
||||
# ---
|
||||
|
||||
H2O = Molecule(
|
||||
name="H2O",
|
||||
multiplicity=1,
|
||||
geometry=[
|
||||
{"element": "O", "x": 0.0000, "y": 0.0000, "z": 0.0000},
|
||||
{"element": "H", "x": 0.7571, "y": 0.0000, "z": 0.5861},
|
||||
{"element": "H", "x": -0.7571, "y": 0.0000, "z": 0.5861}
|
||||
],
|
||||
properties={
|
||||
"properties_rhf":{
|
||||
"cc-pvdz": {
|
||||
"RHF energy": -85.21935817501823,
|
||||
"RHF HOMO energy": -0.493132793449897,
|
||||
"RHF LUMO energy": 0.185534869842355,
|
||||
"RHF dipole moment": 0.233813698748474,
|
||||
"MP2 correlation energy": -0.203978216774657,
|
||||
"CCD correlation energy": -0.212571260121257,
|
||||
"CCSD correlation energy": -0.213302190845899,
|
||||
"drCCD correlation energy": -0.231281853419338,
|
||||
"rCCD correlation energy": -0.277238348710547,
|
||||
"crCCD correlation energy": 0.18014617422324,
|
||||
"lCCD correlation energy": -0.15128653432796,
|
||||
"CIS singlet excitation energy": 0.338828950934568,
|
||||
"CIS triplet excitation energy": 0.304873339484139,
|
||||
"phRPA correlation energy": -0.231281866582435,
|
||||
"phRPAx correlation energy": -0.310796738307943,
|
||||
"crRPA correlation energy": -0.246289801609294,
|
||||
"ppRPA correlation energy": -0.151286536255888,
|
||||
"G0F2 correlation energy": -0.217807591229668,
|
||||
"G0F2 HOMO energy": -0.404541451101377,
|
||||
"G0F2 LUMO energy": 0.16650398400197,
|
||||
"G0W0 correlation energy": -0.23853664665404,
|
||||
"G0W0 HOMO energy": -0.446828623007469,
|
||||
"G0W0 LUMO energy": 0.173026609033024,
|
||||
"evGW correlation energy": -0.239414217281308,
|
||||
"evGW HOMO energy": -0.443076613314424,
|
||||
"evGW LUMO energy": 0.172691758111392,
|
||||
"G0T0pp correlation energy": -0.156214864467344,
|
||||
"G0T0pp HOMO energy": -0.452117482732615,
|
||||
"G0T0pp LUMO energy": 0.16679206983464,
|
||||
}
|
||||
}
|
||||
}
|
||||
)
|
||||
|
||||
# ---
|
||||
|
||||
FeatherBench = [
|
||||
He,
|
||||
H2O
|
||||
]
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
22
tests/inp/methods.RHF
Normal file
22
tests/inp/methods.RHF
Normal file
@ -0,0 +1,22 @@
|
||||
# RHF UHF GHF ROHF
|
||||
T F F F
|
||||
# MP2 MP3
|
||||
T T
|
||||
# CCD pCCD DCD CCSD CCSD(T)
|
||||
T F F T F
|
||||
# drCCD rCCD crCCD lCCD
|
||||
T T T T
|
||||
# CIS CIS(D) CID CISD FCI
|
||||
T F F F F
|
||||
# phRPA phRPAx crRPA ppRPA
|
||||
T T T T
|
||||
# G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3
|
||||
T F F F F F
|
||||
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
|
||||
T T F F F F
|
||||
# G0T0pp evGTpp qsGTpp ufG0T0pp
|
||||
T F F F
|
||||
# G0T0eh evGTeh qsGTeh
|
||||
F F F
|
||||
# Rtest Utest Gtest
|
||||
T F F
|
||||
18
tests/inp/options.RHF
Normal file
18
tests/inp/options.RHF
Normal file
@ -0,0 +1,18 @@
|
||||
# HF: maxSCF thresh DIIS guess mix shift stab search
|
||||
10000 0.0000001 5 1 0.0 0.0 F F
|
||||
# MP: reg
|
||||
F
|
||||
# CC: maxSCF thresh DIIS
|
||||
64 0.0000001 5
|
||||
# spin: TDA singlet triplet
|
||||
F T T
|
||||
# GF: maxSCF thresh DIIS lin eta renorm reg
|
||||
256 0.00001 5 F 0.0 0 F
|
||||
# GW: maxSCF thresh DIIS lin eta TDA_W reg
|
||||
256 0.00001 5 F 0.0 F F
|
||||
# GT: maxSCF thresh DIIS lin eta TDA_T reg
|
||||
256 0.00001 5 F 0.0 F F
|
||||
# ACFDT: AC Kx XBS
|
||||
F F T
|
||||
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
|
||||
F F F F T
|
||||
246
tests/lunch_bench.py
Normal file
246
tests/lunch_bench.py
Normal file
@ -0,0 +1,246 @@
|
||||
|
||||
import time
|
||||
import threading
|
||||
import sys
|
||||
import os
|
||||
import shutil
|
||||
from pathlib import Path
|
||||
import subprocess
|
||||
import platform
|
||||
from datetime import datetime
|
||||
import argparse
|
||||
|
||||
from molecule import get_molecules_from_db
|
||||
from molecule import generate_xyz
|
||||
from utils import print_col, stdout_col
|
||||
|
||||
|
||||
current_date = datetime.now()
|
||||
|
||||
quack_root = os.getenv('QUACK_ROOT')
|
||||
|
||||
# User Name
|
||||
user_name = os.getlogin()
|
||||
|
||||
# Operating System
|
||||
os_name = platform.system()
|
||||
os_release = platform.release()
|
||||
os_version = platform.version()
|
||||
|
||||
# CPU Information
|
||||
machine = platform.machine()
|
||||
processor = platform.processor()
|
||||
|
||||
# System Architecture
|
||||
architecture = platform.architecture()[0]
|
||||
|
||||
# Python Version
|
||||
python_version_full = platform.python_version_tuple()
|
||||
PYTHON_VERSION = "{}.{}".format(python_version_full[0], python_version_full[1])
|
||||
|
||||
|
||||
|
||||
print(f"The current date and time is {current_date.strftime('%Y-%m-%d %H:%M:%S')}")
|
||||
print(f"User Name: {user_name}")
|
||||
print(f"Operating System: {os_name} {os_release} ({os_version})")
|
||||
print(f"CPU: {processor} ({machine})")
|
||||
print(f"System Architecture: {architecture}")
|
||||
print(f"QUACK_ROOT: {quack_root}")
|
||||
print(f"Python version: {python_version_full}\n\n")
|
||||
|
||||
|
||||
parser = argparse.ArgumentParser(description="Benchmark Data Sets")
|
||||
|
||||
parser.add_argument(
|
||||
'-s', '--set_type',
|
||||
choices=['light', 'medium', 'heavy'],
|
||||
default='light',
|
||||
help="Specify the type of data set: light (default), medium, or heavy."
|
||||
)
|
||||
parser.add_argument(
|
||||
'-t', '--thresh',
|
||||
type=float,
|
||||
default=1e-7,
|
||||
help='Threshold for acceptable difference (default: 1e-8)'
|
||||
)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
args = parser.parse_args()
|
||||
|
||||
THRESH = args.thresh
|
||||
|
||||
if args.set_type == 'light':
|
||||
bench = 'FeatherBench'
|
||||
bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
|
||||
elif args.set_type == 'medium':
|
||||
bench = 'BalanceBench'
|
||||
bench_title = "\n\nSelected Medium Benchmark: {}\n\n".format(bench)
|
||||
elif args.set_type == 'heavy':
|
||||
bench = 'TitanBench'
|
||||
bench_title = "\n\nSelected Heavy Benchmark: {}\n\n".format(bench)
|
||||
else:
|
||||
bench_title = "\n\nSelected Light Benchmark: {}\n\n".format(bench)
|
||||
|
||||
print(bench_title.center(150, '-'))
|
||||
print("\n\n")
|
||||
|
||||
# ---
|
||||
|
||||
class Quack_Job:
|
||||
|
||||
def __init__(self, mol, multip, basis, geom, methd):
|
||||
self.mol = mol
|
||||
self.multip = multip
|
||||
self.basis = basis
|
||||
self.geom = geom
|
||||
self.methd = methd
|
||||
|
||||
def prep_inp(self):
|
||||
|
||||
# geometry
|
||||
generate_xyz(self.geom, filename="{}/mol/{}.xyz".format(quack_root, self.mol))
|
||||
|
||||
# input files
|
||||
for inp in ["methods", "options"]:
|
||||
inp_file = "{}.{}".format(inp, self.methd.upper())
|
||||
if os.path.exists("inp/{}".format(inp_file)):
|
||||
shutil.copy("{}/tests/inp/{}".format(quack_root, inp_file),
|
||||
"{}/input/{}".format(quack_root, inp))
|
||||
else:
|
||||
print_col("File 'inp/{}' does not exist.".format(inp_file), "red")
|
||||
sys.exit(1)
|
||||
|
||||
def run(self, work_path):
|
||||
|
||||
def display_spinner():
|
||||
spinner = ['|', '/', '-', '\\']
|
||||
idx = 0
|
||||
while not done_event.is_set():
|
||||
stdout_col(f'\r Testing {self.methd} ({self.basis}) {spinner[idx]}', "cyan")
|
||||
sys.stdout.flush()
|
||||
idx = (idx + 1) % len(spinner)
|
||||
time.sleep(0.05)
|
||||
stdout_col(f'\r Testing {self.methd} ({self.basis}) \n\n', "cyan")
|
||||
|
||||
done_event = threading.Event()
|
||||
spinner_thread = threading.Thread(target=display_spinner)
|
||||
spinner_thread.start()
|
||||
|
||||
try:
|
||||
|
||||
os.chdir('..')
|
||||
#print_col(f" Starting QuAck..", "magenta")
|
||||
#print_col(f" $ cd ..", "magenta")
|
||||
|
||||
command = [
|
||||
'python{}'.format(PYTHON_VERSION), 'PyDuck.py',
|
||||
'-x', '{}'.format(self.mol),
|
||||
'-b', '{}'.format(self.basis),
|
||||
'-m', '{}'.format(self.multip)
|
||||
]
|
||||
#print_col(f" $ {' '.join(command)}", "magenta")
|
||||
|
||||
file_out = "{}/{}/{}_{}_{}.out".format(work_path, self.methd, self.mol, self.multip, self.basis)
|
||||
with open(file_out, 'w') as fobj:
|
||||
result = subprocess.run(command, stdout=fobj, stderr=subprocess.PIPE, text=True)
|
||||
if result.stderr:
|
||||
print("Error output:", result.stderr)
|
||||
|
||||
os.chdir('tests')
|
||||
#print_col(f" $ cd tests", "magenta")
|
||||
|
||||
except Exception as e:
|
||||
|
||||
print_col(f"An error occurred: {str(e)}", "red")
|
||||
|
||||
finally:
|
||||
|
||||
done_event.set()
|
||||
spinner_thread.join()
|
||||
|
||||
def check_data(self, data_ref):
|
||||
filepath = '../test/Rtest.dat'
|
||||
data_new = {}
|
||||
try:
|
||||
# read data_new
|
||||
with open(filepath, 'r') as f:
|
||||
lines = f.readlines()
|
||||
for i in range(0, len(lines) - 1, 2):
|
||||
key = lines[i].strip()
|
||||
value = lines[i + 1].strip()
|
||||
data_new[key] = float(value) # Convert value to float
|
||||
|
||||
# Compare with data_ref
|
||||
for key in data_ref:
|
||||
if key not in data_new:
|
||||
print_col(f" 😐 {key} missing ⚠️ ", "yellow")
|
||||
else:
|
||||
diff = abs(data_new[key] - data_ref[key]) / (1e-15 + abs(data_ref[key]))
|
||||
if(diff <= THRESH):
|
||||
print_col(f" 🙂 {key}", "green")
|
||||
else:
|
||||
print_col(f" ☹️ {key}: ❌ {data_ref[key]} ≠ {data_new[key]}", "red")
|
||||
except FileNotFoundError:
|
||||
print_col(f"Error: The file '{filepath}' does not exist.", "red")
|
||||
sys.exit(1)
|
||||
except Exception as e:
|
||||
print_col(f"An error occurred: {str(e)}", "red")
|
||||
sys.exit(1)
|
||||
|
||||
|
||||
# ---
|
||||
|
||||
|
||||
def main():
|
||||
|
||||
work_path = Path('{}/tests/work'.format(quack_root))
|
||||
if not work_path.exists():
|
||||
work_path.mkdir(parents=True, exist_ok=True)
|
||||
print(f"Directory '{work_path}' created.\n")
|
||||
|
||||
for mol in molecules:
|
||||
|
||||
mol_name = mol.name
|
||||
mol_mult = mol.multiplicity
|
||||
mol_geom = mol.geometry
|
||||
mol_data = mol.properties
|
||||
|
||||
print_col(" Molecule: {} (2S+1 = {})".format(mol_name, mol_mult), "blue")
|
||||
|
||||
for mol_prop_name, mol_prop_data in mol_data.items():
|
||||
|
||||
methd = mol_prop_name[len('properties_'):]
|
||||
|
||||
if(len(mol_prop_data) == 0):
|
||||
continue
|
||||
|
||||
for basis_name, basis_data in mol_prop_data.items():
|
||||
|
||||
if(len(basis_data) == 0):
|
||||
continue
|
||||
|
||||
work_methd = Path('{}/{}'.format(work_path, methd))
|
||||
if not work_methd.exists():
|
||||
work_methd.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
New_Quack_Job = Quack_Job(mol_name, mol_mult, basis_name, mol_geom, methd)
|
||||
New_Quack_Job.prep_inp()
|
||||
New_Quack_Job.run(work_path)
|
||||
New_Quack_Job.check_data(basis_data)
|
||||
|
||||
print()
|
||||
print()
|
||||
print()
|
||||
|
||||
quit()
|
||||
|
||||
|
||||
db_name = '{}.db'.format(bench)
|
||||
|
||||
molecules = get_molecules_from_db(db_name)
|
||||
|
||||
main()
|
||||
|
||||
150
tests/molecule.py
Normal file
150
tests/molecule.py
Normal file
@ -0,0 +1,150 @@
|
||||
|
||||
import os
|
||||
import json
|
||||
import sqlite3
|
||||
|
||||
from utils import print_col
|
||||
|
||||
|
||||
class Molecule:
|
||||
def __init__(self, name, multiplicity, geometry, properties):
|
||||
self.name = name
|
||||
self.multiplicity = multiplicity
|
||||
self.geometry = geometry
|
||||
self.properties = properties
|
||||
|
||||
def to_dict(self):
|
||||
return {
|
||||
"name": self.name,
|
||||
"multiplicity": self.multiplicity,
|
||||
"geometry": self.geometry,
|
||||
"properties": self.properties,
|
||||
}
|
||||
|
||||
@staticmethod
|
||||
def from_dict(data):
|
||||
return Molecule(
|
||||
name=data["name"],
|
||||
multiplicity=data["multiplicity"],
|
||||
geometry=data["geometry"],
|
||||
properties=data["properties"]
|
||||
)
|
||||
|
||||
def save_molecules_to_json(molecules, filename):
|
||||
with open(filename, 'w') as f:
|
||||
json_data = [molecule.to_dict() for molecule in molecules]
|
||||
json.dump(json_data, f, indent=4)
|
||||
|
||||
def load_molecules_from_json(filename):
|
||||
with open(filename, 'r') as f:
|
||||
json_data = json.load(f)
|
||||
return [Molecule.from_dict(data) for data in json_data]
|
||||
|
||||
|
||||
def create_database(db_name):
|
||||
if os.path.exists(db_name):
|
||||
conn = sqlite3.connect(db_name)
|
||||
cursor = conn.cursor()
|
||||
# Check if the table already exists
|
||||
cursor.execute("SELECT name FROM sqlite_master WHERE type='table' AND name='molecules';")
|
||||
table_exists = cursor.fetchone()
|
||||
|
||||
if table_exists:
|
||||
print_col(f"Database '{db_name}' already exists and table 'molecules' is already created.", "yellow")
|
||||
else:
|
||||
# Create the table if it does not exist
|
||||
cursor.execute('''CREATE TABLE molecules
|
||||
(name TEXT, multiplicity INTEGER, geometry TEXT, properties TEXT)''')
|
||||
conn.commit()
|
||||
print_col(f"Table 'molecules' created in existing database '{db_name}' successfully.", "green")
|
||||
conn.close()
|
||||
else:
|
||||
# Create the database and table
|
||||
conn = sqlite3.connect(db_name)
|
||||
cursor = conn.cursor()
|
||||
cursor.execute('''CREATE TABLE molecules
|
||||
(name TEXT, multiplicity INTEGER, geometry TEXT, properties TEXT)''')
|
||||
conn.commit()
|
||||
conn.close()
|
||||
print_col(f"Database '{db_name}' created and table 'molecules' added successfully.", "green")
|
||||
|
||||
def add_molecule_to_db(db_name, molecule):
|
||||
|
||||
conn = sqlite3.connect(db_name)
|
||||
cursor = conn.cursor()
|
||||
|
||||
# Convert geometry and properties to JSON strings
|
||||
geometry_str = json.dumps(molecule.geometry)
|
||||
energies_str = json.dumps(molecule.properties)
|
||||
|
||||
# Check if the molecule already exists
|
||||
cursor.execute("SELECT COUNT(*) FROM molecules WHERE name = ?", (molecule.name,))
|
||||
count = cursor.fetchone()[0]
|
||||
|
||||
if count > 0:
|
||||
print_col(f"Molecule '{molecule.name}' already exists in {db_name}.", "yellow")
|
||||
else:
|
||||
# Insert the molecule if it does not exist
|
||||
cursor.execute("INSERT INTO molecules (name, multiplicity, geometry, properties) VALUES (?, ?, ?, ?)",
|
||||
(molecule.name, molecule.multiplicity, geometry_str, energies_str))
|
||||
conn.commit()
|
||||
print_col(f"'{molecule.name}' added to {db_name} successfully.", "green")
|
||||
|
||||
conn.close()
|
||||
|
||||
|
||||
def remove_database(db_name):
|
||||
if os.path.exists(db_name):
|
||||
os.remove(db_name)
|
||||
print_col(f"Database '{db_name}' removed successfully.", "red")
|
||||
else:
|
||||
print_col(f"Database '{db_name}' does not exist.", "red")
|
||||
|
||||
def get_molecules_from_db(db_name):
|
||||
conn = sqlite3.connect(db_name)
|
||||
cursor = conn.cursor()
|
||||
cursor.execute("SELECT name, multiplicity, geometry, properties FROM molecules")
|
||||
rows = cursor.fetchall()
|
||||
molecules = []
|
||||
for row in rows:
|
||||
name, multiplicity, geometry_str, energies_str = row
|
||||
geometry = json.loads(geometry_str)
|
||||
properties = json.loads(energies_str)
|
||||
molecules.append(Molecule(name, multiplicity, geometry, properties))
|
||||
conn.close()
|
||||
return molecules
|
||||
|
||||
|
||||
def generate_xyz(elements, filename="output.xyz", verbose=False):
|
||||
"""
|
||||
Generate an XYZ file from a list of elements.
|
||||
|
||||
Parameters:
|
||||
elements (list): A list of dictionaries, where each dictionary represents
|
||||
an atom with its element and x, y, z coordinates.
|
||||
filename (str): The name of the output XYZ file. Default is 'output.xyz'.
|
||||
"""
|
||||
|
||||
# Get the number of atoms
|
||||
num_atoms = len(elements)
|
||||
|
||||
# Open the file in write mode
|
||||
with open(filename, 'w') as f:
|
||||
# Write the number of atoms
|
||||
f.write(f"{num_atoms}\n")
|
||||
|
||||
# Write a comment line (can be left blank or customized)
|
||||
f.write("XYZ file generated by generate_xyz function\n")
|
||||
|
||||
# Write the element and coordinates
|
||||
for atom in elements:
|
||||
element = atom['element']
|
||||
x = atom['x']
|
||||
y = atom['y']
|
||||
z = atom['z']
|
||||
f.write(f"{element} {x:.6f} {y:.6f} {z:.6f}\n")
|
||||
|
||||
if(verbose):
|
||||
print(f"XYZ file '{filename}' generated successfully!")
|
||||
|
||||
|
||||
0
tests/titan_bench.py
Normal file
0
tests/titan_bench.py
Normal file
88
tests/utils.py
Normal file
88
tests/utils.py
Normal file
@ -0,0 +1,88 @@
|
||||
|
||||
import sys
|
||||
|
||||
|
||||
|
||||
|
||||
def print_col(text, color):
|
||||
|
||||
if(color == "black"):
|
||||
|
||||
print("\033[30m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "red"):
|
||||
|
||||
print("\033[31m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "green"):
|
||||
|
||||
print("\033[32m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "yellow"):
|
||||
|
||||
print("\033[33m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "blue"):
|
||||
|
||||
print("\033[34m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "magenta"):
|
||||
|
||||
print("\033[35m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "cyan"):
|
||||
|
||||
print("\033[36m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "white"):
|
||||
|
||||
print("\033[37m{}\033[0m".format(text))
|
||||
|
||||
else:
|
||||
|
||||
print("{}".format(text))
|
||||
|
||||
|
||||
# ---
|
||||
|
||||
def stdout_col(text, color):
|
||||
|
||||
if(color == "black"):
|
||||
|
||||
sys.stdout.write("\033[30m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "red"):
|
||||
|
||||
sys.stdout.write("\033[31m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "green"):
|
||||
|
||||
sys.stdout.write("\033[32m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "yellow"):
|
||||
|
||||
sys.stdout.write("\033[33m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "blue"):
|
||||
|
||||
sys.stdout.write("\033[34m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "magenta"):
|
||||
|
||||
sys.stdout.write("\033[35m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "cyan"):
|
||||
|
||||
sys.stdout.write("\033[36m{}\033[0m".format(text))
|
||||
|
||||
elif(color == "white"):
|
||||
|
||||
sys.stdout.write("\033[37m{}\033[0m".format(text))
|
||||
|
||||
else:
|
||||
|
||||
sys.stdout.write("{}".format(text))
|
||||
|
||||
# ---
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user