LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-25 11:57:23 +02:00
Code Releases Activity

Swapped indices in CCSD(T)

Browse Source
This commit is contained in:
Anthony Scemama 2023-05-13 13:32:52 +02:00
parent b8804f058a
commit 6289508c1e
4 changed files with 79 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
scripts/qp_import_trexio.py
View File

@ -95,14 +95,15 @@ def write_ezfio(trexio_filename, filename):
p = re.compile(r'(\d*)$') p = re.compile(r'(\d*)$')
label = [p.sub("", x).capitalize() for x in label] label = [p.sub("", x).capitalize() for x in label]
ezfio.set_nuclei_nucl_label(label) ezfio.set_nuclei_nucl_label(label)
print("OK")
else: else:
ezfio.set_nuclei_nucl_num(1) ezfio.set_nuclei_nucl_num(1)
ezfio.set_nuclei_nucl_charge([0.]) ezfio.set_nuclei_nucl_charge([0.])
ezfio.set_nuclei_nucl_coord([0.,0.,0.]) ezfio.set_nuclei_nucl_coord([0.,0.,0.])
ezfio.set_nuclei_nucl_label(["X"]) ezfio.set_nuclei_nucl_label(["X"])
print("None")
print("OK")
print("Electrons\t...\t", end=' ') print("Electrons\t...\t", end=' ')
@ -110,12 +111,12 @@ def write_ezfio(trexio_filename, filename):
try: try:
num_beta = trexio.read_electron_dn_num(trexio_file) num_beta = trexio.read_electron_dn_num(trexio_file)
except: except:
num_beta = sum(charge)//2 num_beta = int(sum(charge))//2
try: try:
num_alpha = trexio.read_electron_up_num(trexio_file) num_alpha = trexio.read_electron_up_num(trexio_file)
except: except:
num_alpha = sum(charge) - num_beta num_alpha = int(sum(charge)) - num_beta
if num_alpha == 0: if num_alpha == 0:
print("\n\nError: There are zero electrons in the TREXIO file.\n\n") print("\n\nError: There are zero electrons in the TREXIO file.\n\n")
@ -123,7 +124,7 @@ def write_ezfio(trexio_filename, filename):
ezfio.set_electrons_elec_alpha_num(num_alpha) ezfio.set_electrons_elec_alpha_num(num_alpha)
ezfio.set_electrons_elec_beta_num(num_beta) ezfio.set_electrons_elec_beta_num(num_beta)
print("OK") print(f"{num_alpha} {num_beta}")
print("Basis\t\t...\t", end=' ') print("Basis\t\t...\t", end=' ')
@ -265,6 +266,9 @@ def write_ezfio(trexio_filename, filename):
print("OK") print("OK")
else:
print("None")
# _ # _
# |\/| _ _ |_) _. _ o _ # |\/| _ _ |_) _. _ o _
@ -308,10 +312,10 @@ def write_ezfio(trexio_filename, filename):
for i in range(num_beta): for i in range(num_beta):
mo_occ[i] += 1. mo_occ[i] += 1.
ezfio.set_mo_basis_mo_occ(mo_occ) ezfio.set_mo_basis_mo_occ(mo_occ)
except:
pass
print("OK") print("OK")
except:
print("None")
print("Pseudos\t\t...\t", end=' ') print("Pseudos\t\t...\t", end=' ')
@ -391,10 +395,11 @@ def write_ezfio(trexio_filename, filename):
ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl) ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl)
ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl) ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl)
ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl) ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl)
print("OK") print("OK")
else:
print("None")

2
src/ccsd/ccsd_space_orb_sub.irp.f
View File

@ -169,7 +169,7 @@ subroutine run_ccsd_space_orb
! New ! New
print*,'Computing (T) correction...' print*,'Computing (T) correction...'
call wall_time(ta) call wall_time(ta)
call ccsd_par_t_space_v2(nO,nV,t1,t2,cc_space_f_o,cc_space_f_v & call ccsd_par_t_space_v3(nO,nV,t1,t2,cc_space_f_o,cc_space_f_v &
,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t) ,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t)
call wall_time(tb) call wall_time(tb)
print*,'Time: ',tb-ta, ' s' print*,'Time: ',tb-ta, ' s'

100
src/ccsd/ccsd_t_space_orb_abc.irp.f
View File

@ -15,8 +15,8 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
double precision, allocatable :: W_abc(:,:,:), V_abc(:,:,:) double precision, allocatable :: W_abc(:,:,:), V_abc(:,:,:)
double precision, allocatable :: W_cab(:,:,:), W_cba(:,:,:) double precision, allocatable :: W_cab(:,:,:), W_cba(:,:,:)
double precision, allocatable :: W_bca(:,:,:), V_cba(:,:,:) double precision, allocatable :: W_bca(:,:,:), V_cba(:,:,:)
double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:) double precision, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:)
double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:) double precision, allocatable :: T_voov(:,:,:,:), T_oovv(:,:,:,:)
integer :: i,j,k,l,a,b,c,d integer :: i,j,k,l,a,b,c,d
double precision :: e,ta,tb, delta, delta_abc double precision :: e,ta,tb, delta, delta_abc
@ -24,25 +24,25 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
!allocate(V(nV,nV,nV,nO,nO,nO)) !allocate(V(nV,nV,nV,nO,nO,nO))
allocate(W_abc(nO,nO,nO), V_abc(nO,nO,nO), W_cab(nO,nO,nO)) allocate(W_abc(nO,nO,nO), V_abc(nO,nO,nO), W_cab(nO,nO,nO))
allocate(W_bca(nO,nO,nO), V_cba(nO,nO,nO), W_cba(nO,nO,nO)) allocate(W_bca(nO,nO,nO), V_cba(nO,nO,nO), W_cba(nO,nO,nO))
allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO)) allocate(X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV), X_oovv(nO,nO,nV,nV))
allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO)) allocate(T_voov(nV,nO,nO,nV),T_oovv(nO,nO,nV,nV))
! Temporary arrays ! Temporary arrays
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, & !$OMP SHARED(nO,nV,T_voov,T_oovv,X_vovv,X_ooov,X_oovv, &
!$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) & !$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
!$OMP PRIVATE(a,b,c,d,i,j,k,l) & !$OMP PRIVATE(a,b,c,d,i,j,k,l) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!v_vvvo(b,a,d,i) * t2(k,j,c,d) & !v_vvvo(b,a,d,i) * t2(k,j,c,d) &
!X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) !X_vovv(d,i,b,a,i) * T_voov(d,j,c,k)
!$OMP DO collapse(3) !$OMP DO collapse(3)
do i = 1, nO do i = 1, nO
do a = 1, nV do a = 1, nV
do b = 1, nV do b = 1, nV
do d = 1, nV do d = 1, nV
X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i) X_vovv(d,i,b,a) = v_vvvo(b,a,d,i)
enddo enddo
enddo enddo
enddo enddo
@ -54,7 +54,7 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
do k = 1, nO do k = 1, nO
do c = 1, nV do c = 1, nV
do d = 1, nV do d = 1, nV
T_vvoo(d,c,k,j) = t2(k,j,c,d) T_voov(d,k,j,c) = t2(k,j,c,d)
enddo enddo
enddo enddo
enddo enddo
@ -62,14 +62,14 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
!$OMP END DO nowait !$OMP END DO nowait
!v_vooo(c,j,k,l) * t2(i,l,a,b) & !v_vooo(c,j,k,l) * t2(i,l,a,b) &
!X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) & !X_ooov(l,j,k,c) * T_oovv(l,i,a,b) &
!$OMP DO collapse(3) !$OMP DO collapse(3)
do c = 1, nV
do k = 1, nO do k = 1, nO
do j = 1, nO do j = 1, nO
do c = 1, nV
do l = 1, nO do l = 1, nO
X_ovoo(l,c,j,k) = v_vooo(c,j,k,l) X_ooov(l,j,k,c) = v_vooo(c,j,k,l)
enddo enddo
enddo enddo
enddo enddo
@ -81,35 +81,27 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
do b = 1, nV do b = 1, nV
do a = 1, nV do a = 1, nV
do l = 1, nO do l = 1, nO
T_ovvo(l,a,b,i) = t2(i,l,a,b) T_oovv(l,i,a,b) = t2(i,l,a,b)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
!$OMP END DO nowait !$OMP END DO nowait
!v_vvoo(b,c,j,k) * t1(i,a) & !X_oovv(j,k,b,c) * T1_vo(a,i) &
!X_vvoo(b,c,k,j) * T1_vo(a,i) &
!$OMP DO collapse(3) !$OMP DO collapse(3)
do j = 1, nO
do k = 1, nO
do c = 1, nV do c = 1, nV
do b = 1, nV do b = 1, nV
X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k) do j = 1, nO
do k = 1, nO
X_oovv(j,k,b,c) = v_vvoo(b,c,j,k)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
!$OMP END DO nowait !$OMP END DO nowait
!$OMP DO collapse(1)
do i = 1, nO
do a = 1, nV
T_vo(a,i) = t1(i,a)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
call wall_time(ta) call wall_time(ta)
@ -118,13 +110,13 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
do b = 1, nV do b = 1, nV
do a = 1, nV do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c) delta_abc = f_v(a) + f_v(b) + f_v(c)
call form_w_abc(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc) call form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc)
call form_w_abc(nO,nV,b,c,a,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_bca) call form_w_abc(nO,nV,b,c,a,T_voov,T_oovv,X_vovv,X_ooov,W_bca)
call form_w_abc(nO,nV,c,a,b,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_cab) call form_w_abc(nO,nV,c,a,b,T_voov,T_oovv,X_vovv,X_ooov,W_cab)
call form_w_abc(nO,nV,c,b,a,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_cba) call form_w_abc(nO,nV,c,b,a,T_voov,T_oovv,X_vovv,X_ooov,W_cba)
call form_v_abc(nO,nV,a,b,c,T_vo,X_vvoo,W_abc,V_abc) call form_v_abc(nO,nV,a,b,c,t1,X_oovv,W_abc,V_abc)
call form_v_abc(nO,nV,c,b,a,T_vo,X_vvoo,W_cba,V_cba) call form_v_abc(nO,nV,c,b,a,t1,X_oovv,W_cba,V_cba)
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(energy,nO,a,b,c,W_abc,W_cab,W_bca,V_abc,V_cba,f_o,f_v,delta_abc)& !$OMP SHARED(energy,nO,a,b,c,W_abc,W_cab,W_bca,V_abc,V_cba,f_o,f_v,delta_abc)&
!$OMP PRIVATE(i,j,k,e,delta) & !$OMP PRIVATE(i,j,k,e,delta) &
@ -154,26 +146,26 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
energy = energy / 3d0 energy = energy / 3d0
deallocate(W_abc,V_abc,W_cab,V_cba,W_bca,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo) deallocate(W_abc,V_abc,W_cab,V_cba,W_bca,X_vovv,X_ooov,T_voov,T_oovv)
!deallocate(V,W) !deallocate(V,W)
end end
subroutine form_w_abc(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc) subroutine form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc)
implicit none implicit none
integer, intent(in) :: nO,nV,a,b,c integer, intent(in) :: nO,nV,a,b,c
!double precision, intent(in) :: t2(nO,nO,nV,nV) !double precision, intent(in) :: t2(nO,nO,nV,nV)
double precision, intent(in) :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO) double precision, intent(in) :: T_voov(nV,nO,nO,nV), T_oovv(nO,nO,nV,nV)
double precision, intent(in) :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO) double precision, intent(in) :: X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV)
double precision, intent(out) :: W_abc(nO,nO,nO) double precision, intent(out) :: W_abc(nO,nO,nO)
integer :: l,i,j,k,d integer :: l,i,j,k,d
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc) & !$OMP SHARED(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc) &
!$OMP PRIVATE(i,j,k,d,l) & !$OMP PRIVATE(i,j,k,d,l) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
@ -185,23 +177,23 @@ subroutine form_w_abc(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc)
do d = 1, nV do d = 1, nV
W_abc(i,j,k) = W_abc(i,j,k) & W_abc(i,j,k) = W_abc(i,j,k) &
+ X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) & + X_vovv(d,i,b,a) * T_voov(d,k,j,c) &
+ X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) & + X_vovv(d,i,c,a) * T_voov(d,j,k,b) &
+ X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) & + X_vovv(d,k,a,c) * T_voov(d,j,i,b) &
+ X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) & + X_vovv(d,k,b,c) * T_voov(d,i,j,a) &
+ X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) & + X_vovv(d,j,c,b) * T_voov(d,i,k,a) &
+ X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i) + X_vovv(d,j,a,b) * T_voov(d,k,i,c)
enddo enddo
do l = 1, nO do l = 1, nO
W_abc(i,j,k) = W_abc(i,j,k) & W_abc(i,j,k) = W_abc(i,j,k) &
- T_ovvo(l,a,b,i) * X_ovoo(l,c,j,k) & - T_oovv(l,i,a,b) * X_ooov(l,j,k,c) &
- T_ovvo(l,a,c,i) * X_ovoo(l,b,k,j) & ! bc kj - T_oovv(l,i,a,c) * X_ooov(l,k,j,b) & ! bc kj
- T_ovvo(l,c,a,k) * X_ovoo(l,b,i,j) & ! prev ac ik - T_oovv(l,k,c,a) * X_ooov(l,i,j,b) & ! prev ac ik
- T_ovvo(l,c,b,k) * X_ovoo(l,a,j,i) & ! prev ab ij - T_oovv(l,k,c,b) * X_ooov(l,j,i,a) & ! prev ab ij
- T_ovvo(l,b,c,j) * X_ovoo(l,a,k,i) & ! prev bc kj - T_oovv(l,j,b,c) * X_ooov(l,k,i,a) & ! prev bc kj
- T_ovvo(l,b,a,j) * X_ovoo(l,c,i,k) ! prev ac ik - T_oovv(l,j,b,a) * X_ooov(l,i,k,c) ! prev ac ik
enddo enddo
enddo enddo
@ -216,21 +208,21 @@ end
! V_abc ! V_abc
subroutine form_v_abc(nO,nV,a,b,c,T_vo,X_vvoo,W,V) subroutine form_v_abc(nO,nV,a,b,c,T_ov,X_oovv,W,V)
implicit none implicit none
integer, intent(in) :: nO,nV,a,b,c integer, intent(in) :: nO,nV,a,b,c
!double precision, intent(in) :: t1(nO,nV) !double precision, intent(in) :: t1(nO,nV)
double precision, intent(in) :: T_vo(nV,nO) double precision, intent(in) :: T_ov(nO,nV)
double precision, intent(in) :: X_vvoo(nV,nV,nO,nO) double precision, intent(in) :: X_oovv(nO,nO,nV,nV)
double precision, intent(in) :: W(nO,nO,nO) double precision, intent(in) :: W(nO,nO,nO)
double precision, intent(out) :: V(nO,nO,nO) double precision, intent(out) :: V(nO,nO,nO)
integer :: i,j,k integer :: i,j,k
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(nO,nV,a,b,c,T_vo,X_vvoo,W,V) & !$OMP SHARED(nO,nV,a,b,c,T_ov,X_oovv,W,V) &
!$OMP PRIVATE(i,j,k) & !$OMP PRIVATE(i,j,k) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO collapse(2) !$OMP DO collapse(2)
@ -239,9 +231,9 @@ implicit none
do i = 1, nO do i = 1, nO
!V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) & !V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
V(i,j,k) = W(i,j,k) & V(i,j,k) = W(i,j,k) &
+ X_vvoo(b,c,k,j) * T_vo(a,i) & + X_oovv(j,k,b,c) * T_ov(i,a) &
+ X_vvoo(a,c,k,i) * T_vo(b,j) & + X_oovv(i,k,a,c) * T_ov(j,b) &
+ X_vvoo(a,b,j,i) * T_vo(c,k) + X_oovv(i,j,a,b) * T_ov(k,c)
enddo enddo
enddo enddo
enddo enddo

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

@ -1825,39 +1825,37 @@ subroutine pivoted_cholesky( A, rank, tol, ndim, U)
! !
integer :: ndim integer :: ndim
integer, intent(inout) :: rank integer, intent(inout) :: rank
double precision, dimension(ndim, ndim), intent(inout) :: A double precision, intent(inout) :: A(ndim, ndim)
double precision, dimension(ndim, rank), intent(out) :: U double precision, intent(out) :: U(ndim, rank)
double precision, intent(in) :: tol double precision, intent(in) :: tol
integer, dimension(:), allocatable :: piv integer, dimension(:), allocatable :: piv
double precision, dimension(:), allocatable :: work double precision, dimension(:), allocatable :: work
character, parameter :: uplo = "U" character, parameter :: uplo = "U"
integer :: N, LDA integer :: LDA
integer :: info integer :: info
integer :: k, l, rank0 integer :: k, l, rank0
external :: dpstrf
rank0 = rank rank0 = rank
N = size(A, dim=1) LDA = ndim
LDA = N allocate(piv(ndim))
allocate(piv(N)) allocate(work(2*ndim))
allocate(work(2*N)) call dpstrf(uplo, ndim, A, LDA, piv, rank, tol, work, info)
call dpstrf(uplo, N, A, LDA, piv, rank, tol, work, info)
if (rank > rank0) then if (rank > rank0) then
print *, 'Bug: rank > rank0 in pivoted cholesky. Increase rank before calling' print *, 'Bug: rank > rank0 in pivoted cholesky. Increase rank before calling'
stop stop
end if end if
do k = 1, N do k = 1, ndim
A(k+1:, k) = 0.00D+0 A(k+1:ndim, k) = 0.00D+0
end do end do
! TODO: It should be possible to use only one vector of size (1:rank) as a buffer ! TODO: It should be possible to use only one vector of size (1:rank) as a buffer
! to do the swapping in-place ! to do the swapping in-place
U(:,:) = 0.00D+0 U(:,:) = 0.00D+0
do k = 1, N do k = 1, ndim
l = piv(k) l = piv(k)
U(l, :) = A(1:rank, k) U(l, 1:rank) = A(1:rank, k)
end do end do
end subroutine pivoted_cholesky end subroutine pivoted_cholesky