diff --git a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f index fe25e9c0..ef43a78b 100644 --- a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f @@ -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 diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index 50b4bf96..6e3cfed4 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -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 diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f index 552e7069..9486f153 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f @@ -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 diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f index 89f46a05..31cf0624 100644 --- a/src/bi_ort_ints/total_twoe_pot.irp.f +++ b/src/bi_ort_ints/total_twoe_pot.irp.f @@ -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 diff --git a/src/bi_ortho_mos/overlap.irp.f b/src/bi_ortho_mos/overlap.irp.f index 1c88d16f..d7f45c94 100644 --- a/src/bi_ortho_mos/overlap.irp.f +++ b/src/bi_ortho_mos/overlap.irp.f @@ -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 ' diff --git a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f index 8d47f124..4430109d 100644 --- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -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) ! --- diff --git a/src/tc_bi_ortho/print_he_tc_energy.irp.f b/src/tc_bi_ortho/print_he_tc_energy.irp.f new file mode 100644 index 00000000..84d34bcb --- /dev/null +++ b/src/tc_bi_ortho/print_he_tc_energy.irp.f @@ -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) = + 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 + +! --- + + diff --git a/src/tc_bi_ortho/tc_hmat.irp.f b/src/tc_bi_ortho/tc_hmat.irp.f index bf1388e5..49ea6f05 100644 --- a/src/tc_bi_ortho/tc_hmat.irp.f +++ b/src/tc_bi_ortho/tc_hmat.irp.f @@ -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 diff --git a/src/tc_scf/tc_petermann_factor.irp.f b/src/tc_scf/tc_petermann_factor.irp.f new file mode 100644 index 00000000..d3722098 --- /dev/null +++ b/src/tc_scf/tc_petermann_factor.irp.f @@ -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 + +! --- + diff --git a/src/tools/print_he_energy.irp.f b/src/tools/print_he_energy.irp.f index 321b6a92..87488fba 100644 --- a/src/tools/print_he_energy.irp.f +++ b/src/tools/print_he_energy.irp.f @@ -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