mirror of
https://github.com/QuantumPackage/qp2.git
synced 20240723 02:57:24 +02:00
312 lines
8.7 KiB
Fortran
312 lines
8.7 KiB
Fortran


! 




BEGIN_PROVIDER [ double precision, ao_overlap , (ao_num, ao_num) ]


&BEGIN_PROVIDER [ double precision, ao_overlap_x, (ao_num, ao_num) ]


&BEGIN_PROVIDER [ double precision, ao_overlap_y, (ao_num, ao_num) ]


&BEGIN_PROVIDER [ double precision, ao_overlap_z, (ao_num, ao_num) ]




BEGIN_DOC


! Overlap between atomic basis functions:


!


! :math:`\int \chi_i(r) \chi_j(r) dr`


END_DOC




implicit none


integer :: i, j, n, l, dim1, power_A(3), power_B(3)


double precision :: overlap, overlap_x, overlap_y, overlap_z


double precision :: alpha, beta, c


double precision :: A_center(3), B_center(3)




ao_overlap = 0.d0


ao_overlap_x = 0.d0


ao_overlap_y = 0.d0


ao_overlap_z = 0.d0




if(read_ao_integrals_overlap) then




call ezfio_get_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num))


print *, 'AO overlap integrals read from disk'




else




if(use_cosgtos) then


!print*, ' use_cosgtos for ao_overlap ?', use_cosgtos




do j = 1, ao_num


do i = 1, ao_num


ao_overlap (i,j) = ao_overlap_cosgtos (i,j)


ao_overlap_x(i,j) = ao_overlap_cosgtos_x(i,j)


ao_overlap_y(i,j) = ao_overlap_cosgtos_y(i,j)


ao_overlap_z(i,j) = ao_overlap_cosgtos_z(i,j)


enddo


enddo




else




dim1=100


!$OMP PARALLEL DO SCHEDULE(GUIDED) &


!$OMP DEFAULT(NONE) &


!$OMP PRIVATE(A_center,B_center,power_A,power_B,&


!$OMP overlap_x,overlap_y, overlap_z, overlap, &


!$OMP alpha, beta,i,j,c) &


!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &


!$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &


!$OMP ao_expo_ordered_transp,dim1)


do j=1,ao_num


A_center(1) = nucl_coord( ao_nucl(j), 1 )


A_center(2) = nucl_coord( ao_nucl(j), 2 )


A_center(3) = nucl_coord( ao_nucl(j), 3 )


power_A(1) = ao_power( j, 1 )


power_A(2) = ao_power( j, 2 )


power_A(3) = ao_power( j, 3 )


do i= 1,ao_num


B_center(1) = nucl_coord( ao_nucl(i), 1 )


B_center(2) = nucl_coord( ao_nucl(i), 2 )


B_center(3) = nucl_coord( ao_nucl(i), 3 )


power_B(1) = ao_power( i, 1 )


power_B(2) = ao_power( i, 2 )


power_B(3) = ao_power( i, 3 )


do n = 1,ao_prim_num(j)


alpha = ao_expo_ordered_transp(n,j)


do l = 1, ao_prim_num(i)


beta = ao_expo_ordered_transp(l,i)


call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)


c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)


ao_overlap(i,j) += c * overlap


if(isnan(ao_overlap(i,j)))then


print*,'i,j',i,j


print*,'l,n',l,n


print*,'c,overlap',c,overlap


print*,overlap_x,overlap_y,overlap_z


stop


endif


ao_overlap_x(i,j) += c * overlap_x


ao_overlap_y(i,j) += c * overlap_y


ao_overlap_z(i,j) += c * overlap_z


enddo


enddo


enddo


enddo


!$OMP END PARALLEL DO




endif




endif




if (write_ao_integrals_overlap) then


call ezfio_set_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num))


print *, 'AO overlap integrals written to disk'


endif




END_PROVIDER




! 




BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]


implicit none


BEGIN_DOC


! Imaginary part of the overlap


END_DOC


ao_overlap_imag = 0.d0


END_PROVIDER




! 




BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ]


implicit none


BEGIN_DOC


! Overlap for complex AOs


END_DOC


integer :: i,j


do j=1,ao_num


do i=1,ao_num


ao_overlap_complex(i,j) = dcmplx( ao_overlap(i,j), ao_overlap_imag(i,j) )


enddo


enddo


END_PROVIDER




! 




BEGIN_PROVIDER [ double precision, ao_overlap_abs, (ao_num, ao_num) ]




BEGIN_DOC


! Overlap between absolute values of atomic basis functions:


!


! :math:`\int \chi_i(r) \chi_j(r) dr`


END_DOC




implicit none


integer :: i, j, n, l, dim1, power_A(3), power_B(3)


double precision :: overlap_x, overlap_y, overlap_z


double precision :: alpha, beta


double precision :: A_center(3), B_center(3)


double precision :: lower_exp_val, dx




if(is_periodic) then




do j = 1, ao_num


do i = 1, ao_num


ao_overlap_abs(i,j) = cdabs(ao_overlap_complex(i,j))


enddo


enddo




else




dim1=100


lower_exp_val = 40.d0


!$OMP PARALLEL DO SCHEDULE(GUIDED) &


!$OMP DEFAULT(NONE) &


!$OMP PRIVATE(A_center,B_center,power_A,power_B, &


!$OMP overlap_x,overlap_y, overlap_z, &


!$OMP alpha, beta,i,j,dx) &


!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &


!$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl,&


!$OMP ao_expo_ordered_transp,dim1,lower_exp_val)


do j=1,ao_num


A_center(1) = nucl_coord( ao_nucl(j), 1 )


A_center(2) = nucl_coord( ao_nucl(j), 2 )


A_center(3) = nucl_coord( ao_nucl(j), 3 )


power_A(1) = ao_power( j, 1 )


power_A(2) = ao_power( j, 2 )


power_A(3) = ao_power( j, 3 )


do i= 1,ao_num


ao_overlap_abs(i,j)= 0.d0


B_center(1) = nucl_coord( ao_nucl(i), 1 )


B_center(2) = nucl_coord( ao_nucl(i), 2 )


B_center(3) = nucl_coord( ao_nucl(i), 3 )


power_B(1) = ao_power( i, 1 )


power_B(2) = ao_power( i, 2 )


power_B(3) = ao_power( i, 3 )


do n = 1,ao_prim_num(j)


alpha = ao_expo_ordered_transp(n,j)


do l = 1, ao_prim_num(i)


beta = ao_expo_ordered_transp(l,i)


call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)


call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1)


call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1)


ao_overlap_abs(i,j) += abs(ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)) * overlap_x * overlap_y * overlap_z


enddo


enddo


enddo


enddo


!$OMP END PARALLEL DO




endif




END_PROVIDER




! 




BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]


implicit none


BEGIN_DOC


! Inverse of the overlap matrix


END_DOC


call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv, &


size(S_inv,1),lin_dep_cutoff)


END_PROVIDER




BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ]


implicit none


BEGIN_DOC


! Inverse of the overlap matrix


END_DOC


call get_pseudo_inverse_complex(ao_overlap_complex, size(ao_overlap_complex,1),&


ao_num,ao_num,S_inv_complex,size(S_inv_complex,1),lin_dep_cutoff)


END_PROVIDER




BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]




BEGIN_DOC


! :math:`X = S^{1/2}` obtained by SVD


END_DOC




implicit none




integer :: num_linear_dependencies


integer :: LDA, LDC


double precision, allocatable :: U(:,:),Vt(:,:), D(:)


integer :: info, i, j, k


double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d6




LDA = size(AO_overlap,1)


LDC = size(S_half_inv,1)




allocate( &


U(LDC,AO_num), &


Vt(LDA,AO_num), &


D(AO_num))




call svd( &


AO_overlap,LDA, &


U,LDC, &


D, &


Vt,LDA, &


AO_num,AO_num)




num_linear_dependencies = 0


do i=1,AO_num


print*,D(i)


if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then


D(i) = 0.d0


num_linear_dependencies += 1


else


ASSERT (D(i) > 0.d0)


D(i) = 1.d0/sqrt(D(i))


endif


do j=1,AO_num


S_half_inv(j,i) = 0.d0


enddo


enddo


write(*,*) 'linear dependencies',num_linear_dependencies




do k=1,AO_num


if(D(k) /= 0.d0) then


do j=1,AO_num


do i=1,AO_num


S_half_inv(i,j) = S_half_inv(i,j) + U(i,k)*D(k)*Vt(k,j)


enddo


enddo


endif


enddo






END_PROVIDER






BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ]


implicit none


BEGIN_DOC


! :math:`S^{1/2}`


END_DOC




integer :: i,j,k


double precision, allocatable :: U(:,:)


double precision, allocatable :: Vt(:,:)


double precision, allocatable :: D(:)




allocate(U(ao_num,ao_num),Vt(ao_num,ao_num),D(ao_num))




call svd(ao_overlap,size(ao_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_num,ao_num)




do i=1,ao_num


D(i) = dsqrt(D(i))


do j=1,ao_num


S_half(j,i) = 0.d0


enddo


enddo




do k=1,ao_num


do j=1,ao_num


do i=1,ao_num


S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j)


enddo


enddo


enddo




deallocate(U,Vt,D)




END_PROVIDER


