4
1
mirror of https://github.com/pfloos/quack synced 2024-10-19 22:41:50 +02:00

fixed small bugs

This commit is contained in:
Abdallah Ammar 2024-08-29 14:05:21 +02:00
parent 0262e73353
commit 0b28512edc
2 changed files with 239 additions and 235 deletions

View File

@ -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

View File

@ -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