diff --git a/src/ao_basis/spherical_to_cartesian.irp.f b/src/ao_basis/spherical_to_cartesian.irp.f index 91588a52..3564c71d 100644 --- a/src/ao_basis/spherical_to_cartesian.irp.f +++ b/src/ao_basis/spherical_to_cartesian.irp.f @@ -846,70 +846,70 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_cart_num,ao_sphe_num)] - implicit none - BEGIN_DOC -! Coefficients to go from cartesian to spherical coordinates in the current -! basis set +!BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_cart_num,ao_sphe_num)] +! implicit none +! BEGIN_DOC +!!Coefficients to go from cartesian to spherical coordinates in the current +!!basis set ! -! S_cart^-1 - END_DOC - integer :: i - integer, external :: ao_power_index - integer :: ibegin,j,k - integer :: prev, ao_sphe_count - prev = 0 - ao_cart_to_sphe_coef(:,:) = 0.d0 +!! S_cart^-1 +! END_DOC +! integer :: i +! integer, external :: ao_power_index +! integer :: ibegin,j,k +! integer :: prev, ao_sphe_count +! prev = 0 +! ao_cart_to_sphe_coef(:,:) = 0.d0 - if (ao_cartesian) then - ! Identity matrix - do i=1,ao_sphe_num - ao_cart_to_sphe_coef(i,i) = 1.d0 - enddo +! if (ao_cartesian) then +! ! Identity matrix +! do i=1,ao_sphe_num +! ao_cart_to_sphe_coef(i,i) = 1.d0 +! enddo - else - ! Assume order provided by ao_power_index - i = 1 - ao_sphe_count = 0 - do while (i <= ao_num) - select case ( ao_l(i) ) - case (0) - ao_sphe_count += 1 - ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0 - i += 1 - BEGIN_TEMPLATE - case ($SHELL) - if (ao_power(i,1) == $SHELL) then - do k=1,size(cart_to_sphe_$SHELL,2) - do j=1,size(cart_to_sphe_$SHELL,1) - ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k) - enddo - enddo - i += size(cart_to_sphe_$SHELL,1) - ao_sphe_count += size(cart_to_sphe_$SHELL,2) - endif - SUBST [ SHELL ] - 1;; - 2;; - 3;; - 4;; - 5;; - 6;; - 7;; - 8;; - 9;; - END_TEMPLATE - case default - stop 'Error in ao_cart_to_sphe : angular momentum too high' - end select - enddo +! else +! ! Assume order provided by ao_power_index +! i = 1 +! ao_sphe_count = 0 +! do while (i <= ao_num) +! select case ( ao_l(i) ) +! case (0) +! ao_sphe_count += 1 +! ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0 +! i += 1 +! BEGIN_TEMPLATE +! case ($SHELL) +! if (ao_power(i,1) == $SHELL) then +! do k=1,size(cart_to_sphe_$SHELL,2) +! do j=1,size(cart_to_sphe_$SHELL,1) +! ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k) +! enddo +! enddo +! i += size(cart_to_sphe_$SHELL,1) +! ao_sphe_count += size(cart_to_sphe_$SHELL,2) +! endif +! SUBST [ SHELL ] +! 1;; +! 2;; +! 3;; +! 4;; +! 5;; +! 6;; +! 7;; +! 8;; +! 9;; +! END_TEMPLATE +! case default +! stop 'Error in ao_cart_to_sphe : angular momentum too high' +! end select +! enddo - endif +! endif - if (ao_sphe_count /= ao_sphe_num) then - call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num") - endif -END_PROVIDER +! if (ao_sphe_count /= ao_sphe_num) then +! call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num") +! endif +!ND_PROVIDER diff --git a/src/ao_cart_two_e_ints/cholesky.irp.f b/src/ao_cart_two_e_ints/cholesky.irp.f index d565b810..be89c766 100644 --- a/src/ao_cart_two_e_ints/cholesky.irp.f +++ b/src/ao_cart_two_e_ints/cholesky.irp.f @@ -12,7 +12,7 @@ double precision function get_ao_cart$_erf_integ_chol(i,j,k,l) end -BEGIN_PROVIDER [ double precision, cholesky_ao_cart_transp, (cholesky_ao$_erf_cart_num, ao_cart_num, ao_cart_num) ] +BEGIN_PROVIDER [ double precision, cholesky_ao_cart$_erf_transp, (cholesky_ao$_erf_cart_num, ao_cart_num, ao_cart_num) ] implicit none BEGIN_DOC ! Transposed of the Cholesky vectors in AO basis set @@ -21,7 +21,7 @@ BEGIN_PROVIDER [ double precision, cholesky_ao_cart_transp, (cholesky_ao$_erf_ca do j=1,ao_cart_num do i=1,ao_cart_num do k=1,cholesky_ao_cart$_erf_num - cholesky_ao_cart$_erf_transp(k,i,j) = cholesky_ao(i,j,k) + cholesky_ao_cart$_erf_transp(k,i,j) = cholesky_ao_cart$_erf(i,j,k) enddo enddo enddo diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 7d6a7da4..e429047f 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -16,265 +16,6 @@ BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] print*, 'AO map initialized : ', sze END_PROVIDER -subroutine two_e_integrals_index(i,j,k,l,i1) - use map_module - implicit none - BEGIN_DOC -! Gives a unique index for i,j,k,l using permtuation symmetry. -! i <-> k, j <-> l, and (i,k) <-> (j,l) for non-periodic systems - END_DOC - integer, intent(in) :: i,j,k,l - integer(key_kind), intent(out) :: i1 - integer(key_kind) :: p,q,r,s,i2 - p = min(i,k) - r = max(i,k) - p = p+shiftr(r*r-r,1) - q = min(j,l) - s = max(j,l) - q = q+shiftr(s*s-s,1) - i1 = min(p,q) - i2 = max(p,q) - i1 = i1+shiftr(i2*i2-i2,1) -end - - - -subroutine two_e_integrals_index_reverse(i,j,k,l,i1) - use map_module - implicit none - BEGIN_DOC -! Computes the 4 indices $i,j,k,l$ from a unique index $i_1$. -! For 2 indices $i,j$ and $i \le j$, we have -! $p = i(i-1)/2 + j$. -! The key point is that because $j < i$, -! $i(i-1)/2 < p \le i(i+1)/2$. So $i$ can be found by solving -! $i^2 - i - 2p=0$. One obtains $i=1 + \sqrt{1+8p}/2$ -! and $j = p - i(i-1)/2$. -! This rule is applied 3 times. First for the symmetry of the -! pairs (i,k) and (j,l), and then for the symmetry within each pair. - END_DOC - integer, intent(out) :: i(8),j(8),k(8),l(8) - integer(key_kind), intent(in) :: i1 - integer(key_kind) :: i2,i3 - i = 0 - i2 = ceiling(0.5d0*(dsqrt(dble(shiftl(i1,3)+1))-1.d0)) - l(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i2,3)+1))-1.d0)) - i3 = i1 - shiftr(i2*i2-i2,1) - k(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i3,3)+1))-1.d0)) - j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4) - i(1) = int(i3 - shiftr(k(1)*k(1)-k(1),1),4) - - !ijkl - i(2) = i(1) !ilkj - j(2) = l(1) - k(2) = k(1) - l(2) = j(1) - - i(3) = k(1) !kjil - j(3) = j(1) - k(3) = i(1) - l(3) = l(1) - - i(4) = k(1) !klij - j(4) = l(1) - k(4) = i(1) - l(4) = j(1) - - i(5) = j(1) !jilk - j(5) = i(1) - k(5) = l(1) - l(5) = k(1) - - i(6) = j(1) !jkli - j(6) = k(1) - k(6) = l(1) - l(6) = i(1) - - i(7) = l(1) !lijk - j(7) = i(1) - k(7) = j(1) - l(7) = k(1) - - i(8) = l(1) !lkji - j(8) = k(1) - k(8) = j(1) - l(8) = i(1) - - integer :: ii, jj - do ii=2,8 - do jj=1,ii-1 - if ( (i(ii) == i(jj)).and. & - (j(ii) == j(jj)).and. & - (k(ii) == k(jj)).and. & - (l(ii) == l(jj)) ) then - i(ii) = 0 - exit - endif - enddo - enddo -! This has been tested with up to 1000 AOs, and all the reverse indices are -! correct ! We can remove the test -! do ii=1,8 -! if (i(ii) /= 0) then -! call two_e_integrals_index(i(ii),j(ii),k(ii),l(ii),i2) -! if (i1 /= i2) then -! print *, i1, i2 -! print *, i(ii), j(ii), k(ii), l(ii) -! stop 'two_e_integrals_index_reverse failed' -! endif -! endif -! enddo - - -end - - - -subroutine ao_idx2_sq(i,j,ij) - implicit none - integer, intent(in) :: i,j - integer, intent(out) :: ij - if (ij) then - ij=(i-1)*(i-1)+2*j-mod(i,2) - else - ij=i*i - endif -end - -subroutine idx2_tri_int(i,j,ij) - implicit none - integer, intent(in) :: i,j - integer, intent(out) :: ij - integer :: p,q - p = max(i,j) - q = min(i,j) - ij = q+ishft(p*p-p,-1) -end - -subroutine ao_idx2_tri_key(i,j,ij) - use map_module - implicit none - integer, intent(in) :: i,j - integer(key_kind), intent(out) :: ij - integer(key_kind) :: p,q - p = max(i,j) - q = min(i,j) - ij = q+ishft(p*p-p,-1) -end - -subroutine two_e_integrals_index_2fold(i,j,k,l,i1) - use map_module - implicit none - integer, intent(in) :: i,j,k,l - integer(key_kind), intent(out) :: i1 - integer :: ik,jl - - call ao_idx2_sq(i,k,ik) - call ao_idx2_sq(j,l,jl) - call ao_idx2_tri_key(ik,jl,i1) -end - -subroutine ao_idx2_sq_rev(i,k,ik) - BEGIN_DOC - ! reverse square compound index - END_DOC -! p = ceiling(dsqrt(dble(ik))) -! q = ceiling(0.5d0*(dble(ik)-dble((p-1)*(p-1)))) -! if (mod(ik,2)==0) then -! k=p -! i=q -! else -! i=p -! k=q -! endif - integer, intent(in) :: ik - integer, intent(out) :: i,k - integer :: pq(0:1),i1,i2 - pq(0) = ceiling(dsqrt(dble(ik))) - pq(1) = ceiling(0.5d0*(dble(ik)-dble((pq(0)-1)*(pq(0)-1)))) - i1=mod(ik,2) - i2=mod(ik+1,2) - - k=pq(i1) - i=pq(i2) -end - -subroutine ao_idx2_tri_rev_key(i,k,ik) - use map_module - BEGIN_DOC - !return i<=k - END_DOC - integer(key_kind), intent(in) :: ik - integer, intent(out) :: i,k - integer(key_kind) :: tmp_k - k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) - tmp_k = k - i = int(ik - ishft(tmp_k*tmp_k-tmp_k,-1)) -end - -subroutine idx2_tri_rev_int(i,k,ik) - BEGIN_DOC - !return i<=k - END_DOC - integer, intent(in) :: ik - integer, intent(out) :: i,k - k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) - i = int(ik - ishft(k*k-k,-1)) -end - -subroutine two_e_integrals_index_reverse_2fold(i,j,k,l,i1) - use map_module - implicit none - integer, intent(out) :: i(2),j(2),k(2),l(2) - integer(key_kind), intent(in) :: i1 - integer(key_kind) :: i0 - integer :: i2,i3 - i = 0 - call ao_idx2_tri_rev_key(i3,i2,i1) - - call ao_idx2_sq_rev(j(1),l(1),i2) - call ao_idx2_sq_rev(i(1),k(1),i3) - - !ijkl - i(2) = j(1) !jilk - j(2) = i(1) - k(2) = l(1) - l(2) = k(1) - -! i(3) = k(1) !klij complex conjugate -! j(3) = l(1) -! k(3) = i(1) -! l(3) = j(1) -! -! i(4) = l(1) !lkji complex conjugate -! j(4) = k(1) -! k(4) = j(1) -! l(4) = i(1) - - integer :: ii - if ( (i(1)==i(2)).and. & - (j(1)==j(2)).and. & - (k(1)==k(2)).and. & - (l(1)==l(2)) ) then - i(2) = 0 - endif -! This has been tested with up to 1000 AOs, and all the reverse indices are -! correct ! We can remove the test -! do ii=1,2 -! if (i(ii) /= 0) then -! call two_e_integrals_index_2fold(i(ii),j(ii),k(ii),l(ii),i0) -! if (i1 /= i0) then -! print *, i1, i0 -! print *, i(ii), j(ii), k(ii), l(ii) -! stop 'two_e_integrals_index_reverse_2fold failed' -! endif -! endif -! enddo -end - - BEGIN_PROVIDER [ integer, ao_integrals_cache_min ] diff --git a/src/ao_two_e_ints/routines_index.irp.f b/src/ao_two_e_ints/routines_index.irp.f new file mode 100644 index 00000000..665575d6 --- /dev/null +++ b/src/ao_two_e_ints/routines_index.irp.f @@ -0,0 +1,260 @@ + +subroutine two_e_integrals_index(i,j,k,l,i1) + use map_module + implicit none + BEGIN_DOC +! Gives a unique index for i,j,k,l using permtuation symmetry. +! i <-> k, j <-> l, and (i,k) <-> (j,l) for non-periodic systems + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer(key_kind) :: p,q,r,s,i2 + p = min(i,k) + r = max(i,k) + p = p+shiftr(r*r-r,1) + q = min(j,l) + s = max(j,l) + q = q+shiftr(s*s-s,1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+shiftr(i2*i2-i2,1) +end + + + +subroutine two_e_integrals_index_reverse(i,j,k,l,i1) + use map_module + implicit none + BEGIN_DOC +! Computes the 4 indices $i,j,k,l$ from a unique index $i_1$. +! For 2 indices $i,j$ and $i \le j$, we have +! $p = i(i-1)/2 + j$. +! The key point is that because $j < i$, +! $i(i-1)/2 < p \le i(i+1)/2$. So $i$ can be found by solving +! $i^2 - i - 2p=0$. One obtains $i=1 + \sqrt{1+8p}/2$ +! and $j = p - i(i-1)/2$. +! This rule is applied 3 times. First for the symmetry of the +! pairs (i,k) and (j,l), and then for the symmetry within each pair. + END_DOC + integer, intent(out) :: i(8),j(8),k(8),l(8) + integer(key_kind), intent(in) :: i1 + integer(key_kind) :: i2,i3 + i = 0 + i2 = ceiling(0.5d0*(dsqrt(dble(shiftl(i1,3)+1))-1.d0)) + l(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i2,3)+1))-1.d0)) + i3 = i1 - shiftr(i2*i2-i2,1) + k(1) = ceiling(0.5d0*(dsqrt(dble(shiftl(i3,3)+1))-1.d0)) + j(1) = int(i2 - shiftr(l(1)*l(1)-l(1),1),4) + i(1) = int(i3 - shiftr(k(1)*k(1)-k(1),1),4) + + !ijkl + i(2) = i(1) !ilkj + j(2) = l(1) + k(2) = k(1) + l(2) = j(1) + + i(3) = k(1) !kjil + j(3) = j(1) + k(3) = i(1) + l(3) = l(1) + + i(4) = k(1) !klij + j(4) = l(1) + k(4) = i(1) + l(4) = j(1) + + i(5) = j(1) !jilk + j(5) = i(1) + k(5) = l(1) + l(5) = k(1) + + i(6) = j(1) !jkli + j(6) = k(1) + k(6) = l(1) + l(6) = i(1) + + i(7) = l(1) !lijk + j(7) = i(1) + k(7) = j(1) + l(7) = k(1) + + i(8) = l(1) !lkji + j(8) = k(1) + k(8) = j(1) + l(8) = i(1) + + integer :: ii, jj + do ii=2,8 + do jj=1,ii-1 + if ( (i(ii) == i(jj)).and. & + (j(ii) == j(jj)).and. & + (k(ii) == k(jj)).and. & + (l(ii) == l(jj)) ) then + i(ii) = 0 + exit + endif + enddo + enddo +! This has been tested with up to 1000 AOs, and all the reverse indices are +! correct ! We can remove the test +! do ii=1,8 +! if (i(ii) /= 0) then +! call two_e_integrals_index(i(ii),j(ii),k(ii),l(ii),i2) +! if (i1 /= i2) then +! print *, i1, i2 +! print *, i(ii), j(ii), k(ii), l(ii) +! stop 'two_e_integrals_index_reverse failed' +! endif +! endif +! enddo + + +end + + + +subroutine ao_idx2_sq(i,j,ij) + implicit none + integer, intent(in) :: i,j + integer, intent(out) :: ij + if (ij) then + ij=(i-1)*(i-1)+2*j-mod(i,2) + else + ij=i*i + endif +end + +subroutine idx2_tri_int(i,j,ij) + implicit none + integer, intent(in) :: i,j + integer, intent(out) :: ij + integer :: p,q + p = max(i,j) + q = min(i,j) + ij = q+ishft(p*p-p,-1) +end + +subroutine ao_idx2_tri_key(i,j,ij) + use map_module + implicit none + integer, intent(in) :: i,j + integer(key_kind), intent(out) :: ij + integer(key_kind) :: p,q + p = max(i,j) + q = min(i,j) + ij = q+ishft(p*p-p,-1) +end + +subroutine two_e_integrals_index_2fold(i,j,k,l,i1) + use map_module + implicit none + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer :: ik,jl + + call ao_idx2_sq(i,k,ik) + call ao_idx2_sq(j,l,jl) + call ao_idx2_tri_key(ik,jl,i1) +end + +subroutine ao_idx2_sq_rev(i,k,ik) + BEGIN_DOC + ! reverse square compound index + END_DOC +! p = ceiling(dsqrt(dble(ik))) +! q = ceiling(0.5d0*(dble(ik)-dble((p-1)*(p-1)))) +! if (mod(ik,2)==0) then +! k=p +! i=q +! else +! i=p +! k=q +! endif + integer, intent(in) :: ik + integer, intent(out) :: i,k + integer :: pq(0:1),i1,i2 + pq(0) = ceiling(dsqrt(dble(ik))) + pq(1) = ceiling(0.5d0*(dble(ik)-dble((pq(0)-1)*(pq(0)-1)))) + i1=mod(ik,2) + i2=mod(ik+1,2) + + k=pq(i1) + i=pq(i2) +end + +subroutine ao_idx2_tri_rev_key(i,k,ik) + use map_module + BEGIN_DOC + !return i<=k + END_DOC + integer(key_kind), intent(in) :: ik + integer, intent(out) :: i,k + integer(key_kind) :: tmp_k + k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) + tmp_k = k + i = int(ik - ishft(tmp_k*tmp_k-tmp_k,-1)) +end + +subroutine idx2_tri_rev_int(i,k,ik) + BEGIN_DOC + !return i<=k + END_DOC + integer, intent(in) :: ik + integer, intent(out) :: i,k + k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0)) + i = int(ik - ishft(k*k-k,-1)) +end + +subroutine two_e_integrals_index_reverse_2fold(i,j,k,l,i1) + use map_module + implicit none + integer, intent(out) :: i(2),j(2),k(2),l(2) + integer(key_kind), intent(in) :: i1 + integer(key_kind) :: i0 + integer :: i2,i3 + i = 0 + call ao_idx2_tri_rev_key(i3,i2,i1) + + call ao_idx2_sq_rev(j(1),l(1),i2) + call ao_idx2_sq_rev(i(1),k(1),i3) + + !ijkl + i(2) = j(1) !jilk + j(2) = i(1) + k(2) = l(1) + l(2) = k(1) + +! i(3) = k(1) !klij complex conjugate +! j(3) = l(1) +! k(3) = i(1) +! l(3) = j(1) +! +! i(4) = l(1) !lkji complex conjugate +! j(4) = k(1) +! k(4) = j(1) +! l(4) = i(1) + + integer :: ii + if ( (i(1)==i(2)).and. & + (j(1)==j(2)).and. & + (k(1)==k(2)).and. & + (l(1)==l(2)) ) then + i(2) = 0 + endif +! This has been tested with up to 1000 AOs, and all the reverse indices are +! correct ! We can remove the test +! do ii=1,2 +! if (i(ii) /= 0) then +! call two_e_integrals_index_2fold(i(ii),j(ii),k(ii),l(ii),i0) +! if (i1 /= i0) then +! print *, i1, i0 +! print *, i(ii), j(ii), k(ii), l(ii) +! stop 'two_e_integrals_index_reverse_2fold failed' +! endif +! endif +! enddo +end + +