From 6e6a8ac82ac3294ec255dfaca08cfece8bc5e5d5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 Feb 2016 10:56:39 +0100 Subject: [PATCH] Corrected bug in spherical guess MOs --- src/Determinants/connected_to_ref.irp.f | 5 +- src/MO_Basis/ao_ortho_canonical.irp.f | 150 ++++++++++++++++++------ src/Utils/LinearAlgebra.irp.f | 1 - 3 files changed, 116 insertions(+), 40 deletions(-) diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index dc7698b5..b93b18b6 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -95,12 +95,13 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) exit endif enddo - i += 1 - if (i > N_det) then + if (i >= N_det) then return endif + i += 1 + !DIR$ FORCEINLINE do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref) if ( (key(1,1) /= psi_det_sorted_bit(1,1,i)).or. & diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index 62a72149..1916bfe4 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -1,3 +1,87 @@ + BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num_align,ao_num)] +&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ] + implicit none + BEGIN_DOC +! matrix of the coefficients of the mos generated by the +! orthonormalization by the S^{-1/2} canonical transformation of the aos +! ao_cart_to_sphe_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital + END_DOC + integer :: i + integer, external :: ao_power_index + integer :: ibegin,j,k + ao_cart_to_sphe_coef(:,:) = 0.d0 + ! Assume order provided by ao_power_index + i = 1 + ao_cart_to_sphe_num = 0 + do while (i <= ao_num) + select case ( ao_l(i) ) + case (0) + ao_cart_to_sphe_coef(i,i) = 1.d0 + i += 1 + ao_cart_to_sphe_num += 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,i+k-1) = cart_to_sphe_$SHELL(j,k) + enddo + enddo + i += size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_num += 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' + end select + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num,ao_cart_to_sphe_num) ] + implicit none + BEGIN_DOC + ! AO overlap matrix in the spherical basis set + END_DOC + double precision, allocatable :: S(:,:) + allocate (S(ao_cart_to_sphe_num,ao_num)) + + call dgemm('T','N',ao_cart_to_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_cart_to_sphe_num,ao_cart_to_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_cart_to_sphe_num,ao_num) ] + implicit none + BEGIN_DOC + ! AO_cart_to_sphe_coef^(-1) + END_DOC + + call get_pseudo_inverse(ao_cart_to_sphe_coef,ao_num,ao_cart_to_sphe_num, & + ao_cart_to_sphe_inv, size(ao_cart_to_sphe_coef,1)) +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)] &BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ] implicit none @@ -8,47 +92,39 @@ END_DOC integer :: i ao_ortho_canonical_coef(:,:) = 0.d0 + do i=1,ao_num + ao_ortho_canonical_coef(i,i) = 1.d0 + enddo + if (ao_cartesian) then - do i=1,ao_num - ao_ortho_canonical_coef(i,i) = 1.d0 - enddo + + ao_ortho_canonical_num = ao_num + call ortho_canonical(ao_overlap,size(ao_overlap,1), & + ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), & + ao_ortho_canonical_num) + + else - integer, external :: ao_power_index - integer :: ibegin,j,k - ! Assume order provided by ao_power_index - i = 1 - do while (i <= ao_num) - select case ( ao_l(i) ) - case (0) - ao_ortho_canonical_coef(i,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_ortho_canonical_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k) - enddo - enddo - i += size(cart_to_sphe_$SHELL,1) - endif - SUBST [ SHELL ] - 1;; - 2;; - 3;; - 4;; - 5;; - 6;; - 7;; - 8;; - 9;; - END_TEMPLATE - case default - stop 'Error in sphe_to_cart' - end select + + double precision, allocatable :: S(:,:) + + allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num)) + S = 0.d0 + do i=1,ao_cart_to_sphe_num + S(i,i) = 1.d0 enddo + + ao_ortho_canonical_num = ao_cart_to_sphe_num + call ortho_canonical (ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & + ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num) + + call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, & + ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & + S, size(S,1), & + 0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1)) + + deallocate(S) endif - call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,ao_ortho_canonical_num) END_PROVIDER BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)] diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index bcfd43b8..dfa5d982 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -88,7 +88,6 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) endif enddo do i=m+1,n - print *, D(i) D(i) = 0.d0 enddo