mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-19 04:22:32 +01:00
Implementing stochastic (T)
This commit is contained in:
parent
adbc61a218
commit
e3c0df574e
@ -169,7 +169,9 @@ 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_v3(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)
|
||||||
|
call ccsd_par_t_space_stoch(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'
|
||||||
|
|||||||
@ -19,14 +19,13 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
|
|||||||
double precision, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:)
|
double precision, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:)
|
||||||
double precision, allocatable :: T_voov(:,:,:,:), T_oovv(:,:,:,:)
|
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, x1, x2, x3
|
double precision :: e,ta,tb
|
||||||
|
|
||||||
|
call set_multiple_levels_omp(.False.)
|
||||||
|
|
||||||
allocate(X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV), X_oovv(nO,nO,nV,nV))
|
allocate(X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV), X_oovv(nO,nO,nV,nV))
|
||||||
allocate(T_voov(nV,nO,nO,nV),T_oovv(nO,nO,nV,nV))
|
allocate(T_voov(nV,nO,nO,nV),T_oovv(nO,nO,nV,nV))
|
||||||
|
|
||||||
call set_multiple_levels_omp(.False.)
|
|
||||||
|
|
||||||
! Temporary arrays
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP SHARED(nO,nV,T_voov,T_oovv,X_vovv,X_ooov,X_oovv, &
|
!$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) &
|
||||||
@ -103,12 +102,13 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
|
|||||||
|
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
energy = 0d0
|
double precision, external :: ccsd_t_task_aba
|
||||||
|
double precision, external :: ccsd_t_task_abc
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP PRIVATE(a,b,c,x1) &
|
!$OMP PRIVATE(a,b,c,e) &
|
||||||
!$OMP PRIVATE(W_abc, W_cab, W_bca, W_bac, W_cba, W_acb, &
|
!$OMP PRIVATE(W_abc, W_cab, W_bca, W_bac, W_cba, W_acb, &
|
||||||
!$OMP V_abc, V_cab, V_bca, V_bac, V_cba, V_acb ) &
|
!$OMP V_abc, V_cab, V_bca, V_bac, V_cba, V_acb ) &
|
||||||
!$OMP PRIVATE(i,j,k,e,delta,delta_abc) &
|
|
||||||
!$OMP DEFAULT(SHARED)
|
!$OMP DEFAULT(SHARED)
|
||||||
allocate( W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO), &
|
allocate( W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO), &
|
||||||
W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO), &
|
W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO), &
|
||||||
@ -119,46 +119,18 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
|
|||||||
do a = 1, nV
|
do a = 1, nV
|
||||||
do b = a+1, nV
|
do b = a+1, nV
|
||||||
do c = b+1, nV
|
do c = b+1, nV
|
||||||
delta_abc = f_v(a) + f_v(b) + f_v(c)
|
e = e + ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov,V_abc, &
|
||||||
call form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
|
V_acb,V_bac,V_bca,V_cab,V_cba,W_abc,W_acb,W_bac, &
|
||||||
call form_v_abc(nO,nV,a,b,c,t1,X_oovv,W_abc,V_abc,W_cba,V_cba,W_bca,V_bca,W_cab,V_cab,W_bac,V_bac,W_acb,V_acb)
|
W_bca,W_cab,W_cba,X_ooov,X_oovv,X_vovv,f_o,f_v)
|
||||||
do k = 1, nO
|
|
||||||
do j = 1, nO
|
|
||||||
do i = 1, nO
|
|
||||||
delta = 1.d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
|
|
||||||
e = e + delta * ( &
|
|
||||||
(4d0 * (W_abc(i,j,k) - W_cba(i,j,k)) + &
|
|
||||||
W_bca(i,j,k) - W_bac(i,j,k) + &
|
|
||||||
W_cab(i,j,k) - W_acb(i,j,k) ) * (V_abc(i,j,k) - V_cba(i,j,k)) + &
|
|
||||||
(4d0 * (W_acb(i,j,k) - W_bca(i,j,k)) + &
|
|
||||||
W_cba(i,j,k) - W_cab(i,j,k) + &
|
|
||||||
W_bac(i,j,k) - W_abc(i,j,k) ) * (V_acb(i,j,k) - V_bca(i,j,k)) + &
|
|
||||||
(4d0 * (W_bac(i,j,k) - W_cab(i,j,k)) + &
|
|
||||||
W_acb(i,j,k) - W_abc(i,j,k) + &
|
|
||||||
W_cba(i,j,k) - W_bca(i,j,k) ) * (V_bac(i,j,k) - V_cab(i,j,k)) )
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
c = a
|
e = e + ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov,V_abc, &
|
||||||
do b = 1, nV
|
V_acb,V_bac,V_bca,V_cab,V_cba,W_abc,W_acb,W_bac, &
|
||||||
if (b == c) cycle
|
W_bca,W_cab,W_cba,X_ooov,X_oovv,X_vovv,f_o,f_v)
|
||||||
delta_abc = f_v(a) + f_v(b) + f_v(c)
|
|
||||||
call form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
|
e = e + ccsd_t_task_aba(b,a,nO,nV,t1,T_oovv,T_voov,V_abc, &
|
||||||
call form_v_abc(nO,nV,a,b,c,t1,X_oovv,W_abc,V_abc,W_cba,V_cba,W_bca,V_bca,W_cab,V_cab,W_bac,V_bac,W_acb,V_acb)
|
V_acb,V_bac,V_bca,V_cab,V_cba,W_abc,W_acb,W_bac, &
|
||||||
do k = 1, nO
|
W_bca,W_cab,W_cba,X_ooov,X_oovv,X_vovv,f_o,f_v)
|
||||||
do j = 1, nO
|
|
||||||
do i = 1, nO
|
|
||||||
delta = 1.0d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
|
|
||||||
e = e + delta * ( &
|
|
||||||
(4d0 * W_abc(i,j,k) + W_bca(i,j,k) + W_cab(i,j,k)) * (V_abc(i,j,k) - V_cba(i,j,k)) + &
|
|
||||||
(4d0 * W_acb(i,j,k) + W_cba(i,j,k) + W_bac(i,j,k)) * (V_acb(i,j,k) - V_bca(i,j,k)) + &
|
|
||||||
(4d0 * W_bac(i,j,k) + W_acb(i,j,k) + W_cba(i,j,k)) * (V_bac(i,j,k) - V_cab(i,j,k)) )
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO NOWAIT
|
!$OMP END DO NOWAIT
|
||||||
@ -178,6 +150,91 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
|
|||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
double precision function ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov,&
|
||||||
|
V_abc,V_acb,V_bac,V_bca,V_cab,V_cba, &
|
||||||
|
W_abc,W_acb,W_bac,W_bca,W_cab,W_cba, &
|
||||||
|
X_ooov,X_oovv,X_vovv,f_o,f_v) result(e)
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: nO,nV,a,b,c
|
||||||
|
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
|
||||||
|
double precision, intent(in) :: X_oovv(nO,nO,nV,nV)
|
||||||
|
double precision, intent(in) :: T_voov(nV,nO,nO,nV), T_oovv(nO,nO,nV,nV)
|
||||||
|
double precision, intent(in) :: X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV)
|
||||||
|
double precision, intent(in) :: W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO)
|
||||||
|
double precision, intent(in) :: W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO)
|
||||||
|
double precision, intent(in) :: V_abc(nO,nO,nO), V_cab(nO,nO,nO), V_bca(nO,nO,nO)
|
||||||
|
double precision, intent(in) :: V_bac(nO,nO,nO), V_cba(nO,nO,nO), V_acb(nO,nO,nO)
|
||||||
|
|
||||||
|
double precision :: delta, delta_abc
|
||||||
|
integer :: i,j,k
|
||||||
|
|
||||||
|
delta_abc = f_v(a) + f_v(b) + f_v(c)
|
||||||
|
e = 0.d0
|
||||||
|
|
||||||
|
call form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
|
||||||
|
|
||||||
|
call form_v_abc(nO,nV,a,b,c,t1,X_oovv,W_abc,V_abc,W_cba,V_cba,W_bca,V_bca,W_cab,V_cab,W_bac,V_bac,W_acb,V_acb)
|
||||||
|
|
||||||
|
do k = 1, nO
|
||||||
|
do j = 1, nO
|
||||||
|
do i = 1, nO
|
||||||
|
delta = 1.d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
|
||||||
|
e = e + delta * ( &
|
||||||
|
(4d0 * (W_abc(i,j,k) - W_cba(i,j,k)) + &
|
||||||
|
W_bca(i,j,k) - W_bac(i,j,k) + &
|
||||||
|
W_cab(i,j,k) - W_acb(i,j,k) ) * (V_abc(i,j,k) - V_cba(i,j,k)) +&
|
||||||
|
(4d0 * (W_acb(i,j,k) - W_bca(i,j,k)) + &
|
||||||
|
W_cba(i,j,k) - W_cab(i,j,k) + &
|
||||||
|
W_bac(i,j,k) - W_abc(i,j,k) ) * (V_acb(i,j,k) - V_bca(i,j,k)) +&
|
||||||
|
(4d0 * (W_bac(i,j,k) - W_cab(i,j,k)) + &
|
||||||
|
W_acb(i,j,k) - W_abc(i,j,k) + &
|
||||||
|
W_cba(i,j,k) - W_bca(i,j,k) ) * (V_bac(i,j,k) - V_cab(i,j,k)) )
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
double precision function ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov,&
|
||||||
|
V_abc,V_acb,V_bac,V_bca,V_cab,V_cba, &
|
||||||
|
W_abc,W_acb,W_bac,W_bca,W_cab,W_cba, &
|
||||||
|
X_ooov,X_oovv,X_vovv,f_o,f_v) result(e)
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: nO,nV,a,b
|
||||||
|
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
|
||||||
|
double precision, intent(in) :: X_oovv(nO,nO,nV,nV)
|
||||||
|
double precision, intent(in) :: T_voov(nV,nO,nO,nV), T_oovv(nO,nO,nV,nV)
|
||||||
|
double precision, intent(in) :: X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV)
|
||||||
|
double precision, intent(in) :: W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO)
|
||||||
|
double precision, intent(in) :: W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO)
|
||||||
|
double precision, intent(in) :: V_abc(nO,nO,nO), V_cab(nO,nO,nO), V_bca(nO,nO,nO)
|
||||||
|
double precision, intent(in) :: V_bac(nO,nO,nO), V_cba(nO,nO,nO), V_acb(nO,nO,nO)
|
||||||
|
|
||||||
|
double precision :: delta, delta_abc
|
||||||
|
integer :: i,j,k
|
||||||
|
|
||||||
|
delta_abc = f_v(a) + f_v(b) + f_v(a)
|
||||||
|
e = 0.d0
|
||||||
|
|
||||||
|
call form_w_abc(nO,nV,a,b,a,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
|
||||||
|
|
||||||
|
call form_v_abc(nO,nV,a,b,a,t1,X_oovv,W_abc,V_abc,W_cba,V_cba,W_bca,V_bca,W_cab,V_cab,W_bac,V_bac,W_acb,V_acb)
|
||||||
|
|
||||||
|
do k = 1, nO
|
||||||
|
do j = 1, nO
|
||||||
|
do i = 1, nO
|
||||||
|
delta = 1.d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
|
||||||
|
e = e + delta * ( &
|
||||||
|
(4d0 * W_abc(i,j,k) + W_bca(i,j,k) + W_cab(i,j,k)) * (V_abc(i,j,k) - V_cba(i,j,k)) + &
|
||||||
|
(4d0 * W_acb(i,j,k) + W_cba(i,j,k) + W_bac(i,j,k)) * (V_acb(i,j,k) - V_bca(i,j,k)) + &
|
||||||
|
(4d0 * W_bac(i,j,k) + W_acb(i,j,k) + W_cba(i,j,k)) * (V_bac(i,j,k) - V_cab(i,j,k)) )
|
||||||
|
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
subroutine form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
|
subroutine form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|||||||
320
src/ccsd/ccsd_t_space_orb_stoch.irp.f
Normal file
320
src/ccsd/ccsd_t_space_orb_stoch.irp.f
Normal file
@ -0,0 +1,320 @@
|
|||||||
|
! Main
|
||||||
|
|
||||||
|
subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer, intent(in) :: nO,nV
|
||||||
|
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
|
||||||
|
double precision, intent(in) :: t2(nO,nO,nV,nV)
|
||||||
|
double precision, intent(in) :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
|
||||||
|
double precision, intent(out) :: energy
|
||||||
|
|
||||||
|
double precision, allocatable :: W(:,:,:,:,:,:)
|
||||||
|
double precision, allocatable :: V(:,:,:,:,:,:)
|
||||||
|
double precision, allocatable :: W_abc(:,:,:), W_cab(:,:,:), W_bca(:,:,:)
|
||||||
|
double precision, allocatable :: W_bac(:,:,:), W_cba(:,:,:), W_acb(:,:,:)
|
||||||
|
double precision, allocatable :: V_abc(:,:,:), V_cab(:,:,:), V_bca(:,:,:)
|
||||||
|
double precision, allocatable :: V_bac(:,:,:), V_cba(:,:,:), V_acb(:,:,:)
|
||||||
|
double precision, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:)
|
||||||
|
double precision, allocatable :: T_voov(:,:,:,:), T_oovv(:,:,:,:)
|
||||||
|
integer :: i,j,k,l,a,b,c,d
|
||||||
|
double precision :: e,ta,tb
|
||||||
|
|
||||||
|
call set_multiple_levels_omp(.False.)
|
||||||
|
|
||||||
|
allocate(X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV), X_oovv(nO,nO,nV,nV))
|
||||||
|
allocate(T_voov(nV,nO,nO,nV),T_oovv(nO,nO,nV,nV))
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP SHARED(nO,nV,T_voov,T_oovv,X_vovv,X_ooov,X_oovv, &
|
||||||
|
!$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
|
||||||
|
!$OMP PRIVATE(a,b,c,d,i,j,k,l) &
|
||||||
|
!$OMP DEFAULT(NONE)
|
||||||
|
|
||||||
|
!v_vvvo(b,a,d,i) * t2(k,j,c,d) &
|
||||||
|
!X_vovv(d,i,b,a,i) * T_voov(d,j,c,k)
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do a = 1, nV
|
||||||
|
do b = 1, nV
|
||||||
|
do i = 1, nO
|
||||||
|
do d = 1, nV
|
||||||
|
X_vovv(d,i,b,a) = v_vvvo(b,a,d,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO nowait
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do c = 1, nV
|
||||||
|
do j = 1, nO
|
||||||
|
do k = 1, nO
|
||||||
|
do d = 1, nV
|
||||||
|
T_voov(d,k,j,c) = t2(k,j,c,d)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO nowait
|
||||||
|
|
||||||
|
!v_vooo(c,j,k,l) * t2(i,l,a,b) &
|
||||||
|
!X_ooov(l,j,k,c) * T_oovv(l,i,a,b) &
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do c = 1, nV
|
||||||
|
do k = 1, nO
|
||||||
|
do j = 1, nO
|
||||||
|
do l = 1, nO
|
||||||
|
X_ooov(l,j,k,c) = v_vooo(c,j,k,l)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO nowait
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do b = 1, nV
|
||||||
|
do a = 1, nV
|
||||||
|
do i = 1, nO
|
||||||
|
do l = 1, nO
|
||||||
|
T_oovv(l,i,a,b) = t2(i,l,a,b)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO nowait
|
||||||
|
|
||||||
|
!X_oovv(j,k,b,c) * T1_vo(a,i) &
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do c = 1, nV
|
||||||
|
do b = 1, nV
|
||||||
|
do k = 1, nO
|
||||||
|
do j = 1, nO
|
||||||
|
X_oovv(j,k,b,c) = v_vvoo(b,c,j,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO nowait
|
||||||
|
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
double precision, external :: ccsd_t_task_aba
|
||||||
|
double precision, external :: ccsd_t_task_abc
|
||||||
|
|
||||||
|
double precision, allocatable :: memo(:), Pabc(:), waccu(:)
|
||||||
|
logical , allocatable :: computed(:)
|
||||||
|
integer*2 , allocatable :: abc(:,:)
|
||||||
|
integer*8 :: Nabc, i8
|
||||||
|
integer*8, allocatable :: iorder(:)
|
||||||
|
double precision :: eocc
|
||||||
|
double precision :: Pabc_norm, sum_w
|
||||||
|
|
||||||
|
|
||||||
|
! Prepare table of triplets (a,b,c)
|
||||||
|
|
||||||
|
Nabc = (int(nV,8) * int(nV+1,8) * int(nV+2,8))/6_8 - nV
|
||||||
|
allocate (memo(Nabc), computed(Nabc), Pabc(Nabc), waccu(0:Nabc))
|
||||||
|
allocate (abc(4,Nabc), iorder(Nabc))
|
||||||
|
|
||||||
|
! eocc = 3.d0/dble(nO) * sum(f_o(1:nO))
|
||||||
|
memo(:) = 0.d0
|
||||||
|
computed(:) = .False.
|
||||||
|
Nabc = 0_8
|
||||||
|
do a = 1, nV
|
||||||
|
do b = a+1, nV
|
||||||
|
do c = b+1, nV
|
||||||
|
Nabc = Nabc + 1_8
|
||||||
|
! Pabc(Nabc) = 1.d0/((f_v(a) + f_v(b) + f_v(c))*(f_v(a)*f_v(b)*f_v(c))**(1.d0/2.d0))
|
||||||
|
! Pabc(Nabc) = 1.d0/((f_v(a) + f_v(b) + f_v(c))**2)
|
||||||
|
Pabc(Nabc) = 1.d0/(f_v(a) + f_v(b) + f_v(c))
|
||||||
|
abc(1,Nabc) = a
|
||||||
|
abc(2,Nabc) = b
|
||||||
|
abc(3,Nabc) = c
|
||||||
|
enddo
|
||||||
|
|
||||||
|
Nabc = Nabc + 1_8
|
||||||
|
abc(1,Nabc) = a
|
||||||
|
abc(2,Nabc) = b
|
||||||
|
abc(3,Nabc) = a
|
||||||
|
! Pabc(Nabc) = 1.d0/((f_v(a) + f_v(b) + f_v(a))*(f_v(a)*f_v(b)*f_v(a))**(1.d0/2.d0))
|
||||||
|
! Pabc(Nabc) = 1.d0/((f_v(a) + f_v(b) + f_v(a))**2)
|
||||||
|
Pabc(Nabc) = 1.d0/(2.d0*f_v(a) + f_v(b))
|
||||||
|
|
||||||
|
Nabc = Nabc + 1_8
|
||||||
|
abc(1,Nabc) = b
|
||||||
|
abc(2,Nabc) = a
|
||||||
|
abc(3,Nabc) = b
|
||||||
|
! Pabc(Nabc) = 1.d0/((f_v(a) + f_v(b) + f_v(b))*(f_v(b)*f_v(a)*f_v(b))**(1.d0/2.d0))
|
||||||
|
! Pabc(Nabc) = 1.d0/((f_v(a) + f_v(b) + f_v(b))**2)
|
||||||
|
Pabc(Nabc) = 1.d0/(f_v(a) + 2.d0*f_v(b))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do i8=1,Nabc
|
||||||
|
iorder(i8) = i8
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! Sort triplets in decreasing Pabc
|
||||||
|
call dsort_big(Pabc, iorder, Nabc)
|
||||||
|
|
||||||
|
! Normalize
|
||||||
|
Pabc_norm = 0.d0
|
||||||
|
do i8=Nabc,1,-1
|
||||||
|
Pabc_norm = Pabc_norm + Pabc(i8)
|
||||||
|
enddo
|
||||||
|
Pabc_norm = 1.d0/Pabc_norm
|
||||||
|
do i8=Nabc,1,-1
|
||||||
|
Pabc(i8) = Pabc(i8) * Pabc_norm
|
||||||
|
enddo
|
||||||
|
|
||||||
|
call i8set_order_big(abc, iorder, Nabc)
|
||||||
|
|
||||||
|
|
||||||
|
! Cumulative distribution for sampling
|
||||||
|
waccu(Nabc) = 0.d0
|
||||||
|
sum_w = 0.d0
|
||||||
|
do i8=Nabc-1,1,-1
|
||||||
|
waccu(i8) = waccu(i8+1) - Pabc(i8)
|
||||||
|
enddo
|
||||||
|
waccu(:) = waccu(:) + 1.d0
|
||||||
|
waccu(0) = 0.d0
|
||||||
|
|
||||||
|
Pabc(:) = 1.d0/Pabc(:) * (1.d0/3.d0)
|
||||||
|
|
||||||
|
logical :: converged
|
||||||
|
double precision :: ET, ET2, eta, variance, average, error, sample
|
||||||
|
integer*8 :: isample, ieta, Ncomputed
|
||||||
|
integer*8, external :: find_sample
|
||||||
|
|
||||||
|
allocate( W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO), &
|
||||||
|
W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO), &
|
||||||
|
V_abc(nO,nO,nO), V_cab(nO,nO,nO), V_bca(nO,nO,nO), &
|
||||||
|
V_bac(nO,nO,nO), V_cba(nO,nO,nO), V_acb(nO,nO,nO) )
|
||||||
|
|
||||||
|
converged = .False.
|
||||||
|
ET = 0.d0
|
||||||
|
ET2 = 0.d0
|
||||||
|
Ncomputed = 0_8
|
||||||
|
isample = 0_8
|
||||||
|
|
||||||
|
average = 0.d0
|
||||||
|
variance = 0.d0
|
||||||
|
double precision :: t00, t01
|
||||||
|
call wall_time(t00)
|
||||||
|
! do ieta=1,Nabc
|
||||||
|
do while (.not.converged)
|
||||||
|
call random_number(eta)
|
||||||
|
! eta = eta/dble(1000)
|
||||||
|
! do k=0,1000-1
|
||||||
|
! ieta = find_sample(eta+dble(k)/dble(1000),waccu,Nabc)
|
||||||
|
ieta = find_sample(eta,waccu,Nabc)
|
||||||
|
isample = isample+1_8
|
||||||
|
|
||||||
|
if (.not.computed(ieta)) then
|
||||||
|
a = abc(1,ieta)
|
||||||
|
b = abc(2,ieta)
|
||||||
|
c = abc(3,ieta)
|
||||||
|
if (a/=c) then
|
||||||
|
memo(ieta) = ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov,V_abc, &
|
||||||
|
V_acb,V_bac,V_bca,V_cab,V_cba,W_abc,W_acb,W_bac, &
|
||||||
|
W_bca,W_cab,W_cba,X_ooov,X_oovv,X_vovv,f_o,f_v)
|
||||||
|
else
|
||||||
|
memo(ieta) = ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov,V_abc, &
|
||||||
|
V_acb,V_bac,V_bca,V_cab,V_cba,W_abc,W_acb,W_bac, &
|
||||||
|
W_bca,W_cab,W_cba,X_ooov,X_oovv,X_vovv,f_o,f_v)
|
||||||
|
endif
|
||||||
|
computed(ieta) = .True.
|
||||||
|
Ncomputed += 1_8
|
||||||
|
call wall_time(t01)
|
||||||
|
if (t01-t00 > 1.d0) then
|
||||||
|
t00 = t01
|
||||||
|
print *, average, dsqrt(variance/dble(isample)), real(Ncomputed)/real(Nabc), real(isample)/real(Nabc)
|
||||||
|
endif
|
||||||
|
! print *, memo(ieta), Pabc(ieta), memo(ieta) * Pabc(ieta)
|
||||||
|
endif
|
||||||
|
sample = memo(ieta) * Pabc(ieta)
|
||||||
|
ET = ET + sample
|
||||||
|
ET2 = ET2 + sample*sample
|
||||||
|
average = ET/dble(isample)
|
||||||
|
variance = ET2/dble(isample) - average*average
|
||||||
|
converged = (Ncomputed >= (Nabc*90_8)/100_8) .or. (isample>=1000*Nabc)
|
||||||
|
! enddo
|
||||||
|
enddo
|
||||||
|
print *, average, dsqrt(variance/dble(isample)), real(Ncomputed)/real(Nabc), real(isample)/real(Nabc)
|
||||||
|
energy = average
|
||||||
|
|
||||||
|
! !$OMP PARALLEL &
|
||||||
|
! !$OMP PRIVATE(a,b,c,e) &
|
||||||
|
! !$OMP PRIVATE(W_abc, W_cab, W_bca, W_bac, W_cba, W_acb, &
|
||||||
|
! !$OMP V_abc, V_cab, V_bca, V_bac, V_cba, V_acb ) &
|
||||||
|
! !$OMP DEFAULT(SHARED)
|
||||||
|
! allocate( W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO), &
|
||||||
|
! W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO), &
|
||||||
|
! V_abc(nO,nO,nO), V_cab(nO,nO,nO), V_bca(nO,nO,nO), &
|
||||||
|
! V_bac(nO,nO,nO), V_cba(nO,nO,nO), V_acb(nO,nO,nO) )
|
||||||
|
! e = 0d0
|
||||||
|
! !$OMP DO SCHEDULE(dynamic)
|
||||||
|
! do a = 1, nV
|
||||||
|
! do b = a+1, nV
|
||||||
|
! do c = b+1, nV
|
||||||
|
! e = e + ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov,V_abc, &
|
||||||
|
! V_acb,V_bac,V_bca,V_cab,V_cba,W_abc,W_acb,W_bac, &
|
||||||
|
! W_bca,W_cab,W_cba,X_ooov,X_oovv,X_vovv,f_o,f_v)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
!
|
||||||
|
! do b = 1, nV
|
||||||
|
! if (b == a) cycle
|
||||||
|
! e = e + ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov,V_abc, &
|
||||||
|
! V_acb,V_bac,V_bca,V_cab,V_cba,W_abc,W_acb,W_bac, &
|
||||||
|
! W_bca,W_cab,W_cba,X_ooov,X_oovv,X_vovv,f_o,f_v)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! !$OMP END DO NOWAIT
|
||||||
|
!
|
||||||
|
! !$OMP CRITICAL
|
||||||
|
! energy = energy + e
|
||||||
|
! !$OMP END CRITICAL
|
||||||
|
!
|
||||||
|
! deallocate(W_abc, W_cab, W_bca, W_bac, W_cba, W_acb, &
|
||||||
|
! V_abc, V_cab, V_bca, V_bac, V_cba, V_acb )
|
||||||
|
!
|
||||||
|
! !$OMP END PARALLEL
|
||||||
|
|
||||||
|
deallocate(X_vovv,X_ooov,T_voov,T_oovv)
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
integer*8 function find_sample(v, w, n)
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Finds sample v in weights w
|
||||||
|
END_DOC
|
||||||
|
integer*8, intent(in) :: n
|
||||||
|
double precision, intent(in) :: v, w(0:n)
|
||||||
|
integer*8 :: i,l,r
|
||||||
|
|
||||||
|
l=0
|
||||||
|
r=n
|
||||||
|
|
||||||
|
do while(r-l > 1)
|
||||||
|
i = shiftr(r+l,1)
|
||||||
|
if(w(i) < v) then
|
||||||
|
l = i
|
||||||
|
else
|
||||||
|
r = i
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
i = r
|
||||||
|
do r=i+1,n
|
||||||
|
if (w(r) /= w(i)) then
|
||||||
|
exit
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
find_sample = r-1
|
||||||
|
end function
|
||||||
|
|
||||||
Loading…
Reference in New Issue
Block a user