9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-09 06:53:38 +01:00

cycle in integ added

This commit is contained in:
AbdAmmar 2022-11-21 10:38:42 +01:00
parent d70fb23c9a
commit eb37b83763
10 changed files with 876 additions and 61 deletions

View File

@ -275,6 +275,393 @@ end subroutine NAI_pol_x_mult_erf_ao
! ---
double precision function NAI_pol_x_mult_erf_ao_x(i_ao, j_ao, mu_in, C_center)
BEGIN_DOC
!
! Computes the following integral :
!
! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: mu_in, C_center(3)
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, power_xA(3)
double precision :: A_center(3), B_center(3), integral, alpha, beta, coef
double precision :: NAI_pol_mult_erf
NAI_pol_x_mult_erf_ao_x = 0.d0
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) return
num_A = ao_nucl(i_ao)
power_A(1:3) = ao_power(i_ao,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j_ao)
power_B(1:3) = ao_power(j_ao,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
power_xA = power_A
power_xA(1) += 1
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alpha = ao_expo_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
beta = ao_expo_ordered_transp(j,j_ao)
coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
! First term = (x-Ax)**(ax+1)
integral = NAI_pol_mult_erf(A_center, B_center, power_xA, power_B, alpha, beta, C_center, n_pt_in, mu_in)
NAI_pol_x_mult_erf_ao_x += integral * coef
! Second term = Ax * (x-Ax)**(ax)
integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
NAI_pol_x_mult_erf_ao_x += A_center(1) * integral * coef
enddo
enddo
end function NAI_pol_x_mult_erf_ao_x
! ---
double precision function NAI_pol_x_mult_erf_ao_y(i_ao, j_ao, mu_in, C_center)
BEGIN_DOC
!
! Computes the following integral :
!
! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: mu_in, C_center(3)
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, power_xA(3)
double precision :: A_center(3), B_center(3), integral, alpha, beta, coef
double precision :: NAI_pol_mult_erf
NAI_pol_x_mult_erf_ao_y = 0.d0
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) return
num_A = ao_nucl(i_ao)
power_A(1:3) = ao_power(i_ao,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j_ao)
power_B(1:3) = ao_power(j_ao,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
power_xA = power_A
power_xA(2) += 1
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alpha = ao_expo_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
beta = ao_expo_ordered_transp(j,j_ao)
coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
! First term = (x-Ax)**(ax+1)
integral = NAI_pol_mult_erf(A_center, B_center, power_xA, power_B, alpha, beta, C_center, n_pt_in, mu_in)
NAI_pol_x_mult_erf_ao_y += integral * coef
! Second term = Ax * (x-Ax)**(ax)
integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
NAI_pol_x_mult_erf_ao_y += A_center(2) * integral * coef
enddo
enddo
end function NAI_pol_x_mult_erf_ao_y
! ---
double precision function NAI_pol_x_mult_erf_ao_z(i_ao, j_ao, mu_in, C_center)
BEGIN_DOC
!
! Computes the following integral :
!
! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: mu_in, C_center(3)
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, power_xA(3)
double precision :: A_center(3), B_center(3), integral, alpha, beta, coef
double precision :: NAI_pol_mult_erf
NAI_pol_x_mult_erf_ao_z = 0.d0
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) return
num_A = ao_nucl(i_ao)
power_A(1:3) = ao_power(i_ao,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j_ao)
power_B(1:3) = ao_power(j_ao,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
power_xA = power_A
power_xA(3) += 1
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alpha = ao_expo_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
beta = ao_expo_ordered_transp(j,j_ao)
coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
! First term = (x-Ax)**(ax+1)
integral = NAI_pol_mult_erf(A_center, B_center, power_xA, power_B, alpha, beta, C_center, n_pt_in, mu_in)
NAI_pol_x_mult_erf_ao_z += integral * coef
! Second term = Ax * (x-Ax)**(ax)
integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
NAI_pol_x_mult_erf_ao_z += A_center(3) * integral * coef
enddo
enddo
end function NAI_pol_x_mult_erf_ao_z
! ---
double precision function NAI_pol_x_mult_erf_ao_with1s_x(i_ao, j_ao, beta, B_center, mu_in, C_center)
BEGIN_DOC
!
! Computes the following integral :
!
! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3)
integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, power_xA(3)
double precision :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef, coefi
double precision, external :: NAI_pol_mult_erf_with1s
double precision, external :: NAI_pol_x_mult_erf_ao_x
ASSERT(beta .ge. 0.d0)
if(beta .lt. 1d-10) then
NAI_pol_x_mult_erf_ao_with1s_x = NAI_pol_x_mult_erf_ao_x(i_ao, j_ao, mu_in, C_center)
return
endif
NAI_pol_x_mult_erf_ao_with1s_x = 0.d0
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) then
return
endif
power_Ai(1:3) = ao_power(i_ao,1:3)
power_Aj(1:3) = ao_power(j_ao,1:3)
Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
power_xA = power_Ai
power_xA(1) += 1
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alphai = ao_expo_ordered_transp (i,i_ao)
coefi = ao_coef_normalized_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
alphaj = ao_expo_ordered_transp (j,j_ao)
coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao)
! 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 )
NAI_pol_x_mult_erf_ao_with1s_x += 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 )
NAI_pol_x_mult_erf_ao_with1s_x += Ai_center(1) * integral * coef
enddo
enddo
end function NAI_pol_x_mult_erf_ao_with1s_x
! ---
double precision function NAI_pol_x_mult_erf_ao_with1s_y(i_ao, j_ao, beta, B_center, mu_in, C_center)
BEGIN_DOC
!
! Computes the following integral :
!
! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3)
integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, power_xA(3)
double precision :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef, coefi
double precision, external :: NAI_pol_mult_erf_with1s
double precision, external :: NAI_pol_x_mult_erf_ao_y
ASSERT(beta .ge. 0.d0)
if(beta .lt. 1d-10) then
NAI_pol_x_mult_erf_ao_with1s_y = NAI_pol_x_mult_erf_ao_y(i_ao, j_ao, mu_in, C_center)
return
endif
NAI_pol_x_mult_erf_ao_with1s_y = 0.d0
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) then
return
endif
power_Ai(1:3) = ao_power(i_ao,1:3)
power_Aj(1:3) = ao_power(j_ao,1:3)
Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
power_xA = power_Ai
power_xA(2) += 1
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alphai = ao_expo_ordered_transp (i,i_ao)
coefi = ao_coef_normalized_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
alphaj = ao_expo_ordered_transp (j,j_ao)
coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao)
! 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 )
NAI_pol_x_mult_erf_ao_with1s_y += 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 )
NAI_pol_x_mult_erf_ao_with1s_y += Ai_center(2) * integral * coef
enddo
enddo
end function NAI_pol_x_mult_erf_ao_with1s_y
! ---
double precision function NAI_pol_x_mult_erf_ao_with1s_z(i_ao, j_ao, beta, B_center, mu_in, C_center)
BEGIN_DOC
!
! Computes the following integral :
!
! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3)
integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, power_xA(3)
double precision :: Ai_center(3), Aj_center(3), integral, alphai, alphaj, coef, coefi
double precision, external :: NAI_pol_mult_erf_with1s
double precision, external :: NAI_pol_x_mult_erf_ao_z
ASSERT(beta .ge. 0.d0)
if(beta .lt. 1d-10) then
NAI_pol_x_mult_erf_ao_with1s_z = NAI_pol_x_mult_erf_ao_z(i_ao, j_ao, mu_in, C_center)
return
endif
NAI_pol_x_mult_erf_ao_with1s_z = 0.d0
if(ao_overlap_abs(j_ao,i_ao) .lt. 1.d-12) then
return
endif
power_Ai(1:3) = ao_power(i_ao,1:3)
power_Aj(1:3) = ao_power(j_ao,1:3)
Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
power_xA = power_Ai
power_xA(3) += 1
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alphai = ao_expo_ordered_transp (i,i_ao)
coefi = ao_coef_normalized_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
alphaj = ao_expo_ordered_transp (j,j_ao)
coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao)
! 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 )
NAI_pol_x_mult_erf_ao_with1s_z += 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 )
NAI_pol_x_mult_erf_ao_with1s_z += Ai_center(3) * integral * coef
enddo
enddo
end function NAI_pol_x_mult_erf_ao_with1s_z
! ---
subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints)
BEGIN_DOC

View File

@ -32,7 +32,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2)
!$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)
@ -42,22 +41,41 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
do j = i, ao_num
tmp = 0.d0
do i_1s = 1, List_all_comb_b3_size
do i_fit = 1, n_max_fit_slat
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = coef_gauss_1_erf_x_2(i_fit)
do i_fit = 1, n_max_fit_slat
! ---
expo_fit = expo_gauss_1_erf_x_2(i_fit)
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)
coef = List_all_comb_b3_coef (1)
beta = List_all_comb_b3_expo (1)
B_center(1) = List_all_comb_b3_cent(1,1)
B_center(2) = List_all_comb_b3_cent(2,1)
B_center(3) = List_all_comb_b3_cent(3,1)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
if(dabs(int_fit) .lt. 1d-10) cycle
tmp += -0.25d0 * coef * coef_fit * int_fit
! ---
do i_1s = 2, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += -0.25d0 * coef * coef_fit * int_fit
enddo
! ---
enddo
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = tmp
@ -112,7 +130,6 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u2_j1b2)
!$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)
@ -122,22 +139,41 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final
do j = i, ao_num
tmp = 0.d0
do i_1s = 1, List_all_comb_b3_size
do i_fit = 1, n_max_fit_slat
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
expo_fit = expo_gauss_j_mu_x_2(i_fit)
coef_fit = coef_gauss_j_mu_x_2(i_fit)
do i_fit = 1, n_max_fit_slat
! ---
expo_fit = expo_gauss_j_mu_x_2(i_fit)
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)
coef = List_all_comb_b3_coef (1)
beta = List_all_comb_b3_expo (1)
B_center(1) = List_all_comb_b3_cent(1,1)
B_center(2) = List_all_comb_b3_cent(2,1)
B_center(3) = List_all_comb_b3_cent(3,1)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
if(dabs(int_fit) .lt. 1d-10) cycle
tmp += coef * coef_fit * int_fit
! ---
do i_1s = 2, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit
enddo
! ---
enddo
int2_u2_j1b2(j,i,ipoint) = tmp
@ -205,21 +241,52 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_p
tmp_x = 0.d0
tmp_y = 0.d0
tmp_z = 0.d0
do i_1s = 1, List_all_comb_b3_size
do i_fit = 1, n_max_fit_slat
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
! ---
coef = List_all_comb_b3_coef (1)
beta = List_all_comb_b3_expo (1)
B_center(1) = List_all_comb_b3_cent(1,1)
B_center(2) = List_all_comb_b3_cent(2,1)
B_center(3) = List_all_comb_b3_cent(3,1)
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
alpha_1s = beta + expo_fit
alpha_1s_inv = 1.d0 / alpha_1s
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
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 * dist
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)
if( (dabs(int_fit(1)) + dabs(int_fit(2)) + dabs(int_fit(3))) .lt. 3d-10 ) cycle
tmp_x += coef_tmp * int_fit(1)
tmp_y += coef_tmp * int_fit(2)
tmp_z += coef_tmp * int_fit(3)
! ---
do i_1s = 2, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
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))
alpha_1s = beta + expo_fit
alpha_1s_inv = 1.d0 / alpha_1s
@ -229,9 +296,8 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_p
centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
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
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)
@ -239,6 +305,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (3, ao_num, ao_num, n_p
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
@ -306,21 +375,49 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
r(3) = final_grid_points(3,ipoint)
tmp = 0.d0
do i_1s = 1, List_all_comb_b3_size
do i_fit = 1, n_max_fit_slat
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
! ---
coef = List_all_comb_b3_coef (1)
beta = List_all_comb_b3_expo (1)
B_center(1) = List_all_comb_b3_cent(1,1)
B_center(2) = List_all_comb_b3_cent(2,1)
B_center(3) = List_all_comb_b3_cent(3,1)
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
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_fit = expo_gauss_j_mu_1_erf(i_fit)
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
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)
if(dabs(int_fit) .lt. 1d-10) cycle
tmp += coef_tmp * int_fit
! ---
do i_1s = 2, List_all_comb_b3_size
coef = List_all_comb_b3_coef (i_1s)
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
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))
alpha_1s = beta + expo_fit
alpha_1s_inv = 1.d0 / alpha_1s
@ -329,14 +426,16 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
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
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 += coef_tmp * int_fit
enddo
! ---
enddo
int2_u_grad1u_j1b2(j,i,ipoint) = tmp

View File

@ -38,7 +38,24 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
do j = i, ao_num
tmp = 0.d0
do i_1s = 1, List_all_comb_b2_size
! ---
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
int_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)
if(dabs(int_mu - int_coulomb) .lt. 1d-10) cycle
tmp += coef * (int_mu - int_coulomb)
! ---
do i_1s = 2, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
beta = List_all_comb_b2_expo (i_1s)
@ -52,6 +69,8 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
tmp += coef * (int_mu - int_coulomb)
enddo
! ---
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = tmp
enddo
enddo
@ -138,7 +157,27 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_
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 (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
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)
if( (dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle
tmp_x += coef * (ints(1) - ints_coulomb(1))
tmp_y += coef * (ints(2) - ints_coulomb(2))
tmp_z += coef * (ints(3) - ints_coulomb(3))
! ---
do i_1s = 2, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
beta = List_all_comb_b2_expo (i_1s)
@ -154,6 +193,8 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_
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
@ -222,22 +263,41 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
do j = i, ao_num
tmp = 0.d0
do i_1s = 1, List_all_comb_b2_size
do i_fit = 1, n_max_fit_slat
coef = List_all_comb_b2_coef (i_1s)
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
expo_fit = expo_gauss_j_mu_x(i_fit)
coef_fit = coef_gauss_j_mu_x(i_fit)
do i_fit = 1, n_max_fit_slat
! ---
expo_fit = expo_gauss_j_mu_x(i_fit)
coef_fit = coef_gauss_j_mu_x(i_fit)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
if(dabs(int_fit) .lt. 1d-10) cycle
tmp += coef * coef_fit * int_fit
! ---
do i_1s = 2, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit
enddo
! ---
enddo
v_ij_u_cst_mu_j1b(j,i,ipoint) = tmp

View File

@ -45,6 +45,9 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
! 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)
!print *, ' sym integ = ', integral_sym
!print *, ' non-sym integ = ', integral_nsym
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

View File

@ -58,7 +58,7 @@
accu_nd = accu_nd/dble(mo_num**2-mo_num)
if(dabs(accu_d-1.d0).gt.1.d-10.or.dabs(accu_nd).gt.1.d-10)then
print*,'Warning !!!'
print*,'Average trace of overlap_bi_ortho is different from 1 by ', accu_d
print*,'Average trace of overlap_bi_ortho is different from 1 by ', dabs(accu_d-1.d0)
print*,'And bi orthogonality is off by an average of ',accu_nd
print*,'****************'
print*,'Overlap matrix betwee mo_l_coef and mo_r_coef '

View File

@ -356,6 +356,7 @@ subroutine non_hrmt_real_diag(n, A, leigvec, reigvec, n_real_eigv, eigval)
! Eigvalue(n) = WR(n) + i * WI(n)
allocate(WR(n), WI(n), VL(n,n), VR(n,n), Aw(n,n))
Aw = A
!print *, ' matrix to diagonalize', Aw
call lapack_diag_non_sym(n, Aw, WR, WI, VL, VR)
! ---
@ -1212,6 +1213,7 @@ subroutine impose_orthog_svd(n, m, C)
num_linear_dependencies = 0
do i = 1, m
if(abs(D(i)) <= threshold) then
write(*,*) ' D(i) = ', D(i)
D(i) = 0.d0
num_linear_dependencies += 1
else
@ -1690,11 +1692,13 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
! S(i,i) = -1
do i = 1, m
if( (S(i,i) + 1.d0) .lt. thr_d ) then
if(S(i,i) .lt. 0.d0) then
!if( (S(i,i) + 1.d0) .lt. thr_d ) then
do j = 1, n
Vl(j,i) = -1.d0 * Vl(j,i)
enddo
S(i,i) = 1.d0
!S(i,i) = 1.d0
S(i,i) = -S(i,i)
endif
enddo
@ -1726,6 +1730,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
Vr(j,i) = Vr(j,i) * s_tmp
enddo
endif
enddo
endif
@ -1978,7 +1983,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
do i = 1, n
if(deg_num(i).gt.1) then
print *, ' degen on', i, deg_num(i)
print *, ' degen on', i, deg_num(i), e0(i)
endif
enddo
@ -2001,6 +2006,8 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
call impose_orthog_svd(n, m, L)
call impose_orthog_svd(n, m, R)
!call impose_orthog_GramSchmidt(n, m, L)
!call impose_orthog_GramSchmidt(n, m, R)
! ---

View File

@ -0,0 +1,142 @@
! ---
program print_he_tc_energy
implicit none
call print_overlap()
call print_energy1()
end
! ---
subroutine print_overlap()
implicit none
integer :: i, j, k, l
double precision :: S_ij
print *, ' ao_overlap:'
do i = 1, ao_num
do j = 1, ao_num
print *, j, i, ao_overlap(j,i)
enddo
enddo
print *, ' mo_overlap:'
do i = 1, mo_num
do j = 1, mo_num
S_ij = 0.d0
do k = 1, ao_num
do l = 1, ao_num
S_ij += mo_l_coef(k,i) * ao_overlap(k,l) * mo_r_coef(l,j)
enddo
enddo
print *, i, j, S_ij
enddo
enddo
end subroutine print_overlap
! ---
subroutine print_energy1()
implicit none
integer :: i, j, k, l
double precision :: e, n, e_tmp, n_tmp, e_ns
double precision, external :: ao_two_e_integral
e = 0.d0
n = 0.d0
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! < phi_1 phi_1 | h1 | phi_1 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_l_coef(i,1) * ao_one_e_integrals(i,j) * mo_r_coef(j,1)
n_tmp += mo_l_coef(i,1) * ao_overlap(i,j) * mo_r_coef(j,1)
enddo
enddo
e += e_tmp * n_tmp
! ---
! < phi_1 phi_1 | h2 | phi_1 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
n_tmp += mo_l_coef(i,1) * ao_overlap(i,j) * mo_r_coef(j,1)
e_tmp += mo_l_coef(i,1) * ao_one_e_integrals(i,j) * mo_r_coef(j,1)
enddo
enddo
e += e_tmp * n_tmp
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! ---
e_ns = 0.d0
do i = 1, ao_num
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
! ao_two_e_tc_tot(i,j,k,l) = <k i| V^TC(r_12) |l j>
e += mo_l_coef(i,1) * mo_l_coef(k,1) * ao_two_e_tc_tot(i,j,k,l) * mo_r_coef(j,1) * mo_r_coef(l,1)
e_ns += mo_l_coef(i,1) * mo_l_coef(k,1) * ao_non_hermit_term_chemist(i,j,k,l) * mo_r_coef(j,1) * mo_r_coef(l,1)
enddo
enddo
enddo
enddo
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! ---
! < phi_1 phi_1 | phi_1 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_l_coef(i,1) * ao_overlap(i,j) * mo_r_coef(j,1)
n_tmp += mo_l_coef(i,1) * ao_overlap(i,j) * mo_r_coef(j,1)
enddo
enddo
n += e_tmp * n_tmp
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
e = e / n
e_ns = e_ns / n
print *, ' tc energy = ', e
print *, ' non-sym energy = ', e_ns
end subroutine print_energy1
! ---

View File

@ -18,6 +18,10 @@
do j = 1, N_det
! < J | Htilde | I >
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
print *, ' hmono = ', hmono
print *, ' htwoe = ', htwoe
print *, ' hthree = ', hthree
htilde_matrix_elmt_bi_ortho(j,i) = htot
enddo
enddo

View File

@ -0,0 +1,78 @@
! ---
program tc_petermann_factor
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
! my_n_pt_r_grid = 10 ! small grid for quick debug
! my_n_pt_a_grid = 26 ! small grid for quick debug
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call main()
end
! ---
subroutine main()
implicit none
integer :: i, j
double precision :: Pf_diag_av
double precision, allocatable :: Sl(:,:), Sr(:,:), Pf(:,:)
allocate(Sl(mo_num,mo_num), Sr(mo_num,mo_num), Pf(mo_num,mo_num))
call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
, mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
, 0.d0, Sl, size(Sl, 1) )
print *, ''
print *, ' left-orthog matrix:'
do i = 1, mo_num
write(*,'(100(F8.4,X))') Sl(:,i)
enddo
call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
, mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
, 0.d0, Sr, size(Sr, 1) )
print *, ''
print *, ' right-orthog matrix:'
do i = 1, mo_num
write(*,'(100(F8.4,X))') Sr(:,i)
enddo
print *, ''
print *, ' Petermann matrix:'
do i = 1, mo_num
do j = 1, mo_num
Pf(j,i) = Sl(j,i) * Sr(j,i)
enddo
write(*,'(100(F8.4,X))') Pf(:,i)
enddo
Pf_diag_av = 0.d0
do i = 1, mo_num
Pf_diag_av = Pf_diag_av + Pf(i,i)
enddo
Pf_diag_av = Pf_diag_av / dble(mo_num)
print *, ''
print *, ' mean of the diagonal Petermann factor = ', Pf_diag_av
deallocate(Sl, Sr, Pf)
return
end subroutine
! ---

View File

@ -5,13 +5,48 @@ program print_he_energy
implicit none
!call print_energy1()
call print_overlap()
call print_energy1()
call print_energy2()
end
! ---
subroutine print_overlap()
implicit none
integer :: i, j, k, l
double precision :: S_ij
print *, ' ao_overlap:'
do i = 1, ao_num
do j = 1, ao_num
print *, j, i, ao_overlap(j,i)
enddo
enddo
print *, ' mo_overlap:'
do i = 1, mo_num
do j = 1, mo_num
S_ij = 0.d0
do k = 1, ao_num
do l = 1, ao_num
S_ij += mo_coef(k,i) * ao_overlap(k,l) * mo_coef(l,j)
enddo
enddo
print *, i, j, S_ij
enddo
enddo
end subroutine print_overlap
! ---
subroutine print_energy1()
implicit none