diff --git a/src/Dets/s2.irp.f b/src/Dets/s2.irp.f index 588ab4cc..cd1d9fda 100644 --- a/src/Dets/s2.irp.f +++ b/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 diff --git a/src/MOs/utils.irp.f b/src/MOs/utils.irp.f index 7e23f744..f13982ac 100644 --- a/src/MOs/utils.irp.f +++ b/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 diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index b29654b1..7d173246 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/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 - diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index 273b4f48..79db5473 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/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)