PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-11 13:53:49 +02:00
Code Releases Activity

Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

Browse Source
This commit is contained in:
eginer 2023-08-21 17:29:17 +02:00
commit 99d856570a
2 changed files with 106 additions and 111 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

164
src/ao_two_e_ints/cholesky.irp.f
View File

@ -23,38 +23,38 @@ END_PROVIDER
! !
! Last dimension of cholesky_ao is cholesky_ao_num ! Last dimension of cholesky_ao is cholesky_ao_num
END_DOC END_DOC
integer :: rank, ndim integer :: rank, ndim
double precision :: tau double precision :: tau
double precision, pointer :: L(:,:), L_old(:,:) double precision, pointer :: L(:,:), L_old(:,:)
double precision :: s double precision :: s
double precision, parameter :: dscale = 1.d0 double precision, parameter :: dscale = 1.d0
double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:) double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:)
integer, allocatable :: Lset(:), Dset(:), addr(:,:) integer, allocatable :: Lset(:), Dset(:), addr(:,:)
logical, allocatable :: computed(:) logical, allocatable :: computed(:)
integer :: i,j,k,m,p,q, qj, dj, p2, q2 integer :: i,j,k,m,p,q, qj, dj, p2, q2
integer :: N, np, nq integer :: N, np, nq
double precision :: Dmax, Dmin, Qmax, f double precision :: Dmax, Dmin, Qmax, f
double precision, external :: get_ao_two_e_integral double precision, external :: get_ao_two_e_integral
logical, external :: ao_two_e_integral_zero logical, external :: ao_two_e_integral_zero
double precision, external :: ao_two_e_integral double precision, external :: ao_two_e_integral
integer :: block_size, iblock, ierr integer :: block_size, iblock, ierr
double precision :: mem double precision :: mem
double precision, external :: memory_of_double, memory_of_int double precision, external :: memory_of_double, memory_of_int
integer, external :: getUnitAndOpen integer, external :: getUnitAndOpen
integer :: iunit integer :: iunit
ndim = ao_num*ao_num ndim = ao_num*ao_num
deallocate(cholesky_ao) deallocate(cholesky_ao)
if (read_ao_cholesky) then if (read_ao_cholesky) then
print *, 'Reading Cholesky vectors from disk...' print *, 'Reading Cholesky vectors from disk...'
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R') iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R')
@ -63,11 +63,11 @@ END_PROVIDER
read(iunit) cholesky_ao read(iunit) cholesky_ao
close(iunit) close(iunit)
cholesky_ao_num = rank cholesky_ao_num = rank
else else
PROVIDE nucl_coord PROVIDE nucl_coord
if (do_direct_integrals) then if (do_direct_integrals) then
if (ao_two_e_integral(1,1,1,1) < huge(1.d0)) then if (ao_two_e_integral(1,1,1,1) < huge(1.d0)) then
! Trigger providers inside ao_two_e_integral ! Trigger providers inside ao_two_e_integral
@ -76,16 +76,16 @@ END_PROVIDER
else else
PROVIDE ao_two_e_integrals_in_map PROVIDE ao_two_e_integrals_in_map
endif endif
tau = ao_cholesky_threshold tau = ao_cholesky_threshold
mem = 6.d0 * memory_of_double(ndim) + 6.d0 * memory_of_int(ndim) mem = 6.d0 * memory_of_double(ndim) + 6.d0 * memory_of_int(ndim)
call check_mem(mem, irp_here) call check_mem(mem, irp_here)
call print_memory_usage() call print_memory_usage()
allocate(L(ndim,1)) allocate(L(ndim,1))
print *, '' print *, ''
print *, 'Cholesky decomposition of AO integrals' print *, 'Cholesky decomposition of AO integrals'
print *, '======================================' print *, '======================================'
@ -93,13 +93,13 @@ END_PROVIDER
print *, '============ =============' print *, '============ ============='
print *, ' Rank Threshold' print *, ' Rank Threshold'
print *, '============ =============' print *, '============ ============='
rank = 0 rank = 0
allocate( D(ndim), Lset(ndim), Dset(ndim) ) allocate( D(ndim), Lset(ndim), Dset(ndim) )
allocate( addr(3,ndim) ) allocate( addr(3,ndim) )
! 1. ! 1.
k=0 k=0
do j=1,ao_num do j=1,ao_num
@ -110,7 +110,7 @@ END_PROVIDER
addr(3,k) = (i-1)*ao_num + j addr(3,k) = (i-1)*ao_num + j
enddo enddo
enddo enddo
if (do_direct_integrals) then if (do_direct_integrals) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided)
do i=1,ndim do i=1,ndim
@ -127,9 +127,9 @@ END_PROVIDER
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif endif
Dmax = maxval(D) Dmax = maxval(D)
! 2. ! 2.
np=0 np=0
do p=1,ndim do p=1,ndim
@ -138,26 +138,26 @@ END_PROVIDER
Lset(np) = p Lset(np) = p
endif endif
enddo enddo
! 3. ! 3.
N = 0 N = 0
! 4. ! 4.
i = 0 i = 0
! 5. ! 5.
do while ( (Dmax > tau).and.(rank < ndim) ) do while ( (Dmax > tau).and.(rank < ndim) )
! a. ! a.
i = i+1 i = i+1
s = 0.01d0 s = 0.01d0
! Inrease s until the arrays fit in memory ! Inrease s until the arrays fit in memory
do while (.True.) do while (.True.)
! b. ! b.
Dmin = max(s*Dmax,tau) Dmin = max(s*Dmax,tau)
! c. ! c.
nq=0 nq=0
do p=1,np do p=1,np
@ -166,30 +166,30 @@ END_PROVIDER
Dset(nq) = Lset(p) Dset(nq) = Lset(p)
endif endif
enddo enddo
call total_memory(mem) call total_memory(mem)
mem = mem & mem = mem &
+ np*memory_of_double(nq) &! Delta(np,nq) + np*memory_of_double(nq) &! Delta(np,nq)
+ (rank+nq)* memory_of_double(ndim) &! L(ndim,rank+nq) + (rank+nq)* memory_of_double(ndim) &! L(ndim,rank+nq)
+ (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size) + (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
if (mem > qp_max_mem) then if (mem > qp_max_mem) then
s = s*2.d0 s = s*2.d0
else else
exit exit
endif endif
if ((s > 1.d0).or.(nq == 0)) then if ((s > 1.d0).or.(nq == 0)) then
call print_memory_usage() call print_memory_usage()
print *, 'Not enough memory. Reduce cholesky threshold' print *, 'Not enough memory. Reduce cholesky threshold'
stop -1 stop -1
endif endif
enddo enddo
! d., e. ! d., e.
block_size = max(N,24) block_size = max(N,24)
L_old => L L_old => L
allocate(L(ndim,rank+nq), stat=ierr) allocate(L(ndim,rank+nq), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
@ -197,48 +197,52 @@ END_PROVIDER
print *, irp_here, ': allocation failed : (L(ndim,rank+nq))' print *, irp_here, ': allocation failed : (L(ndim,rank+nq))'
stop -1 stop -1
endif endif
!$OMP PARALLEL DO PRIVATE(k) !$OMP PARALLEL DO PRIVATE(k,j)
do k=1,rank do k=1,rank
L(:,k) = L_old(:,k) do j=1,ndim
L(j,k) = L_old(j,k)
enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
deallocate(L_old) deallocate(L_old)
allocate(Delta(np,nq), stat=ierr) allocate(Delta(np,nq), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
print *, irp_here, ': allocation failed : (Delta(np,nq))' print *, irp_here, ': allocation failed : (Delta(np,nq))'
stop -1 stop -1
endif endif
allocate(Ltmp_p(np,block_size), stat=ierr) allocate(Ltmp_p(np,block_size), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
print *, irp_here, ': allocation failed : (Ltmp_p(np,block_size))' print *, irp_here, ': allocation failed : (Ltmp_p(np,block_size))'
stop -1 stop -1
endif endif
allocate(Ltmp_q(nq,block_size), stat=ierr) allocate(Ltmp_q(nq,block_size), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
print *, irp_here, ': allocation failed : (Ltmp_q(nq,block_size))' print *, irp_here, ': allocation failed : (Ltmp_q(nq,block_size))'
stop -1 stop -1
endif endif
allocate(computed(nq)) allocate(computed(nq))
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j) !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j)
!$OMP DO !$OMP DO
do q=1,nq do q=1,nq
Delta(:,q) = 0.d0 do j=1,np
Delta(j,q) = 0.d0
enddo
computed(q) = .False. computed(q) = .False.
enddo enddo
!$OMP ENDDO NOWAIT !$OMP ENDDO NOWAIT
!$OMP DO !$OMP DO
do k=1,N do k=1,N
do p=1,np do p=1,np
@ -249,36 +253,36 @@ END_PROVIDER
enddo enddo
enddo enddo
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
!$OMP BARRIER !$OMP BARRIER
!$OMP END PARALLEL !$OMP END PARALLEL
if (N>0) then if (N>0) then
call dgemm('N','T', np, nq, N, -1.d0, & call dgemm('N','T', np, nq, N, -1.d0, &
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
endif endif
! f. ! f.
Qmax = D(Dset(1)) Qmax = D(Dset(1))
do q=1,nq do q=1,nq
Qmax = max(Qmax, D(Dset(q))) Qmax = max(Qmax, D(Dset(q)))
enddo enddo
! g. ! g.
iblock = 0 iblock = 0
do j=1,nq do j=1,nq
if ( (Qmax <= Dmin).or.(N+j > ndim) ) exit if ( (Qmax <= Dmin).or.(N+j > ndim) ) exit
! i. ! i.
rank = N+j rank = N+j
if (iblock == block_size) then if (iblock == block_size) then
call dgemm('N','T',np,nq,block_size,-1.d0, & call dgemm('N','T',np,nq,block_size,-1.d0, &
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
iblock = 0 iblock = 0
endif endif
! ii. ! ii.
do dj=1,nq do dj=1,nq
qj = Dset(dj) qj = Dset(dj)
@ -286,9 +290,9 @@ END_PROVIDER
exit exit
endif endif
enddo enddo
L(1:ndim, rank) = 0.d0 L(1:ndim, rank) = 0.d0
if (.not.computed(dj)) then if (.not.computed(dj)) then
m = dj m = dj
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided) !$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided)
@ -314,16 +318,16 @@ END_PROVIDER
do p=1,np do p=1,np
Ltmp_p(p,iblock) = Delta(p,dj) Ltmp_p(p,iblock) = Delta(p,dj)
enddo enddo
! iv. ! iv.
if (iblock > 1) then if (iblock > 1) then
call dgemv('N', np, iblock-1, -1.d0, Ltmp_p, np, Ltmp_q(dj,1), nq, 1.d0,& call dgemv('N', np, iblock-1, -1.d0, Ltmp_p, np, Ltmp_q(dj,1), nq, 1.d0,&
Ltmp_p(1,iblock), 1) Ltmp_p(1,iblock), 1)
endif endif
! iii. ! iii.
f = 1.d0/dsqrt(Qmax) f = 1.d0/dsqrt(Qmax)
!$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared) !$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared)
!$OMP DO !$OMP DO
do p=1,np do p=1,np
@ -332,38 +336,38 @@ END_PROVIDER
D(Lset(p)) = D(Lset(p)) - Ltmp_p(p,iblock) * Ltmp_p(p,iblock) D(Lset(p)) = D(Lset(p)) - Ltmp_p(p,iblock) * Ltmp_p(p,iblock)
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP DO !$OMP DO
do q=1,nq do q=1,nq
Ltmp_q(q,iblock) = L(Dset(q), rank) Ltmp_q(q,iblock) = L(Dset(q), rank)
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
Qmax = D(Dset(1)) Qmax = D(Dset(1))
do q=1,nq do q=1,nq
Qmax = max(Qmax, D(Dset(q))) Qmax = max(Qmax, D(Dset(q)))
enddo enddo
enddo enddo
print '(I10, 4X, ES12.3)', rank, Qmax print '(I10, 4X, ES12.3)', rank, Qmax
deallocate(computed) deallocate(computed)
deallocate(Delta) deallocate(Delta)
deallocate(Ltmp_p) deallocate(Ltmp_p)
deallocate(Ltmp_q) deallocate(Ltmp_q)
! i. ! i.
N = rank N = rank
! j. ! j.
Dmax = D(Lset(1)) Dmax = D(Lset(1))
do p=1,np do p=1,np
Dmax = max(Dmax, D(Lset(p))) Dmax = max(Dmax, D(Lset(p)))
enddo enddo
np=0 np=0
do p=1,ndim do p=1,ndim
if ( dscale*dscale*Dmax*D(p) > tau*tau ) then if ( dscale*dscale*Dmax*D(p) > tau*tau ) then
@ -371,9 +375,9 @@ END_PROVIDER
Lset(np) = p Lset(np) = p
endif endif
enddo enddo
enddo enddo
allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
@ -387,10 +391,10 @@ END_PROVIDER
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
deallocate(L) deallocate(L)
cholesky_ao_num = rank cholesky_ao_num = rank
print *, '============ =============' print *, '============ ============='
print *, '' print *, ''
if (write_ao_cholesky) then if (write_ao_cholesky) then
print *, 'Writing Cholesky vectors to disk...' print *, 'Writing Cholesky vectors to disk...'
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W') iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W')
@ -401,9 +405,9 @@ END_PROVIDER
endif endif
endif endif
print *, 'Rank : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)' print *, 'Rank : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
print *, '' print *, ''
END_PROVIDER END_PROVIDER

53
src/ccsd/ccsd_space_orb_sub_chol.irp.f
View File

@ -244,7 +244,7 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
allocate(W_oovo(nO,nO,nV,nO)) allocate(W_oovo(nO,nO,nV,nO))
!$omp parallel & !$omp parallel &
!$omp shared(nO,nV,cc_space_v_vooo,W_oovo) & !$omp shared(nO,nV,cc_space_v_oovo,W_oovo) &
!$omp private(u,a,i,j) & !$omp private(u,a,i,j) &
!$omp default(none) !$omp default(none)
do u = 1, nO do u = 1, nO
@ -252,7 +252,8 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
do a = 1, nV do a = 1, nV
do j = 1, nO do j = 1, nO
do i = 1, nO do i = 1, nO
W_oovo(i,j,a,u) = 2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i) ! W_oovo(i,j,a,u) = 2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)
W_oovo(i,j,a,u) = 2d0 * cc_space_v_oovo(i,j,a,u) - cc_space_v_oovo(j,i,a,u)
enddo enddo
enddo enddo
enddo enddo
@ -514,10 +515,10 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
cc_space_v_vv_chol(1,1,gam), cholesky_mo_num, & cc_space_v_vv_chol(1,1,gam), cholesky_mo_num, &
0.d0, tmpB1, nV*block_size) 0.d0, tmpB1, nV*block_size)
call dgemm('T','N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, 1.d0, & call dgemm('T','N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, &
cc_space_v_vv_chol(1,1,iblock), cholesky_mo_num, & 1.d0, cc_space_v_vv_chol(1,1,iblock), cholesky_mo_num, &
tmp_cc2, cholesky_mo_num, & tmp_cc2, cholesky_mo_num, &
1.d0, tmpB1, nV*block_size) 1.d0, tmpB1, nV*block_size)
do beta = iblock, min(nV, iblock+block_size-1) do beta = iblock, min(nV, iblock+block_size-1)
do b = 1, nV do b = 1, nV
@ -1071,37 +1072,27 @@ subroutine compute_g_occ_chol(nO,nV,t1,t2,H_oo,g_occ)
double precision, intent(in) :: t2(nO, nO, nV, nV) double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: g_occ(nO, nO) double precision, intent(out) :: g_occ(nO, nO)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam g_occ = H_oo
call dgemm('N','N',nO,nO,nV, & call dgemm('N','N',nO,nO,nV, &
1d0, t1, size(t1,1), & 1d0, t1, size(t1,1), &
cc_space_f_vo, size(cc_space_f_vo,1), & cc_space_f_vo, size(cc_space_f_vo,1), &
0d0, g_occ, size(g_occ,1)) 1d0, g_occ, size(g_occ,1))
!$omp parallel & double precision, allocatable :: X(:)
!$omp shared(nO,nV,g_occ,H_oo, cc_space_v_ovoo,t1) & allocate(X(cholesky_mo_num))
!$omp private(i,j,a,u) & call dgemv('N',cholesky_mo_num,nO*nV,2.d0, &
!$omp default(none) cc_space_v_ov_chol, cholesky_mo_num, &
!$omp do t1, 1, 0.d0, X, 1)
do i = 1, nO
do u = 1, nO
g_occ(u,i) = g_occ(u,i) + H_oo(u,i)
enddo
enddo
!$omp end do
!$omp do call dgemv('T',cholesky_mo_num,nO*nO,1.d0, &
do i = 1, nO cc_space_v_oo_chol, cholesky_mo_num, &
do j = 1, nO X, 1, 1.d0, g_occ, 1)
do a = 1, nV deallocate(X)
do u = 1, nO
g_occ(u,i) = g_occ(u,i) + (2d0 * cc_space_v_ovoo(u,a,i,j) - cc_space_v_ovoo(u,a,j,i)) * t1(j,a) call dgemv('T',nO*nV,nO*nO,-1.d0, &
enddo cc_space_v_ovoo, nO*nV, &
enddo t1, 1, 1.d0, g_occ, 1)
enddo
enddo
!$omp end do
!$omp end parallel
end end