LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
Code Releases Activity

new jast added in QP

Browse Source
This commit is contained in:
AbdAmmar 2022-10-21 23:27:51 +02:00
parent 813f2c3601
commit fccf2fe438
33 changed files with 3076 additions and 1834 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
src/ao_many_one_e_ints/ao_erf_gauss.irp.f
View File

@ -298,10 +298,16 @@ subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_cen
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 :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef, coefi
double precision, external :: NAI_pol_mult_erf_with1s
ASSERT(beta .ge. 0.d0)
if(beta .lt. 1d-10) then
call NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
return
endif
ints = 0.d0
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) then
return
@ -316,26 +322,27 @@ subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_cen
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alphai = ao_expo_ordered_transp(i,i_ao)
alphai = ao_expo_ordered_transp (i,i_ao)
coefi = ao_coef_normalized_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 = power_Ai
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)
alphaj = ao_expo_ordered_transp (j,j_ao)
coef = coefi * ao_coef_normalized_ordered_transp(j,j_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 )
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 )
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

26
src/ao_many_one_e_ints/ao_gaus_gauss.irp.f
View File

@ -116,7 +116,7 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
double precision, intent(in) :: D_center(3), delta
integer :: power_A(3), power_B(3), l, k
double precision :: A_center(3), B_center(3), alpha, beta, coef, analytical_j
double precision :: A_center(3), B_center(3), alpha, beta, coef, coef1, analytical_j
double precision, external :: overlap_gauss_r12
@ -133,10 +133,12 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
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)
alpha = ao_expo_ordered_transp (l,i)
coef1 = ao_coef_normalized_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)
coef = coef1 * ao_coef_normalized_ordered_transp(k,j)
if(dabs(coef) .lt. 1d-12) cycle
@ -153,7 +155,9 @@ 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
@ -161,7 +165,7 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center,
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 :: A1_center(3), A2_center(3), alpha1, alpha2, coef1, coef12, analytical_j
double precision :: G_center(3), gama, fact_g, gama_inv
double precision, external :: overlap_gauss_r12, overlap_gauss_r12_ao
@ -188,8 +192,8 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center,
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
if(fact_g .gt. 80d0) return
fact_g = dexp(-fact_g)
! ---
@ -200,11 +204,13 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center,
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)
alpha1 = ao_expo_ordered_transp (l,i)
coef1 = fact_g * ao_coef_normalized_ordered_transp(l,i)
!if(dabs(coef1) .lt. 1d-12) cycle
do k = 1, ao_prim_num(j)
alpha2 = ao_expo_ordered_transp (k,j)
coef12 = coef1 * 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)

204
src/ao_many_one_e_ints/grad2_jmu_modif.irp.f
View File

@ -13,8 +13,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
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 :: tmp
double precision :: wall0, wall1
double precision, allocatable :: tmp(:,:,:)
double precision, external :: overlap_gauss_r12_ao_with1s
@ -31,19 +31,17 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2)
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
tmp = 0.d0
!$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)
r(3) = final_grid_points(3,ipoint)
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)
tmp = 0.d0
do i_1s = 1, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
@ -58,29 +56,19 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
coef_fit = coef_gauss_1_erf_x_2(i_fit)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp(j,i,ipoint) += -0.25d0 * coef * coef_fit * int_fit
tmp += -0.25d0 * coef * coef_fit * int_fit
enddo
enddo
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = tmp
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
int2_grad1u2_grad2u2_j1b2(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 i = 2, ao_num
do j = 1, i-1
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
enddo
@ -105,9 +93,8 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
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 :: coef, beta, B_center(3), tmp
double precision :: wall0, wall1
double precision, allocatable :: tmp(:,:,:)
double precision, external :: overlap_gauss_r12_ao_with1s
@ -124,19 +111,17 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u2_j1b2)
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
tmp = 0.d0
!$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)
r(3) = final_grid_points(3,ipoint)
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)
tmp = 0.d0
do i_1s = 1, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
@ -151,29 +136,19 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
coef_fit = coef_gauss_j_mu_x_2(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
tmp += coef * coef_fit * int_fit
enddo
enddo
int2_u2_j1b2(j,i,ipoint) = tmp
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
int2_u2_j1b2(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 i = 2, ao_num
do j = 1, i-1
int2_u2_j1b2(j,i,ipoint) = int2_u2_j1b2(i,j,ipoint)
enddo
@ -187,7 +162,7 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
@ -196,39 +171,40 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_po
END_DOC
implicit none
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), int_fit(3), expo_fit, coef_fit
double precision :: coef, beta, B_center(3)
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coeff_1s
double precision :: wall0, wall1
double precision, allocatable :: tmp(:,:,:,:)
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), int_fit(3), expo_fit, coef_fit
double precision :: coef, beta, B_center(3), dist
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp
double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
int2_u_grad1u_x_j1b = 0.d0
int2_u_grad1u_x_j1b2 = 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, alpha_1s, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coeff_1s) &
!$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, &
!$OMP tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP final_grid_points, n_max_fit_slat, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b)
allocate( tmp(3,ao_num,ao_num,n_points_final_grid) )
tmp = 0.d0
!$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b2)
!$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)
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)
tmp_x = 0.d0
tmp_y = 0.d0
tmp_z = 0.d0
do i_1s = 1, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
@ -236,6 +212,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_po
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3))
do i_fit = 1, n_max_fit_slat
@ -244,56 +223,45 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b, (3, ao_num, ao_num, n_po
alpha_1s = beta + expo_fit
alpha_1s_inv = 1.d0 / alpha_1s
centr_1s(1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
centr_1s(2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
expo_coef_1s = -beta * expo_fit * alpha_1s_inv &
* ( (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3)) )
if(expo_coef_1s .gt. 80.d0) cycle
coeff_1s = dexp(-expo_coef_1s)
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
!if(expo_coef_1s .gt. 80.d0) cycle
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
!if(dabs(coef_tmp) .lt. 1d-10) cycle
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
tmp(1,j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit(1)
tmp(2,j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit(2)
tmp(3,j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit(3)
tmp_x += coef_tmp * int_fit(1)
tmp_y += coef_tmp * int_fit(2)
tmp_z += coef_tmp * int_fit(3)
enddo
enddo
int2_u_grad1u_x_j1b2(1,j,i,ipoint) = tmp_x
int2_u_grad1u_x_j1b2(2,j,i,ipoint) = tmp_y
int2_u_grad1u_x_j1b2(3,j,i,ipoint) = tmp_z
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
int2_u_grad1u_x_j1b(1,j,i,ipoint) += tmp(1,j,i,ipoint)
int2_u_grad1u_x_j1b(2,j,i,ipoint) += tmp(2,j,i,ipoint)
int2_u_grad1u_x_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 i = 2, ao_num
do j = 1, i-1
int2_u_grad1u_x_j1b(1,j,i,ipoint) = int2_u_grad1u_x_j1b(1,i,j,ipoint)
int2_u_grad1u_x_j1b(2,j,i,ipoint) = int2_u_grad1u_x_j1b(2,i,j,ipoint)
int2_u_grad1u_x_j1b(3,j,i,ipoint) = int2_u_grad1u_x_j1b(3,i,j,ipoint)
int2_u_grad1u_x_j1b2(1,j,i,ipoint) = int2_u_grad1u_x_j1b2(1,i,j,ipoint)
int2_u_grad1u_x_j1b2(2,j,i,ipoint) = int2_u_grad1u_x_j1b2(2,i,j,ipoint)
int2_u_grad1u_x_j1b2(3,j,i,ipoint) = int2_u_grad1u_x_j1b2(3,i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_x_j1b', wall1 - wall0
print*, ' wall time for int2_u_grad1u_x_j1b2', wall1 - wall0
END_PROVIDER
@ -309,11 +277,10 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
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 :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coeff_1s
double precision :: r(3), int_fit, expo_fit, coef_fit, coef_tmp
double precision :: coef, beta, B_center(3), dist
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp
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
@ -323,17 +290,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
!$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, alpha_1s, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coeff_1s) &
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP final_grid_points, n_max_fit_slat, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2)
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
@ -342,6 +305,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
tmp = 0.d0
do i_1s = 1, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
@ -349,6 +313,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3))
do i_fit = 1, n_max_fit_slat
@ -360,39 +327,27 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
centr_1s(1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
centr_1s(2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
expo_coef_1s = -beta * expo_fit * alpha_1s_inv &
* ( (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3)) )
if(expo_coef_1s .gt. 80.d0) cycle
coeff_1s = dexp(-expo_coef_1s)
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
!if(expo_coef_1s .gt. 80.d0) cycle
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
!if(dabs(coef_tmp) .lt. 1d-10) cycle
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r)
tmp(j,i,ipoint) += coef * coef_fit * coeff_1s * int_fit
tmp += coef_tmp * int_fit
enddo
enddo
int2_u_grad1u_j1b2(j,i,ipoint) = tmp
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
int2_u_grad1u_j1b2(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 i = 2, ao_num
do j = 1, i-1
int2_u_grad1u_j1b2(j,i,ipoint) = int2_u_grad1u_j1b2(i,j,ipoint)
enddo
@ -405,3 +360,4 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
END_PROVIDER
! ---

138
src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f
View File

@ -10,13 +10,12 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
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
integer :: i, j, ipoint, i_1s
double precision :: r(3), int_mu, int_coulomb
double precision :: coef, beta, B_center(3)
double precision :: tmp
double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
@ -28,19 +27,17 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_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, 10
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 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)
tmp = 0.d0
do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
@ -52,28 +49,18 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
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)
tmp += coef * (int_mu - int_coulomb)
enddo
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = tmp
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 i = 2, 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
@ -123,33 +110,34 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_
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(:,:,:,:)
integer :: i, j, ipoint, i_1s
double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1
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 PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, &
!$OMP tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,&
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_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, 10
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 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)
tmp_x = 0.d0
tmp_y = 0.d0
tmp_z = 0.d0
do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
@ -161,32 +149,22 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_
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))
tmp_x += coef * (ints(1) - ints_coulomb(1))
tmp_y += coef * (ints(2) - ints_coulomb(2))
tmp_z += coef * (ints(3) - ints_coulomb(3))
enddo
x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) = tmp_x
x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) = tmp_y
x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint) = tmp_z
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 i = 2, ao_num
do j = 1, i-1
x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b(1,i,j,ipoint)
x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b(2,i,j,ipoint)
@ -202,6 +180,7 @@ END_PROVIDER
! ---
! TODO analytically
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
@ -211,13 +190,13 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
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(:,:,:)
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 :: tmp
double precision :: wall0, wall1
double precision, external :: overlap_gauss_r12_ao_with1s
double precision, external :: overlap_gauss_r12_ao_with1s
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
@ -232,19 +211,17 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
!$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)
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
tmp = 0.d0
!$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)
r(3) = final_grid_points(3,ipoint)
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)
tmp = 0.d0
do i_1s = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
@ -259,29 +236,19 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
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
tmp += coef * coef_fit * int_fit
enddo
enddo
v_ij_u_cst_mu_j1b(j,i,ipoint) = tmp
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 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)
enddo
@ -294,3 +261,4 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
END_PROVIDER
! ---

450
src/ao_many_one_e_ints/grad_related_ints.irp.f
View File

@ -28,11 +28,12 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points
!$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
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
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)
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)
@ -45,7 +46,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do i = 2, ao_num
do j = 1, i-1
v_ij_erf_rk_cst_mu(j,i,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint)
enddo
@ -53,54 +54,61 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu, (ao_num, ao_num, n_points
enddo
call wall_time(wall1)
print*, 'wall time for v_ij_erf_rk_cst_mu ', wall1 - wall0
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
! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R| - 1)/|r - R|
END_DOC
integer :: i,j,ipoint
double precision :: mu,r(3),NAI_pol_mult_erf_ao
double precision :: int_mu, int_coulomb
provide mu_erf final_grid_points
double precision :: wall0, wall1
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,ipoint,mu,r,int_mu,int_coulomb) &
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 :: r(3)
double precision :: int_mu, int_coulomb
double precision :: wall0, wall1
double precision :: NAI_pol_mult_erf_ao
provide mu_erf final_grid_points
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,ipoint,r,int_mu,int_coulomb) &
!$OMP SHARED (ao_num,n_points_final_grid,v_ij_erf_rk_cst_mu_transp,final_grid_points,mu_erf)
!$OMP DO SCHEDULE (dynamic)
do i = 1, ao_num
do j = i, ao_num
do ipoint = 1, n_points_final_grid
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_coulomb = NAI_pol_mult_erf_ao(i,j,1.d+9,r)
v_ij_erf_rk_cst_mu_transp(ipoint,j,i)= (int_mu - int_coulomb )
enddo
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 i = 1, ao_num
do j = i, ao_num
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_transp(ipoint,j,i) = int_mu - int_coulomb
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, ao_num
do j = 1, i-1
do ipoint = 1, n_points_final_grid
v_ij_erf_rk_cst_mu_transp(ipoint,j,i)= v_ij_erf_rk_cst_mu_transp(ipoint,i,j)
do i = 2, ao_num
do j = 1, i-1
do ipoint = 1, n_points_final_grid
v_ij_erf_rk_cst_mu_transp(ipoint,j,i) = v_ij_erf_rk_cst_mu_transp(ipoint,i,j)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for v_ij_erf_rk_cst_mu_transp ',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for v_ij_erf_rk_cst_mu_transp ', wall1 - wall0
END_PROVIDER
! ---
@ -112,30 +120,31 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3, ao_num, ao_num,
END_DOC
implicit none
integer :: i, j, ipoint, m
integer :: i, j, ipoint
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,r,ints,m,ints_coulomb) &
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,ipoint,r,ints,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
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
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)
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
x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint) = ints(1) - ints_coulomb(1)
x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint) = ints(2) - ints_coulomb(2)
x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint) = ints(3) - ints_coulomb(3)
enddo
enddo
enddo
@ -143,11 +152,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp, (3, ao_num, ao_num,
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do i = 2, ao_num
do j = 1, i-1
do m = 1, 3
x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(m,i,j,ipoint)
enddo
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
@ -160,208 +169,249 @@ 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
! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R|
END_DOC
integer :: i,j,ipoint,m
double precision :: mu,r(3),ints,ints_coulomb
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
do m = 1, 3
x_v_ij_erf_rk_cst_mu(j,i,ipoint,m)= x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for x_v_ij_erf_rk_cst_mu',wall1 - wall0
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
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(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint)
x_v_ij_erf_rk_cst_mu(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint)
x_v_ij_erf_rk_cst_mu(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for x_v_ij_erf_rk_cst_mu', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_transp, (ao_num, ao_num,3,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
integer :: i,j,ipoint,m
double precision :: mu,r(3),ints,ints_coulomb
double precision :: wall0, wall1
call wall_time(wall0)
do ipoint = 1, n_points_final_grid
do m = 1, 3
do i = 1, ao_num
do j = 1, ao_num
x_v_ij_erf_rk_cst_mu_transp(j,i,m,ipoint)= x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for x_v_ij_erf_rk_cst_mu_transp',wall1 - wall0
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
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_transp(j,i,1,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_transp(j,i,2,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_transp(j,i,3,ipoint) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for x_v_ij_erf_rk_cst_mu_transp', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_transp_bis, (n_points_final_grid,ao_num, ao_num,3)]
implicit none
BEGIN_DOC
! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/|r - R|
END_DOC
integer :: i,j,ipoint,m
double precision :: mu,r(3),ints,ints_coulomb
double precision :: wall0, wall1
call wall_time(wall0)
do m = 1, 3
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,m)= x_v_ij_erf_rk_cst_mu_tmp(m,j,i,ipoint)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for x_v_ij_erf_rk_cst_mu_transp',wall1 - wall0
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
double precision :: wall0, wall1
call wall_time(wall0)
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,1) = x_v_ij_erf_rk_cst_mu_tmp(1,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,2) = x_v_ij_erf_rk_cst_mu_tmp(2,j,i,ipoint)
x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,j,i,3) = x_v_ij_erf_rk_cst_mu_tmp(3,j,i,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for x_v_ij_erf_rk_cst_mu_transp_bis', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu_tmp, (3, n_points_final_grid, ao_num, ao_num)]
BEGIN_DOC
! d_dx_v_ij_erf_rk_cst_mu_tmp(m,R,j,i) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
END_DOC
BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu_tmp, (3,n_points_final_grid,ao_num, ao_num)]
implicit none
BEGIN_DOC
! d_dx_v_ij_erf_rk_cst_mu_tmp(m,R,j,i) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
END_DOC
integer :: i,j,ipoint,m
double precision :: mu,r(3),ints(3),ints_coulomb(3)
integer :: i, j, ipoint
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 PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,ipoint,r,ints,ints_coulomb) &
!$OMP SHARED (ao_num,n_points_final_grid,d_dx_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf)
!$OMP DO SCHEDULE (dynamic)
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
mu = mu_erf
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)
call phi_j_erf_mu_r_dxyz_phi(j,i,mu, r, ints)
call phi_j_erf_mu_r_dxyz_phi(j,i,1.d+9, r, ints_coulomb)
do m = 1, 3
d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i) = ( ints(m) - ints_coulomb(m))
do i = 1, ao_num
do j = 1, ao_num
call phi_j_erf_mu_r_dxyz_phi(j, i, mu_erf, r, ints)
call phi_j_erf_mu_r_dxyz_phi(j, i, 1.d+9, r, ints_coulomb)
d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) = ints(1) - ints_coulomb(1)
d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) = ints(2) - ints_coulomb(2)
d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) = ints(3) - ints_coulomb(3)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for d_dx_v_ij_erf_rk_cst_mu_tmp',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for d_dx_v_ij_erf_rk_cst_mu_tmp', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid,ao_num, ao_num,3)]
implicit none
BEGIN_DOC
! d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
END_DOC
integer :: i,j,ipoint,m
double precision :: mu,r(3),ints,ints_coulomb
double precision :: wall0, wall1
call wall_time(wall0)
do i = 1, ao_num
do j = 1, ao_num
do m = 1, 3
do ipoint = 1, n_points_final_grid
d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,m)= d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i)
! ---
BEGIN_PROVIDER [ double precision, d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid, ao_num, ao_num, 3)]
BEGIN_DOC
!
! d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
!
END_DOC
implicit none
integer :: i, j, ipoint
double precision :: wall0, wall1
call wall_time(wall0)
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,1) = d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i)
d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,2) = d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i)
d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,3) = d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i)
enddo
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for d_dx_v_ij_erf_rk_cst_mu',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for d_dx_v_ij_erf_rk_cst_mu', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu_tmp, (3,n_points_final_grid,ao_num, ao_num)]
implicit none
BEGIN_DOC
! x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,j,i,R) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
END_DOC
integer :: i,j,ipoint,m
double precision :: mu,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) &
! ---
BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu_tmp, (3, n_points_final_grid, ao_num, ao_num)]
BEGIN_DOC
!
! x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,j,i,R) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
!
END_DOC
implicit none
integer :: i, j, ipoint
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,r,ints,ints_coulomb) &
!$OMP SHARED (ao_num,n_points_final_grid,x_d_dx_v_ij_erf_rk_cst_mu_tmp,final_grid_points,mu_erf)
!$OMP DO SCHEDULE (dynamic)
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
mu = mu_erf
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)
call phi_j_erf_mu_r_xyz_dxyz_phi(j,i,mu, r, ints)
call phi_j_erf_mu_r_xyz_dxyz_phi(j,i,1.d+9, r, ints_coulomb)
do m = 1, 3
x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i) = ( ints(m) - ints_coulomb(m))
do i = 1, ao_num
do j = 1, ao_num
call phi_j_erf_mu_r_xyz_dxyz_phi(j, i, mu_erf, r, ints)
call phi_j_erf_mu_r_xyz_dxyz_phi(j, i, 1.d+9, r, ints_coulomb)
x_d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i) = ints(1) - ints_coulomb(1)
x_d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i) = ints(2) - ints_coulomb(2)
x_d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i) = ints(3) - ints_coulomb(3)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for x_d_dx_v_ij_erf_rk_cst_mu_tmp',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for x_d_dx_v_ij_erf_rk_cst_mu_tmp', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, x_d_dx_v_ij_erf_rk_cst_mu, (n_points_final_grid,ao_num, ao_num,3)]
implicit none
BEGIN_DOC
! x_d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
END_DOC
integer :: i,j,ipoint,m
double precision :: mu,r(3),ints,ints_coulomb
double precision :: wall0, wall1
call wall_time(wall0)
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
do m = 1, 3
x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,m)= x_d_dx_v_ij_erf_rk_cst_mu_tmp(m,ipoint,j,i)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for x_d_dx_v_ij_erf_rk_cst_mu',wall1 - wall0
BEGIN_DOC
!
! x_d_dx_v_ij_erf_rk_cst_mu_tmp(j,i,R,m) = int dr x phi_j(r)) (erf(mu(R) |r - R|) - 1)/|r - R| d/dx (phi_i(r)
!
! with m == 1 -> d/dx , m == 2 -> d/dy , m == 3 -> d/dz
!
END_DOC
implicit none
integer :: i, j, ipoint
double precision :: wall0, wall1
call wall_time(wall0)
do i = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,1) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(1,ipoint,j,i)
x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,2) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(2,ipoint,j,i)
x_d_dx_v_ij_erf_rk_cst_mu(ipoint,j,i,3) = x_d_dx_v_ij_erf_rk_cst_mu_tmp(3,ipoint,j,i)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for x_d_dx_v_ij_erf_rk_cst_mu', wall1 - wall0
END_PROVIDER
! ---

29
src/ao_one_e_ints/pot_ao_erf_ints.irp.f
View File

@ -78,7 +78,7 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center,
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 :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao
@ -98,11 +98,12 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center,
NAI_pol_mult_erf_ao_with1s = 0.d0
do i = 1, ao_prim_num(i_ao)
alpha1 = ao_expo_ordered_transp(i,i_ao)
alpha1 = ao_expo_ordered_transp (i,i_ao)
coef1 = ao_coef_normalized_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)
coef12 = coef1 * ao_coef_normalized_ordered_transp(j,j_ao)
if(dabs(coef12) .lt. 1d-14) cycle
integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
@ -242,9 +243,9 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A
A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv
A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv
A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv
dist12 = ( (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) &
+ (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) &
+ (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) )
dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) &
+ (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) &
+ (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3))
const_factor12 = dist12 * rho12
if(const_factor12 > 80.d0) then
@ -262,9 +263,9 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv
dist = ( (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) &
+ (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) &
+ (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3)) )
dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) &
+ (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) &
+ (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3))
const_factor = const_factor12 + dist * rho
if(const_factor > 80.d0) then
@ -272,11 +273,9 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A
return
endif
dist_integral = 0.d0
do i = 1, 3
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
enddo
dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) &
+ (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) &
+ (P_center(3) - C_center(3)) * (P_center(3) - C_center(3))
! ---

1
src/ao_tc_eff_map/compute_ints_eff_pot.irp.f
View File

@ -61,7 +61,6 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va
integral = integral + j1b_gauss_2e_j2(i, k, j, l)
endif
if(abs(integral) < thr) then
cycle
endif

0
src/non_h_ints_mu/fit_j.irp.f → src/ao_tc_eff_map/fit_j.irp.f
View File

1
src/ao_tc_eff_map/two_e_ints_gauss.irp.f
View File

@ -254,6 +254,7 @@ double precision function general_primitive_integral_gauss(dim, &
rho_old = (p*q)/(p+q)
prefactor = pi_3 * inv_pq_3_2 * fact_p * fact_q
do i = 1, n_gauss_eff_pot ! browse the gaussians with different expo/coef
!do i = 1, n_gauss_eff_pot-1
aa = expo_gauss_eff_pot(i)
c_a = coef_gauss_eff_pot(i)
t_a = dsqrt( aa /(rho_old + aa) )

3
src/ao_two_e_ints/map_integrals.irp.f
View File

@ -321,8 +321,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
!$OMP END PARALLEL DO
END_PROVIDER
! ---
double precision function get_ao_two_e_integral(i,j,k,l,map) result(result)
double precision function get_ao_two_e_integral(i, j, k, l, map) result(result)
use map_module
implicit none
BEGIN_DOC

40
src/bi_ort_ints/biorthog_mo_for_h.irp.f
View File

@ -1,37 +1,6 @@
! ---
BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num) ]
BEGIN_DOC
!
! ao_two_e_coul(k,i,l,j) = ( k i | 1/r12 | l j ) = < l k | 1/r12 | j i >
!
END_DOC
integer :: i, j, k, l
double precision :: integral
double precision, external :: get_ao_two_e_integral
PROVIDE ao_integrals_map
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
integral = get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
ao_two_e_coul(k,i,l,j) = integral
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
double precision function bi_ortho_mo_coul_ints(l, k, j, i)
BEGIN_DOC
@ -155,7 +124,7 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_coul_e, (mo_num, mo_num, mo_num, m
do i = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! < k l | V12 | i j > (k i|l j)
! < k l | V12 | i j > (k i|l j)
mo_bi_ortho_coul_e(k,l,i,j) = mo_bi_ortho_coul_e_chemist(k,i,l,j)
enddo
enddo
@ -169,13 +138,14 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bi_ortho_one_e, (mo_num, mo_num)]
BEGIN_DOC
! mo_bi_ortho_one_e(k,i) = <MO^L_k | h_c | MO^R_i>
!
! mo_bi_ortho_one_e(k,i) = < MO^L_k | h_c | MO^R_i >
!
END_DOC
implicit none
call ao_to_mo_bi_ortho( ao_one_e_integrals, ao_num &
, mo_bi_ortho_one_e , mo_num )
call ao_to_mo_bi_ortho(ao_one_e_integrals, ao_num, mo_bi_ortho_one_e , mo_num)
END_PROVIDER

2
src/bi_ort_ints/one_e_bi_ort.irp.f
View File

@ -10,7 +10,7 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
provide j1b_type
if(j1b_type .ne. 0) then
if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
do i = 1, ao_num
do j = 1, ao_num

66
src/bi_ort_ints/semi_num_ints_mo.irp.f
View File

@ -86,10 +86,8 @@ BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo
call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp (1,1,1,ipoint), size(x_v_ij_erf_rk_cst_mu_transp, 1) &
, mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,1,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp (1,1,2,ipoint), size(x_v_ij_erf_rk_cst_mu_transp, 1) &
, mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,2,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp (1,1,3,ipoint), size(x_v_ij_erf_rk_cst_mu_transp, 1) &
, mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,3,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
@ -103,7 +101,55 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, n_points_final_grid)]
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
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(1,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(2,j,i,ipoint)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(3,j,i,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)]
implicit none
integer :: ipoint
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint) &
!$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp)
!$OMP DO SCHEDULE (dynamic)
do ipoint = 1, n_points_final_grid
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, mo_num, mo_num, n_points_final_grid)]
BEGIN_DOC
!
@ -121,14 +167,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, ao_num, ao_num, n_po
!$OMP DO SCHEDULE (dynamic)
do ipoint = 1, n_points_final_grid
call ao_to_mo_bi_ortho( int2_grad1_u12_ao (1,1,1,ipoint), size(int2_grad1_u12_ao , 1) &
, int2_grad1_u12_bimo(1,1,1,ipoint), size(int2_grad1_u12_bimo, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao (2,1,1,ipoint), size(int2_grad1_u12_ao , 1) &
, int2_grad1_u12_bimo(2,1,1,ipoint), size(int2_grad1_u12_bimo, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao (3,1,1,ipoint), size(int2_grad1_u12_ao , 1) &
, int2_grad1_u12_bimo(3,1,1,ipoint), size(int2_grad1_u12_bimo, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao (1,1,1,ipoint), size(int2_grad1_u12_ao , 2) &
, int2_grad1_u12_bimo(1,1,1,ipoint), size(int2_grad1_u12_bimo, 2) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao (2,1,1,ipoint), size(int2_grad1_u12_ao , 2) &
, int2_grad1_u12_bimo(2,1,1,ipoint), size(int2_grad1_u12_bimo, 2) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao (3,1,1,ipoint), size(int2_grad1_u12_ao , 2) &
, int2_grad1_u12_bimo(3,1,1,ipoint), size(int2_grad1_u12_bimo, 2) )
enddo
!$OMP END DO

460
src/bi_ort_ints/three_body_ijm.irp.f
View File

@ -1,304 +1,366 @@
! ---
BEGIN_PROVIDER [ double precision, three_e_3_idx_direct_bi_ort, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the direct terms
!
! three_e_3_idx_direct_bi_ort(m,j,i) = <mji|-L|mji>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_3_idx_direct_bi_ort = 0.d0
print*,'Providing the three_e_3_idx_direct_bi_ort ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the direct terms
!
! three_e_3_idx_direct_bi_ort(m,j,i) = <mji|-L|mji>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0
three_e_3_idx_direct_bi_ort = 0.d0
print *, ' Providing the three_e_3_idx_direct_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_3_idx_direct_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m,j,i,m,j,i,integral)
three_e_3_idx_direct_bi_ort(m,j,i) = -1.d0 * integral
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m, j, i, m, j, i, integral)
three_e_3_idx_direct_bi_ort(m,j,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_3_idx_direct_bi_ort',wall1 - wall0
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_3_idx_direct_bi_ort(m,j,i) = three_e_3_idx_direct_bi_ort(j,m,i)
do j = 1, mo_num
do m = 1, j
three_e_3_idx_direct_bi_ort(m,j,i) = three_e_3_idx_direct_bi_ort(j,m,i)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_direct_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the first cyclic permutation
!
! three_e_3_idx_direct_bi_ort(m,j,i) = <mji|-L|jim>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_3_idx_cycle_1_bi_ort = 0.d0
print*,'Providing the three_e_3_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the first cyclic permutation
!
! three_e_3_idx_direct_bi_ort(m,j,i) = <mji|-L|jim>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0
three_e_3_idx_cycle_1_bi_ort = 0.d0
print *, ' Providing the three_e_3_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_3_idx_cycle_1_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m,j,i,j,i,m,integral)
three_e_3_idx_cycle_1_bi_ort(m,j,i) = -1.d0 * integral
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m, j, i, j, i, m, integral)
three_e_3_idx_cycle_1_bi_ort(m,j,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_3_idx_cycle_1_bi_ort(m,j,i) = three_e_3_idx_cycle_1_bi_ort(j,m,i)
do j = 1, mo_num
do m = 1, j
three_e_3_idx_cycle_1_bi_ort(m,j,i) = three_e_3_idx_cycle_1_bi_ort(j,m,i)
enddo
enddo
enddo
enddo
print*,'wall time for three_e_3_idx_cycle_1_bi_ort',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_cycle_1_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the second cyclic permutation
!
! three_e_3_idx_direct_bi_ort(m,j,i) = <mji|-L|imj>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
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)
name_file = 'six_index_tensor'
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the second cyclic permutation
!
! three_e_3_idx_direct_bi_ort(m,j,i) = <mji|-L|imj>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0
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)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_3_idx_cycle_2_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m,j,i,i,m,j,integral)
three_e_3_idx_cycle_2_bi_ort(m,j,i) = -1.d0 * integral
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m, j, i, i, m, j, integral)
three_e_3_idx_cycle_2_bi_ort(m,j,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_3_idx_cycle_2_bi_ort(m,j,i) = three_e_3_idx_cycle_2_bi_ort(j,m,i)
do j = 1, mo_num
do m = 1, j
three_e_3_idx_cycle_2_bi_ort(m,j,i) = three_e_3_idx_cycle_2_bi_ort(j,m,i)
enddo
enddo
enddo
enddo
print*,'wall time for three_e_3_idx_cycle_2_bi_ort',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_cycle_2_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_3_idx_exch23_bi_ort, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 2 and 3
!
! three_e_3_idx_exch23_bi_ort(m,j,i) = <mji|-L|jmi>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_3_idx_exch23_bi_ort = 0.d0
print*,'Providing the three_e_3_idx_exch23_bi_ort ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 2 and 3
!
! three_e_3_idx_exch23_bi_ort(m,j,i) = <mji|-L|jmi>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0
three_e_3_idx_exch23_bi_ort = 0.d0
print*,'Providing the three_e_3_idx_exch23_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_3_idx_exch23_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m,j,i,j,m,i,integral)
three_e_3_idx_exch23_bi_ort(m,j,i) = -1.d0 * integral
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m, j, i, j, m, i, integral)
three_e_3_idx_exch23_bi_ort(m,j,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_3_idx_exch23_bi_ort(m,j,i) = three_e_3_idx_exch23_bi_ort(j,m,i)
do j = 1, mo_num
do m = 1, j
three_e_3_idx_exch23_bi_ort(m,j,i) = three_e_3_idx_exch23_bi_ort(j,m,i)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for three_e_3_idx_exch23_bi_ort',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch23_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_3_idx_exch13_bi_ort, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 3
!
! three_e_3_idx_exch13_bi_ort(m,j,i) = <mji|-L|ijm>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_3_idx_exch13_bi_ort = 0.d0
print*,'Providing the three_e_3_idx_exch13_bi_ort ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 3
!
! three_e_3_idx_exch13_bi_ort(m,j,i) = <mji|-L|ijm>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i,j,m
double precision :: integral, wall1, wall0
three_e_3_idx_exch13_bi_ort = 0.d0
print *, ' Providing the three_e_3_idx_exch13_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_3_idx_exch13_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m,j,i,i,j,m,integral)
three_e_3_idx_exch13_bi_ort(m,j,i) = -1.d0 * integral
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m, j, i, i, j, m,integral)
three_e_3_idx_exch13_bi_ort(m,j,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_3_idx_exch13_bi_ort(m,j,i) = three_e_3_idx_exch13_bi_ort(j,m,i)
do j = 1, mo_num
do m = 1, j
three_e_3_idx_exch13_bi_ort(m,j,i) = three_e_3_idx_exch13_bi_ort(j,m,i)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for three_e_3_idx_exch13_bi_ort',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch13_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2
!
! three_e_3_idx_exch12_bi_ort(m,j,i) = <mji|-L|mij>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_3_idx_exch12_bi_ort = 0.d0
print*,'Providing the three_e_3_idx_exch12_bi_ort ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2
!
! three_e_3_idx_exch12_bi_ort(m,j,i) = <mji|-L|mij>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0
three_e_3_idx_exch12_bi_ort = 0.d0
print *, ' Providing the three_e_3_idx_exch12_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_3_idx_exch12_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,j,i,m,i,j,integral)
three_e_3_idx_exch12_bi_ort(m,j,i) = -1.d0 * integral
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, j, i, m, i, j, integral)
three_e_3_idx_exch12_bi_ort(m,j,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_3_idx_exch12_bi_ort',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch12_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort_new, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2
!
! three_e_3_idx_exch12_bi_ort_new(m,j,i) = <mji|-L|mij>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_3_idx_exch12_bi_ort_new = 0.d0
print*,'Providing the three_e_3_idx_exch12_bi_ort_new ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the permutations of particle 1 and 2
!
! three_e_3_idx_exch12_bi_ort_new(m,j,i) = <mji|-L|mij>
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0
three_e_3_idx_exch12_bi_ort_new = 0.d0
print *, ' Providing the three_e_3_idx_exch12_bi_ort_new ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_3_idx_exch12_bi_ort_new)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m,j,i,m,i,j,integral)
three_e_3_idx_exch12_bi_ort_new(m,j,i) = -1.d0 * integral
do j = 1, mo_num
do m = j, mo_num
call give_integrals_3_body_bi_ort(m, j, i, m, i, j, integral)
three_e_3_idx_exch12_bi_ort_new(m,j,i) = -1.d0 * integral
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_3_idx_exch12_bi_ort_new(m,j,i) = three_e_3_idx_exch12_bi_ort_new(j,m,i)
do j = 1, mo_num
do m = 1, j
three_e_3_idx_exch12_bi_ort_new(m,j,i) = three_e_3_idx_exch12_bi_ort_new(j,m,i)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print*,'wall time for three_e_3_idx_exch12_bi_ort_new',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_e_3_idx_exch12_bi_ort_new', wall1 - wall0
END_PROVIDER
! ---

428
src/bi_ort_ints/three_body_ijmk.irp.f
View File

@ -1,228 +1,284 @@
! ---
BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_direct_bi_ort(m,j,k,i) = <mjk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_4_idx_direct_bi_ort(m,j,k,i) = <mjk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m
integer :: i, j, k, m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_4_idx_direct_bi_ort = 0.d0
print*,'Providing the three_e_4_idx_direct_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_direct_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
three_e_4_idx_direct_bi_ort = 0.d0
print *, ' Providing the three_e_4_idx_direct_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_direct_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,j,k,m,j,i,integral)
three_e_4_idx_direct_bi_ort(m,j,k,i) = -1.d0 * integral
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, j, k, m, j, i, integral)
three_e_4_idx_direct_bi_ort(m,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_4_idx_direct_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_direct_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = <mjk|-L|jim> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_4_idx_cycle_1_bi_ort = 0.d0
print*,'Providing the three_e_4_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_cycle_1_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = <mjk|-L|jim> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0
three_e_4_idx_cycle_1_bi_ort = 0.d0
print *, ' Providing the three_e_4_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_cycle_1_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,j,k,j,i,m,integral)
three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = -1.d0 * integral
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, j, k, j, i, m, integral)
three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_4_idx_cycle_1_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_cycle_1_bi_ort', wall1 - wall0
END_PROVIDER
! --
BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = <mjk|-L|imj> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_4_idx_cycle_2_bi_ort = 0.d0
print*,'Providing the three_e_4_idx_cycle_2_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_cycle_2_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = <mjk|-L|imj> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0
three_e_4_idx_cycle_2_bi_ort = 0.d0
print *, ' Providing the three_e_4_idx_cycle_2_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_cycle_2_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,j,k,i,m,j,integral)
three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = -1.d0 * integral
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, j, k, i, m, j, integral)
three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_4_idx_cycle_2_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_cycle_2_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_4_idx_exch23_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_4_idx_exch23_bi_ort = 0.d0
print*,'Providing the three_e_4_idx_exch23_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_exch23_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,j,k,j,m,i,integral)
three_e_4_idx_exch23_bi_ort(m,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_4_idx_exch23_bi_ort',wall1 - wall0
END_PROVIDER
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_exch23_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0
three_e_4_idx_exch23_bi_ort = 0.d0
print *, ' Providing the three_e_4_idx_exch23_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_exch23_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, j, k, j, m, i, integral)
three_e_4_idx_exch23_bi_ort(m,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_exch23_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_4_idx_exch13_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_4_idx_exch13_bi_ort = 0.d0
print*,'Providing the three_e_4_idx_exch13_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_exch13_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_exch13_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0
three_e_4_idx_exch13_bi_ort = 0.d0
print *, ' Providing the three_e_4_idx_exch13_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_exch13_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,j,k,i,j,m,integral)
three_e_4_idx_exch13_bi_ort(m,j,k,i) = -1.d0 * integral
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, j, k, i, j, m, integral)
three_e_4_idx_exch13_bi_ort(m,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_4_idx_exch13_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_exch13_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_4_idx_exch12_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_4_idx_exch12_bi_ort = 0.d0
print*,'Providing the three_e_4_idx_exch12_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_exch12_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_exch12_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0
three_e_4_idx_exch12_bi_ort = 0.d0
print *, ' Providing the three_e_4_idx_exch12_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_4_idx_exch12_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,j,k,m,i,j,integral)
three_e_4_idx_exch12_bi_ort(m,j,k,i) = -1.d0 * integral
do j = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, j, k, m, i, j, integral)
three_e_4_idx_exch12_bi_ort(m,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_4_idx_exch12_bi_ort',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_exch12_bi_ort', wall1 - wall0
END_PROVIDER
! ---

452
src/bi_ort_ints/three_body_ijmkl.irp.f
View File

@ -1,240 +1,296 @@
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mjk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m,l
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_5_idx_direct_bi_ort = 0.d0
print*,'Providing the three_e_5_idx_direct_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mjk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0
three_e_5_idx_direct_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_direct_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,l,k,m,j,i,integral)
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, m, j, i, integral)
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_5_idx_direct_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_direct_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = <mlk|-L|jim> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m,l
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_5_idx_cycle_1_bi_ort = 0.d0
print*,'Providing the three_e_5_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = <mlk|-L|jim> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0
three_e_5_idx_cycle_1_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,l,k,j,i,m,integral)
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, j, i, m, integral)
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_5_idx_cycle_1_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_1_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = <mlk|-L|imj> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m,l
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_5_idx_cycle_2_bi_ort = 0.d0
print*,'Providing the three_e_5_idx_cycle_2_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = <mlk|-L|imj> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0
three_e_5_idx_cycle_2_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_cycle_2_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
do l = 1, mo_num
call give_integrals_3_body_bi_ort(m,l,k,i,m,j,integral)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
do j = 1, mo_num
do m = 1, mo_num
do l = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, i, m, j, integral)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_5_idx_cycle_2_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_2_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m,l
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_5_idx_exch23_bi_ort = 0.d0
print*,'Providing the three_e_5_idx_exch23_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0
three_e_5_idx_exch23_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_exch23_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,l,k,j,m,i,integral)
three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, j, m, i, integral)
three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_5_idx_exch23_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch23_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m,l
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_5_idx_exch13_bi_ort = 0.d0
print*,'Providing the three_e_5_idx_exch13_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0
three_e_5_idx_exch13_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_exch13_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,l,k,i,j,m,integral)
three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, i, j, m, integral)
three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_5_idx_exch13_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch13_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
!three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
integer :: i,j,k,m,l
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_e_5_idx_exch12_bi_ort = 0.d0
print*,'Providing the three_e_5_idx_exch12_bi_ort ...'
call wall_time(wall0)
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0
three_e_5_idx_exch12_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_exch12_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m,l,k,m,i,j,integral)
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, m, i, j, integral)
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'wall time for three_e_5_idx_exch12_bi_ort',wall1 - wall0
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch12_bi_ort', wall1 - wall0
END_PROVIDER
! ---

105
src/bi_ort_ints/three_body_ints_bi_ort.irp.f
View File

@ -1,17 +1,24 @@
! ---
BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator
!
! notice the -1 sign: in this way three_body_ints_bi_ort can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
implicit none
integer :: i, j, k, l, m, n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_ints_bi_ort = 0.d0
print*,'Providing the three_body_ints_bi_ort ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
character*(128) :: name_file
three_body_ints_bi_ort = 0.d0
print*,'Providing the three_body_ints_bi_ort ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
! if(read_three_body_ints_bi_ort)then
! call read_fcidump_3_tc(three_body_ints_bi_ort)
! else
@ -19,32 +26,37 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n
! print*,'Reading three_body_ints_bi_ort from disk ...'
! call read_array_6_index_tensor(mo_num,three_body_ints_bi_ort,name_file)
! else
provide x_W_ki_bi_ortho_erf_rk mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_ints_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
!provide x_W_ki_bi_ortho_erf_rk
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_ints_bi_ort)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
do n = 1, mo_num
call give_integrals_3_body_bi_ort(n,l,k,m,j,i,integral)
three_body_ints_bi_ort(n,l,k,m,j,i) = -1.d0 * integral
do m = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
do n = 1, mo_num
call give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
three_body_ints_bi_ort(n,l,k,m,j,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
enddo
!$OMP END DO
!$OMP END PARALLEL
! endif
! endif
call wall_time(wall1)
print*,'wall time for three_body_ints_bi_ort',wall1 - wall0
call wall_time(wall1)
print *, ' wall time for three_body_ints_bi_ort', wall1 - wall0
! if(write_three_body_ints_bi_ort)then
! print*,'Writing three_body_ints_bi_ort on disk ...'
! call write_array_6_index_tensor(mo_num,three_body_ints_bi_ort,name_file)
@ -64,7 +76,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
END_DOC
implicit none
integer, intent(in) :: n,l,k,m,j,i
integer, intent(in) :: n, l, k, m, j, i
double precision, intent(out) :: integral
integer :: ipoint
double precision :: weight
@ -86,18 +98,31 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
! + x_W_ki_bi_ortho_erf_rk(ipoint,2,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,2,k,i) &
! + x_W_ki_bi_ortho_erf_rk(ipoint,3,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,3,k,i) )
integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
* ( int2_grad1_u12_bimo(1,ipoint,n,m) * int2_grad1_u12_bimo(1,ipoint,l,j) &
+ int2_grad1_u12_bimo(2,ipoint,n,m) * int2_grad1_u12_bimo(2,ipoint,l,j) &
+ int2_grad1_u12_bimo(3,ipoint,n,m) * int2_grad1_u12_bimo(3,ipoint,l,j) )
integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
* ( int2_grad1_u12_bimo(1,ipoint,n,m) * int2_grad1_u12_bimo(1,ipoint,k,i) &
+ int2_grad1_u12_bimo(2,ipoint,n,m) * int2_grad1_u12_bimo(2,ipoint,k,i) &
+ int2_grad1_u12_bimo(3,ipoint,n,m) * int2_grad1_u12_bimo(3,ipoint,k,i) )
integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
* ( int2_grad1_u12_bimo(1,ipoint,l,j) * int2_grad1_u12_bimo(1,ipoint,k,i) &
+ int2_grad1_u12_bimo(2,ipoint,l,j) * int2_grad1_u12_bimo(2,ipoint,k,i) &
+ int2_grad1_u12_bimo(3,ipoint,l,j) * int2_grad1_u12_bimo(3,ipoint,k,i) )
! integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
! * ( int2_grad1_u12_bimo(1,n,m,ipoint) * int2_grad1_u12_bimo(1,l,j,ipoint) &
! + int2_grad1_u12_bimo(2,n,m,ipoint) * int2_grad1_u12_bimo(2,l,j,ipoint) &
! + int2_grad1_u12_bimo(3,n,m,ipoint) * int2_grad1_u12_bimo(3,l,j,ipoint) )
! integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
! * ( int2_grad1_u12_bimo(1,n,m,ipoint) * int2_grad1_u12_bimo(1,k,i,ipoint) &
! + int2_grad1_u12_bimo(2,n,m,ipoint) * int2_grad1_u12_bimo(2,k,i,ipoint) &
! + int2_grad1_u12_bimo(3,n,m,ipoint) * int2_grad1_u12_bimo(3,k,i,ipoint) )
! integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
! * ( int2_grad1_u12_bimo(1,l,j,ipoint) * int2_grad1_u12_bimo(1,k,i,ipoint) &
! + int2_grad1_u12_bimo(2,l,j,ipoint) * int2_grad1_u12_bimo(2,k,i,ipoint) &
! + int2_grad1_u12_bimo(3,l,j,ipoint) * int2_grad1_u12_bimo(3,k,i,ipoint) )
integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
* ( int2_grad1_u12_bimo_transp(n,m,1,ipoint) * int2_grad1_u12_bimo_transp(l,j,1,ipoint) &
+ int2_grad1_u12_bimo_transp(n,m,2,ipoint) * int2_grad1_u12_bimo_transp(l,j,2,ipoint) &
+ int2_grad1_u12_bimo_transp(n,m,3,ipoint) * int2_grad1_u12_bimo_transp(l,j,3,ipoint) )
integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
* ( int2_grad1_u12_bimo_transp(n,m,1,ipoint) * int2_grad1_u12_bimo_transp(k,i,1,ipoint) &
+ int2_grad1_u12_bimo_transp(n,m,2,ipoint) * int2_grad1_u12_bimo_transp(k,i,2,ipoint) &
+ int2_grad1_u12_bimo_transp(n,m,3,ipoint) * int2_grad1_u12_bimo_transp(k,i,3,ipoint) )
integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
* ( int2_grad1_u12_bimo_transp(l,j,1,ipoint) * int2_grad1_u12_bimo_transp(k,i,1,ipoint) &
+ int2_grad1_u12_bimo_transp(l,j,2,ipoint) * int2_grad1_u12_bimo_transp(k,i,2,ipoint) &
+ int2_grad1_u12_bimo_transp(l,j,3,ipoint) * int2_grad1_u12_bimo_transp(k,i,3,ipoint) )
enddo

45
src/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -17,23 +17,42 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
double precision :: integral_sym, integral_nsym
double precision, external :: get_ao_tc_sym_two_e_pot
PROVIDE ao_tc_sym_two_e_pot_in_map
provide j1b_type
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
if(j1b_type .eq. 3) then
integral_sym = get_ao_tc_sym_two_e_pot(i,j,k,l,ao_tc_sym_two_e_pot_map)
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
integral_nsym = ao_non_hermit_term_chemist(k,i,l,j)
ao_two_e_tc_tot(k,i,l,j) = integral_sym + integral_nsym
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_two_e_tc_tot(k,i,l,j) = ao_tc_int_chemist(k,i,l,j)
!write(222,*) ao_two_e_tc_tot(k,i,l,j)
enddo
enddo
enddo
enddo
enddo
else
PROVIDE ao_tc_sym_two_e_pot_in_map
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
integral_sym = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map)
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
integral_nsym = ao_non_hermit_term_chemist(k,i,l,j)
ao_two_e_tc_tot(k,i,l,j) = integral_sym + integral_nsym
!write(111,*) ao_two_e_tc_tot(k,i,l,j)
enddo
enddo
enddo
enddo
endif
END_PROVIDER
@ -42,9 +61,11 @@ END_PROVIDER
double precision function bi_ortho_mo_ints(l, k, j, i)
BEGIN_DOC
!
! <mo^L_k mo^L_l | V^TC(r_12) | mo^R_i mo^R_j>
!
! WARNING :: very naive, super slow, only used to DEBUG.
!
END_DOC
implicit none

28
src/bi_ortho_mos/mos_rl.irp.f
View File

@ -1,33 +1,37 @@
! ---
subroutine ao_to_mo_bi_ortho(A_ao, LDA_ao, A_mo, LDA_mo)
BEGIN_DOC
!
! Transform A from the |AO| basis to the BI ORTHONORMAL MOS
!
! $C_L^\dagger.A_{ao}.C_R$ where C_L and C_R are the LEFT and RIGHT MO coefs
!
END_DOC
implicit none
integer, intent(in) :: LDA_ao,LDA_mo
integer, intent(in) :: LDA_ao, LDA_mo
double precision, intent(in) :: A_ao(LDA_ao,ao_num)
double precision, intent(out) :: A_mo(LDA_mo,mo_num)
double precision, allocatable :: T(:,:)
allocate ( T(ao_num,mo_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
integer :: i,j,p,q
call dgemm('N', 'N', ao_num, mo_num, ao_num, &
1.d0, A_ao, LDA_ao, &
mo_r_coef, size(mo_r_coef, 1), &
0.d0, T, size(T, 1))
! T = A_ao x mo_r_coef
call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0 &
, A_ao, LDA_ao, mo_r_coef, size(mo_r_coef, 1) &
, 0.d0, T, size(T, 1) )
call dgemm('T', 'N', mo_num, mo_num, ao_num, &
1.d0, mo_l_coef, size(mo_l_coef, 1), &
T, ao_num, &
0.d0, A_mo, size(A_mo, 1))
! A_mo = mo_l_coef.T x T
call dgemm( 'T', 'N', mo_num, mo_num, ao_num, 1.d0 &
, mo_l_coef, size(mo_l_coef, 1), T, size(T, 1) &
, 0.d0, A_mo, LDA_mo )
! call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-12)
deallocate(T)
deallocate(T)
end subroutine ao_to_mo_bi_ortho
@ -131,7 +135,7 @@ BEGIN_PROVIDER [ double precision, mo_l_coef, (ao_num, mo_num) ]
IRP_ENDIF
else
print*, 'mo_r_coef are mo_coef'
print*, 'mo_l_coef are mo_coef'
do i = 1, mo_num
do j = 1, ao_num
mo_l_coef(j,i) = mo_coef(j,i)

512
src/non_h_ints_mu/debug_fit.irp.f Normal file
View File

@ -0,0 +1,512 @@
! --
program debug_fit
implicit none
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
!my_n_pt_r_grid = 100
!my_n_pt_a_grid = 170
!my_n_pt_r_grid = 150
!my_n_pt_a_grid = 194
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
PROVIDE mu_erf j1b_pen
!call test_j1b_nucl()
call test_grad_j1b_nucl()
!call test_lapl_j1b_nucl()
!call test_list_b2()
!call test_list_b3()
call test_fit_u()
!call test_fit_u2()
!call test_fit_ugradu()
end
! ---
subroutine test_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_j1b_nucl ...'
PROVIDE v_1b
eps_ij = 1d-7
acc_tot = 0.d0
normalz = 0.d0
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)
i_exc = v_1b(ipoint)
i_num = j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in v_1b on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_j1b_nucl
! ---
subroutine test_grad_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
print*, ' test_grad_j1b_nucl ...'
PROVIDE v_1b_grad
eps_ij = 1d-7
acc_tot = 0.d0
normalz = 0.d0
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)
i_exc = v_1b_grad(1,ipoint)
i_num = grad_x_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in x of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
i_exc = v_1b_grad(2,ipoint)
i_num = grad_y_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in y of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
i_exc = v_1b_grad(3,ipoint)
i_num = grad_z_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in z of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_grad_j1b_nucl
! ---
subroutine test_lapl_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: lapl_j1b_nucl
print*, ' test_lapl_j1b_nucl ...'
PROVIDE v_1b_lapl
eps_ij = 1d-5
acc_tot = 0.d0
normalz = 0.d0
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)
i_exc = v_1b_lapl(ipoint)
i_num = lapl_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in v_1b_lapl on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_lapl_j1b_nucl
! ---
subroutine test_list_b2()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_list_b2 ...'
PROVIDE v_1b_list_b2
eps_ij = 1d-7
acc_tot = 0.d0
normalz = 0.d0
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)
i_exc = v_1b_list_b2(ipoint)
i_num = j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in list_b2 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_list_b2
! ---
subroutine test_list_b3()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_list_b3 ...'
PROVIDE v_1b_list_b3
eps_ij = 1d-7
acc_tot = 0.d0
normalz = 0.d0
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)
i_exc = v_1b_list_b3(ipoint)
i_tmp = j1b_nucl(r)
i_num = i_tmp * i_tmp
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in list_b3 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_list_b3
! ---
subroutine test_fit_ugradu()
implicit none
integer :: jpoint, ipoint, i
double precision :: i_exc, i_fit, i_num, x2, tmp, dx, dy, dz
double precision :: r1(3), r2(3), grad(3)
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
double precision, external :: j12_mu
print*, ' test_fit_ugradu ...'
eps_ij = 1d-3
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
x2 = dx * dx + dy * dy + dz * dz
if(x2 .lt. 1d-10) cycle
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_1_erf(i)
coef = coef_gauss_j_mu_1_erf(i)
i_fit += coef * dexp(-expo*x2)
enddo
i_fit = i_fit / dsqrt(x2)
tmp = j12_mu(r1, r2)
call grad1_j12_mu_exc(r1, r2, grad)
! ---
i_exc = tmp * grad(1)
i_num = i_fit * dx
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on x in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
i_exc = tmp * grad(2)
i_num = i_fit * dy
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on y in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
i_exc = tmp * grad(3)
i_num = i_fit * dz
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on z in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
enddo
if( (acc_tot/normalz) .gt. 1d-3 ) then
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
endif
enddo
return
end subroutine test_fit_ugradu
! ---
subroutine test_fit_u()
implicit none
integer :: jpoint, ipoint, i
double precision :: i_exc, i_fit, i_num, x2
double precision :: r1(3), r2(3), dx, dy, dz
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
double precision, external :: j12_mu
print*, ' test_fit_u ...'
eps_ij = 1d-3
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
x2 = dx * dx + dy * dy + dz * dz
if(x2 .lt. 1d-10) cycle
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_x(i)
coef = coef_gauss_j_mu_x(i)
i_fit += coef * dexp(-expo*x2)
enddo
i_exc = j12_mu(r1, r2)
i_num = i_fit
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in test_fit_u on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
enddo
if( (acc_tot/normalz) .gt. 1d-3 ) then
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
endif
enddo
return
end subroutine test_fit_u
! ---
subroutine test_fit_u2()
implicit none
integer :: jpoint, ipoint, i
double precision :: i_exc, i_fit, i_num, x2
double precision :: r1(3), r2(3), dx, dy, dz, tmp
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
double precision, external :: j12_mu
print*, ' test_fit_u2 ...'
eps_ij = 1d-3
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
x2 = dx * dx + dy * dy + dz * dz
if(x2 .lt. 1d-10) cycle
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_x_2(i)
coef = coef_gauss_j_mu_x_2(i)
i_fit += coef * dexp(-expo*x2)
enddo
tmp = j12_mu(r1, r2)
i_exc = tmp * tmp
i_num = i_fit
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in test_fit_u2 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
enddo
if( (acc_tot/normalz) .gt. 1d-3 ) then
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
endif
enddo
return
end subroutine test_fit_u2
! ---

632
src/non_h_ints_mu/debug_integ_jmu_modif.irp.f
View File

@ -17,25 +17,19 @@ program debug_integ_jmu_modif
PROVIDE mu_erf j1b_pen
!call test_j1b_nucl()
!call test_grad_j1b_nucl()
!call test_lapl_j1b_nucl()
!call test_list_b2()
!call test_list_b3()
!call test_fit_u()
call test_fit_u2()
!call test_fit_ugradu()
!call test_v_ij_u_cst_mu_j1b()
!call test_v_ij_erf_rk_cst_mu_j1b()
!call test_x_v_ij_erf_rk_cst_mu_j1b()
!call test_int2_u2_j1b2()
!call test_int2_grad1u2_grad2u2_j1b2()
!call test_int2_grad1_u12_ao()
!call test_gradu_squared_u_ij_mu()
call test_v_ij_u_cst_mu_j1b()
! call test_v_ij_erf_rk_cst_mu_j1b()
! call test_x_v_ij_erf_rk_cst_mu_j1b()
! call test_int2_u2_j1b2()
! call test_int2_grad1u2_grad2u2_j1b2()
! call test_int2_u_grad1u_total_j1b2()
!
! call test_int2_grad1_u12_ao()
!
! call test_grad12_j12()
! call test_u12sq_j1bsq()
! call test_u12_grad1_u12_j1b_grad1_j1b()
! !call test_gradu_squared_u_ij_mu()
end
@ -52,8 +46,9 @@ subroutine test_v_ij_u_cst_mu_j1b()
PROVIDE v_ij_u_cst_mu_j1b
eps_ij = 1d-8
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
@ -76,9 +71,8 @@ subroutine test_v_ij_u_cst_mu_j1b()
enddo
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_v_ij_u_cst_mu_j1b
@ -96,8 +90,9 @@ subroutine test_v_ij_erf_rk_cst_mu_j1b()
PROVIDE v_ij_erf_rk_cst_mu_j1b
eps_ij = 1d-8
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
@ -120,9 +115,8 @@ subroutine test_v_ij_erf_rk_cst_mu_j1b()
enddo
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_v_ij_erf_rk_cst_mu_j1b
@ -140,8 +134,9 @@ subroutine test_x_v_ij_erf_rk_cst_mu_j1b()
PROVIDE x_v_ij_erf_rk_cst_mu_j1b
eps_ij = 1d-8
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
@ -190,9 +185,8 @@ subroutine test_x_v_ij_erf_rk_cst_mu_j1b()
enddo
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_x_v_ij_erf_rk_cst_mu_j1b
@ -210,8 +204,9 @@ subroutine test_int2_u2_j1b2()
PROVIDE int2_u2_j1b2
eps_ij = 1d-8
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
@ -235,8 +230,8 @@ subroutine test_int2_u2_j1b2()
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_int2_u2_j1b2
@ -254,8 +249,9 @@ subroutine test_int2_grad1u2_grad2u2_j1b2()
PROVIDE int2_grad1u2_grad2u2_j1b2
eps_ij = 1d-8
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
@ -278,9 +274,8 @@ subroutine test_int2_grad1u2_grad2u2_j1b2()
enddo
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_int2_grad1u2_grad2u2_j1b2
@ -298,8 +293,9 @@ subroutine test_int2_grad1_u12_ao()
PROVIDE int2_grad1_u12_ao
eps_ij = 1d-6
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
@ -347,15 +343,90 @@ subroutine test_int2_grad1_u12_ao()
enddo
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_int2_grad1_u12_ao
! ---
subroutine test_int2_u_grad1u_total_j1b2()
implicit none
integer :: i, j, ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: x, y, z
double precision :: integ(3)
print*, ' test_int2_u_grad1u_total_j1b2 ...'
PROVIDE int2_u_grad1u_j1b2
PROVIDE int2_u_grad1u_x_j1b2
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
do j = 1, ao_num
do i = 1, ao_num
call num_int2_u_grad1u_total_j1b2(i, j, ipoint, integ)
i_exc = x * int2_u_grad1u_j1b2(i,j,ipoint) - int2_u_grad1u_x_j1b2(1,i,j,ipoint)
i_num = integ(1)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in x part of int2_u_grad1u_total_j1b2 on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
i_exc = y * int2_u_grad1u_j1b2(i,j,ipoint) - int2_u_grad1u_x_j1b2(2,i,j,ipoint)
i_num = integ(2)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in y part of int2_u_grad1u_total_j1b2 on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
i_exc = z * int2_u_grad1u_j1b2(i,j,ipoint) - int2_u_grad1u_x_j1b2(3,i,j,ipoint)
i_num = integ(3)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in z part of int2_u_grad1u_total_j1b2 on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
enddo
enddo
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_int2_u_grad1u_total_j1b2
! ---
subroutine test_gradu_squared_u_ij_mu()
implicit none
@ -367,8 +438,9 @@ subroutine test_gradu_squared_u_ij_mu()
PROVIDE gradu_squared_u_ij_mu
eps_ij = 1d-6
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
@ -390,458 +462,140 @@ subroutine test_gradu_squared_u_ij_mu()
enddo
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_gradu_squared_u_ij_mu
! ---
subroutine test_j1b_nucl()
subroutine test_grad12_j12()
implicit none
integer :: ipoint
integer :: i, j, ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
double precision, external :: num_grad12_j12
print*, ' test_j1b_nucl ...'
print*, ' test_grad12_j12 ...'
PROVIDE v_1b
PROVIDE grad12_j12
eps_ij = 1d-7
eps_ij = 1d-3
acc_tot = 0.d0
normalz = 0.d0
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)
i_exc = v_1b(ipoint)
i_num = j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in v_1b on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_j1b_nucl
! ---
subroutine test_grad_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
print*, ' test_grad_j1b_nucl ...'
PROVIDE v_1b_grad
eps_ij = 1d-6
acc_tot = 0.d0
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)
i_exc = v_1b_grad(1,ipoint)
i_num = grad_x_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in x of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
i_exc = v_1b_grad(2,ipoint)
i_num = grad_y_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in y of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
i_exc = v_1b_grad(3,ipoint)
i_num = grad_z_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in z of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_grad_j1b_nucl
! ---
subroutine test_lapl_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: lapl_j1b_nucl
print*, ' test_lapl_j1b_nucl ...'
PROVIDE v_1b_lapl
eps_ij = 1d-5
acc_tot = 0.d0
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)
i_exc = v_1b_lapl(ipoint)
i_num = lapl_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in v_1b_lapl on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_lapl_j1b_nucl
! ---
subroutine test_list_b2()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_list_b2 ...'
PROVIDE v_1b_list_b2
eps_ij = 1d-7
acc_tot = 0.d0
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)
i_exc = v_1b_list_b2(ipoint)
i_num = j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in list_b2 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_list_b2
! ---
subroutine test_list_b3()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_list_b3 ...'
PROVIDE v_1b_list_b3
eps_ij = 1d-7
acc_tot = 0.d0
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)
i_exc = v_1b_list_b3(ipoint)
i_tmp = j1b_nucl(r)
i_num = i_tmp * i_tmp
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in list_b3 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_list_b3
! ---
subroutine test_fit_ugradu()
implicit none
integer :: ipoint, i
double precision :: i_exc, i_fit, i_num, x2
double precision :: r1(3), r2(3), grad(3)
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
double precision, external :: j12_mu
print*, ' test_fit_ugradu ...'
eps_ij = 1d-7
acc_tot = 0.d0
r2 = 0.d0
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3)
if(x2 .lt. 1d-10) cycle
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_1_erf(i)
coef = coef_gauss_j_mu_1_erf(i)
i_fit += coef * dexp(-expo*x2)
do j = 1, ao_num
do i = 1, ao_num
i_exc = grad12_j12(i,j,ipoint)
i_num = num_grad12_j12(i, j, ipoint)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in grad12_j12 on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
enddo
i_fit = i_fit / dsqrt(x2)
call grad1_j12_mu_exc(r1, r2, grad)
! ---
i_exc = j12_mu(r1, r2) * grad(1)
i_num = i_fit * r1(1)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on x in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
i_exc = j12_mu(r1, r2) * grad(2)
i_num = i_fit * r1(2)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on y in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
i_exc = j12_mu(r1, r2) * grad(3)
i_num = i_fit * r1(3)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem on z in test_fit_ugradu on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
! ---
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_fit_ugradu
end subroutine test_grad12_j12
! ---
subroutine test_fit_u()
subroutine test_u12sq_j1bsq()
implicit none
integer :: i, j, ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision, external :: num_u12sq_j1bsq
integer :: ipoint, i
double precision :: i_exc, i_fit, i_num, x2
double precision :: r1(3), r2(3)
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
print*, ' test_u12sq_j1bsq ...'
double precision, external :: j12_mu
PROVIDE u12sq_j1bsq
print*, ' test_fit_u ...'
eps_ij = 1d-7
eps_ij = 1d-3
acc_tot = 0.d0
r2 = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
do i = 1, ao_num
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3)
if(x2 .lt. 1d-10) cycle
i_exc = u12sq_j1bsq(i,j,ipoint)
i_num = num_u12sq_j1bsq(i, j, ipoint)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in u12sq_j1bsq on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_x(i)
coef = coef_gauss_j_mu_x(i)
i_fit += coef * dexp(-expo*x2)
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
enddo
i_exc = j12_mu(r1, r2)
i_num = i_fit
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in test_fit_u on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_fit_u
end subroutine test_u12sq_j1bsq
! ---
subroutine test_fit_u2()
subroutine test_u12_grad1_u12_j1b_grad1_j1b()
implicit none
integer :: i, j, ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision, external :: num_u12_grad1_u12_j1b_grad1_j1b
integer :: ipoint, i
double precision :: i_exc, i_fit, i_num, x2
double precision :: r1(3), r2(3)
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
print*, ' test_u12_grad1_u12_j1b_grad1_j1b ...'
double precision, external :: j12_mu
PROVIDE u12_grad1_u12_j1b_grad1_j1b
print*, ' test_fit_u2 ...'
eps_ij = 1d-7
eps_ij = 1d-3
acc_tot = 0.d0
r2 = 0.d0
normalz = 0.d0
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
do i = 1, ao_num
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3)
if(x2 .lt. 1d-10) cycle
i_exc = u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint)
i_num = num_u12_grad1_u12_j1b_grad1_j1b(i, j, ipoint)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in u12_grad1_u12_j1b_grad1_j1b on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
i_fit = 0.d0
do i = 1, n_max_fit_slat
expo = expo_gauss_j_mu_x_2(i)
coef = coef_gauss_j_mu_x_2(i)
i_fit += coef * dexp(-expo*x2)
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
enddo
i_exc = j12_mu(r1, r2) * j12_mu(r1, r2)
i_num = i_fit
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in test_fit_u2 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_exc)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
print*, ' acc_tot = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_fit_u2
end subroutine test_u12_grad1_u12_j1b_grad1_j1b,
! ---

238
src/non_h_ints_mu/grad_squared.irp.f
View File

@ -23,52 +23,63 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
! + -1.00 x v1 (grad_1 v1) \int r2 \phi_i(2) \phi_j(2) (grad_1 u12) v2^2
! = v1^2 x int2_grad1u2_grad2u2_j1b2
! + -0.5 x (grad_1 v1)^2 x int2_u2_j1b2
! + -1.0 X V1 x (grad_1 v1) \cdot int2_u_grad1u_x_j1b
! + -1.0 X V1 x (grad_1 v1) \cdot [ int2_u_grad1u_j1b2 x r - int2_u_grad1u_x_j1b ]
!
!
END_DOC
implicit none
integer :: ipoint, i, j, m, igauss
double precision :: r(3), delta, coef
double precision :: tmp_v, tmp_x, tmp_y, tmp_z, tmp1, tmp2, tmp3, tmp4, tmp5
double precision :: x, y, z, r(3), delta, coef
double precision :: tmp_v, tmp_x, tmp_y, tmp_z
double precision :: tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9
double precision :: time0, time1
double precision, external :: overlap_gauss_r12_ao
print*, ' providing gradu_squared_u_ij_mu ...'
call wall_time(time0)
PROVIDE j1b_type j1b_pen
PROVIDE j1b_type
if(j1b_type .eq. 3) then
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
tmp_v = v_1b (ipoint)
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint)
tmp1 = tmp_v * tmp_v
tmp2 = 0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
tmp2 = -0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
tmp3 = tmp_v * tmp_x
tmp4 = tmp_v * tmp_y
tmp5 = tmp_v * tmp_z
tmp6 = -x * tmp3
tmp7 = -y * tmp4
tmp8 = -z * tmp5
do j = 1, ao_num
do i = 1, ao_num
gradu_squared_u_ij_mu(j,i,ipoint) += tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint) &
- tmp2 * int2_u2_j1b2 (i,j,ipoint) &
- tmp3 * int2_u_grad1u_x_j1b (1,i,j,ipoint) &
- tmp4 * int2_u_grad1u_x_j1b (2,i,j,ipoint) &
- tmp5 * int2_u_grad1u_x_j1b (3,i,j,ipoint)
tmp9 = int2_u_grad1u_j1b2(i,j,ipoint)
gradu_squared_u_ij_mu(i,j,ipoint) = tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint) &
+ tmp2 * int2_u2_j1b2 (i,j,ipoint) &
+ tmp6 * tmp9 + tmp3 * int2_u_grad1u_x_j1b2(1,i,j,ipoint) &
+ tmp7 * tmp9 + tmp4 * int2_u_grad1u_x_j1b2(2,i,j,ipoint) &
+ tmp8 * tmp9 + tmp5 * int2_u_grad1u_x_j1b2(3,i,j,ipoint)
enddo
enddo
enddo
else
gradu_squared_u_ij_mu = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
@ -78,7 +89,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
do igauss = 1, n_max_fit_slat
delta = expo_gauss_1_erf_x_2(igauss)
coef = coef_gauss_1_erf_x_2(igauss)
gradu_squared_u_ij_mu(j,i,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j)
gradu_squared_u_ij_mu(i,j,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j)
enddo
enddo
enddo
@ -93,6 +104,57 @@ END_PROVIDER
! ---
!BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)]
!
! BEGIN_DOC
! !
! ! tc_grad_square_ao(k,i,l,j) = -1/2 <kl | |\grad_1 u(r1,r2)|^2 + |\grad_1 u(r1,r2)|^2 | ij>
! !
! END_DOC
!
! implicit none
! integer :: ipoint, i, j, k, l
! double precision :: weight1, ao_ik_r, ao_i_r
! double precision, allocatable :: ac_mat(:,:,:,:)
!
! allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
! ac_mat = 0.d0
!
! do ipoint = 1, n_points_final_grid
! weight1 = final_weight_at_r_vector(ipoint)
!
! do i = 1, ao_num
! ao_i_r = weight1 * aos_in_r_array_transp(ipoint,i)
!
! do k = 1, ao_num
! ao_ik_r = ao_i_r * aos_in_r_array_transp(ipoint,k)
!
! do j = 1, ao_num
! do l = 1, ao_num
! ac_mat(k,i,l,j) += ao_ik_r * gradu_squared_u_ij_mu(l,j,ipoint)
! enddo
! enddo
! enddo
! enddo
! enddo
!
! do j = 1, ao_num
! do l = 1, ao_num
! do i = 1, ao_num
! do k = 1, ao_num
! tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
! !write(11,*) tc_grad_square_ao(k,i,l,j)
! enddo
! enddo
! enddo
! enddo
!
! deallocate(ac_mat)
!
!END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)]
BEGIN_DOC
@ -103,27 +165,27 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
implicit none
integer :: ipoint, i, j, k, l
double precision :: contrib, weight1, ao_k_r, ao_i_r
double precision, allocatable :: ac_mat(:,:,:,:)
double precision :: weight1, ao_ik_r, ao_i_r
double precision, allocatable :: ac_mat(:,:,:,:), bc_mat(:,:,:,:)
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
ac_mat = 0.d0
allocate(bc_mat(ao_num,ao_num,ao_num,ao_num))
bc_mat = 0.d0
do ipoint = 1, n_points_final_grid
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
weight1 = final_weight_at_r_vector(ipoint)
do i = 1, ao_num
ao_i_r = aos_in_r_array_transp(ipoint,i)
ao_i_r = weight1 * aos_in_r_array_transp(ipoint,i)
do k = 1, ao_num
ao_k_r = aos_in_r_array_transp(ipoint,k)
ao_ik_r = ao_i_r * aos_in_r_array_transp(ipoint,k)
do j = 1, ao_num
do l = 1, ao_num
contrib = gradu_squared_u_ij_mu(l,j,ipoint) * ao_k_r * ao_i_r
ac_mat(k,i,l,j) += weight1 * contrib
ac_mat(k,i,l,j) += ao_ik_r * ( u12sq_j1bsq(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(l,j,ipoint) )
bc_mat(k,i,l,j) += ao_ik_r * grad12_j12(l,j,ipoint)
enddo
enddo
enddo
@ -134,13 +196,147 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + bc_mat(k,i,l,j)
enddo
enddo
enddo
enddo
deallocate(ac_mat)
deallocate(bc_mat)
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_grid) ]
implicit none
integer :: ipoint, i, j, m, igauss
double precision :: r(3), delta, coef
double precision :: tmp1
double precision :: time0, time1
double precision, external :: overlap_gauss_r12_ao
print*, ' providing grad12_j12 ...'
call wall_time(time0)
PROVIDE j1b_type
if(j1b_type .eq. 3) then
do ipoint = 1, n_points_final_grid
tmp1 = v_1b(ipoint)
tmp1 = tmp1 * tmp1
do j = 1, ao_num
do i = 1, ao_num
grad12_j12(i,j,ipoint) = tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
enddo
enddo
enddo
else
grad12_j12 = 0.d0
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
do igauss = 1, n_max_fit_slat
delta = expo_gauss_1_erf_x_2(igauss)
coef = coef_gauss_1_erf_x_2(igauss)
grad12_j12(i,j,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j)
enddo
enddo
enddo
enddo
endif
call wall_time(time1)
print*, ' Wall time for grad12_j12 = ', time1 - time0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid) ]
implicit none
integer :: ipoint, i, j
double precision :: tmp_x, tmp_y, tmp_z
double precision :: tmp1
double precision :: time0, time1
print*, ' providing u12sq_j1bsq ...'
call wall_time(time0)
do ipoint = 1, n_points_final_grid
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint)
tmp1 = -0.5d0 * (tmp_x * tmp_x + tmp_y * tmp_y + tmp_z * tmp_z)
do j = 1, ao_num
do i = 1, ao_num
u12sq_j1bsq(i,j,ipoint) = tmp1 * int2_u2_j1b2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(time1)
print*, ' Wall time for u12sq_j1bsq = ', time1 - time0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ]
implicit none
integer :: ipoint, i, j, m, igauss
double precision :: x, y, z
double precision :: tmp_v, tmp_x, tmp_y, tmp_z
double precision :: tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9
double precision :: time0, time1
double precision, external :: overlap_gauss_r12_ao
print*, ' providing u12_grad1_u12_j1b_grad1_j1b ...'
call wall_time(time0)
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
tmp_v = v_1b (ipoint)
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint)
tmp3 = tmp_v * tmp_x
tmp4 = tmp_v * tmp_y
tmp5 = tmp_v * tmp_z
tmp6 = -x * tmp3
tmp7 = -y * tmp4
tmp8 = -z * tmp5
do j = 1, ao_num
do i = 1, ao_num
tmp9 = int2_u_grad1u_j1b2(i,j,ipoint)
u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) = tmp6 * tmp9 + tmp3 * int2_u_grad1u_x_j1b2(1,i,j,ipoint) &
+ tmp7 * tmp9 + tmp4 * int2_u_grad1u_x_j1b2(2,i,j,ipoint) &
+ tmp8 * tmp9 + tmp5 * int2_u_grad1u_x_j1b2(3,i,j,ipoint)
enddo
enddo
enddo
call wall_time(time1)
print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b = ', time1 - time0
END_PROVIDER

290
src/non_h_ints_mu/grad_tc_int.irp.f
View File

@ -1,67 +1,75 @@
! ---
BEGIN_PROVIDER [double precision, ao_non_hermit_term_chemist, (ao_num, ao_num, ao_num, ao_num)]
implicit none
BEGIN_DOC
! 1 1 2 2 1 2 1 2
!
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
END_DOC
integer :: i,j,k,l,ipoint,m
double precision :: weight1,thr,r(3)
thr = 1.d-8
double precision, allocatable :: b_mat(:,:,:,:),ac_mat(:,:,:,:)
! provide v_ij_erf_rk_cst_mu
provide v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu
! ao_non_hermit_term_chemist = non_h_ints
! return
BEGIN_DOC
! 1 1 2 2 1 2 1 2
!
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
!
END_DOC
implicit none
integer :: i, j, k, l, ipoint, m
double precision :: weight1, r(3)
double precision :: wall1, wall0
double precision, allocatable :: b_mat(:,:,:,:), ac_mat(:,:,:,:)
provide v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu
call wall_time(wall0)
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3),ac_mat(ao_num, ao_num, ao_num, ao_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,k,m,ipoint,r,weight1) &
!$OMP SHARED (aos_in_r_array_transp,aos_grad_in_r_array_transp_bis,b_mat)&
!$OMP SHARED (ao_num,n_points_final_grid,final_grid_points,final_weight_at_r_vector)
!$OMP DO SCHEDULE (static)
do m = 1, 3
do i = 1, ao_num
do k = 1, ao_num
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)
weight1 = final_weight_at_r_vector(ipoint)
b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * r(m) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m)
do m = 1, 3
do i = 1, ao_num
do k = 1, ao_num
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)
weight1 = final_weight_at_r_vector(ipoint)
b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * r(m) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! (A) b_mat(ipoint,k,i,m) X v_ij_erf_rk_cst_mu(j,l,r1)
! 1/2 \int dr1 x1 phi_k(1) d/dx1 phi_i(1) \int dr2 (1 - erf(mu_r12))/r12 phi_j(2) phi_l(2)
! (A) b_mat(ipoint,k,i,m) X v_ij_erf_rk_cst_mu(j,l,r1)
! 1/2 \int dr1 x1 phi_k(1) d/dx1 phi_i(1) \int dr2 (1 - erf(mu_r12))/r12 phi_j(2) phi_l(2)
ac_mat = 0.d0
do m = 1, 3
! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA
call dgemm("N","N",ao_num*ao_num,ao_num*ao_num,n_points_final_grid,1.d0,v_ij_erf_rk_cst_mu(1,1,1),ao_num*ao_num &
,b_mat(1,1,1,m),n_points_final_grid,1.d0,ac_mat,ao_num*ao_num)
! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, v_ij_erf_rk_cst_mu(1,1,1), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
, 1.d0, ac_mat, ao_num*ao_num)
enddo
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,k,m,ipoint,weight1) &
!$OMP SHARED (aos_in_r_array_transp,aos_grad_in_r_array_transp_bis,b_mat,ao_num,n_points_final_grid,final_weight_at_r_vector)
!$OMP DO SCHEDULE (static)
do m = 1, 3
do i = 1, ao_num
do k = 1, ao_num
do ipoint = 1, n_points_final_grid
weight1 = final_weight_at_r_vector(ipoint)
b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m)
do m = 1, 3
do i = 1, ao_num
do k = 1, ao_num
do ipoint = 1, n_points_final_grid
weight1 = final_weight_at_r_vector(ipoint)
b_mat(ipoint,k,i,m) = 0.5d0 * aos_in_r_array_transp(ipoint,k) * weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,m)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
@ -69,117 +77,141 @@ END_DOC
! 1/2 \int dr1 phi_k(1) d/dx1 phi_i(1) \int dr2 x2(1 - erf(mu_r12))/r12 phi_j(2) phi_l(2)
do m = 1, 3
! A B^T dim(A,1) dim(B,2) dim(A,2) alpha * A LDA
call dgemm("N","N",ao_num*ao_num,ao_num*ao_num,n_points_final_grid,-1.d0,x_v_ij_erf_rk_cst_mu(1,1,1,m),ao_num*ao_num &
,b_mat(1,1,1,m),n_points_final_grid,1.d0,ac_mat,ao_num*ao_num)
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, -1.d0 &
, x_v_ij_erf_rk_cst_mu(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
, 1.d0, ac_mat, ao_num*ao_num)
enddo
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,k,j,l) &
!$OMP SHARED (ac_mat,ao_non_hermit_term_chemist,ao_num)
!$OMP DO SCHEDULE (static)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! (ki|lj) (ki|lj) (lj|ki)
ao_non_hermit_term_chemist(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! (ki|lj) (ki|lj) (lj|ki)
ao_non_hermit_term_chemist(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
double precision :: wall1, wall0
call wall_time(wall1)
print*,'wall time dgemm ',wall1 - wall0
print *, ' wall time dgemm ', wall1 - wall0
END_PROVIDER
! ---
! TODO :: optimization :: transform into DGEM
BEGIN_PROVIDER [double precision, mo_non_hermit_term_chemist, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 1 1 2 2 1 2 1 2
!
! mo_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis
END_DOC
integer :: i,j,k,l,m,n,p,q
double precision, allocatable :: mo_tmp_1(:,:,:,:),mo_tmp_2(:,:,:,:),mo_tmp_3(:,:,:,:)
BEGIN_DOC
! 1 1 2 2 1 2 1 2
!
! mo_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis
END_DOC
implicit none
integer :: i, j, k, l, m, n, p, q
double precision, allocatable :: mo_tmp_1(:,:,:,:), mo_tmp_2(:,:,:,:)
allocate(mo_tmp_1(mo_num,ao_num,ao_num,ao_num))
! TODO :: optimization :: transform into DGEM
mo_tmp_1 = 0.d0
do m = 1, ao_num
do p = 1, ao_num
do n = 1, ao_num
do q = 1, ao_num
do k = 1, mo_num
! (k n|p m) = sum_q c_qk * (q n|p m)
mo_tmp_1(k,n,p,m) += mo_coef_transp(k,q) * ao_non_hermit_term_chemist(q,n,p,m)
enddo
allocate(mo_tmp_1(mo_num,ao_num,ao_num,ao_num))
mo_tmp_1 = 0.d0
do m = 1, ao_num
do p = 1, ao_num
do n = 1, ao_num
do q = 1, ao_num
do k = 1, mo_num
! (k n|p m) = sum_q c_qk * (q n|p m)
mo_tmp_1(k,n,p,m) += mo_coef_transp(k,q) * ao_non_hermit_term_chemist(q,n,p,m)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
free ao_non_hermit_term_chemist
allocate(mo_tmp_2(mo_num,mo_num,ao_num,ao_num))
mo_tmp_2 = 0.d0
do m = 1, ao_num
do p = 1, ao_num
do n = 1, ao_num
do i = 1, mo_num
do k = 1, mo_num
! (k i|p m) = sum_n c_ni * (k n|p m)
mo_tmp_2(k,i,p,m) += mo_coef_transp(i,n) * mo_tmp_1(k,n,p,m)
enddo
free ao_non_hermit_term_chemist
allocate(mo_tmp_2(mo_num,mo_num,ao_num,ao_num))
mo_tmp_2 = 0.d0
do m = 1, ao_num
do p = 1, ao_num
do n = 1, ao_num
do i = 1, mo_num
do k = 1, mo_num
! (k i|p m) = sum_n c_ni * (k n|p m)
mo_tmp_2(k,i,p,m) += mo_coef_transp(i,n) * mo_tmp_1(k,n,p,m)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(mo_tmp_1)
allocate(mo_tmp_1(mo_num,mo_num,mo_num,ao_num))
mo_tmp_1 = 0.d0
do m = 1, ao_num
do p = 1, ao_num
do l = 1, mo_num
do i = 1, mo_num
do k = 1, mo_num
mo_tmp_1(k,i,l,m) += mo_coef_transp(l,p) * mo_tmp_2(k,i,p,m)
enddo
deallocate(mo_tmp_1)
allocate(mo_tmp_1(mo_num,mo_num,mo_num,ao_num))
mo_tmp_1 = 0.d0
do m = 1, ao_num
do p = 1, ao_num
do l = 1, mo_num
do i = 1, mo_num
do k = 1, mo_num
mo_tmp_1(k,i,l,m) += mo_coef_transp(l,p) * mo_tmp_2(k,i,p,m)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(mo_tmp_2)
mo_non_hermit_term_chemist = 0.d0
do m = 1, ao_num
do j = 1, mo_num
do l = 1, mo_num
do i = 1, mo_num
do k = 1, mo_num
mo_non_hermit_term_chemist(k,i,l,j) += mo_coef_transp(j,m) * mo_tmp_1(k,i,l,m)
enddo
deallocate(mo_tmp_2)
mo_non_hermit_term_chemist = 0.d0
do m = 1, ao_num
do j = 1, mo_num
do l = 1, mo_num
do i = 1, mo_num
do k = 1, mo_num
mo_non_hermit_term_chemist(k,i,l,j) += mo_coef_transp(j,m) * mo_tmp_1(k,i,l,m)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(mo_tmp_1)
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, mo_non_hermit_term, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 1 2 1 2 1 2 1 2
!
! mo_non_hermit_term(k,l,i,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis
END_DOC
integer :: i,j,k,l
do j = 1, mo_num
do i = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
mo_non_hermit_term(k,l,i,j) = mo_non_hermit_term_chemist(k,i,l,j)
BEGIN_DOC
! 1 2 1 2 1 2 1 2
!
! mo_non_hermit_term(k,l,i,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the MO basis
END_DOC
implicit none
integer :: i, j, k, l
do j = 1, mo_num
do i = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
mo_non_hermit_term(k,l,i,j) = mo_non_hermit_term_chemist(k,i,l,j)
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
! ---

34
src/non_h_ints_mu/j12_nucl_utils.irp.f
View File

@ -586,4 +586,38 @@ end subroutine grad1_j12_mu_exc
! ---
subroutine grad1_jmu_modif_num(r1, r2, grad)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision, intent(out) :: grad(3)
double precision :: tmp0, tmp1, tmp2, tmp3, tmp4, grad_u12(3)
double precision, external :: j12_mu
double precision, external :: j1b_nucl
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
call grad1_j12_mu_exc(r1, r2, grad_u12)
tmp0 = j1b_nucl(r1)
tmp1 = j1b_nucl(r2)
tmp2 = j12_mu(r1, r2)
tmp3 = tmp0 * tmp1
tmp4 = tmp2 * tmp1
grad(1) = tmp3 * grad_u12(1) + tmp4 * grad_x_j1b_nucl(r1)
grad(2) = tmp3 * grad_u12(2) + tmp4 * grad_y_j1b_nucl(r1)
grad(3) = tmp3 * grad_u12(3) + tmp4 * grad_z_j1b_nucl(r1)
return
end subroutine grad1_jmu_modif_num
! ---

57
src/non_h_ints_mu/new_grad_tc.irp.f
View File

@ -17,10 +17,10 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin
! if J(r1,r2) = u12 x v1 x v2
!
! int2_grad1_u12_ao(:,i,j,ipoint) = v1 x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ]
! + \grad_1 v1 x [ \int dr2 u12 v2 \phi_i(r2) \phi_j(r2) ]
! - \grad_1 v1 x [ \int dr2 u12 v2 \phi_i(r2) \phi_j(r2) ]
! = 0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:)
! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:)
! + v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
!
!
END_DOC
@ -29,7 +29,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin
integer :: ipoint, i, j
double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
PROVIDE j1b_type j1b_pen
PROVIDE j1b_type
if(j1b_type .eq. 3) then
@ -46,12 +46,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin
do j = 1, ao_num
do i = 1, ao_num
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
int2_grad1_u12_ao(1,i,j,ipoint) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) + tmp_x * tmp2
int2_grad1_u12_ao(2,i,j,ipoint) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) + tmp_y * tmp2
int2_grad1_u12_ao(3,i,j,ipoint) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) + tmp_z * tmp2
int2_grad1_u12_ao(1,i,j,ipoint) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(1,i,j,ipoint) - tmp2 * tmp_x
int2_grad1_u12_ao(2,i,j,ipoint) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(2,i,j,ipoint) - tmp2 * tmp_y
int2_grad1_u12_ao(3,i,j,ipoint) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(3,i,j,ipoint) - tmp2 * tmp_z
enddo
enddo
enddo
@ -62,11 +62,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_poin
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
do j = 1, ao_num
do i = 1, ao_num
int2_grad1_u12_ao(1,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)
int2_grad1_u12_ao(2,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)
int2_grad1_u12_ao(3,i,j,ipoint) = v_ij_erf_rk_cst_mu(i,j,ipoint) * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)
tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
int2_grad1_u12_ao(1,i,j,ipoint) = tmp1 * x - x_v_ij_erf_rk_cst_mu_tmp(1,i,j,ipoint)
int2_grad1_u12_ao(2,i,j,ipoint) = tmp1 * y - x_v_ij_erf_rk_cst_mu_tmp(2,i,j,ipoint)
int2_grad1_u12_ao(3,i,j,ipoint) = tmp1 * z - x_v_ij_erf_rk_cst_mu_tmp(3,i,j,ipoint)
enddo
enddo
enddo
@ -93,9 +96,8 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
implicit none
integer :: ipoint, i, j, k, l
double precision :: contrib, weight1, contrib_x, contrib_y, contrib_z
double precision :: ao_k_r, ao_k_dx, ao_k_dy, ao_k_dz
double precision :: ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz
double precision :: weight1, contrib_x, contrib_y, contrib_z, tmp_x, tmp_y, tmp_z
double precision :: ao_k_r, ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz
double precision, allocatable :: ac_mat(:,:,:,:)
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num))
@ -105,27 +107,26 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
do i = 1, ao_num
ao_i_r = aos_in_r_array_transp (ipoint,i)
ao_i_dx = aos_grad_in_r_array_transp_bis(ipoint,i,1)
ao_i_dy = aos_grad_in_r_array_transp_bis(ipoint,i,2)
ao_i_dz = aos_grad_in_r_array_transp_bis(ipoint,i,3)
ao_i_r = weight1 * aos_in_r_array_transp (ipoint,i)
ao_i_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,1)
ao_i_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,2)
ao_i_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,3)
do k = 1, ao_num
ao_k_r = aos_in_r_array_transp (ipoint,k)
ao_k_dx = aos_grad_in_r_array_transp_bis(ipoint,k,1)
ao_k_dy = aos_grad_in_r_array_transp_bis(ipoint,k,2)
ao_k_dz = aos_grad_in_r_array_transp_bis(ipoint,k,3)
ao_k_r = aos_in_r_array_transp(ipoint,k)
tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)
tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)
tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)
do j = 1, ao_num
do l = 1, ao_num
contrib_x = int2_grad1_u12_ao(1,l,j,ipoint) * ( ao_k_r * ao_i_dx - ao_i_r * ao_k_dx )
contrib_y = int2_grad1_u12_ao(2,l,j,ipoint) * ( ao_k_r * ao_i_dy - ao_i_r * ao_k_dy )
contrib_z = int2_grad1_u12_ao(3,l,j,ipoint) * ( ao_k_r * ao_i_dz - ao_i_r * ao_k_dz )
contrib_x = int2_grad1_u12_ao(1,l,j,ipoint) * tmp_x
contrib_y = int2_grad1_u12_ao(2,l,j,ipoint) * tmp_y
contrib_z = int2_grad1_u12_ao(3,l,j,ipoint) * tmp_z
contrib = weight1 * ( contrib_x + contrib_y + contrib_z )
ac_mat(k,i,l,j) += contrib
ac_mat(k,i,l,j) += contrib_x + contrib_y + contrib_z
enddo
enddo
enddo

341
src/non_h_ints_mu/numerical_integ.irp.f
View File

@ -55,6 +55,7 @@ double precision function num_int2_u2_j1b2(i, j, ipoint)
double precision, external :: ao_value
double precision, external :: j1b_nucl
double precision, external :: j12_mu
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
@ -74,13 +75,14 @@ double precision function num_int2_u2_j1b2(i, j, ipoint)
tmp1 = j1b_nucl(r2)
tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
tmp3 = 0.d0
do i_fit = 1, n_max_fit_slat
expo = expo_gauss_j_mu_x_2(i_fit)
coef = coef_gauss_j_mu_x_2(i_fit)
tmp3 += coef * dexp(-expo*x2)
enddo
!tmp3 = 0.d0
!do i_fit = 1, n_max_fit_slat
! expo = expo_gauss_j_mu_x_2(i_fit)
! coef = coef_gauss_j_mu_x_2(i_fit)
! tmp3 += coef * dexp(-expo*x2)
!enddo
tmp3 = j12_mu(r1, r2)
tmp3 = tmp3 * tmp3
num_int2_u2_j1b2 += tmp2 * tmp3
enddo
@ -127,13 +129,15 @@ double precision function num_int2_grad1u2_grad2u2_j1b2(i, j, ipoint)
tmp1 = j1b_nucl(r2)
tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
tmp3 = 0.d0
do i_fit = 1, n_max_fit_slat
expo = expo_gauss_1_erf_x_2(i_fit)
coef = coef_gauss_1_erf_x_2(i_fit)
!tmp3 = 0.d0
!do i_fit = 1, n_max_fit_slat
! expo = expo_gauss_1_erf_x_2(i_fit)
! coef = coef_gauss_1_erf_x_2(i_fit)
! tmp3 += coef * dexp(-expo*x2)
!enddo
tmp3 = derf(mu_erf*r12) - 1.d0
tmp3 = tmp3 * tmp3
tmp3 += coef * dexp(-expo*x2)
enddo
tmp3 = -0.25d0 * tmp3
num_int2_grad1u2_grad2u2_j1b2 += tmp2 * tmp3
@ -246,6 +250,12 @@ end subroutine num_x_v_ij_erf_rk_cst_mu_j1b
subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ)
BEGIN_DOC
!
! \int dr2 [-grad_1 u12] \phi_i(r2) \phi_j(r2) x v12_1b(r1, r2)
!
END_DOC
implicit none
integer, intent(in) :: i, j, ipoint
@ -256,7 +266,6 @@ subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ)
double precision :: tmp_x, tmp_y, tmp_z
double precision, external :: ao_value
double precision, external :: j12_nucl
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
@ -269,9 +278,9 @@ subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ)
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint)
tmp = ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
call grad1_j12_mu_exc(r1, r2, grad)
call grad1_jmu_modif_num(r1, r2, grad)
tmp_x += tmp * (-1.d0 * grad(1))
tmp_y += tmp * (-1.d0 * grad(2))
@ -289,16 +298,33 @@ end subroutine num_int2_grad1_u12_ao
double precision function num_gradu_squared_u_ij_mu(i, j, ipoint)
BEGIN_DOC
!
! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2
! [ v1^2 ((grad_1 u12)^2 + (grad_2 u12^2)])
! + u12^2 (grad_1 v1)^2
! + 2 u12 v1 (grad_1 u12) . (grad_1 v1)
!
END_DOC
implicit none
integer, intent(in) :: i, j, ipoint
integer :: jpoint
double precision :: tmp, r1(3), r2(3), r12
double precision :: tmp_x, tmp_y, tmp_z, tmp1, tmp2
double precision :: r1(3), r2(3)
double precision :: tmp_x, tmp_y, tmp_z, r12
double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3)
double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp
double precision :: fst_term, scd_term, thd_term, tmp
double precision, external :: ao_value
double precision, external :: j12_nucl
double precision, external :: j1b_nucl
double precision, external :: j12_mu
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
@ -306,16 +332,32 @@ double precision function num_gradu_squared_u_ij_mu(i, j, ipoint)
num_gradu_squared_u_ij_mu = 0.d0
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp_x = r1(1) - r2(1)
tmp_y = r1(2) - r2(2)
tmp_z = r1(3) - r2(3)
r12 = dsqrt( tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z )
tmp1 = 1.d0 - derf(mu_erf * r12)
tmp2 = j12_nucl(r1, r2)
tmp = -0.25d0 * tmp1 * tmp1 * tmp2 * tmp2 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
dx1_v1 = grad_x_j1b_nucl(r1)
dy1_v1 = grad_y_j1b_nucl(r1)
dz1_v1 = grad_z_j1b_nucl(r1)
call grad1_j12_mu_exc(r1, r2, grad_u12)
tmp1 = 1.d0 - derf(mu_erf * r12)
v1_tmp = j1b_nucl(r1)
v2_tmp = j1b_nucl(r2)
u12_tmp = j12_mu(r1, r2)
fst_term = 0.5d0 * tmp1 * tmp1 * v1_tmp * v1_tmp
scd_term = u12_tmp * u12_tmp * (dx1_v1*dx1_v1 + dy1_v1*dy1_v1 + dz1_v1*dz1_v1)
thd_term = 2.d0 * v1_tmp * u12_tmp * (dx1_v1*grad_u12(1) + dy1_v1*grad_u12(2) + dz1_v1*grad_u12(3))
tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * (fst_term + scd_term + thd_term) * v2_tmp * v2_tmp
num_gradu_squared_u_ij_mu += tmp
enddo
@ -325,4 +367,257 @@ end function num_gradu_squared_u_ij_mu
! ---
double precision function num_grad12_j12(i, j, ipoint)
BEGIN_DOC
!
! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [v1^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) ]
!
END_DOC
implicit none
integer, intent(in) :: i, j, ipoint
integer :: jpoint
double precision :: r1(3), r2(3)
double precision :: tmp_x, tmp_y, tmp_z, r12
double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3)
double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp
double precision :: fst_term, scd_term, thd_term, tmp
double precision, external :: ao_value
double precision, external :: j1b_nucl
double precision, external :: j12_mu
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
num_grad12_j12 = 0.d0
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp_x = r1(1) - r2(1)
tmp_y = r1(2) - r2(2)
tmp_z = r1(3) - r2(3)
r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
dx1_v1 = grad_x_j1b_nucl(r1)
dy1_v1 = grad_y_j1b_nucl(r1)
dz1_v1 = grad_z_j1b_nucl(r1)
call grad1_j12_mu_exc(r1, r2, grad_u12)
tmp1 = 1.d0 - derf(mu_erf * r12)
v1_tmp = j1b_nucl(r1)
v2_tmp = j1b_nucl(r2)
u12_tmp = j12_mu(r1, r2)
fst_term = 0.5d0 * tmp1 * tmp1 * v1_tmp * v1_tmp
tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * fst_term * v2_tmp * v2_tmp
num_grad12_j12 += tmp
enddo
return
end function num_grad12_j12
! ---
double precision function num_u12sq_j1bsq(i, j, ipoint)
BEGIN_DOC
!
! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [ u12^2 (grad_1 v1)^2 ]
!
END_DOC
implicit none
integer, intent(in) :: i, j, ipoint
integer :: jpoint
double precision :: r1(3), r2(3)
double precision :: tmp_x, tmp_y, tmp_z, r12
double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3)
double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp
double precision :: fst_term, scd_term, thd_term, tmp
double precision, external :: ao_value
double precision, external :: j1b_nucl
double precision, external :: j12_mu
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
num_u12sq_j1bsq = 0.d0
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp_x = r1(1) - r2(1)
tmp_y = r1(2) - r2(2)
tmp_z = r1(3) - r2(3)
r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
dx1_v1 = grad_x_j1b_nucl(r1)
dy1_v1 = grad_y_j1b_nucl(r1)
dz1_v1 = grad_z_j1b_nucl(r1)
call grad1_j12_mu_exc(r1, r2, grad_u12)
tmp1 = 1.d0 - derf(mu_erf * r12)
v1_tmp = j1b_nucl(r1)
v2_tmp = j1b_nucl(r2)
u12_tmp = j12_mu(r1, r2)
scd_term = u12_tmp * u12_tmp * (dx1_v1*dx1_v1 + dy1_v1*dy1_v1 + dz1_v1*dz1_v1)
tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * scd_term * v2_tmp * v2_tmp
num_u12sq_j1bsq += tmp
enddo
return
end function num_u12sq_j1bsq
! ---
double precision function num_u12_grad1_u12_j1b_grad1_j1b(i, j, ipoint)
BEGIN_DOC
!
! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [ 2 u12 v1 (grad_1 u12) . (grad_1 v1) ]
!
END_DOC
implicit none
integer, intent(in) :: i, j, ipoint
integer :: jpoint
double precision :: r1(3), r2(3)
double precision :: tmp_x, tmp_y, tmp_z, r12
double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3)
double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp
double precision :: fst_term, scd_term, thd_term, tmp
double precision, external :: ao_value
double precision, external :: j1b_nucl
double precision, external :: j12_mu
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
num_u12_grad1_u12_j1b_grad1_j1b = 0.d0
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp_x = r1(1) - r2(1)
tmp_y = r1(2) - r2(2)
tmp_z = r1(3) - r2(3)
r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
dx1_v1 = grad_x_j1b_nucl(r1)
dy1_v1 = grad_y_j1b_nucl(r1)
dz1_v1 = grad_z_j1b_nucl(r1)
call grad1_j12_mu_exc(r1, r2, grad_u12)
tmp1 = 1.d0 - derf(mu_erf * r12)
v1_tmp = j1b_nucl(r1)
v2_tmp = j1b_nucl(r2)
u12_tmp = j12_mu(r1, r2)
thd_term = 2.d0 * v1_tmp * u12_tmp * (dx1_v1*grad_u12(1) + dy1_v1*grad_u12(2) + dz1_v1*grad_u12(3))
tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * thd_term * v2_tmp * v2_tmp
num_u12_grad1_u12_j1b_grad1_j1b += tmp
enddo
return
end function num_u12_grad1_u12_j1b_grad1_j1b
! ---
subroutine num_int2_u_grad1u_total_j1b2(i, j, ipoint, integ)
BEGIN_DOC
!
! \int dr2 u12 (grad_1 u12) \phi_i(r2) \phi_j(r2) x v_1b(r2)^2
!
END_DOC
implicit none
integer, intent(in) :: i, j, ipoint
double precision, intent(out) :: integ(3)
integer :: jpoint
double precision :: r1(3), r2(3), grad(3)
double precision :: dx, dy, dz, r12, tmp0, tmp1, tmp2
double precision :: tmp_x, tmp_y, tmp_z
double precision, external :: ao_value
double precision, external :: j1b_nucl
double precision, external :: j12_mu
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
tmp_x = 0.d0
tmp_y = 0.d0
tmp_z = 0.d0
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
dx = r1(1) - r2(1)
dy = r1(2) - r2(2)
dz = r1(3) - r2(3)
r12 = dsqrt( dx * dx + dy * dy + dz * dz )
if(r12 .lt. 1d-10) cycle
tmp0 = j1b_nucl(r2)
tmp1 = 0.5d0 * j12_mu(r1, r2) * (1.d0 - derf(mu_erf * r12)) / r12
tmp2 = tmp0 * tmp0 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
tmp_x += tmp2 * dx
tmp_y += tmp2 * dy
tmp_z += tmp2 * dz
enddo
integ(1) = tmp_x
integ(2) = tmp_y
integ(3) = tmp_z
return
end subroutine num_int2_u_grad1u_total_j1b2
! ---

60
src/non_h_ints_mu/total_tc_int.irp.f Normal file
View File

@ -0,0 +1,60 @@
! ---
BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao_num)]
implicit none
integer :: i, j, k, l
double precision :: wall1, wall0
call wall_time(wall0)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist ', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num) ]
BEGIN_DOC
!
! ao_two_e_coul(k,i,l,j) = ( k i | 1/r12 | l j ) = < l k | 1/r12 | j i >
!
END_DOC
integer :: i, j, k, l
double precision :: integral
double precision, external :: get_ao_two_e_integral
PROVIDE ao_integrals_map
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! < 1:k, 2:l | 1:i, 2:j >
integral = get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
ao_two_e_coul(k,i,l,j) = integral
enddo
enddo
enddo
enddo
END_PROVIDER
! ---

1
src/tc_bi_ortho/compute_deltamu_right.irp.f
View File

@ -34,6 +34,7 @@ subroutine delta_right()
!do k = 1, 1
! get < I_left | H_mu - H | psi_right >
!call get_h_bitc_right(psi_det, psi_r_coef_bi_ortho(:,k), N_det, N_int, delta(:,k))
call get_delta_bitc_right(psi_det, psi_r_coef_bi_ortho(:,k), N_det, N_int, delta(:,k))
! order as QMCCHEM

105
src/tc_bi_ortho/dressing_vectors_lr.irp.f
View File

@ -23,10 +23,12 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
double precision :: htc_mono, htc_twoe, htc_three, htc_tot
double precision :: delta_mat
print *, ' get_delta_bitc_right ...'
i = 1
j = 1
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono , h_twoe , h_tot)
call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
delta = 0.d0
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
@ -39,7 +41,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
! < I | Htilde | J >
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
! < I | H | J >
call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono , h_twoe , h_tot)
call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
delta_mat = htc_tot - h_tot
@ -52,3 +54,102 @@ end subroutine get_delta_bitc_right
! ---
subroutine get_htc_bitc_right(psidet, psicoef, ndet, Nint, delta)
BEGIN_DOC
!
! delta(I) = < I_left | H_TC | Psi_right >
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: ndet, Nint
double precision, intent(in) :: psicoef(ndet)
integer(bit_kind), intent(in) :: psidet(Nint,2,ndet)
double precision, intent(out) :: delta(ndet)
integer :: i, j
double precision :: htc_mono, htc_twoe, htc_three, htc_tot
print *, ' get_htc_bitc_right ...'
i = 1
j = 1
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
delta = 0.d0
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
!$OMP SHARED(delta, ndet, psidet, psicoef, Nint) &
!$OMP PRIVATE(i, j, htc_mono, htc_twoe, htc_three, htc_tot)
do i = 1, ndet
do j = 1, ndet
! < I | Htilde | J >
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
delta(i) = delta(i) + psicoef(j) * htc_tot
enddo
enddo
!$OMP END PARALLEL DO
end subroutine get_htc_bitc_right
! ---
subroutine get_h_bitc_right(psidet, psicoef, ndet, Nint, delta)
BEGIN_DOC
!
! delta(I) = < I_left | H | Psi_right >
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: ndet, Nint
double precision, intent(in) :: psicoef(ndet)
integer(bit_kind), intent(in) :: psidet(Nint,2,ndet)
double precision, intent(out) :: delta(ndet)
integer :: i, j
double precision :: h_mono, h_twoe, h_tot
print *, ' get_h_bitc_right ...'
i = 1
j = 1
call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
!double precision :: norm
!norm = 0.d0
!do i = 1, ndet
! norm += psicoef(i) * psicoef(i)
!enddo
!print*, ' norm = ', norm
call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
delta = 0.d0
! !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
! !$OMP SHARED(delta, ndet, psidet, psicoef, Nint) &
! !$OMP PRIVATE(i, j, h_mono, h_twoe, h_tot)
do i = 1, ndet
do j = 1, ndet
! < I | H | J >
call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
delta(i) = delta(i) + psicoef(j) * h_tot
enddo
enddo
! !$OMP END PARALLEL DO
end subroutine get_h_bitc_right
! ---

32
src/tc_bi_ortho/h_biortho.irp.f
View File

@ -5,7 +5,7 @@ subroutine hmat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, htot)
BEGIN_DOC
!
! <key_j | H | key_i > where | key_j > is developed on the LEFT basis and | key_i > is developed on the RIGHT basis
! < key_j | H | key_i > where | key_j > is developed on the LEFT basis and | key_i > is developed on the RIGHT basis
!
END_DOC
@ -13,11 +13,11 @@ subroutine hmat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, htot)
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, htot
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, htot
integer :: degree
integer :: degree
hmono = 0.d0
htwoe = 0.d0
@ -31,11 +31,11 @@ subroutine hmat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, htot)
call diag_hmat_bi_ortho(Nint, key_i, hmono, htwoe)
htot = htot + nuclear_repulsion
else if (degree == 1)then
else if (degree == 1) then
call single_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe)
else if(degree == 2)then
else if(degree == 2) then
call double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe)
@ -59,8 +59,7 @@ subroutine diag_hmat_bi_ortho(Nint, key_i, hmono, htwoe)
double precision, intent(out) :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer(bit_kind) :: key_i_core(Nint,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin
hmono = 0.d0
htwoe = 0.d0
@ -125,13 +124,11 @@ subroutine single_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe)
double precision, intent(out) :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: Ne(2), i, j, ii, ispin, jspin
integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
integer :: other_spin(2)
integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2)
double precision :: phase
double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13
other_spin(1) = 2
other_spin(2) = 1
@ -201,11 +198,10 @@ subroutine double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe)
double precision, intent(out) :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: Ne(2), i, j, ii, ispin, jspin
integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
integer :: other_spin(2)
integer(bit_kind) :: key_i_core(Nint,2)
double precision :: phase
other_spin(1) = 2
@ -225,7 +221,7 @@ subroutine double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe)
call get_double_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
if(s1.ne.s2) then
if(s1 .ne. s2) then
htwoe = mo_bi_ortho_coul_e(p2,p1,h2,h1)
@ -233,10 +229,8 @@ subroutine double_hmat_bi_ortho(Nint, key_j, key_i, hmono, htwoe)
! same spin two-body
! direct terms
htwoe = mo_bi_ortho_coul_e(p2,p1,h2,h1)
! exchange terms
htwoe -= mo_bi_ortho_coul_e(p1,p2,h2,h1)
! direct terms exchange terms
htwoe = mo_bi_ortho_coul_e(p2,p1,h2,h1) - mo_bi_ortho_coul_e(p1,p2,h2,h1)
endif

18
src/tc_bi_ortho/psi_r_l_prov.irp.f
View File

@ -41,6 +41,15 @@ BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states)
enddo
deallocate(psi_l_coef_bi_ortho_read)
else
print*, 'psi_l_coef_bi_ortho are psi_coef'
do k=1,N_states
do i=1,N_det
psi_l_coef_bi_ortho(i,k) = psi_coef(i,k)
enddo
enddo
endif
endif
endif
@ -100,6 +109,15 @@ BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states)
enddo
deallocate(psi_r_coef_bi_ortho_read)
else
print*, 'psi_r_coef_bi_ortho are psi_coef'
do k=1,N_states
do i=1,N_det
psi_r_coef_bi_ortho(i,k) = psi_coef(i,k)
enddo
enddo
endif
endif
endif

80
src/tc_bi_ortho/slater_tc.irp.f
View File

@ -1,4 +1,6 @@
!!!!!!
! ---
subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot)
BEGIN_DOC
@ -15,13 +17,14 @@ subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot)
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
double precision, intent(out) :: htot
double precision :: hmono,htwoe,hthree
double precision :: hmono, htwoe, hthree
integer :: degree
call get_excitation_degree(key_j, key_i, degree, Nint)
if(degree.gt.2)then
htot = 0.d0
htot = 0.d0
else
call htilde_mu_mat_bi_ortho(key_j,key_i, Nint, hmono,htwoe,hthree,htot)
call htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
endif
end subroutine htilde_mu_mat_bi_ortho_tot
@ -29,52 +32,63 @@ end subroutine htilde_mu_mat_bi_ortho_tot
! --
subroutine htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
implicit none
use bitmasks
BEGIN_DOC
!
! <key_j | H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
!!
! Returns the detail of the matrix element in terms of single, two and three electron contribution.
!! WARNING !!
!
! Non hermitian !!
!
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2),key_j(Nint,2)
double precision, intent(out) :: hmono,htwoe,hthree,htot
integer :: degree
hmono = 0.d0
htwoe= 0.d0
htot = 0.d0
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: degree
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
hthree = 0.D0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.gt.2)return
if(degree.gt.2) return
if(degree == 0)then
call diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
call diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
else if (degree == 1)then
call single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
call single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
else if(degree == 2)then
call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
endif
if(three_body_h_tc) then
if(degree == 2) then
if(.not.double_normal_ord) then
call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
endif
else if(degree == 1)then
call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
else if(degree == 0)then
call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
endif
endif
htot = hmono + htwoe + hthree
if(degree==0)then
htot += nuclear_repulsion
endif
call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
endif
if(three_body_h_tc) then
if(degree == 2) then
if(.not.double_normal_ord) then
call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
endif
else if(degree == 1) then
call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
else if(degree == 0) then
call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
endif
endif
htot = hmono + htwoe + hthree
if(degree==0) then
htot += nuclear_repulsion
endif
end
! ---
subroutine diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
BEGIN_DOC

7
src/tc_bi_ortho/slater_tc_3e.irp.f
View File

@ -207,6 +207,8 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
end
! ---
subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
BEGIN_DOC
@ -244,7 +246,7 @@ subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
return
endif
if(core_tc_op)then
if(core_tc_op) then
do i = 1, Nint
key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
@ -291,3 +293,6 @@ subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
endif
hthree *= phase
end
! ---