From 0bcbdb44c5108c30a82cecb2d437f0273dcaa60c Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 15 Apr 2025 15:43:22 +0200 Subject: [PATCH 1/2] minor modif in ao_spherical and ao_ortho_canonical --- src/ao_one_e_ints/ao_ortho_canonical.irp.f | 42 +------ src/ao_one_e_ints/ao_spherical.irp.f | 130 +++++++++++++++++++++ 2 files changed, 131 insertions(+), 41 deletions(-) create mode 100644 src/ao_one_e_ints/ao_spherical.irp.f diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f index e015c89e..654633b3 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f @@ -1,5 +1,4 @@ - BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_num)] -&BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_normalization, (ao_num)] + BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_normalization, (ao_num)] implicit none BEGIN_DOC ! Coefficients to go from cartesian to spherical coordinates in the current @@ -12,39 +11,23 @@ integer :: ibegin,j,k integer :: prev, ao_sphe_count prev = 0 - ao_cart_to_sphe_coef(:,:) = 0.d0 ao_cart_to_sphe_normalization(:) = 1.d0 - 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 ao_cart_to_sphe_normalization(i) = 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 do j=1,size(cart_to_sphe_$SHELL,1) ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j) enddo i += size(cart_to_sphe_$SHELL,1) - ao_sphe_count += size(cart_to_sphe_$SHELL,2) endif SUBST [ SHELL ] 1;; @@ -62,35 +45,12 @@ end select enddo - 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 -BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_sphe_num,ao_sphe_num) ] - implicit none - BEGIN_DOC - ! T^T . S . T - END_DOC - double precision, allocatable :: S(:,:) - allocate (S(ao_sphe_num,ao_num)) - - call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, & - ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), & - ao_overlap,size(ao_overlap,1), 0.d0, & - S, size(S,1)) - - call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, & - S, size(S,1), & - ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), 0.d0, & - ao_cart_to_sphe_overlap,size(ao_cart_to_sphe_overlap,1)) - - deallocate(S) - -END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_sphe_num,ao_num) ] implicit none BEGIN_DOC diff --git a/src/ao_one_e_ints/ao_spherical.irp.f b/src/ao_one_e_ints/ao_spherical.irp.f new file mode 100644 index 00000000..30110e14 --- /dev/null +++ b/src/ao_one_e_ints/ao_spherical.irp.f @@ -0,0 +1,130 @@ +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_sphe_num)] + implicit none + BEGIN_DOC + ! Coefficients to go from current cartesian AO basis set to spherical AO basis set + ! + ! ao_cart_to_sphe_coef(i,j) = coefficient of the i-th cartesian |AO| on the j-th spherical |AO| + ! + ! :math:`\chi^s_{\nu} = \sum_{\nu}^{N_{\text{cart}}} B^{sc}_{\mu\nu} \chi^c_{\mu}` + ! + ! where :math:`\chi^s_{\nu}` is an element of the spherical AO basis, + ! :math:`\chi^c_{\mu}` is an element of the cartesian AO basis, + ! and :math:`B^{sc}_{\mu\nu}` is an change of basis matrix + END_DOC + integer :: row,col,k,j + ! + ao_cart_to_sphe_coef(:,:) = 0.d0 + row = 1 + col = 1 + do while (row <= ao_num) + ! Select case based on azimuthal quantum number of i-th AO orbitals + select case ( ao_l(row) ) + case (0) + ! S orbital + ao_cart_to_sphe_coef(row,col) = 1.d0 + row += 1 + col += 1 + BEGIN_TEMPLATE + case ($SHELL) + ! P,D,F,... orbitals + do k=1,size(cart_to_sphe_$SHELL,2) + do j=1,size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_coef(row+j-1,col+k-1) = cart_to_sphe_$SHELL(j,k) + enddo + enddo + row += size(cart_to_sphe_$SHELL,1) + col += size(cart_to_sphe_$SHELL,2) + 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 +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef_transp, (ao_sphe_num,ao_num)] + implicit none + BEGIN_DOC + ! Transpose of :c:data:`ao_cart_to_sphe_coef` + END_DOC + integer :: i,j + do i = 1, ao_num + do j = 1, ao_sphe_num + ao_cart_to_sphe_coef_transp(j,i) = ao_cart_to_sphe_coef(i,j) + enddo + enddo +END_PROVIDER + + +subroutine ao_cart_to_ao_sphe(A_ao_cart,LDA_ao_cart,A_ao_sphe,LDA_ao_sphe) + implicit none + BEGIN_DOC + ! Transform matrix A from the |AO| cartesian basis to the |AO| spherical basis + ! + ! :math:`(B^{sc})^T.A^c.B^{sc}` + ! + ! where :math:`B^{sc}` is :c:data:`ao_cart_to_sphe_coef`, + ! the matrix of coefficients from the cartesian AO basis to spherical one, + ! and :math:`B^{sc}` is :c:data:`ao_cart_to_sphe_coef_transp`, its transpose. + END_DOC + integer, intent(in) :: LDA_ao_cart,LDA_ao_sphe + double precision, intent(in) :: A_ao_cart(LDA_ao_cart,ao_num) + double precision, intent(out) :: A_ao_sphe(LDA_ao_sphe,ao_sphe_num) + double precision, allocatable :: T(:,:) + ! + allocate (T(ao_num,ao_sphe_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + call dgemm('N','N', ao_num, ao_sphe_num, ao_num, & + 1.d0, A_ao_cart, LDA_ao_cart, & + ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & + 0.d0, T, size(T,1)) + ! Notice that for the following dgemm we could have used + ! ao_cart_to_sphe_coef_transp, but instead we transposed with the 'T' argument + call dgemm('T','N', ao_sphe_num, ao_sphe_num, ao_num, & + 1.d0, ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & + T, ao_num, & + 0.d0, A_ao_sphe, size(A_ao_sphe,1)) + ! + ! call restore_symmetry(mo_num,mo_num,A_ao_sphe,size(A_ao_sphe,1),1.d-15) + deallocate(T) +end + + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_sphe_num,ao_sphe_num) ] + implicit none + BEGIN_DOC + ! |AO| overlap matrix in the spherical basis set obtained as + ! + ! :math:`(B^{sc})^T.S^c.B^{sc}` + ! + ! where :math:`S^c` is the overlap matrix in the cartesian AO basis + END_DOC + ! + call ao_cart_to_ao_sphe(ao_overlap,ao_num,ao_cart_to_sphe_overlap,ao_sphe_num) + ! +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap_inv, (ao_sphe_num,ao_sphe_num) ] + implicit none + BEGIN_DOC + ! Inverse of :c:data:`ao_cart_to_sphe_overlap` + END_DOC + ! + call get_pseudo_inverse( & + ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & + ao_sphe_num,ao_sphe_num, & + ao_cart_to_sphe_overlap_inv, size(ao_cart_to_sphe_overlap_inv,1), & + lin_dep_cutoff) +END_PROVIDER From 5f1ab387cbdcce371aea3deb011734382086c306 Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 16 Apr 2025 11:21:10 +0200 Subject: [PATCH 2/2] added the identity matrix in ao_one_e_ints/ao_spherical.irp.f --- src/ao_one_e_ints/ao_spherical.irp.f | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ao_one_e_ints/ao_spherical.irp.f b/src/ao_one_e_ints/ao_spherical.irp.f index 30110e14..b82b8f33 100644 --- a/src/ao_one_e_ints/ao_spherical.irp.f +++ b/src/ao_one_e_ints/ao_spherical.irp.f @@ -13,6 +13,13 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_sphe_num)] END_DOC integer :: row,col,k,j ! + if (ao_cartesian) then + ! Identity matrix + integer :: i + do i=1,ao_sphe_num + ao_cart_to_sphe_coef(i,i) = 1.d0 + enddo + else ao_cart_to_sphe_coef(:,:) = 0.d0 row = 1 col = 1 @@ -49,6 +56,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_sphe_num)] stop 'Error in ao_cart_to_sphe : angular momentum too high' end select enddo + endif END_PROVIDER