Merge pull request #133 from scemama/master

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,36 +1,34 @@
- BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
+ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num_align,ao_num)]
-&BEGIN_PROVIDER [ integer, ao_ortho_canonical_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_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
+! ao_cart_to_sphe_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
   END_DOC
   integer :: i
-  ao_ortho_canonical_coef(:,:) = 0.d0
-  if (ao_cartesian) then
-    do i=1,ao_num
-      ao_ortho_canonical_coef(i,i) = 1.d0
-    enddo
-  else
   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_ortho_canonical_coef(i,i) = 1.d0
+        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_ortho_canonical_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k)
+              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;;
@@ -44,11 +42,89 @@
         9;;
       END_TEMPLATE
       case default
-          stop 'Error in sphe_to_cart'
+        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
+  BEGIN_DOC
+! matrix of the coefficients of the mos generated by the 
+! orthonormalization by the S^{-1/2} canonical transformation of the aos
+! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
+  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
+
+    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
+
+    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)]

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
 

tests/bats/qp.bats

@@ -100,11 +100,11 @@ function run_FCI() {
 }
 
 @test "SCF H2O cc-pVDZ" {
-  run_HF  h2o.ezfio  -76.0270218692377 
+  run_HF  h2o.ezfio  -76.0230273511649
 }
 
 @test "FCI H2O cc-pVDZ" {
-  run_FCI h2o.ezfio 10000  -76.2382562245107  -76.2433933430971
+  run_FCI h2o.ezfio 10000  -76.2315960087713  -76.237428352199
 }
 
 @test "CAS_SD H2O cc-pVDZ" {
@@ -113,10 +113,10 @@ function run_FCI() {
   ezfio set_file $INPUT
   ezfio set perturbation do_pt2_end False
   ezfio set determinants n_det_max 1000
-  qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-25]"
+  qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-23]"
   qp_run cas_sd_selected $INPUT 
   energy="$(ezfio get cas_sd energy)"
-  eq $energy -76.2219816183745 1.E-6
+  eq $energy -76.2085511296300 1.E-6
 }
 
 @test "MRCC H2O cc-pVDZ" {
@@ -128,7 +128,8 @@ function run_FCI() {
   ezfio set determinants read_wf True
   qp_run mrcc_cassd $INPUT 
   energy="$(ezfio get mrcc_cassd energy)"
-  eq $energy -76.2313853324571 1.E-3
+  eq $energy -76.2165731870755 1.E-3
+  
 }
 
 
@@ -138,11 +139,11 @@ function run_FCI() {
 }
 
 @test "SCF H2O VDZ pseudo" {
-  run_HF  h2o_pseudo.ezfio  -16.9483703904632
+  run_HF  h2o_pseudo.ezfio  -16.9457263818675
 }
 
 @test "FCI H2O VDZ pseudo" {
-  run_FCI h2o_pseudo.ezfio 2000    -17.1550015533975  -17.1645044211228
+  run_FCI h2o_pseudo.ezfio 2000    -17.1476897854369  -17.1598005211929
 }
 
 #=== Convert