mirror of
https://github.com/LCPQ/quantum_package
synced 2024-11-06 22:24:00 +01:00
423 lines
15 KiB
Fortran
423 lines
15 KiB
Fortran
BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num_align,ao_num)]
|
|
&BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num_align,ao_num)]
|
|
&BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num_align,ao_num)]
|
|
BEGIN_DOC
|
|
! array of the integrals of AO_i * x^2 AO_j
|
|
! array of the integrals of AO_i * y^2 AO_j
|
|
! array of the integrals of AO_i * z^2 AO_j
|
|
END_DOC
|
|
implicit none
|
|
integer :: i,j,n,l
|
|
double precision :: f, tmp
|
|
integer :: dim1
|
|
double precision :: overlap, overlap_x, overlap_y, overlap_z
|
|
double precision :: alpha, beta
|
|
double precision :: A_center(3), B_center(3)
|
|
integer :: power_A(3), power_B(3)
|
|
double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z
|
|
dim1=500
|
|
lower_exp_val = 40.d0
|
|
ao_spread_x= 0.d0
|
|
ao_spread_y= 0.d0
|
|
ao_spread_z= 0.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, overlap, &
|
|
!$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) &
|
|
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
|
|
!$OMP ao_spread_x,ao_spread_y,ao_spread_z,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 )
|
|
!DEC$ VECTOR ALIGNED
|
|
!DEC$ VECTOR ALWAYS
|
|
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 )
|
|
accu_x = 0.d0
|
|
accu_y = 0.d0
|
|
accu_z = 0.d0
|
|
do n = 1,ao_prim_num(j)
|
|
alpha = ao_expo_ordered_transp(n,j)
|
|
!DEC$ VECTOR ALIGNED
|
|
do l = 1, ao_prim_num(i)
|
|
c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,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)
|
|
call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1)
|
|
accu_x += c*(tmp*overlap_y*overlap_z)
|
|
call overlap_bourrin_spread(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1)
|
|
accu_y += c*(tmp*overlap_x*overlap_z)
|
|
call overlap_bourrin_spread(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1)
|
|
accu_z += c*(tmp*overlap_y*overlap_x)
|
|
enddo
|
|
enddo
|
|
ao_spread_x(i,j) = accu_x
|
|
ao_spread_y(i,j) = accu_y
|
|
ao_spread_z(i,j) = accu_z
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num_align,ao_num)]
|
|
&BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num_align,ao_num)]
|
|
&BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num_align,ao_num)]
|
|
BEGIN_DOC
|
|
! array of the integrals of AO_i * x AO_j
|
|
! array of the integrals of AO_i * y AO_j
|
|
! array of the integrals of AO_i * z AO_j
|
|
END_DOC
|
|
implicit none
|
|
integer :: i,j,n,l
|
|
double precision :: f, tmp
|
|
integer :: dim1
|
|
double precision :: overlap, overlap_x, overlap_y, overlap_z,accu_x,accu_y,accu_z
|
|
double precision :: alpha, beta
|
|
double precision :: A_center(3), B_center(3)
|
|
integer :: power_A(3), power_B(3)
|
|
double precision :: lower_exp_val, dx, c
|
|
dim1=500
|
|
lower_exp_val = 40.d0
|
|
ao_dipole_x= 0.d0
|
|
ao_dipole_y= 0.d0
|
|
ao_dipole_z= 0.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, overlap, &
|
|
!$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) &
|
|
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
|
|
!$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,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 )
|
|
!DEC$ VECTOR ALIGNED
|
|
!DEC$ VECTOR ALWAYS
|
|
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 )
|
|
accu_x = 0.d0
|
|
accu_y = 0.d0
|
|
accu_z = 0.d0
|
|
do n = 1,ao_prim_num(j)
|
|
alpha = ao_expo_ordered_transp(n,j)
|
|
!DEC$ VECTOR ALIGNED
|
|
do l = 1, ao_prim_num(i)
|
|
beta = ao_expo_ordered_transp(l,i)
|
|
c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j)
|
|
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
|
|
|
|
call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1)
|
|
accu_x = accu_x + c*(tmp*overlap_y*overlap_z)
|
|
call overlap_bourrin_dipole(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1)
|
|
accu_y = accu_y + c*(tmp*overlap_x*overlap_z)
|
|
call overlap_bourrin_dipole(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1)
|
|
accu_z = accu_z + c*(tmp*overlap_y*overlap_x)
|
|
enddo
|
|
enddo
|
|
ao_dipole_x(i,j) = accu_x
|
|
ao_dipole_y(i,j) = accu_y
|
|
ao_dipole_z(i,j) = accu_z
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num_align,ao_num)]
|
|
&BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num_align,ao_num)]
|
|
&BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num_align,ao_num)]
|
|
BEGIN_DOC
|
|
! array of the integrals of AO_i * d/dx AO_j
|
|
! array of the integrals of AO_i * d/dy AO_j
|
|
! array of the integrals of AO_i * d/dz AO_j
|
|
END_DOC
|
|
implicit none
|
|
integer :: i,j,n,l
|
|
double precision :: f, tmp
|
|
integer :: dim1
|
|
double precision :: overlap, overlap_x, overlap_y, overlap_z
|
|
double precision :: alpha, beta
|
|
double precision :: A_center(3), B_center(3)
|
|
integer :: power_A(3), power_B(3)
|
|
double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z
|
|
integer :: i_component
|
|
dim1=500
|
|
lower_exp_val = 40.d0
|
|
ao_deriv_1_x= 0.d0
|
|
ao_deriv_1_y= 0.d0
|
|
ao_deriv_1_z= 0.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, overlap, &
|
|
!$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) &
|
|
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
|
|
!$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,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 )
|
|
!DEC$ VECTOR ALIGNED
|
|
!DEC$ VECTOR ALWAYS
|
|
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 )
|
|
accu_x = 0.d0
|
|
accu_y = 0.d0
|
|
accu_z = 0.d0
|
|
do n = 1,ao_prim_num(j)
|
|
alpha = ao_expo_ordered_transp(n,j)
|
|
!DEC$ VECTOR ALIGNED
|
|
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(l,i) * ao_coef_normalized_ordered_transp(n,j)
|
|
i_component = 1
|
|
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
|
|
accu_x += c*(tmp*overlap_y*overlap_z)
|
|
i_component = 2
|
|
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
|
|
accu_y += c*(tmp*overlap_x*overlap_z)
|
|
i_component = 3
|
|
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
|
|
accu_z += c*(tmp*overlap_y*overlap_x)
|
|
enddo
|
|
enddo
|
|
ao_deriv_1_x(i,j) = accu_x
|
|
ao_deriv_1_y(i,j) = accu_y
|
|
ao_deriv_1_z(i,j) = accu_z
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
subroutine overlap_bourrin_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
|
implicit none
|
|
! compute the following integral :
|
|
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ]
|
|
integer :: i,j,k,l
|
|
integer,intent(in) :: power_A,power_B
|
|
double precision, intent(in) :: lower_exp_val
|
|
double precision,intent(in) :: A_center, B_center,alpha,beta
|
|
double precision, intent(out) :: overlap_x,dx
|
|
integer, intent(in) :: nx
|
|
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
|
double precision :: P_center,pouet_timy
|
|
if(power_A.lt.0.or.power_B.lt.0)then
|
|
overlap_x = 0.d0
|
|
dx = 0.d0
|
|
return
|
|
endif
|
|
p = alpha + beta
|
|
p_inv= 1.d0/p
|
|
rho = alpha * beta * p_inv
|
|
dist = (A_center - B_center)*(A_center - B_center)
|
|
P_center = (alpha * A_center + beta * B_center) * p_inv
|
|
factor = dexp(-rho * dist)
|
|
pouet_timy = dsqrt(lower_exp_val/p)
|
|
x_min = P_center - pouet_timy
|
|
x_max = P_center + pouet_timy
|
|
!print*,'xmin = ',x_min
|
|
!print*,'xmax = ',x_max
|
|
domain = x_max-x_min
|
|
dx = domain/dble(nx)
|
|
overlap_x = 0.d0
|
|
x = x_min
|
|
do i = 1, nx
|
|
x += dx
|
|
overlap_x += (power(power_A,x-A_center) * power(power_B,x-B_center)) * dexp(-p * (x-P_center)*(x-P_center))
|
|
enddo
|
|
overlap_x *= factor * dx
|
|
end
|
|
|
|
subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
|
! compute the following integral :
|
|
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ]
|
|
! needed for the dipole and those things
|
|
implicit none
|
|
integer :: i,j,k,l
|
|
integer,intent(in) :: power_A,power_B
|
|
double precision, intent(in) :: lower_exp_val
|
|
double precision,intent(in) :: A_center, B_center,alpha,beta
|
|
double precision, intent(out) :: overlap_x,dx
|
|
integer, intent(in) :: nx
|
|
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
|
double precision :: P_center,pouet_timy
|
|
if(power_A.lt.0.or.power_B.lt.0)then
|
|
overlap_x = 0.d0
|
|
dx = 0.d0
|
|
return
|
|
endif
|
|
p = alpha + beta
|
|
p_inv= 1.d0/p
|
|
rho = alpha * beta * p_inv
|
|
dist = (A_center - B_center)*(A_center - B_center)
|
|
P_center = (alpha * A_center + beta * B_center) * p_inv
|
|
factor = dexp(-rho * dist)
|
|
if(factor.lt.0.000001d0)then
|
|
! print*,'factor = ',factor
|
|
dx = 0.d0
|
|
overlap_x = 0.d0
|
|
return
|
|
endif
|
|
pouet_timy = dsqrt(lower_exp_val/p)
|
|
x_min = P_center - pouet_timy
|
|
x_max = P_center + pouet_timy
|
|
domain = x_max-x_min
|
|
dx = domain/dble(nx)
|
|
overlap_x = 0.d0
|
|
x = x_min
|
|
do i = 1, nx
|
|
x += dx
|
|
overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x * x
|
|
enddo
|
|
overlap_x *= factor * dx
|
|
|
|
end
|
|
|
|
double precision function power(n,x)
|
|
implicit none
|
|
integer :: i,n
|
|
double precision :: x,accu
|
|
power = x**n
|
|
return
|
|
end
|
|
|
|
subroutine overlap_bourrin_dipole(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
|
! compute the following integral :
|
|
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ]
|
|
! needed for the dipole and those things
|
|
implicit none
|
|
integer :: i,j,k,l
|
|
integer,intent(in) :: power_A,power_B
|
|
double precision, intent(in) :: lower_exp_val
|
|
double precision,intent(in) :: A_center, B_center,alpha,beta
|
|
double precision, intent(out) :: overlap_x,dx
|
|
integer, intent(in) :: nx
|
|
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
|
double precision :: P_center
|
|
if(power_A.lt.0.or.power_B.lt.0)then
|
|
overlap_x = 0.d0
|
|
dx = 0.d0
|
|
return
|
|
endif
|
|
p = alpha + beta
|
|
p_inv= 1.d0/p
|
|
rho = alpha * beta * p_inv
|
|
dist = (A_center - B_center)*(A_center - B_center)
|
|
P_center = (alpha * A_center + beta * B_center) * p_inv
|
|
factor = dexp(-rho * dist)
|
|
if(power_B == 0 .and. power_A ==0)then
|
|
double precision :: F_integral
|
|
overlap_x = P_center * F_integral(0,p) * factor
|
|
dx = 0.d0
|
|
return
|
|
endif
|
|
double precision :: pouet_timy
|
|
|
|
pouet_timy = dsqrt(lower_exp_val/p)
|
|
x_min = P_center - pouet_timy
|
|
x_max = P_center + pouet_timy
|
|
domain = x_max-x_min
|
|
dx = domain/dble(nx)
|
|
overlap_x = 0.d0
|
|
x = x_min
|
|
do i = 1, nx
|
|
x += dx
|
|
overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x
|
|
enddo
|
|
overlap_x *= factor * dx
|
|
|
|
end
|
|
|
|
subroutine overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,overlap_x,nx)
|
|
implicit none
|
|
integer :: i,j,k,l
|
|
integer,intent(in) :: power_A(3),power_B(3),i_component
|
|
double precision,intent(in) :: A_center(3), B_center(3),alpha,beta,lower_exp_val
|
|
double precision, intent(out) :: overlap_x,dx
|
|
integer, intent(in) :: nx
|
|
double precision :: overlap_first, overlap_second
|
|
! computes : <phi_i|d/dx|phi_j> = (a_x_i <phi_i_x|phi_j_x(a_x_j-1)> - 2 alpha <phi_i_x|phi_j_w(a_x_j+1)>)
|
|
|
|
call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)-1,power_B(i_component),overlap_first,lower_exp_val,dx,nx)
|
|
call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)+1,power_B(i_component),overlap_second,lower_exp_val,dx,nx)
|
|
overlap_x = (power_A(i_component) * overlap_first - 2.d0 * alpha * overlap_second)
|
|
end
|
|
|
|
subroutine overlap_bourrin_x(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
|
implicit none
|
|
! compute the following integral :
|
|
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ]
|
|
integer :: i,j,k,l
|
|
integer,intent(in) :: power_A,power_B
|
|
double precision, intent(in) :: lower_exp_val
|
|
double precision,intent(in) :: A_center, B_center,alpha,beta
|
|
double precision, intent(out) :: overlap_x,dx
|
|
integer, intent(in) :: nx
|
|
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
|
double precision :: P_center,pouet_timy
|
|
if(power_A.lt.0.or.power_B.lt.0)then
|
|
overlap_x = 0.d0
|
|
dx = 0.d0
|
|
return
|
|
endif
|
|
p = alpha + beta
|
|
p_inv= 1.d0/p
|
|
rho = alpha * beta * p_inv
|
|
dist = (A_center - B_center)*(A_center - B_center)
|
|
P_center = (alpha * A_center + beta * B_center) * p_inv
|
|
factor = dexp(-rho * dist)
|
|
if(factor.lt.0.000001d0)then
|
|
dx = 0.d0
|
|
overlap_x = 0.d0
|
|
return
|
|
endif
|
|
|
|
pouet_timy = dsqrt(lower_exp_val/p)
|
|
x_min = P_center - pouet_timy
|
|
x_max = P_center + pouet_timy
|
|
domain = x_max-x_min
|
|
dx = domain/dble(nx)
|
|
overlap_x = 0.d0
|
|
x = x_min
|
|
do i = 1, nx
|
|
x += dx
|
|
overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center))
|
|
enddo
|
|
overlap_x *= factor * dx
|
|
end
|
|
|