diff --git a/.gitignore b/.gitignore index 288e294..899b091 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,5 @@ *.o *. +__pycache__ + +.ninja_deps diff --git a/src/GT/GTpp_excitation_density.f90 b/src/GT/GTpp_excitation_density.f90 index 05f9c2d..6bfe164 100644 --- a/src/GT/GTpp_excitation_density.f90 +++ b/src/GT/GTpp_excitation_density.f90 @@ -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 diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90 index ee92faf..acd44d3 100644 --- a/src/GW/evRGW.f90 +++ b/src/GW/evRGW.f90 @@ -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 diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 64d533c..582e28f 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -72,38 +72,43 @@ subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA) 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)) - !! 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 .eq. 0) .or. (nVV .eq. 0)) then + ! 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) - 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) + else - 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 + 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 - 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 + endif end if diff --git a/src/test/check_test_value.f90 b/src/test/check_test_value.f90 index 2ee9bde..8828f82 100755 --- a/src/test/check_test_value.f90 +++ b/src/test/check_test_value.f90 @@ -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 diff --git a/src/test/dump_test_value.f90 b/src/test/dump_test_value.f90 index ea00afe..e903e28 100755 --- a/src/test/dump_test_value.f90 +++ b/src/test/dump_test_value.f90 @@ -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 diff --git a/src/test/init_test.f90 b/src/test/init_test.f90 index 602ba54..b5ef295 100755 --- a/src/test/init_test.f90 +++ b/src/test/init_test.f90 @@ -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 diff --git a/src/test/stop_test.f90 b/src/test/stop_test.f90 index d41e2b0..185cbc5 100755 --- a/src/test/stop_test.f90 +++ b/src/test/stop_test.f90 @@ -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 diff --git a/src/utils/non_sym_diag.f90 b/src/utils/non_sym_diag.f90 index 04279ac..0f44099 100644 --- a/src/utils/non_sym_diag.f90 +++ b/src/utils/non_sym_diag.f90 @@ -3,7 +3,7 @@ subroutine diagonalize_nonsym_matrix(N, A, L, e_re, thr_d, thr_nd, thr_deg, imp_bio, verbose) - ! Diagonalize a non-symmetric matrix + ! Diagonalize a non-symmetric matrix A ! ! Output ! right-eigenvectors are saved in A @@ -278,7 +278,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) do i = 1, m if(S(i,i) <= 0.d0) then - print *, ' overap negative' + print *, ' negative overlap !' print *, i, S(i,i) exit endif diff --git a/test/export_tobench.py b/test/export_tobench.py new file mode 100644 index 0000000..526716b --- /dev/null +++ b/test/export_tobench.py @@ -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) + diff --git a/tests/.gitignore b/tests/.gitignore new file mode 100644 index 0000000..afb251f --- /dev/null +++ b/tests/.gitignore @@ -0,0 +1,7 @@ + +FeatherBench.db +FeatherBench.json + +*.xyz +work + diff --git a/tests/balance_bench.py b/tests/balance_bench.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/create_database.py b/tests/create_database.py new file mode 100644 index 0000000..14ed3f6 --- /dev/null +++ b/tests/create_database.py @@ -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) + + + diff --git a/tests/feather_bench.py b/tests/feather_bench.py new file mode 100644 index 0000000..b6d8f26 --- /dev/null +++ b/tests/feather_bench.py @@ -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 +] + + + + + diff --git a/tests/inp/methods.RHF b/tests/inp/methods.RHF new file mode 100644 index 0000000..2ddb2bb --- /dev/null +++ b/tests/inp/methods.RHF @@ -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 diff --git a/tests/inp/options.RHF b/tests/inp/options.RHF new file mode 100644 index 0000000..92084cd --- /dev/null +++ b/tests/inp/options.RHF @@ -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 diff --git a/tests/lunch_bench.py b/tests/lunch_bench.py new file mode 100644 index 0000000..02d6db5 --- /dev/null +++ b/tests/lunch_bench.py @@ -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() + diff --git a/tests/molecule.py b/tests/molecule.py new file mode 100644 index 0000000..b5f33eb --- /dev/null +++ b/tests/molecule.py @@ -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!") + + diff --git a/tests/titan_bench.py b/tests/titan_bench.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/utils.py b/tests/utils.py new file mode 100644 index 0000000..a5f090f --- /dev/null +++ b/tests/utils.py @@ -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)) + +# --- + +