Corrected bug in spherical guess MOs

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.                               &

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
+    ao_ortho_canonical_num = ao_num
-    enddo
+    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
+    double precision, allocatable :: S(:,:)
-    ! Assume order provided by ao_power_index
+
-    i = 1
+    allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num))
-    do while (i <= ao_num)
+    S = 0.d0
-      select case ( ao_l(i) )
+    do i=1,ao_cart_to_sphe_num
-        case (0)
+      S(i,i) = 1.d0
-          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
     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)]

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