mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-10-05 15:45:57 +02:00
merged with Abdallah
This commit is contained in:
commit
0b88c8d3c0
2
external/ezfio
vendored
2
external/ezfio
vendored
@ -1 +1 @@
|
||||
Subproject commit 0520b5e2cf70e2451c37ce5b7f2f64f6d2e5e956
|
||||
Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93
|
@ -65,46 +65,60 @@ double precision function primitive_value(i,j,r)
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_all_aos_at_r(r,aos_array)
|
||||
implicit none
|
||||
BEGIN_dOC
|
||||
! input : r == r(1) = x and so on
|
||||
!
|
||||
! output : aos_array(i) = aos(i) evaluated in $\textbf{r}$
|
||||
END_DOC
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out):: aos_array(ao_num)
|
||||
subroutine give_all_aos_at_r(r, tmp_array)
|
||||
|
||||
integer :: power_ao(3)
|
||||
integer :: i,j,k,l,m
|
||||
double precision :: dx,dy,dz,r2
|
||||
double precision :: dx2,dy2,dz2
|
||||
double precision :: center_ao(3)
|
||||
double precision :: beta
|
||||
do i = 1, nucl_num
|
||||
center_ao(1:3) = nucl_coord(i,1:3)
|
||||
dx = (r(1) - center_ao(1))
|
||||
dy = (r(2) - center_ao(2))
|
||||
dz = (r(3) - center_ao(3))
|
||||
r2 = dx*dx + dy*dy + dz*dz
|
||||
do j = 1,Nucl_N_Aos(i)
|
||||
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
|
||||
aos_array(k) = 0.d0
|
||||
power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
|
||||
dx2 = dx**power_ao(1)
|
||||
dy2 = dy**power_ao(2)
|
||||
dz2 = dz**power_ao(3)
|
||||
do l = 1,ao_prim_num(k)
|
||||
beta = ao_expo_ordered_transp_per_nucl(l,j,i)
|
||||
if(dabs(beta*r2).gt.40.d0)cycle
|
||||
aos_array(k)+= ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
|
||||
enddo
|
||||
aos_array(k) = aos_array(k) * dx2 * dy2 * dz2
|
||||
BEGIN_dOC
|
||||
!
|
||||
! input : r == r(1) = x and so on
|
||||
!
|
||||
! output : tmp_array(i) = aos(i) evaluated in $\textbf{r}$
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: tmp_array(ao_num)
|
||||
integer :: p_ao(3)
|
||||
integer :: i, j, k, l, m
|
||||
double precision :: dx, dy, dz, r2
|
||||
double precision :: dx2, dy2, dz2
|
||||
double precision :: c_ao(3)
|
||||
double precision :: beta
|
||||
|
||||
do i = 1, nucl_num
|
||||
|
||||
c_ao(1:3) = nucl_coord(i,1:3)
|
||||
dx = r(1) - c_ao(1)
|
||||
dy = r(2) - c_ao(2)
|
||||
dz = r(3) - c_ao(3)
|
||||
r2 = dx*dx + dy*dy + dz*dz
|
||||
|
||||
do j = 1, Nucl_N_Aos(i)
|
||||
|
||||
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
|
||||
p_ao(1:3) = ao_power_ordered_transp_per_nucl(1:3,j,i)
|
||||
dx2 = dx**p_ao(1)
|
||||
dy2 = dy**p_ao(2)
|
||||
dz2 = dz**p_ao(3)
|
||||
|
||||
tmp_array(k) = 0.d0
|
||||
do l = 1,ao_prim_num(k)
|
||||
beta = ao_expo_ordered_transp_per_nucl(l,j,i)
|
||||
if(dabs(beta*r2).gt.40.d0) cycle
|
||||
|
||||
tmp_array(k) += ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
|
||||
enddo
|
||||
|
||||
tmp_array(k) = tmp_array(k) * dx2 * dy2 * dz2
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
|
||||
implicit none
|
||||
|
@ -1,20 +1,28 @@
|
||||
BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! List of AOs attached on each atom
|
||||
END_DOC
|
||||
integer :: i
|
||||
integer, allocatable :: nucl_tmp(:)
|
||||
allocate(nucl_tmp(nucl_num))
|
||||
nucl_tmp = 0
|
||||
Nucl_Aos = 0
|
||||
do i = 1, ao_num
|
||||
nucl_tmp(ao_nucl(i))+=1
|
||||
Nucl_Aos_transposed(nucl_tmp(ao_nucl(i)),ao_nucl(i)) = i
|
||||
enddo
|
||||
deallocate(nucl_tmp)
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
! List of AOs attached on each atom
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i
|
||||
integer, allocatable :: nucl_tmp(:)
|
||||
|
||||
allocate(nucl_tmp(nucl_num))
|
||||
nucl_tmp = 0
|
||||
do i = 1, ao_num
|
||||
nucl_tmp(ao_nucl(i)) += 1
|
||||
Nucl_Aos_transposed(nucl_tmp(ao_nucl(i)),ao_nucl(i)) = i
|
||||
enddo
|
||||
deallocate(nucl_tmp)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_expo_ordered_transp_per_nucl, (ao_prim_num_max,N_AOs_max,nucl_num) ]
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
|
@ -212,9 +212,7 @@ subroutine NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
|
||||
! 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
|
||||
@ -279,9 +277,7 @@ subroutine NAI_pol_x_mult_erf_ao_v0(i_ao, j_ao, mu_in, C_center, LD_C, ints, LD_
|
||||
! 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
|
||||
@ -1111,3 +1107,141 @@ end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine NAI_pol_x2_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^2 * \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^2 * \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^2 * \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, m
|
||||
integer :: power_A1(3), power_A2(3)
|
||||
double precision :: Ai_center(3), Aj_center(3), alphai, alphaj, coef, coefi
|
||||
double precision :: integral0, integral1, integral2
|
||||
|
||||
double precision, external :: NAI_pol_mult_erf_with1s
|
||||
|
||||
ASSERT(beta .ge. 0.d0)
|
||||
if(beta .lt. 1d-10) then
|
||||
call NAI_pol_x2_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
|
||||
return
|
||||
endif
|
||||
|
||||
ints = 0.d0
|
||||
|
||||
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)
|
||||
coefi = ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
|
||||
do m = 1, 3
|
||||
|
||||
power_A1 = power_Ai
|
||||
power_A1(m) += 1
|
||||
|
||||
power_A2 = power_Ai
|
||||
power_A2(m) += 2
|
||||
|
||||
do j = 1, ao_prim_num(j_ao)
|
||||
alphaj = ao_expo_ordered_transp (j,j_ao)
|
||||
coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao)
|
||||
|
||||
integral0 = 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)
|
||||
integral1 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A1, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in)
|
||||
integral2 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A2, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in)
|
||||
|
||||
ints(m) += coef * (integral2 + Ai_center(m) * (2.d0*integral1 + Ai_center(m)*integral0))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine NAI_pol_x2_mult_erf_ao_with1s
|
||||
|
||||
! ---
|
||||
|
||||
subroutine NAI_pol_x2_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
! $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
! $\int_{-\infty}^{infty} dr z^2 * \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)
|
||||
|
||||
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, m
|
||||
integer :: power_A1(3), power_A2(3)
|
||||
double precision :: A_center(3), B_center(3), alpha, beta, coef
|
||||
double precision :: integral0, integral1, integral2
|
||||
|
||||
double precision :: NAI_pol_mult_erf
|
||||
|
||||
ints = 0.d0
|
||||
|
||||
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
|
||||
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alpha = ao_expo_ordered_transp(i,i_ao)
|
||||
|
||||
do m = 1, 3
|
||||
|
||||
power_A1 = power_A
|
||||
power_A1(m) += 1
|
||||
|
||||
power_A2 = power_A
|
||||
power_A2(m) += 2
|
||||
|
||||
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)
|
||||
|
||||
integral0 = NAI_pol_mult_erf(A_center, B_center, power_A , power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
integral1 = NAI_pol_mult_erf(A_center, B_center, power_A1, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
integral2 = NAI_pol_mult_erf(A_center, B_center, power_A2, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
|
||||
ints(m) += coef * (integral2 + A_center(m) * (2.d0*integral1 + A_center(m)*integral0))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine NAI_pol_x2_mult_erf_ao
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -128,6 +128,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
|
||||
do i_1s = 2, List_all_comb_b3_size
|
||||
|
||||
coef = List_all_comb_b3_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b3_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b3_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b3_cent(2,i_1s)
|
||||
@ -222,6 +223,7 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
|
||||
do i_1s = 2, List_all_comb_b3_size
|
||||
|
||||
coef = List_all_comb_b3_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b3_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b3_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b3_cent(2,i_1s)
|
||||
@ -322,6 +324,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
|
||||
do i_1s = 2, List_all_comb_b3_size
|
||||
|
||||
coef = List_all_comb_b3_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b3_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b3_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b3_cent(2,i_1s)
|
||||
@ -436,6 +439,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
|
||||
do i_1s = 2, List_all_comb_b3_size
|
||||
|
||||
coef = List_all_comb_b3_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b3_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b3_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b3_cent(2,i_1s)
|
||||
|
@ -60,6 +60,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
|
||||
do i_1s = 2, List_all_comb_b2_size
|
||||
|
||||
coef = List_all_comb_b2_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b2_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b2_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b2_cent(2,i_1s)
|
||||
@ -154,6 +155,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
|
||||
do i_1s = 2, List_all_comb_b2_size
|
||||
|
||||
coef = List_all_comb_b2_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b2_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b2_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b2_cent(2,i_1s)
|
||||
@ -195,8 +197,7 @@ END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
! TODO analytically
|
||||
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
|
||||
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
@ -213,12 +214,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
|
||||
double precision, external :: overlap_gauss_r12_ao_with1s
|
||||
|
||||
print*, ' providing v_ij_u_cst_mu_j1b ...'
|
||||
print*, ' providing v_ij_u_cst_mu_j1b_fit ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x
|
||||
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
|
||||
|
||||
v_ij_u_cst_mu_j1b = 0.d0
|
||||
v_ij_u_cst_mu_j1b_fit = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
@ -227,9 +230,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
!$OMP final_grid_points, ng_fit_jast, &
|
||||
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
|
||||
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
|
||||
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b)
|
||||
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_fit)
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
@ -240,7 +242,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
|
||||
tmp = 0.d0
|
||||
do i_fit = 1, ng_fit_jast
|
||||
|
||||
expo_fit = expo_gauss_j_mu_x(i_fit)
|
||||
coef_fit = coef_gauss_j_mu_x(i_fit)
|
||||
|
||||
@ -253,7 +254,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
B_center(3) = List_all_comb_b2_cent(3,1)
|
||||
|
||||
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
! if(dabs(int_fit*coef) .lt. 1d-12) cycle
|
||||
|
||||
tmp += coef * coef_fit * int_fit
|
||||
|
||||
@ -262,6 +262,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
do i_1s = 2, List_all_comb_b2_size
|
||||
|
||||
coef = List_all_comb_b2_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b2_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b2_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b2_cent(2,i_1s)
|
||||
@ -276,7 +277,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
|
||||
enddo
|
||||
|
||||
v_ij_u_cst_mu_j1b(j,i,ipoint) = tmp
|
||||
v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = tmp
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -286,13 +287,149 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 2, 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)
|
||||
v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = v_ij_u_cst_mu_j1b_fit(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for v_ij_u_cst_mu_j1b', wall1 - wall0
|
||||
print*, ' wall time for v_ij_u_cst_mu_j1b_fit', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
|
||||
!
|
||||
! TODO
|
||||
! one subroutine for all integrals
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s
|
||||
double precision :: r(3), r1_2
|
||||
double precision :: int_c1, int_e1, int_o
|
||||
double precision :: int_c2(3), int_e2(3)
|
||||
double precision :: int_c3(3), int_e3(3)
|
||||
double precision :: coef, beta, B_center(3)
|
||||
double precision :: tmp, ct
|
||||
double precision :: wall0, wall1
|
||||
|
||||
double precision, external :: overlap_gauss_r12_ao_with1s
|
||||
double precision, external :: NAI_pol_mult_erf_ao_with1s
|
||||
|
||||
print*, ' providing v_ij_u_cst_mu_j1b_an ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
|
||||
|
||||
ct = inv_sq_pi_2 / mu_erf
|
||||
|
||||
v_ij_u_cst_mu_j1b_an = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
|
||||
!$OMP r1_2, tmp, int_c1, int_e1, int_o, int_c2, &
|
||||
!$OMP int_e2, int_c3, int_e3) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, &
|
||||
!$OMP final_grid_points, mu_erf, ct, &
|
||||
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
|
||||
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an)
|
||||
!$OMP DO
|
||||
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)
|
||||
r1_2 = 0.5d0 * (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
|
||||
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
|
||||
! ---
|
||||
|
||||
coef = List_all_comb_b2_coef (1)
|
||||
beta = List_all_comb_b2_expo (1)
|
||||
B_center(1) = List_all_comb_b2_cent(1,1)
|
||||
B_center(2) = List_all_comb_b2_cent(2,1)
|
||||
B_center(3) = List_all_comb_b2_cent(3,1)
|
||||
|
||||
int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
|
||||
int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
|
||||
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c2)
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e2)
|
||||
|
||||
call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c3)
|
||||
call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e3)
|
||||
|
||||
int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j)
|
||||
|
||||
tmp = coef &
|
||||
* ( r1_2 * (int_c1 - int_e1) &
|
||||
- r(1) * (int_c2(1) - int_e2(1)) - r(2) * (int_c2(2) - int_e2(2)) - r(3) * (int_c2(3) - int_e2(3)) &
|
||||
+ 0.5d0 * (int_c3(1) + int_c3(2) + int_c3(3) - int_e3(1) - int_e3(2) - int_e3(3)) &
|
||||
- ct * int_o &
|
||||
)
|
||||
|
||||
! ---
|
||||
|
||||
do i_1s = 2, List_all_comb_b2_size
|
||||
|
||||
coef = List_all_comb_b2_coef (i_1s)
|
||||
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
|
||||
beta = List_all_comb_b2_expo (i_1s)
|
||||
B_center(1) = List_all_comb_b2_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b2_cent(2,i_1s)
|
||||
B_center(3) = List_all_comb_b2_cent(3,i_1s)
|
||||
|
||||
int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
|
||||
int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
|
||||
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c2)
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e2)
|
||||
|
||||
call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c3)
|
||||
call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e3)
|
||||
|
||||
int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j)
|
||||
|
||||
tmp = tmp + coef &
|
||||
* ( r1_2 * (int_c1 - int_e1) &
|
||||
- r(1) * (int_c2(1) - int_e2(1)) - r(2) * (int_c2(2) - int_e2(2)) - r(3) * (int_c2(3) - int_e2(3)) &
|
||||
+ 0.5d0 * (int_c3(1) + int_c3(2) + int_c3(3) - int_e3(1) - int_e3(2) - int_e3(3)) &
|
||||
- ct * int_o &
|
||||
)
|
||||
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
v_ij_u_cst_mu_j1b_an(j,i,ipoint) = tmp
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 2, ao_num
|
||||
do j = 1, i-1
|
||||
v_ij_u_cst_mu_j1b_an(j,i,ipoint) = v_ij_u_cst_mu_j1b_an(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for v_ij_u_cst_mu_j1b_an', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -62,6 +62,7 @@ END_PROVIDER
|
||||
double precision :: tmp_cent_x, tmp_cent_y, tmp_cent_z
|
||||
|
||||
provide j1b_pen
|
||||
provide j1b_pen_coef
|
||||
|
||||
List_all_comb_b2_coef = 0.d0
|
||||
List_all_comb_b2_expo = 0.d0
|
||||
@ -127,8 +128,8 @@ END_PROVIDER
|
||||
List_all_comb_b2_expo( 1) = 0.d0
|
||||
List_all_comb_b2_cent(1:3,1) = 0.d0
|
||||
do i = 1, nucl_num
|
||||
List_all_comb_b2_coef( i+1) = -1.d0
|
||||
List_all_comb_b2_expo( i+1) = j1b_pen( i)
|
||||
List_all_comb_b2_coef( i+1) = -1.d0 * j1b_pen_coef(i)
|
||||
List_all_comb_b2_expo( i+1) = j1b_pen(i)
|
||||
List_all_comb_b2_cent(1,i+1) = nucl_coord(i,1)
|
||||
List_all_comb_b2_cent(2,i+1) = nucl_coord(i,2)
|
||||
List_all_comb_b2_cent(3,i+1) = nucl_coord(i,3)
|
||||
@ -225,6 +226,7 @@ END_PROVIDER
|
||||
double precision :: dx, dy, dz, r2
|
||||
|
||||
provide j1b_pen
|
||||
provide j1b_pen_coef
|
||||
|
||||
List_all_comb_b3_coef = 0.d0
|
||||
List_all_comb_b3_expo = 0.d0
|
||||
@ -296,8 +298,8 @@ END_PROVIDER
|
||||
|
||||
do i = 1, nucl_num
|
||||
ii = ii + 1
|
||||
List_all_comb_b3_coef( ii) = -2.d0
|
||||
List_all_comb_b3_expo( ii) = j1b_pen( i)
|
||||
List_all_comb_b3_coef( ii) = -2.d0 * j1b_pen_coef(i)
|
||||
List_all_comb_b3_expo( ii) = j1b_pen(i)
|
||||
List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
|
||||
List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
|
||||
List_all_comb_b3_cent(3,ii) = nucl_coord(i,3)
|
||||
@ -305,7 +307,7 @@ END_PROVIDER
|
||||
|
||||
do i = 1, nucl_num
|
||||
ii = ii + 1
|
||||
List_all_comb_b3_coef( ii) = 1.d0
|
||||
List_all_comb_b3_coef( ii) = 1.d0 * j1b_pen_coef(i) * j1b_pen_coef(i)
|
||||
List_all_comb_b3_expo( ii) = 2.d0 * j1b_pen(i)
|
||||
List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
|
||||
List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
|
||||
@ -337,7 +339,7 @@ END_PROVIDER
|
||||
|
||||
ii = ii + 1
|
||||
! x 2 to avoid doing integrals twice
|
||||
List_all_comb_b3_coef( ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2)
|
||||
List_all_comb_b3_coef( ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2) * j1b_pen_coef(i) * j1b_pen_coef(j)
|
||||
List_all_comb_b3_expo( ii) = tmp3
|
||||
List_all_comb_b3_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj)
|
||||
List_all_comb_b3_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj)
|
||||
|
@ -36,16 +36,25 @@ END_PROVIDER
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) is fitted with a gaussian and a Slater
|
||||
!
|
||||
! \approx - 1/sqrt(pi) * exp(-alpha * x ) exp(-beta * x**2)
|
||||
!
|
||||
! where alpha = expo_j_xmu(1) and beta = expo_j_xmu(2)
|
||||
END_DOC
|
||||
expo_j_xmu(1) = 1.7477d0
|
||||
expo_j_xmu(2) = 0.668662d0
|
||||
|
||||
BEGIN_DOC
|
||||
! F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) is fitted with a gaussian and a Slater
|
||||
!
|
||||
! \approx - 1/sqrt(pi) * exp(-alpha * x ) exp(-beta * x**2)
|
||||
!
|
||||
! where alpha = expo_j_xmu(1) and beta = expo_j_xmu(2)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
!expo_j_xmu(1) = 1.7477d0
|
||||
!expo_j_xmu(2) = 0.668662d0
|
||||
|
||||
!expo_j_xmu(1) = 1.74766377595541d0
|
||||
!expo_j_xmu(2) = 0.668719925486403d0
|
||||
|
||||
expo_j_xmu(1) = 1.74770446934522d0
|
||||
expo_j_xmu(2) = 0.668659706559979d0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -9,19 +9,19 @@ program bi_ort_ints
|
||||
implicit none
|
||||
|
||||
my_grid_becke = .True.
|
||||
!my_n_pt_r_grid = 10
|
||||
!my_n_pt_a_grid = 14
|
||||
my_n_pt_r_grid = 30
|
||||
my_n_pt_a_grid = 50
|
||||
PROVIDE tc_grid1_a tc_grid1_r
|
||||
my_n_pt_r_grid = tc_grid1_r
|
||||
my_n_pt_a_grid = tc_grid1_a
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
|
||||
! call test_3e
|
||||
! call test_5idx
|
||||
! call test_5idx2
|
||||
!call test_4idx
|
||||
call test_4idx2()
|
||||
call test_5idx2
|
||||
call test_5idx
|
||||
call test_4idx()
|
||||
call test_4idx_n4()
|
||||
!call test_4idx2()
|
||||
!call test_5idx2
|
||||
!call test_5idx
|
||||
end
|
||||
|
||||
subroutine test_5idx2
|
||||
@ -211,13 +211,138 @@ end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_4idx_n4()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, l
|
||||
double precision :: accu, contrib, new, ref, thr
|
||||
|
||||
thr = 1d-10
|
||||
|
||||
PROVIDE three_e_4_idx_direct_bi_ort_old
|
||||
PROVIDE three_e_4_idx_direct_bi_ort_n4
|
||||
|
||||
accu = 0.d0
|
||||
do i = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
|
||||
new = three_e_4_idx_direct_bi_ort_n4 (l,k,j,i)
|
||||
ref = three_e_4_idx_direct_bi_ort_old(l,k,j,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. thr) then
|
||||
print*, ' problem in three_e_4_idx_direct_bi_ort_n4'
|
||||
print*, l, k, j, i
|
||||
print*, ref, new, contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*, ' accu on three_e_4_idx_direct_bi_ort_n4 = ', accu / dble(mo_num)**4
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE three_e_4_idx_exch13_bi_ort_old
|
||||
PROVIDE three_e_4_idx_exch13_bi_ort_n4
|
||||
|
||||
accu = 0.d0
|
||||
do i = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
|
||||
new = three_e_4_idx_exch13_bi_ort_n4 (l,k,j,i)
|
||||
ref = three_e_4_idx_exch13_bi_ort_old(l,k,j,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. thr) then
|
||||
print*, ' problem in three_e_4_idx_exch13_bi_ort_n4'
|
||||
print*, l, k, j, i
|
||||
print*, ref, new, contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*, ' accu on three_e_4_idx_exch13_bi_ort_n4 = ', accu / dble(mo_num)**4
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE three_e_4_idx_cycle_1_bi_ort_old
|
||||
PROVIDE three_e_4_idx_cycle_1_bi_ort_n4
|
||||
|
||||
accu = 0.d0
|
||||
do i = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
|
||||
new = three_e_4_idx_cycle_1_bi_ort_n4 (l,k,j,i)
|
||||
ref = three_e_4_idx_cycle_1_bi_ort_old(l,k,j,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. thr) then
|
||||
print*, ' problem in three_e_4_idx_cycle_1_bi_ort_n4'
|
||||
print*, l, k, j, i
|
||||
print*, ref, new, contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*, ' accu on three_e_4_idx_cycle_1_bi_ort_n4 = ', accu / dble(mo_num)**4
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE three_e_4_idx_exch23_bi_ort_old
|
||||
PROVIDE three_e_4_idx_exch23_bi_ort_n4
|
||||
|
||||
accu = 0.d0
|
||||
do i = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
|
||||
new = three_e_4_idx_exch23_bi_ort_n4 (l,k,j,i)
|
||||
ref = three_e_4_idx_exch23_bi_ort_old(l,k,j,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. thr) then
|
||||
print*, ' problem in three_e_4_idx_exch23_bi_ort_n4'
|
||||
print*, l, k, j, i
|
||||
print*, ref, new, contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*, ' accu on three_e_4_idx_exch23_bi_ort_n4 = ', accu / dble(mo_num)**4
|
||||
|
||||
! ---
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_4idx()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, l
|
||||
double precision :: accu, contrib, new, ref, thr
|
||||
|
||||
thr = 1d-5
|
||||
thr = 1d-10
|
||||
|
||||
PROVIDE three_e_4_idx_direct_bi_ort_old
|
||||
PROVIDE three_e_4_idx_direct_bi_ort
|
||||
@ -275,34 +400,6 @@ subroutine test_4idx()
|
||||
|
||||
! ---
|
||||
|
||||
! PROVIDE three_e_4_idx_exch12_bi_ort_old
|
||||
! PROVIDE three_e_4_idx_exch12_bi_ort
|
||||
!
|
||||
! accu = 0.d0
|
||||
! do i = 1, mo_num
|
||||
! do j = 1, mo_num
|
||||
! do k = 1, mo_num
|
||||
! do l = 1, mo_num
|
||||
!
|
||||
! new = three_e_4_idx_exch12_bi_ort (l,k,j,i)
|
||||
! ref = three_e_4_idx_exch12_bi_ort_old(l,k,j,i)
|
||||
! contrib = dabs(new - ref)
|
||||
! accu += contrib
|
||||
! if(contrib .gt. thr) then
|
||||
! print*, ' problem in three_e_4_idx_exch12_bi_ort'
|
||||
! print*, l, k, j, i
|
||||
! print*, ref, new, contrib
|
||||
! stop
|
||||
! endif
|
||||
!
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! print*, ' accu on three_e_4_idx_exch12_bi_ort = ', accu / dble(mo_num)**4
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE three_e_4_idx_cycle_1_bi_ort_old
|
||||
PROVIDE three_e_4_idx_cycle_1_bi_ort
|
||||
|
||||
@ -331,34 +428,6 @@ subroutine test_4idx()
|
||||
|
||||
! ---
|
||||
|
||||
! PROVIDE three_e_4_idx_cycle_2_bi_ort_old
|
||||
! PROVIDE three_e_4_idx_cycle_2_bi_ort
|
||||
!
|
||||
! accu = 0.d0
|
||||
! do i = 1, mo_num
|
||||
! do j = 1, mo_num
|
||||
! do k = 1, mo_num
|
||||
! do l = 1, mo_num
|
||||
!
|
||||
! new = three_e_4_idx_cycle_2_bi_ort (l,k,j,i)
|
||||
! ref = three_e_4_idx_cycle_2_bi_ort_old(l,k,j,i)
|
||||
! contrib = dabs(new - ref)
|
||||
! accu += contrib
|
||||
! if(contrib .gt. thr) then
|
||||
! print*, ' problem in three_e_4_idx_cycle_2_bi_ort'
|
||||
! print*, l, k, j, i
|
||||
! print*, ref, new, contrib
|
||||
! stop
|
||||
! endif
|
||||
!
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! print*, ' accu on three_e_4_idx_cycle_2_bi_ort = ', accu / dble(mo_num)**4
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE three_e_4_idx_exch23_bi_ort_old
|
||||
PROVIDE three_e_4_idx_exch23_bi_ort
|
||||
|
||||
|
@ -6,7 +6,7 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
|
||||
implicit none
|
||||
integer :: i, j
|
||||
|
||||
ao_one_e_integrals_tc_tot = ao_one_e_integrals
|
||||
ao_one_e_integrals_tc_tot = ao_one_e_integrals
|
||||
|
||||
!provide j1b_type
|
||||
|
||||
|
@ -140,8 +140,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
|
||||
enddo
|
||||
enddo
|
||||
|
||||
FREE int2_grad1_u12_ao
|
||||
|
||||
endif
|
||||
|
||||
call wall_time(wall1)
|
||||
@ -225,6 +223,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3,
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
|
||||
PROVIDE int2_grad1_u12_ao
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
|
@ -17,6 +17,8 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_direct_bi_ort, (mo_num, mo_num,
|
||||
integer :: i, j, m
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
PROVIDE mo_l_coef mo_r_coef
|
||||
|
||||
three_e_3_idx_direct_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_3_idx_direct_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
@ -125,6 +127,8 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_2_bi_ort, (mo_num, mo_num
|
||||
integer :: i, j, m
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
PROVIDE mo_l_coef mo_r_coef
|
||||
|
||||
three_e_3_idx_cycle_2_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_3_idx_cycle_2_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
@ -179,6 +183,8 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch23_bi_ort, (mo_num, mo_num,
|
||||
integer :: i, j, m
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
PROVIDE mo_l_coef mo_r_coef
|
||||
|
||||
three_e_3_idx_exch23_bi_ort = 0.d0
|
||||
print*,'Providing the three_e_3_idx_exch23_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
@ -233,6 +239,8 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch13_bi_ort, (mo_num, mo_num,
|
||||
integer :: i,j,m
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
PROVIDE mo_l_coef mo_r_coef
|
||||
|
||||
three_e_3_idx_exch13_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_3_idx_exch13_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
@ -287,6 +295,8 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort, (mo_num, mo_num,
|
||||
integer :: i, j, m
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
PROVIDE mo_l_coef mo_r_coef
|
||||
|
||||
three_e_3_idx_exch12_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_3_idx_exch12_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
|
@ -3,9 +3,8 @@
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
|
||||
!&BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
|
||||
!&BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
@ -13,28 +12,25 @@
|
||||
!
|
||||
! three_e_4_idx_direct_bi_ort (m,j,k,i) = < m j k | -L | m j i > ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
! three_e_4_idx_exch13_bi_ort (m,j,k,i) = < m j k | -L | i j m > ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
! three_e_4_idx_exch12_bi_ort (m,j,k,i) = < m j k | -L | m i j > ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
! = three_e_4_idx_exch13_bi_ort (j,m,k,i)
|
||||
! three_e_4_idx_exch23_bi_ort (m,j,k,i) = < m j k | -L | j m i > ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = < m j k | -L | j i m > ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
! three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = < m j k | -L | i m j > ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
! = three_e_4_idx_cycle_1_bi_ort(j,m,k,i)
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_4_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
! three_e_4_idx_direct_bi_ort (m,j,k,i) : Lk Ri Imm Ijj + Lj Rj Imm Iki + Lm Rm Ijj Iki
|
||||
! three_e_4_idx_exch13_bi_ort (m,j,k,i) : Lk Rm Imi Ijj + Lj Rj Imi Ikm + Lm Ri Ijj Ikm
|
||||
! three_e_4_idx_exch23_bi_ort (m,j,k,i) : Lk Ri Imj Ijm + Lj Rm Imj Iki + Lm Rj Ijm Iki
|
||||
! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) : Lk Rm Imj Iji + Lj Ri Imj Ikm + Lm Rj Iji Ikm
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: ipoint, i, j, k, l, m
|
||||
integer :: ipoint, i, j, k, m, n
|
||||
double precision :: wall1, wall0
|
||||
double precision, allocatable :: tmp1(:,:,:,:), tmp2(:,:,:,:), tmp3(:,:,:,:)
|
||||
double precision, allocatable :: tmp_4d(:,:,:,:)
|
||||
double precision, allocatable :: tmp4(:,:,:)
|
||||
double precision, allocatable :: tmp5(:,:)
|
||||
double precision, allocatable :: tmp_3d(:,:,:)
|
||||
double precision :: tmp_loc_1, tmp_loc_2
|
||||
double precision, allocatable :: tmp1(:,:,:), tmp2(:,:,:)
|
||||
double precision, allocatable :: tmp_2d(:,:)
|
||||
double precision, allocatable :: tmp_aux_1(:,:,:), tmp_aux_2(:,:)
|
||||
|
||||
print *, ' Providing the three_e_4_idx_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
@ -42,324 +38,181 @@
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
|
||||
allocate(tmp_4d(mo_num,mo_num,mo_num,mo_num))
|
||||
|
||||
allocate(tmp1(n_points_final_grid,3,mo_num,mo_num))
|
||||
allocate(tmp2(n_points_final_grid,3,mo_num,mo_num))
|
||||
allocate(tmp3(n_points_final_grid,3,mo_num,mo_num))
|
||||
! to reduce the number of operations
|
||||
allocate(tmp_aux_1(n_points_final_grid,4,mo_num))
|
||||
allocate(tmp_aux_2(n_points_final_grid,mo_num))
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, l, ipoint) &
|
||||
!$OMP PRIVATE (n, ipoint) &
|
||||
!$OMP SHARED (mo_num, n_points_final_grid, &
|
||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||
!$OMP tmp1, tmp2, tmp3)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do ipoint = 1, n_points_final_grid
|
||||
|
||||
tmp1(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,l,l) * mos_l_in_r_array_transp(ipoint,i) * final_weight_at_r_vector(ipoint)
|
||||
tmp1(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,l,l) * mos_l_in_r_array_transp(ipoint,i) * final_weight_at_r_vector(ipoint)
|
||||
tmp1(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,l,l) * mos_l_in_r_array_transp(ipoint,i) * final_weight_at_r_vector(ipoint)
|
||||
|
||||
tmp2(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,l,l) * mos_r_in_r_array_transp(ipoint,i)
|
||||
tmp2(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,l,l) * mos_r_in_r_array_transp(ipoint,i)
|
||||
tmp2(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,l,l) * mos_r_in_r_array_transp(ipoint,i)
|
||||
|
||||
tmp3(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,l,i) * mos_r_in_r_array_transp(ipoint,l)
|
||||
tmp3(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,l,i) * mos_r_in_r_array_transp(ipoint,l)
|
||||
tmp3(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,l,i) * mos_r_in_r_array_transp(ipoint,l)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
|
||||
, tmp1, 3*n_points_final_grid, tmp2, 3*n_points_final_grid &
|
||||
, 0.d0, tmp_4d, mo_num*mo_num)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
three_e_4_idx_direct_bi_ort(m,j,k,i) = -tmp_4d(m,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
|
||||
, tmp3, 3*n_points_final_grid, tmp1, 3*n_points_final_grid &
|
||||
, 0.d0, tmp_4d, mo_num*mo_num)
|
||||
|
||||
! deallocate(tmp1)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
three_e_4_idx_exch13_bi_ort(m,j,k,i) = -tmp_4d(m,i,j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
|
||||
! allocate(tmp1(n_points_final_grid, 2, mo_num, mo_num))
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, l, ipoint) &
|
||||
!$OMP SHARED (mo_num, n_points_final_grid, &
|
||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||
!$OMP tmp1)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do ipoint = 1, n_points_final_grid
|
||||
tmp1(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,i,l) * mos_l_in_r_array_transp(ipoint,l) * final_weight_at_r_vector(ipoint)
|
||||
tmp1(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,i,l) * mos_l_in_r_array_transp(ipoint,l) * final_weight_at_r_vector(ipoint)
|
||||
tmp1(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,i,l) * mos_l_in_r_array_transp(ipoint,l) * final_weight_at_r_vector(ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
|
||||
, tmp1, 3*n_points_final_grid, tmp2, 3*n_points_final_grid &
|
||||
, 0.d0, tmp_4d, mo_num*mo_num)
|
||||
|
||||
deallocate(tmp2)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
three_e_4_idx_exch13_bi_ort(m,j,k,i) = three_e_4_idx_exch13_bi_ort(m,j,k,i) - tmp_4d(m,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||