mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-22 19:43:32 +01:00
1st version of grad + lapl of Jmu_modif
This commit is contained in:
parent
6513358da3
commit
7a4f732e15
@ -1,4 +1,6 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine phi_j_erf_mu_r_xyz_phi(i,j,mu_in, C_center, xyz_ints)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -49,26 +51,37 @@ subroutine phi_j_erf_mu_r_xyz_phi(i,j,mu_in, C_center, xyz_ints)
|
||||
enddo
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
double precision function phi_j_erf_mu_r_phi(i, j, mu_in, C_center)
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! phi_j_erf_mu_r_phi = int dr phi_j(r) [erf(mu |r - C|)/|r-C|] phi_i(r)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i,j
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
|
||||
integer :: num_A, power_A(3), num_b, power_B(3)
|
||||
double precision :: alpha, beta, A_center(3), B_center(3),contrib,NAI_pol_mult_erf
|
||||
integer :: n_pt_in, l, m
|
||||
double precision :: alpha, beta, A_center(3), B_center(3), contrib
|
||||
|
||||
double precision :: NAI_pol_mult_erf
|
||||
|
||||
phi_j_erf_mu_r_phi = 0.d0
|
||||
|
||||
if(ao_overlap_abs(j,i).lt.1.d-12) then
|
||||
return
|
||||
endif
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
! j
|
||||
num_A = ao_nucl(j)
|
||||
power_A(1:3) = ao_power(j,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
|
||||
! i
|
||||
num_B = ao_nucl(i)
|
||||
power_B(1:3) = ao_power(i,1:3)
|
||||
@ -78,14 +91,16 @@ double precision function phi_j_erf_mu_r_phi(i,j,mu_in, C_center)
|
||||
alpha = ao_expo_ordered_transp(l,j)
|
||||
do m = 1, ao_prim_num(i)
|
||||
beta = ao_expo_ordered_transp(m,i)
|
||||
|
||||
contrib = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
phi_j_erf_mu_r_phi += contrib * ao_coef_normalized_ordered_transp(l,j) &
|
||||
* ao_coef_normalized_ordered_transp(m,i)
|
||||
|
||||
phi_j_erf_mu_r_phi += contrib * ao_coef_normalized_ordered_transp(l,j) * ao_coef_normalized_ordered_transp(m,i)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
end function phi_j_erf_mu_r_phi
|
||||
|
||||
! ---
|
||||
|
||||
subroutine erfc_mu_gauss_xyz_ij_ao(i, j, mu, C_center, delta, gauss_ints)
|
||||
implicit none
|
||||
@ -132,29 +147,41 @@ subroutine erfc_mu_gauss_xyz_ij_ao(i,j,mu, C_center, delta,gauss_ints)
|
||||
enddo
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine erf_mu_gauss_ij_ao(i, j, mu, C_center, delta, gauss_ints)
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! gauss_ints(m) = \int dr exp(-delta (r - C)^2 ) * erf(mu |r-r'|)/ |r-r'| * AO_i(r') * AO_j(r')
|
||||
!
|
||||
! gauss_ints = \int dr exp(-delta (r - C)^2) * erf(mu |r-C|) / |r-C| * AO_i(r) * AO_j(r)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i, j
|
||||
double precision, intent(in) :: mu, C_center(3), delta
|
||||
double precision, intent(out) :: gauss_ints
|
||||
|
||||
integer :: n_pt_in, l, m
|
||||
integer :: num_A, power_A(3), num_b, power_B(3)
|
||||
double precision :: alpha, beta, A_center(3), B_center(3),contrib,NAI_pol_mult_erf
|
||||
double precision :: integral , erf_mu_gauss
|
||||
integer :: n_pt_in,l,m,mm
|
||||
double precision :: alpha, beta, A_center(3), B_center(3), coef
|
||||
double precision :: integral
|
||||
|
||||
double precision :: erf_mu_gauss
|
||||
|
||||
gauss_ints = 0.d0
|
||||
|
||||
if(ao_overlap_abs(j,i).lt.1.d-12) then
|
||||
return
|
||||
endif
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
! j
|
||||
num_A = ao_nucl(j)
|
||||
power_A(1:3) = ao_power(j,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
|
||||
! i
|
||||
num_B = ao_nucl(i)
|
||||
power_B(1:3) = ao_power(i,1:3)
|
||||
@ -164,63 +191,160 @@ subroutine erf_mu_gauss_ij_ao(i,j,mu, C_center, delta,gauss_ints)
|
||||
alpha = ao_expo_ordered_transp(l,j)
|
||||
do m = 1, ao_prim_num(i)
|
||||
beta = ao_expo_ordered_transp(m,i)
|
||||
if(dabs(ao_coef_normalized_ordered_transp(l,j) * ao_coef_normalized_ordered_transp(m,i)).lt.1.d-12)cycle
|
||||
integral = erf_mu_gauss(C_center,delta,mu,A_center,B_center,power_A,power_B,alpha,beta,n_pt_in)
|
||||
gauss_ints += integral * ao_coef_normalized_ordered_transp(l,j) &
|
||||
* ao_coef_normalized_ordered_transp(m,i)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
coef = ao_coef_normalized_ordered_transp(l,j) * ao_coef_normalized_ordered_transp(m,i)
|
||||
|
||||
if(dabs(coef) .lt. 1.d-12) cycle
|
||||
|
||||
integral = erf_mu_gauss(C_center, delta, mu, A_center, B_center, power_A, power_B, alpha, beta, n_pt_in)
|
||||
|
||||
gauss_ints += integral * coef
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine erf_mu_gauss_ij_ao
|
||||
|
||||
! ---
|
||||
|
||||
subroutine NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
!
|
||||
! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
!
|
||||
! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
double precision, intent(out) :: ints(3)
|
||||
double precision :: A_center(3), B_center(3),integral, alpha,beta
|
||||
double precision :: NAI_pol_mult_erf
|
||||
|
||||
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, power_xA(3), m
|
||||
double precision :: A_center(3), B_center(3), integral, alpha, beta, coef
|
||||
|
||||
double precision :: NAI_pol_mult_erf
|
||||
|
||||
ints = 0.d0
|
||||
if(ao_overlap_abs(j_ao,i_ao).lt.1.d-12) then
|
||||
return
|
||||
endif
|
||||
|
||||
num_A = ao_nucl(i_ao)
|
||||
power_A(1:3) = ao_power(i_ao,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
num_B = ao_nucl(j_ao)
|
||||
power_B(1:3) = ao_power(j_ao,1:3)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alpha = ao_expo_ordered_transp(i,i_ao)
|
||||
|
||||
do m = 1, 3
|
||||
|
||||
power_xA = power_A
|
||||
! x * phi_i(r) = x * (x-Ax)**ax = (x-Ax)**(ax+1) + Ax * (x-Ax)**ax
|
||||
power_xA(m) += 1
|
||||
|
||||
do j = 1, ao_prim_num(j_ao)
|
||||
beta = ao_expo_ordered_transp(j,j_ao)
|
||||
coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
|
||||
! First term = (x-Ax)**(ax+1)
|
||||
integral = NAI_pol_mult_erf(A_center, B_center, power_xA, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
ints(m) += integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
ints(m) += integral * coef
|
||||
|
||||
! Second term = Ax * (x-Ax)**(ax)
|
||||
integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
ints(m) += A_center(m) * integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
ints(m) += A_center(m) * integral * coef
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
end subroutine NAI_pol_x_mult_erf_ao
|
||||
|
||||
! ---
|
||||
|
||||
subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
!
|
||||
! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
!
|
||||
! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3)
|
||||
double precision, intent(out) :: ints(3)
|
||||
|
||||
integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, power_xA(3), m
|
||||
double precision :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef
|
||||
|
||||
double precision, external :: NAI_pol_mult_erf_with1s
|
||||
|
||||
ints = 0.d0
|
||||
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) then
|
||||
return
|
||||
endif
|
||||
|
||||
power_Ai(1:3) = ao_power(i_ao,1:3)
|
||||
power_Aj(1:3) = ao_power(j_ao,1:3)
|
||||
|
||||
Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
|
||||
Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alphai = ao_expo_ordered_transp(i,i_ao)
|
||||
|
||||
do m = 1, 3
|
||||
|
||||
power_xA = power_Ai
|
||||
! x * phi_i(r) = x * (x-Ax)**ax = (x-Ax)**(ax+1) + Ax * (x-Ax)**ax
|
||||
power_xA(m) += 1
|
||||
|
||||
do j = 1, ao_prim_num(j_ao)
|
||||
alphaj = ao_expo_ordered_transp(j,j_ao)
|
||||
coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
|
||||
! First term = (x-Ax)**(ax+1)
|
||||
integral = NAI_pol_mult_erf_with1s( Ai_center, Aj_center, power_xA, power_Aj, alphai, alphaj &
|
||||
, beta, b_center, c_center, n_pt_in, mu_in )
|
||||
ints(m) += integral * coef
|
||||
|
||||
! Second term = Ax * (x-Ax)**(ax)
|
||||
integral = NAI_pol_mult_erf_with1s( Ai_center, Aj_center, power_Ai, power_Aj, alphai, alphaj &
|
||||
, beta, b_center, c_center, n_pt_in, mu_in )
|
||||
ints(m) += Ai_center(m) * integral * coef
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine NAI_pol_x_mult_erf_ao_with1s
|
||||
|
||||
! ---
|
||||
|
||||
subroutine NAI_pol_x_specify_mult_erf_ao(i_ao,j_ao,mu_in,C_center,m,ints)
|
||||
implicit none
|
||||
@ -249,7 +373,6 @@ subroutine NAI_pol_x_specify_mult_erf_ao(i_ao,j_ao,mu_in,C_center,m,ints)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alpha = ao_expo_ordered_transp(i,i_ao)
|
||||
power_xA = power_A
|
||||
@ -267,3 +390,5 @@ subroutine NAI_pol_x_specify_mult_erf_ao(i_ao,j_ao,mu_in,C_center,m,ints)
|
||||
enddo
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -102,36 +102,118 @@ subroutine overlap_gauss_r12_all_ao(D_center,delta,aos_ints)
|
||||
enddo
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
! TODO :: PUT CYCLES IN LOOPS
|
||||
double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! \int dr AO_i(r) AO_j(r) e^{-delta |r-D_center|^2}
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i, j
|
||||
double precision, intent(in) :: D_center(3), delta
|
||||
|
||||
integer :: num_a,num_b,power_A(3), power_B(3),l,k
|
||||
double precision :: A_center(3), B_center(3),overlap_gauss_r12,alpha,beta,analytical_j
|
||||
integer :: power_A(3), power_B(3), l, k
|
||||
double precision :: A_center(3), B_center(3), alpha, beta, coef, analytical_j
|
||||
|
||||
double precision, external :: overlap_gauss_r12
|
||||
|
||||
overlap_gauss_r12_ao = 0.d0
|
||||
|
||||
if(ao_overlap_abs(j,i).lt.1.d-12) then
|
||||
return
|
||||
endif
|
||||
! TODO :: PUT CYCLES IN LOOPS
|
||||
num_A = ao_nucl(i)
|
||||
|
||||
power_A(1:3) = ao_power(i,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
num_B = ao_nucl(j)
|
||||
power_B(1:3) = ao_power(j,1:3)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
|
||||
A_center(1:3) = nucl_coord(ao_nucl(i),1:3)
|
||||
B_center(1:3) = nucl_coord(ao_nucl(j),1:3)
|
||||
|
||||
do l = 1, ao_prim_num(i)
|
||||
alpha = ao_expo_ordered_transp(l,i)
|
||||
do k = 1, ao_prim_num(j)
|
||||
beta = ao_expo_ordered_transp(k,j)
|
||||
coef = ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(k,j)
|
||||
|
||||
if(dabs(coef) .lt. 1d-12) cycle
|
||||
|
||||
analytical_j = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
|
||||
overlap_gauss_r12_ao += analytical_j * ao_coef_normalized_ordered_transp(l,i) &
|
||||
* ao_coef_normalized_ordered_transp(k,j)
|
||||
|
||||
overlap_gauss_r12_ao += coef * analytical_j
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
end function overlap_gauss_r12_ao
|
||||
|
||||
! ---
|
||||
|
||||
double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center, delta, i, j)
|
||||
|
||||
BEGIN_DOC
|
||||
! \int dr AO_i(r) AO_j(r) e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2}
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i, j
|
||||
double precision, intent(in) :: B_center(3), beta, D_center(3), delta
|
||||
|
||||
integer :: power_A1(3), power_A2(3), l, k
|
||||
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, analytical_j
|
||||
double precision :: G_center(3), gama, fact_g, gama_inv
|
||||
|
||||
double precision, external :: overlap_gauss_r12, overlap_gauss_r12_ao
|
||||
|
||||
ASSERT(beta .gt. 0.d0)
|
||||
if(beta .lt. 1d-10) then
|
||||
overlap_gauss_r12_ao_with1s = overlap_gauss_r12_ao(D_center, delta, i, j)
|
||||
return
|
||||
endif
|
||||
|
||||
overlap_gauss_r12_ao_with1s = 0.d0
|
||||
|
||||
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
|
||||
return
|
||||
endif
|
||||
|
||||
! e^{-beta |r-B_center^2|} e^{-delta |r-D_center|^2} = fact_g e^{-gama |r - G|^2}
|
||||
|
||||
gama = beta + delta
|
||||
gama_inv = 1.d0 / gama
|
||||
G_center(1) = (beta * B_center(1) + delta * D_center(1)) * gama_inv
|
||||
G_center(2) = (beta * B_center(2) + delta * D_center(2)) * gama_inv
|
||||
G_center(3) = (beta * B_center(3) + delta * D_center(3)) * gama_inv
|
||||
fact_g = beta * delta * gama_inv * ( (B_center(1) - D_center(1)) * (B_center(1) - D_center(1)) &
|
||||
+ (B_center(2) - D_center(2)) * (B_center(2) - D_center(2)) &
|
||||
+ (B_center(3) - D_center(3)) * (B_center(3) - D_center(3)) )
|
||||
fact_g = dexp(-fact_g)
|
||||
if(fact_g .lt. 1.d-12) return
|
||||
|
||||
! ---
|
||||
|
||||
power_A1(1:3) = ao_power(i,1:3)
|
||||
power_A2(1:3) = ao_power(j,1:3)
|
||||
|
||||
A1_center(1:3) = nucl_coord(ao_nucl(i),1:3)
|
||||
A2_center(1:3) = nucl_coord(ao_nucl(j),1:3)
|
||||
|
||||
do l = 1, ao_prim_num(i)
|
||||
alpha1 = ao_expo_ordered_transp(l,i)
|
||||
do k = 1, ao_prim_num(j)
|
||||
alpha2 = ao_expo_ordered_transp(k,j)
|
||||
coef12 = fact_g * ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(k,j)
|
||||
|
||||
if(dabs(coef12) .lt. 1d-12) cycle
|
||||
|
||||
analytical_j = overlap_gauss_r12(G_center, gama, A1_center, A2_center, power_A1, power_A2, alpha1, alpha2)
|
||||
|
||||
overlap_gauss_r12_ao_with1s += coef12 * analytical_j
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end function overlap_gauss_r12_ao_with1s
|
||||
|
||||
! ---
|
||||
|
||||
|
398
src/ao_many_one_e_ints/grad_J1b_ints.irp.f
Normal file
398
src/ao_many_one_e_ints/grad_J1b_ints.irp.f
Normal file
@ -0,0 +1,398 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ integer, List_all_comb_size]
|
||||
|
||||
implicit none
|
||||
|
||||
List_all_comb_size = 2**nucl_num
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ integer, List_all_comb, (nucl_num, List_all_comb_size)]
|
||||
|
||||
implicit none
|
||||
integer :: i, j
|
||||
|
||||
if(nucl_num .gt. 32) then
|
||||
print *, ' nucl_num = ', nucl_num, '> 32'
|
||||
stop
|
||||
endif
|
||||
|
||||
List_all_comb = 0
|
||||
|
||||
do i = 0, List_all_comb_size-1
|
||||
do j = 0, nucl_num-1
|
||||
if (btest(i,j)) then
|
||||
List_all_comb(j+1,i+1) = 1
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, List_all_j1b1s_coef, ( List_all_comb_size)]
|
||||
&BEGIN_PROVIDER [ double precision, List_all_j1b1s_expo, ( List_all_comb_size)]
|
||||
&BEGIN_PROVIDER [ double precision, List_all_j1b1s_cent, (3, List_all_comb_size)]
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, phase
|
||||
double precision :: tmp_alphaj, tmp_alphak
|
||||
|
||||
provide j1b_pen
|
||||
|
||||
List_all_j1b1s_coef = 0.d0
|
||||
List_all_j1b1s_expo = 0.d0
|
||||
List_all_j1b1s_cent = 0.d0
|
||||
|
||||
do i = 1, List_all_comb_size
|
||||
|
||||
do j = 1, nucl_num
|
||||
tmp_alphaj = dble(List_all_comb(j,i)) * j1b_pen(j)
|
||||
|
||||
List_all_j1b1s_expo(i) += tmp_alphaj
|
||||
List_all_j1b1s_cent(1,i) += tmp_alphaj * nucl_coord(j,1)
|
||||
List_all_j1b1s_cent(2,i) += tmp_alphaj * nucl_coord(j,2)
|
||||
List_all_j1b1s_cent(3,i) += tmp_alphaj * nucl_coord(j,3)
|
||||
|
||||
enddo
|
||||
|
||||
ASSERT(List_all_j1b1s_expo(i) .gt. 0d0)
|
||||
if(List_all_j1b1s_expo(i) .lt. 1d-10) cycle
|
||||
|
||||
List_all_j1b1s_cent(1,i) = List_all_j1b1s_cent(1,i) / List_all_j1b1s_expo(i)
|
||||
List_all_j1b1s_cent(2,i) = List_all_j1b1s_cent(2,i) / List_all_j1b1s_expo(i)
|
||||
List_all_j1b1s_cent(3,i) = List_all_j1b1s_cent(3,i) / List_all_j1b1s_expo(i)
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
do i = 1, List_all_comb_size
|
||||
|
||||
do j = 2, nucl_num, 1
|
||||
tmp_alphaj = dble(List_all_comb(j,i)) * j1b_pen(j)
|
||||
do k = 1, j-1, 1
|
||||
tmp_alphak = dble(List_all_comb(k,i)) * j1b_pen(k)
|
||||
|
||||
List_all_j1b1s_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
|
||||
+ (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
|
||||
+ (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if(List_all_j1b1s_expo(i) .lt. 1d-10) cycle
|
||||
|
||||
List_all_j1b1s_coef(i) = List_all_j1b1s_coef(i) / List_all_j1b1s_expo(i)
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
do i = 1, List_all_comb_size
|
||||
|
||||
phase = 0
|
||||
do j = 1, nucl_num
|
||||
phase += List_all_comb(j,i)
|
||||
enddo
|
||||
|
||||
List_all_j1b1s_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_j1b1s_coef(i))
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R|
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s
|
||||
double precision :: r(3), int_mu, int_coulomb
|
||||
double precision :: coef, beta, B_center(3)
|
||||
double precision :: wall0, wall1
|
||||
double precision, allocatable :: tmp(:,:,:)
|
||||
|
||||
double precision, external :: NAI_pol_mult_erf_ao_with1s
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
call wall_time(wall0)
|
||||
|
||||
v_ij_erf_rk_cst_mu_j1b = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_size, final_grid_points, &
|
||||
!$OMP List_all_j1b1s_coef, List_all_j1b1s_expo, List_all_j1b1s_cent, &
|
||||
!$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf)
|
||||
|
||||
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
|
||||
do i_1s = 1, List_all_comb_size
|
||||
|
||||
coef = List_all_j1b1s_coef (i_1s)
|
||||
beta = List_all_j1b1s_expo (i_1s)
|
||||
B_center(1) = List_all_j1b1s_cent(1,i_1s)
|
||||
B_center(2) = List_all_j1b1s_cent(2,i_1s)
|
||||
B_center(3) = List_all_j1b1s_cent(3,i_1s)
|
||||
|
||||
int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
|
||||
int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
|
||||
|
||||
tmp(j,i,ipoint) += coef * (int_mu - int_coulomb)
|
||||
enddo
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP CRITICAL
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END CRITICAL
|
||||
|
||||
deallocate( tmp )
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for v_ij_erf_rk_cst_mu_j1b', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid, 3)]
|
||||
|
||||
BEGIN_DOC
|
||||
! int dr x * phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: wall0, wall1
|
||||
|
||||
call wall_time(wall0)
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
! int dr x * phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s
|
||||
double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
|
||||
double precision :: wall0, wall1
|
||||
double precision, allocatable :: tmp(:,:,:,:)
|
||||
|
||||
call wall_time(wall0)
|
||||
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, tmp) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_size, final_grid_points, &
|
||||
!$OMP List_all_j1b1s_coef, List_all_j1b1s_expo, List_all_j1b1s_cent, &
|
||||
!$OMP x_v_ij_erf_rk_cst_mu_tmp_j1b, mu_erf)
|
||||
|
||||
allocate( tmp(3,ao_num,ao_num,n_points_final_grid) )
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
|
||||
do i_1s = 1, List_all_comb_size
|
||||
|
||||
coef = List_all_j1b1s_coef (i_1s)
|
||||
beta = List_all_j1b1s_expo (i_1s)
|
||||
B_center(1) = List_all_j1b1s_cent(1,i_1s)
|
||||
B_center(2) = List_all_j1b1s_cent(2,i_1s)
|
||||
B_center(3) = List_all_j1b1s_cent(3,i_1s)
|
||||
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
|
||||
|
||||
tmp(1,j,i,ipoint) += coef * (ints(1) - ints_coulomb(1))
|
||||
tmp(2,j,i,ipoint) += coef * (ints(2) - ints_coulomb(2))
|
||||
tmp(3,j,i,ipoint) += coef * (ints(3) - ints_coulomb(3))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP CRITICAL
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) += tmp(1,j,i,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) += tmp(2,j,i,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint) += tmp(3,j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END CRITICAL
|
||||
|
||||
deallocate( tmp )
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(1,i,j,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(2,i,j,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(3,i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for x_v_ij_erf_rk_cst_mu_tmp', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s, i_fit
|
||||
double precision :: r(3), int_fit, expo_fit, coef_fit
|
||||
double precision :: coef, beta, B_center(3)
|
||||
double precision :: wall0, wall1
|
||||
double precision, allocatable :: tmp(:,:,:)
|
||||
|
||||
double precision, external :: overlap_gauss_r12_ao_with1s
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
call wall_time(wall0)
|
||||
|
||||
v_ij_u_cst_mu_j1b = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_size, &
|
||||
!$OMP final_grid_points, n_max_fit_slat, &
|
||||
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
|
||||
!$OMP List_all_j1b1s_coef, List_all_j1b1s_expo, &
|
||||
!$OMP List_all_j1b1s_cent, v_ij_u_cst_mu_j1b)
|
||||
|
||||
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
|
||||
do i_1s = 1, List_all_comb_size
|
||||
|
||||
coef = List_all_j1b1s_coef (i_1s)
|
||||
beta = List_all_j1b1s_expo (i_1s)
|
||||
B_center(1) = List_all_j1b1s_cent(1,i_1s)
|
||||
B_center(2) = List_all_j1b1s_cent(2,i_1s)
|
||||
B_center(3) = List_all_j1b1s_cent(3,i_1s)
|
||||
|
||||
do i_fit = 1, n_max_fit_slat
|
||||
|
||||
expo_fit = expo_gauss_j_mu_x(i_fit)
|
||||
coef_fit = coef_gauss_j_mu_x(i_fit)
|
||||
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
|
||||
tmp(j,i,ipoint) += coef * coef_fit * int_fit
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP CRITICAL
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
v_ij_u_cst_mu_j1b(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END CRITICAL
|
||||
|
||||
deallocate( tmp )
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
v_ij_u_cst_mu_j1b(j,i,ipoint) = v_ij_u_cst_mu_j1b(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for v_ij_u_cst_mu_j1b', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
@ -1,29 +1,43 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points_final_grid)]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1) / |r - R|
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: mu,r(3),NAI_pol_mult_erf_ao
|
||||
double precision :: r(3)
|
||||
double precision :: int_mu, int_coulomb
|
||||
provide mu_erf final_grid_points
|
||||
double precision :: wall0, wall1
|
||||
|
||||
double precision :: NAI_pol_mult_erf_ao
|
||||
|
||||
provide mu_erf final_grid_points
|
||||
call wall_time(wall0)
|
||||
|
||||
v_ij_erf_rk_cst_mu = 0.d0
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,ipoint,mu,r,int_mu,int_coulomb) &
|
||||
!$OMP PRIVATE (i, j, ipoint, r, int_mu, int_coulomb) &
|
||||
!$OMP SHARED (ao_num, n_points_final_grid, v_ij_erf_rk_cst_mu, final_grid_points, mu_erf)
|
||||
!$OMP DO SCHEDULE (dynamic)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
mu = mu_erf
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
int_mu = NAI_pol_mult_erf_ao(i,j,mu,r)
|
||||
|
||||
int_mu = NAI_pol_mult_erf_ao(i, j, mu_erf, r)
|
||||
int_coulomb = NAI_pol_mult_erf_ao(i, j, 1.d+9, r)
|
||||
v_ij_erf_rk_cst_mu(j,i,ipoint)= (int_mu - int_coulomb )
|
||||
|
||||
v_ij_erf_rk_cst_mu(j,i,ipoint) = int_mu - int_coulomb
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -40,8 +54,11 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, ( ao_num, ao_num,n_points
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, 'wall time for v_ij_erf_rk_cst_mu ', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_transp, (n_points_final_grid, ao_num, ao_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -86,30 +103,36 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_transp, (n_points_final_gr
|
||||
print*,'wall time for v_ij_erf_rk_cst_mu_transp ',wall1 - wall0
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3, ao_num, ao_num, n_points_final_grid)]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R|
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, m
|
||||
double precision :: mu,r(3),ints(3),ints_coulomb(3)
|
||||
double precision :: r(3), ints(3), ints_coulomb(3)
|
||||
double precision :: wall0, wall1
|
||||
|
||||
call wall_time(wall0)
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,ipoint,mu,r,ints,m,ints_coulomb) &
|
||||
!$OMP PRIVATE (i,j,ipoint,r,ints,m,ints_coulomb) &
|
||||
!$OMP SHARED (ao_num,n_points_final_grid,x_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf)
|
||||
!$OMP DO SCHEDULE (dynamic)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
mu = mu_erf
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
call NAI_pol_x_mult_erf_ao(i,j,mu,r,ints)
|
||||
|
||||
call NAI_pol_x_mult_erf_ao(i, j, mu_erf, r, ints )
|
||||
call NAI_pol_x_mult_erf_ao(i, j, 1.d+9 , r, ints_coulomb)
|
||||
|
||||
do m = 1, 3
|
||||
x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) = (ints(m) - ints_coulomb(m))
|
||||
enddo
|
||||
@ -128,12 +151,14 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3,ao_num, ao_num,n
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*,'wall time for x_v_ij_erf_rk_cst_mu_tmp',wall1 - wall0
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu, (ao_num, ao_num,n_points_final_grid,3)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
@ -1,3 +1,6 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -15,41 +18,108 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
|
||||
enddo
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
double precision function NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center)
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
|
||||
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in
|
||||
double precision :: A_center(3), B_center(3), integral, alpha, beta
|
||||
|
||||
double precision :: NAI_pol_mult_erf
|
||||
|
||||
num_A = ao_nucl(i_ao)
|
||||
power_A(1:3) = ao_power(i_ao,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
num_B = ao_nucl(j_ao)
|
||||
power_B(1:3) = ao_power(j_ao,1:3)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
NAI_pol_mult_erf_ao = 0.d0
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alpha = ao_expo_ordered_transp(i,i_ao)
|
||||
do j = 1, ao_prim_num(j_ao)
|
||||
beta = ao_expo_ordered_transp(j,j_ao)
|
||||
|
||||
integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in,mu_in)
|
||||
|
||||
NAI_pol_mult_erf_ao += integral * ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
end function NAI_pol_mult_erf_ao
|
||||
|
||||
! ---
|
||||
|
||||
double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: beta, B_center(3)
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
|
||||
integer :: i, j, power_A1(3), power_A2(3), n_pt_in
|
||||
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, integral
|
||||
|
||||
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao
|
||||
|
||||
ASSERT(beta .lt. 0.d0)
|
||||
if(beta .lt. 1d-10) then
|
||||
NAI_pol_mult_erf_ao_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center)
|
||||
return
|
||||
endif
|
||||
|
||||
power_A1(1:3) = ao_power(i_ao,1:3)
|
||||
power_A2(1:3) = ao_power(j_ao,1:3)
|
||||
|
||||
A1_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
|
||||
A2_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s = 0.d0
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alpha1 = ao_expo_ordered_transp(i,i_ao)
|
||||
do j = 1, ao_prim_num(j_ao)
|
||||
alpha2 = ao_expo_ordered_transp(j,j_ao)
|
||||
|
||||
coef12 = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
if(coef12 .lt. 1d-14) cycle
|
||||
|
||||
integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
|
||||
, beta, B_center, C_center, n_pt_in, mu_in )
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s += integral * coef12
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end function NAI_pol_mult_erf_ao_with1s
|
||||
|
||||
! ---
|
||||
|
||||
double precision function NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math::
|
||||
@ -59,19 +129,22 @@ double precision function NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alp
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta,mu_in
|
||||
integer, intent(in) :: power_A(3), power_B(3)
|
||||
integer :: i,j,k,l,n_pt
|
||||
double precision :: P_center(3)
|
||||
double precision, intent(in) :: C_center(3), A_center(3), B_center(3), alpha, beta, mu_in
|
||||
|
||||
integer :: i, n_pt, n_pt_out
|
||||
double precision :: P_center(3)
|
||||
double precision :: d(0:n_pt_in), coeff, dist, const, factor
|
||||
double precision :: const_factor, dist_integral
|
||||
double precision :: accu, p_inv, p, rho, p_inv_2
|
||||
double precision :: p_new
|
||||
|
||||
double precision :: rint
|
||||
|
||||
double precision :: d(0:n_pt_in),pouet,coeff,dist,const,pouet_2,factor
|
||||
double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi
|
||||
double precision :: V_e_n,const_factor,dist_integral,tmp
|
||||
double precision :: accu,rint,p_inv,p,rho,p_inv_2
|
||||
integer :: n_pt_out,lmax
|
||||
include 'utils/constants.include.F'
|
||||
p = alpha + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
@ -89,68 +162,184 @@ double precision function NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alp
|
||||
NAI_pol_mult_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
double precision :: p_new
|
||||
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
lmax = 20
|
||||
|
||||
! print*, "b"
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) )
|
||||
const = p * dist_integral * p_new * p_new
|
||||
if(n_pt == 0) then
|
||||
pouet = rint(0,const)
|
||||
NAI_pol_mult_erf = coeff * pouet
|
||||
NAI_pol_mult_erf = coeff * rint(0, const)
|
||||
return
|
||||
endif
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
! call give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
|
||||
p_new = p_new * p_new
|
||||
call give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center)
|
||||
|
||||
call give_polynomial_mult_center_one_e_erf_opt( A_center, B_center, power_A, power_B, C_center &
|
||||
, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
NAI_pol_mult_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
accu = 0.d0
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
NAI_pol_mult_erf = accu * coeff
|
||||
|
||||
end
|
||||
end function NAI_pol_mult_erf
|
||||
|
||||
! ---
|
||||
|
||||
double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
|
||||
, beta, B_center, C_center, n_pt_in, mu_in )
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math::
|
||||
!
|
||||
! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2)
|
||||
! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2)
|
||||
! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2)
|
||||
! \exp(-\beta (r - B)^2)
|
||||
! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer, intent(in) :: power_A1(3), power_A2(3)
|
||||
double precision, intent(in) :: C_center(3), A1_center(3), A2_center(3), B_center(3)
|
||||
double precision, intent(in) :: alpha1, alpha2, beta, mu_in
|
||||
|
||||
integer :: i, n_pt, n_pt_out
|
||||
double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12
|
||||
double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor
|
||||
double precision :: dist_integral
|
||||
double precision :: d(0:n_pt_in), coeff, const, factor
|
||||
double precision :: accu
|
||||
double precision :: p_new
|
||||
|
||||
double precision :: rint
|
||||
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,beta,&
|
||||
power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in,p,p_inv,p_inv_2,p_new,P_center)
|
||||
! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{K12} e^{-alpha12 (r - A12)^2}
|
||||
alpha12 = alpha1 + alpha2
|
||||
alpha12_inv = 1.d0 / alpha12
|
||||
alpha12_inv_2 = 0.5d0 * alpha12_inv
|
||||
rho12 = alpha1 * alpha2 * alpha12_inv
|
||||
|
||||
dist12 = 0.d0
|
||||
do i = 1, 3
|
||||
A12_center(i) = (alpha1 * A1_center(i) + alpha2 * A2_center(i)) * alpha12_inv
|
||||
dist12 += (A1_center(i) - A2_center(i)) * (A1_center(i) - A2_center(i))
|
||||
enddo
|
||||
|
||||
const_factor12 = dist12 * rho12
|
||||
if(const_factor12 > 80.d0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
! e^{K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{K} e^{-p (r - P)^2}
|
||||
p = alpha12 + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha12 * beta * p_inv
|
||||
|
||||
dist = 0.d0
|
||||
dist_integral = 0.d0
|
||||
do i = 1, 3
|
||||
P_center(i) = (alpha12 * A12_center(i) + beta * B_center(i)) * p_inv
|
||||
dist += (A12_center(i) - B_center(i)) * (A12_center(i) - B_center(i))
|
||||
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
|
||||
enddo
|
||||
|
||||
const_factor = const_factor12 + dist * rho
|
||||
if(const_factor > 80.d0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
|
||||
n_pt = 2 * ( (power_A1(1) + power_A2(1)) + (power_A1(2) + power_A2(2)) + (power_A1(3) + power_A2(3)) )
|
||||
const = p * dist_integral * p_new * p_new
|
||||
if(n_pt == 0) then
|
||||
NAI_pol_mult_erf_with1s = coeff * rint(0, const)
|
||||
return
|
||||
endif
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
p_new = p_new * p_new
|
||||
|
||||
call give_polynomial_mult_center_one_e_erf_opt( A1_center, A2_center, power_A1, power_A2, C_center &
|
||||
, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
NAI_pol_mult_erf_with1s = accu * coeff
|
||||
|
||||
end function NAI_pol_mult_erf_with1s
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf_opt( A_center, B_center, power_A, power_B, C_center &
|
||||
, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
|
||||
|
||||
BEGIN_DOC
|
||||
! Returns the explicit polynomial in terms of the $t$ variable of the
|
||||
! following polynomial:
|
||||
!
|
||||
! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$.
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer,intent(out) :: n_pt_out
|
||||
double precision, intent(in) :: A_center(3), B_center(3),C_center(3),p,p_inv,p_inv_2,p_new,P_center(3)
|
||||
double precision, intent(in) :: alpha,beta,mu_in
|
||||
integer, intent(in) :: power_A(3), power_B(3)
|
||||
double precision, intent(in) :: A_center(3), B_center(3), C_center(3), p_inv_2, p_new, P_center(3)
|
||||
integer, intent(out) :: n_pt_out
|
||||
double precision, intent(out) :: d(0:n_pt_in)
|
||||
|
||||
integer :: a_x, b_x, a_y, b_y, a_z, b_z
|
||||
double precision :: d(0:n_pt_in)
|
||||
integer :: n_pt1, n_pt2, n_pt3, dim, i
|
||||
integer :: n_pt_tmp
|
||||
double precision :: d1(0:n_pt_in)
|
||||
double precision :: d2(0:n_pt_in)
|
||||
double precision :: d3(0:n_pt_in)
|
||||
double precision :: accu
|
||||
double precision :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2)
|
||||
|
||||
accu = 0.d0
|
||||
ASSERT (n_pt_in > 1)
|
||||
|
||||
double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
|
||||
R1x(0) = (P_center(1) - A_center(1))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(1) - C_center(1))* p_new
|
||||
@ -163,19 +352,14 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet
|
||||
R2x(1) = 0.d0
|
||||
R2x(2) = -p_inv_2 * p_new
|
||||
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d (i) = 0.d0
|
||||
enddo
|
||||
do i = 0,n_pt_in
|
||||
d1(i) = 0.d0
|
||||
enddo
|
||||
do i = 0,n_pt_in
|
||||
d2(i) = 0.d0
|
||||
enddo
|
||||
do i = 0,n_pt_in
|
||||
d3(i) = 0.d0
|
||||
enddo
|
||||
integer :: n_pt1,n_pt2,n_pt3,dim,i
|
||||
|
||||
n_pt1 = n_pt_in
|
||||
n_pt2 = n_pt_in
|
||||
n_pt3 = n_pt_in
|
||||
@ -209,7 +393,6 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet
|
||||
return
|
||||
endif
|
||||
|
||||
|
||||
R1x(0) = (P_center(3) - A_center(3))
|
||||
R1x(1) = 0.d0
|
||||
R1x(2) = -(P_center(3) - C_center(3)) * p_new
|
||||
@ -229,7 +412,7 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet
|
||||
enddo
|
||||
return
|
||||
endif
|
||||
integer :: n_pt_tmp
|
||||
|
||||
n_pt_tmp = 0
|
||||
call multiply_poly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp)
|
||||
do i = 0, n_pt_tmp
|
||||
@ -241,9 +424,9 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center,B_center,alpha,bet
|
||||
d(i) = d1(i)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
end subroutine give_polynomial_mult_center_one_e_erf_opt
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,&
|
||||
|
@ -12,10 +12,13 @@ BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, expo_gauss_j_mu_x, (n_max_fit_slat)]
|
||||
&BEGIN_PROVIDER [double precision, coef_gauss_j_mu_x, (n_max_fit_slat)]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! J(mu,r12) = 1/2 r12 * (1 - erf(mu*r12)) - 1/(2 sqrt(pi)*mu) exp(-(mu*r12)^2) is expressed as
|
||||
!
|
||||
! J(mu,r12) = 0.5/mu * F(r12*mu) where F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)
|
||||
@ -27,18 +30,27 @@ END_PROVIDER
|
||||
! See Appendix 2 of JCP 154, 084119 (2021)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: tmp
|
||||
double precision :: expos(n_max_fit_slat), alpha, beta
|
||||
|
||||
tmp = -0.5d0 / (mu_erf * sqrt(dacos(-1.d0)))
|
||||
|
||||
alpha = expo_j_xmu(1) * mu_erf
|
||||
call expo_fit_slater_gam(alpha, expos)
|
||||
beta = expo_j_xmu(2) * mu_erf**2.d0
|
||||
beta = expo_j_xmu(2) * mu_erf * mu_erf
|
||||
|
||||
do i = 1, n_max_fit_slat
|
||||
expo_gauss_j_mu_x(i) = expos(i) + beta
|
||||
coef_gauss_j_mu_x(i) = coef_fit_slat_gauss(i) / (2.d0 * mu_erf) * (- 1/dsqrt(dacos(-1.d0)))
|
||||
coef_gauss_j_mu_x(i) = tmp * coef_fit_slat_gauss(i)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
double precision function F_x_j(x)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -89,3 +101,6 @@ double precision function j_mu_fit_gauss(x)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -1,30 +1,84 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num,n_points_final_grid, 3)]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! grad_1_u_ij_mu(i,j,ipoint) = -1 * \int dr2 \grad_r1 u(r1,r2) \phi_i(r2) \phi_j(r2)
|
||||
!
|
||||
! grad_1_u_ij_mu(i,j,ipoint) = \int dr2 [-1 * \grad_r1 u(r1,r2)] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2)
|
||||
! = \int dr2 [(r1 - r2) (erf(mu * r12)-1)/2 r_12] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2)
|
||||
!
|
||||
! where r1 = r(ipoint)
|
||||
!
|
||||
! grad_1_u_ij_mu(i,j,ipoint) = \int dr2 (r1 - r2) (erf(mu * r12)-1)/2 r_12 \phi_i(r2) \phi_j(r2)
|
||||
END_DOC
|
||||
integer :: ipoint,i,j,m
|
||||
|
||||
implicit none
|
||||
integer :: ipoint, i, j, i_1s
|
||||
double precision :: r(3)
|
||||
do m = 1, 3
|
||||
double precision :: a, d, e, tmp1, tmp2, fact_r, fact_xyz(3)
|
||||
|
||||
PROVIDE j1b_type j1b_pen
|
||||
|
||||
if(j1b_type .eq. 3) then
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
|
||||
fact_r = 1.d0
|
||||
fact_xyz(1) = 1.d0
|
||||
fact_xyz(2) = 1.d0
|
||||
fact_xyz(3) = 1.d0
|
||||
do i_1s = 1, nucl_num
|
||||
a = j1b_pen(i_1s)
|
||||
d = (r(1) - nucl_coord(i_1s,1)) * (r(1) - nucl_coord(i_1s,1)) &
|
||||
+ (r(2) - nucl_coord(i_1s,2)) * (r(2) - nucl_coord(i_1s,2)) &
|
||||
+ (r(3) - nucl_coord(i_1s,3)) * (r(3) - nucl_coord(i_1s,3))
|
||||
e = dexp(-a*d)
|
||||
|
||||
fact_r = fact_r * (1.d0 - e)
|
||||
fact_xyz(1) = fact_xyz(1) * (2.d0 * a * (r(1) - nucl_coord(i_1s,1)) * e)
|
||||
fact_xyz(2) = fact_xyz(2) * (2.d0 * a * (r(2) - nucl_coord(i_1s,2)) * e)
|
||||
fact_xyz(3) = fact_xyz(3) * (2.d0 * a * (r(3) - nucl_coord(i_1s,3)) * e)
|
||||
enddo
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
tmp1 = fact_r * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
|
||||
tmp2 = v_ij_u_cst_mu_j1b (i,j,ipoint)
|
||||
|
||||
grad_1_u_ij_mu(i,j,ipoint,1) = tmp1 * r(1) - fact_r * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) + fact_xyz(1) * tmp2
|
||||
grad_1_u_ij_mu(i,j,ipoint,2) = tmp1 * r(2) - fact_r * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) + fact_xyz(2) * tmp2
|
||||
grad_1_u_ij_mu(i,j,ipoint,3) = tmp1 * r(3) - fact_r * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) + fact_xyz(3) * tmp2
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
grad_1_u_ij_mu(i,j,ipoint,m) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(m) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,m)
|
||||
enddo
|
||||
grad_1_u_ij_mu(i,j,ipoint,1) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(1) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)
|
||||
grad_1_u_ij_mu(i,j,ipoint,2) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(2) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)
|
||||
grad_1_u_ij_mu(i,j,ipoint,3) = v_ij_erf_rk_cst_mu(i,j,ipoint) * r(3) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
grad_1_u_ij_mu *= 0.5d0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
Loading…
Reference in New Issue
Block a user