LCPQ
/
qp2
mirror of https://github.com/QuantumPackage/qp2.git
9
1
Fork 0
Code Releases Activity
qp2/src/ao_two_e_ints/cholesky.irp.f

296 lines
6.6 KiB
Fortran
Raw Blame History

BEGIN_PROVIDER [ integer, mini_basis_size, (128) ]
implicit none
BEGIN_DOC
! Size of the minimal basis set per element
END_DOC
mini_basis_size(1:2) = 1
mini_basis_size(3:4) = 2
mini_basis_size(5:10) = 5
mini_basis_size(11:12) = 6
mini_basis_size(13:18) = 9
mini_basis_size(19:20) = 13
mini_basis_size(21:36) = 18
mini_basis_size(37:38) = 22
mini_basis_size(39:54) = 27
mini_basis_size(55:) = 36
END_PROVIDER
BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
implicit none
BEGIN_DOC
! Number of Cholesky vectors in AO basis
END_DOC
cholesky_ao_num_guess = ao_num*ao_num !sum(mini_basis_size(int(nucl_charge(:))))
END_PROVIDER
BEGIN_PROVIDER [ integer, cholesky_ao_num ]
&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, cholesky_ao_num_guess) ]
use mmap_module
implicit none
BEGIN_DOC
! Cholesky vectors in AO basis: (ik|a):
! <ij|kl> = (ik|jl) = sum_a (ik|a).(a|jl)
END_DOC
cholesky_ao_num = cholesky_ao_num_guess
call direct_cholesky(cholesky_ao, ao_num*ao_num, cholesky_ao_num, ao_cholesky_threshold)
print *, 'Rank : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
END_PROVIDER
BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num, ao_num) ]
implicit none
BEGIN_DOC
! Transposed of the Cholesky vectors in AO basis set
END_DOC
integer :: i,j,k
do j=1,ao_num
do i=1,ao_num
do k=1,ao_num
cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k)
enddo
enddo
enddo
END_PROVIDER
subroutine direct_cholesky(L, ndim, rank, tau)
implicit none
BEGIN_DOC
! Cholesky-decomposed AOs.
!
! https://www.diva-portal.org/smash/get/diva2:396223/FULLTEXT01.pdf :
! Page 32, section 13.5
END_DOC
integer :: ndim
integer, intent(out) :: rank
double precision, intent(out) :: L(ndim, ndim)
double precision, intent(in) :: tau
double precision, parameter :: s = 1.d-2
double precision, parameter :: dscale = 1.d0
double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:)
integer*8, allocatable :: Lset(:), Dset(:), addr(:,:), LDmap(:), DLmap(:)
integer*8, allocatable :: Lset_rev(:), Dset_rev(:)
integer*8 :: i,j,k,m,p,q, qj, dj, p2, q2
integer*8 :: N, np, nq
double precision :: Dmax, Dmin, Qmax, f
double precision, external :: get_ao_two_e_integral
logical, external :: ao_two_e_integral_zero
print *, 'Entering Cholesky'
rank = 0
allocate( D(ndim), Lset(ndim), LDmap(ndim), DLmap(ndim), Dset(ndim) )
allocate( Lset_rev(ndim), Dset_rev(ndim) )
allocate( addr(3,ndim) )
! 1.
k=0
do j=1,ao_num
do i=1,ao_num
k = k+1
addr(1,k) = i
addr(2,k) = j
addr(3,k) = (i-1)*ao_num + j
enddo
enddo
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
do i=1,ndim
D(i) = get_ao_two_e_integral(addr(1,i), addr(1,i), &
addr(2,i), addr(2,i), &
ao_integrals_map)
enddo
!$OMP END PARALLEL DO
Dmax = maxval(D)
! 2.
np=0
Lset_rev = 0
do p=1,ndim
if ( dscale*dscale*Dmax*D(p) > tau*tau ) then
np = np+1
Lset(np) = p
Lset_rev(p) = np
endif
enddo
! 3.
N = 0
! 4.
i = 0
! 5.
do while (Dmax > tau)
! a.
i = i+1
! b.
Dmin = max(s*Dmax, tau)
! c.
nq=0
LDmap = 0
DLmap = 0
do p=1,np
if ( D(Lset(p)) > Dmin ) then
nq = nq+1
Dset(nq) = Lset(p)
Dset_rev(Dset(nq)) = nq
LDmap(p) = nq
DLmap(nq) = p
endif
enddo
! d., e.
allocate(Delta(np,nq), Ltmp_p(max(np,1),max(N,1)), Ltmp_q(max(nq,1),max(N,1)))
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q)
!$OMP DO
do k=1,N
do p=1,np
Ltmp_p(p,k) = L(Lset(p),k)
enddo
do q=1,nq
Ltmp_q(q,k) = L(Dset(q),k)
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,8)
do m=1,nq
do k=1, nq
! Apply only to (k,m) pairs both in Dset
p = DLmap(k)
q = Lset_rev(addr(3,Dset(k)))
if ((0 < q).and.(q < p)) cycle
if (ao_two_e_integral_zero( addr(1,Dset(k)), addr(1,Dset(m)), &
addr(2,Dset(k)), addr(2,Dset(m)) ) ) then
Delta(p,m) = 0.d0
else
Delta(p,m) = get_ao_two_e_integral( addr(1,Dset(k)), addr(1,Dset(m)), &
addr(2,Dset(k)), addr(2,Dset(m)), ao_integrals_map)
endif
Delta(q,m) = Delta(p,m)
enddo
do k=1,np
! Apply only to (k,m) pairs where k is not in Dset
if (LDmap(k) /= 0) cycle
q = Lset_rev(addr(3,Lset(k)))
if ((0 < q).and.(q < k)) cycle
if (ao_two_e_integral_zero( addr(1,Lset(k)), addr(1,Dset(m)), &
addr(2,Lset(k)), addr(2,Dset(m)) ) ) then
Delta(k,m) = 0.d0
else
Delta(k,m) = get_ao_two_e_integral( addr(1,Lset(k)), addr(1,Dset(m)), &
addr(2,Lset(k)), addr(2,Dset(m)), ao_integrals_map)
endif
Delta(q,m) = Delta(k,m)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm('N','T', int(np,4), int(nq,4), int(N,4), -1.d0, &
Ltmp_p, int(np,4), Ltmp_q, int(nq,4), 1.d0, Delta, int(np,4))
! f.
Qmax = D(Dset(1))
do q=1,nq
Qmax = max(Qmax, D(Dset(q)))
enddo
! g.
do j=1,nq
if ( (Qmax <= Dmin).or.(N+j > ndim) ) exit
! i.
rank = N+j
! ii.
do dj=1,nq
qj = Dset(dj)
if (D(qj) == Qmax) then
exit
endif
enddo
L(:, rank) = 0.d0
! iii.
f = 1.d0/dsqrt(Qmax)
!$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared)
!$OMP DO
do p=1,np
Ltmp_p(p,1) = Delta(p,dj) * f
L(Lset(p), rank) = Ltmp_p(p,1)
D(Lset(p)) = D(Lset(p)) - Ltmp_p(p,1) * Ltmp_p(p,1)
enddo
!$OMP END DO
!$OMP DO
do q=1,nq
Ltmp_q(q,1) = L(Dset(q), rank)
enddo
!$OMP END DO
! iv.
!$OMP DO SCHEDULE(static)
do m=1, nq
do k=1, np
Delta(k,m) = Delta(k,m) - Ltmp_p(k,1) * Ltmp_q(m,1)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
Qmax = D(Dset(1))
do q=1,np
Qmax = max(Qmax, D(Lset(q)))
enddo
enddo
print *, Qmax
deallocate(Delta, Ltmp_p, Ltmp_q)
! i.
N = N+j
! j.
Dmax = D(Lset(1))
do p=1,np
Dmax = max(Dmax, D(Lset(p)))
enddo
np=0
Lset_rev = 0
do p=1,ndim
if ( dscale*dscale*Dmax*D(p) > tau*tau ) then
np = np+1
Lset(np) = p
Lset_rev(p) = np
endif
enddo
enddo
end
View Git Blame Copy Permalink