Minor changes

src/Dets/s2.irp.f

@@ -64,6 +64,20 @@ BEGIN_PROVIDER [ double precision, expected_s2]
 
 END_PROVIDER 
 
+BEGIN_PROVIDER [ double precision, s2_values, (N_states) ]
+ implicit none
+ BEGIN_DOC
+! array of the averaged values of the S^2 operator on the various states
+ END_DOC
+ integer :: i
+ double precision :: s2
+ do i = 1, N_states
+  call get_s2_u0(psi_det,psi_coef(1,i),n_det,psi_det_size,s2)
+  s2_values(i) = s2
+ enddo
+
+END_PROVIDER
+
 
 subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
  implicit none

src/MOs/utils.irp.f

@@ -38,6 +38,9 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label)
   mo_coef_new = mo_coef
   
   call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2))
+  do i = 1, mo_tot_num
+   print*,R(i,4)
+  enddo
   integer :: i
   write (output_mos,'(A)'), 'MOs are now **'//trim(label)//'**'
   write (output_mos,'(A)'), ''
@@ -58,3 +61,83 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label)
   mo_label = label
   SOFT_TOUCH mo_coef mo_label
 end
+subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)
+  implicit none
+  integer,intent(in)             :: n,m
+  character*(64), intent(in)     :: label
+  double precision, intent(in)   :: matrix(n,m),observable(n,n)
+  
+  double precision, allocatable  :: mo_coef_new(:,:), R(:,:),eigvalues(:),value(:)
+  integer,allocatable :: iorder(:)
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R
+  
+  call write_time(output_mos)
+  if (m /= mo_tot_num) then
+    print *, irp_here, ': Error : m/= mo_tot_num'
+    stop 1
+  endif
+  allocate(R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m),value(m),iorder(m))
+  mo_coef_new = mo_coef
+  
+  call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2))
+
+  do i = 1, mo_tot_num
+   value(i) = 0.d0
+   iorder(i) = i
+   do j = 1, mo_tot_num
+    do k = 1, mo_tot_num
+     value(i) += R(k,i) * R(j,i) * observable(k,j)
+    enddo
+   enddo
+!  print*,'value(i) = ',i,value(i)
+  enddo
+  integer :: i,j,k,index
+  double precision :: R_tmp(m,m),obs(m)
+  print*,'sort ....'
+  call dsort(value,iorder,n)
+  do i = 1, mo_tot_num
+   index = iorder(i)
+!  print*,'iorder(i) = ',iorder(i)
+   obs(i) = eigvalues(iorder(i))
+   do j = 1, mo_tot_num
+    R_tmp(j,i) = R(j,index)
+   enddo
+  enddo
+  do i = 1, mo_tot_num
+   do j = 1, mo_tot_num
+    R(j,i) = R_tmp(j,i)
+   enddo
+  enddo
+
+  do i = 1, mo_tot_num
+   value(i) = 0.d0
+   do j = 1, mo_tot_num
+    do k = 1, mo_tot_num
+     value(i) += R(k,i) * R(j,i) * observable(k,j)
+    enddo
+   enddo
+   print*,'i = ',i
+   print*,'value(i) = ',value(i)
+   print*,'obs(i)   = ',obs(i)
+   print*,''
+  enddo
+
+  write (output_mos,'(A)'), 'MOs are now **'//trim(label)//'**'
+  write (output_mos,'(A)'), ''
+  write (output_mos,'(A)'), 'Eigenvalues'
+  write (output_mos,'(A)'), '-----------'
+  write (output_mos,'(A)'), ''
+  write (output_mos,'(A)'), '======== ================'
+  do i = 1, m
+   write (output_mos,'(I8,X,F16.10)'), i,eigvalues(i)
+  enddo
+  write (output_mos,'(A)'), '======== ================'
+  write (output_mos,'(A)'), ''
+  
+  call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
+  deallocate(mo_coef_new,R,eigvalues)
+  call write_time(output_mos)
+  
+  mo_label = label
+  SOFT_TOUCH mo_coef mo_label
+end

src/Utils/LinearAlgebra.irp.f

@@ -430,88 +430,24 @@ subroutine lapack_partial_diag(eigvalues,eigvectors,H,nmax,n,n_st)
 end
 
 
+subroutine set_zero_extra_diag(i1,i2,matrix,lda,m)
+ implicit none
+ integer, intent(in) :: i1,i2,lda,m
+ double precision, intent(inout) :: matrix(lda,m)
+ integer :: i,j
+ do j=i1,i2
+  do i = 1,i1-1
+   matrix(i,j) = 0.d0
+   matrix(j,i) = 0.d0
+  enddo
+ enddo
+
+ do i = i2,i1
+  do j=i2+1,m
+   matrix(i,j) = 0.d0
+   matrix(j,i) = 0.d0
+  enddo
+ enddo
+
 
-subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
-  implicit none
-  BEGIN_DOC
-  ! Transform A from the AO basis to the MO basis
-  END_DOC
-  double precision, intent(in)   :: A_ao(LDA_ao)
-  double precision, intent(out)  :: A_mo(LDA_mo)
-  integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, allocatable  :: T(:,:)
-  
-  allocate ( T(ao_num_align,mo_tot_num) )
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
-  
-  call dgemm('N','N', ao_num, mo_tot_num, ao_num,                    &
-      1.d0, A_ao,LDA_ao,                                             &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0, T, ao_num_align)
-  
-  call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num,                &
-      1.d0, mo_coef,size(mo_coef,1),                                 &
-      T, ao_num_align,                                               &
-      0.d0, A_mo, LDA_mo)
-  
-  deallocate(T)
 end
-
-subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
-  implicit none
-  BEGIN_DOC
-  ! Transform A from the MO basis to the AO basis
-  END_DOC
-  double precision, intent(in)   :: A_mo(LDA_mo)
-  double precision, intent(out)  :: A_ao(LDA_ao)
-  integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, allocatable  :: T(:,:), SC(:,:)
-  
-  allocate ( SC(ao_num_align,mo_tot_num) )
-  allocate ( T(mo_tot_num_align,ao_num) )
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
-  
-  call dgemm('N','N', ao_num, mo_tot_num, ao_num,                    &
-      1.d0, ao_overlap,size(ao_overlap,1),                           &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0, SC, ao_num_align)
-  
-  call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
-      1.d0, A_mo,LDA_mo,                                             &
-      SC, size(SC,1),                                      &
-      0.d0, T, mo_tot_num_align)
-  
-  call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
-      1.d0, SC,size(SC,1),                                 &
-      T, mo_tot_num_align,                                           &
-      0.d0, A_ao, LDA_ao)
-  
-  deallocate(T,SC)
-end
-
-subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
-  implicit none
-  BEGIN_DOC
-  ! Transform A from the MO basis to the S^-1 AO basis
-  END_DOC
-  double precision, intent(in)   :: A_mo(LDA_mo)
-  double precision, intent(out)  :: A_ao(LDA_ao)
-  integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, allocatable  :: T(:,:)
-  
-  allocate ( T(mo_tot_num_align,ao_num) )
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
-  
-  call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
-      1.d0, A_mo,LDA_mo,                                             &
-      mo_coef, size(mo_coef,1),                                      &
-      0.d0, T, mo_tot_num_align)
-  
-  call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
-      1.d0, mo_coef,size(mo_coef,1),                                 &
-      T, mo_tot_num_align,                                           &
-      0.d0, A_ao, LDA_ao)
-  
-  deallocate(T)
-end
-

src/Utils/one_e_integration.irp.f

@@ -31,8 +31,123 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
   overlap_gaussian_x*= fact_p
 end
 
+subroutine test(alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap_x,overlap_y,overlap_z,overlap)
+   implicit none
+  include 'constants.F'
+ integer, intent(in)            :: a(3),b(3)         ! powers : (x-xa)**a_x = (x-A(1))**a(1)
+ double precision, intent(in)   :: alpha, beta, gama ! exponents
+ double precision, intent(in)   :: A_center(3)       ! A center
+ double precision, intent(in)   :: B_center (3)      ! B center
+ double precision, intent(in)   :: Nucl_center(3)    ! B center
+ double precision, intent(out)  :: overlap_x,overlap_y,overlap_z,overlap
+ integer :: i,j
+ double precision :: dx,Lx,nx,x(3)
+ nx = 100000000
+ Lx = 25.d0
+ dx = dble(Lx/nx)
+ overlap_x = 0.d0
+ overlap_y = 0.d0
+ overlap_z = 0.d0
+ x(1) = -12.5d0
+ x(2) = -12.5d0
+ x(3) = -12.5d0
+ do i = 1,nx 
+  overlap_x += (x(1) - A_center(1))**a(1) * (x(1) - B_center(1))**b(1) &
+              * dexp(-alpha*(x(1) - A_center(1))**2) * dexp(-beta*(x(1) - B_center(1))**2) * dexp(-gama*(x(1) - Nucl_center(1))**2) 
+
+  overlap_y += (x(2) - A_center(2))**a(2) * (x(2) - B_center(2))**b(2) &
+              * dexp(-alpha*(x(2) - A_center(2))**2) * dexp(-beta*(x(2) - B_center(2))**2) * dexp(-gama*(x(2) - Nucl_center(2))**2) 
+  overlap_z += (x(3) - A_center(3))**a(3) * (x(3) - B_center(3))**b(3) &
+              * dexp(-alpha*(x(3) - A_center(3))**2) * dexp(-beta*(x(3) - B_center(3))**2) * dexp(-gama*(x(3) - Nucl_center(3))**2) 
+  x(1) += dx
+  x(2) += dx
+  x(3) += dx
+ enddo
+ overlap_x = overlap_x * dx
+ overlap_y = overlap_y * dx
+ overlap_z = overlap_z * dx
+ overlap = overlap_x * overlap_y * overlap_z
+
+end
 
 
+subroutine overlap_A_B_C(dim,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap)
+ implicit none
+  include 'constants.F'
+ integer, intent(in)            :: dim                                                             
+ integer, intent(in)            :: a(3),b(3)         ! powers : (x-xa)**a_x = (x-A(1))**a(1)
+ double precision, intent(in)   :: alpha, beta, gama ! exponents
+ double precision, intent(in)   :: A_center(3)       ! A center
+ double precision, intent(in)   :: B_center (3)      ! B center
+ double precision, intent(in)   :: Nucl_center(3)    ! B center
+ double precision, intent(out)  :: overlap
+
+ double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,p
+ double precision               :: F_integral_tab(0:max_dim)
+ integer                        :: iorder_p(3)
+ double precision :: overlap_x,overlap_z,overlap_y
+
+ call give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_p,iorder_p,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim)
+  if(fact_p.lt.1d-10)then
+!   overlap_x = 0.d0
+!   overlap_y = 0.d0
+!   overlap_z = 0.d0
+    overlap = 0.d0
+    return
+  endif
+
+  integer                        :: nmax
+  double precision               :: F_integral
+  nmax = maxval(iorder_p)
+  do i = 0,nmax
+    F_integral_tab(i) = F_integral(i,p)
+  enddo
+  overlap_x = P_new(0,1) * F_integral_tab(0)
+  overlap_y = P_new(0,2) * F_integral_tab(0)
+  overlap_z = P_new(0,3) * F_integral_tab(0)
+  
+  integer                        :: i
+  do i = 1,iorder_p(1)
+    overlap_x += P_new(i,1) * F_integral_tab(i)
+  enddo
+  do i = 1,iorder_p(2)
+    overlap_y += P_new(i,2) * F_integral_tab(i)
+  enddo
+  do i = 1,iorder_p(3)
+    overlap_z += P_new(i,3) * F_integral_tab(i)
+  enddo
+  overlap = overlap_x * overlap_y * overlap_z * fact_p 
+  
+!double precision :: overlap_x_1,overlap_y_1,overlap_z_1,overlap_1
+!call test(alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap_x_1,overlap_y_1,overlap_z_1,overlap_1)
+!print*,'overlap_1 = ',overlap_1
+!print*,'overlap   = ',overlap
+!if(dabs(overlap - overlap_1).ge.1.d-3)then
+!  print*,'power_A(1) = ',a(1)
+!  print*,'power_A(2) = ',a(2)
+!  print*,'power_A(3) = ',a(3)
+!  print*,'power_B(1) = ',b(1)
+!  print*,'power_B(2) = ',b(2)
+!  print*,'power_B(3) = ',b(3)
+!  print*,'alpha = ',alpha
+!  print*,'beta = ',beta
+!  print*,'gama = ',gama
+!  print*,'A_center(1) = ',A_center(1)
+!  print*,'A_center(2) = ',A_center(2)
+!  print*,'A_center(3) = ',A_center(3)
+!  print*,'B_center(1) = ',B_center(1)
+!  print*,'B_center(2) = ',B_center(2)
+!  print*,'B_center(3) = ',B_center(3)
+!  print*,'Nucl_center(1) = ',Nucl_center(1)
+!  print*,'Nucl_center(2) = ',Nucl_center(2)
+!  print*,'Nucl_center(3) = ',Nucl_center(3)
+!  print*,'overlap = ',overlap
+!  print*,'overlap_1=',overlap_1
+
+!  stop
+!endif
+
+end
 
 subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
       power_B,overlap_x,overlap_y,overlap_z,overlap,dim)