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 b847a630..fab50805 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 @@ -90,7 +90,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid, 3)] BEGIN_DOC - ! int dr x * phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| + ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| END_DOC implicit none @@ -119,7 +119,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_num, n_points_final_grid)] BEGIN_DOC - ! int dr x * phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| + ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| END_DOC implicit none diff --git a/src/bi_ort_ints/semi_num_ints_mo.irp.f b/src/bi_ort_ints/semi_num_ints_mo.irp.f index 6c4b44c0..27fcb7de 100644 --- a/src/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f @@ -1,81 +1,154 @@ -BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu, ( mo_num, mo_num,n_points_final_grid)] - implicit none - BEGIN_DOC -! mo_v_ki_bi_ortho_erf_rk_cst_mu(k,i,ip) = int dr chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1 )/(2|r - R_ip|) on the BI-ORTHO MO basis -! -! where phi_k(r) is a LEFT MOs and phi_i(r) is a RIGHT MO -! -! R_ip = the "ip"-th point of the DFT Grid - END_DOC - integer :: ipoint - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & + +! --- + +! TODO :: optimization : transform into a DGEMM + +BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! mo_v_ki_bi_ortho_erf_rk_cst_mu(k,i,ip) = int dr chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1 )/(2|r - R_ip|) on the BI-ORTHO MO basis + ! + ! where phi_k(r) is a LEFT MOs and phi_i(r) is a RIGHT MO + ! + ! R_ip = the "ip"-th point of the DFT Grid + ! + END_DOC + + implicit none + integer :: ipoint + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint) & !$OMP SHARED (n_points_final_grid,v_ij_erf_rk_cst_mu,mo_v_ki_bi_ortho_erf_rk_cst_mu) !$OMP DO SCHEDULE (dynamic) -! TODO :: optimization : transform into a DGEMM - do ipoint = 1, n_points_final_grid - call ao_to_mo_bi_ortho(v_ij_erf_rk_cst_mu(1,1,ipoint),size(v_ij_erf_rk_cst_mu,1),mo_v_ki_bi_ortho_erf_rk_cst_mu(1,1,ipoint),size(mo_v_ki_bi_ortho_erf_rk_cst_mu,1)) - enddo + do ipoint = 1, n_points_final_grid + call ao_to_mo_bi_ortho( v_ij_erf_rk_cst_mu (1,1,ipoint), size(v_ij_erf_rk_cst_mu, 1) & + , mo_v_ki_bi_ortho_erf_rk_cst_mu(1,1,ipoint), size(mo_v_ki_bi_ortho_erf_rk_cst_mu, 1) ) + enddo !$OMP END DO !$OMP END PARALLEL - mo_v_ki_bi_ortho_erf_rk_cst_mu = mo_v_ki_bi_ortho_erf_rk_cst_mu * 0.5d0 -END_PROVIDER - -BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu_transp, ( n_points_final_grid,mo_num, mo_num)] - implicit none - BEGIN_DOC -! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the BI-ORTHO MO basis - END_DOC - integer :: ipoint,i,j - do i = 1, mo_num - do j = 1, mo_num - do ipoint = 1, n_points_final_grid - mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,j,i) = mo_v_ki_bi_ortho_erf_rk_cst_mu(j,i,ipoint) - enddo - enddo - enddo -! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu -END_PROVIDER - -BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, ( mo_num, mo_num,3,n_points_final_grid)] - implicit none - BEGIN_DOC -! mo_x_v_ki_bi_ortho_erf_rk_cst_mu(k,i,m,ip) = int dr x(m) * chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1)/2|r - R_ip| on the BI-ORTHO MO basis -! -! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => x(m) = x, m=2 => x(m) = y, m=3 => x(m) = z, -! -! R_ip = the "ip"-th point of the DFT Grid - END_DOC - integer :: ipoint,m - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint,m) & - !$OMP SHARED (n_points_final_grid,x_v_ij_erf_rk_cst_mu_transp,mo_x_v_ki_bi_ortho_erf_rk_cst_mu) - !$OMP DO SCHEDULE (dynamic) -! TODO :: optimization : transform into a DGEMM - do ipoint = 1, n_points_final_grid - do m = 1, 3 - call ao_to_mo_bi_ortho(x_v_ij_erf_rk_cst_mu_transp(1,1,m,ipoint),size(x_v_ij_erf_rk_cst_mu_transp,1),mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,m,ipoint),size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu,1)) - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - mo_x_v_ki_bi_ortho_erf_rk_cst_mu = 0.5d0 * mo_x_v_ki_bi_ortho_erf_rk_cst_mu + mo_v_ki_bi_ortho_erf_rk_cst_mu = mo_v_ki_bi_ortho_erf_rk_cst_mu * 0.5d0 END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, 3, mo_num, mo_num)] + +BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, mo_num, mo_num)] + + BEGIN_DOC + ! + ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the BI-ORTHO MO basis + ! + END_DOC + implicit none - integer :: i, j, m, ipoint + integer :: ipoint, i, j + do i = 1, mo_num do j = 1, mo_num - do m = 1, 3 - do ipoint = 1, n_points_final_grid - mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,m,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,m,ipoint) - enddo + do ipoint = 1, n_points_final_grid + mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,j,i) = mo_v_ki_bi_ortho_erf_rk_cst_mu(j,i,ipoint) + enddo + enddo + enddo + +! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu + +END_PROVIDER + +! --- + +! TODO :: optimization : transform into a DGEMM + +BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, 3, n_points_final_grid)] + + BEGIN_DOC + ! + ! mo_x_v_ki_bi_ortho_erf_rk_cst_mu(k,i,m,ip) = int dr x(m) * chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1)/2|r - R_ip| on the BI-ORTHO MO basis + ! + ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => x(m) = x, m=2 => x(m) = y, m=3 => x(m) = z, + ! + ! R_ip = the "ip"-th point of the DFT Grid + ! + END_DOC + + implicit none + integer :: ipoint + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint) & + !$OMP SHARED (n_points_final_grid,x_v_ij_erf_rk_cst_mu_transp,mo_x_v_ki_bi_ortho_erf_rk_cst_mu) + !$OMP DO SCHEDULE (dynamic) + do ipoint = 1, n_points_final_grid + + 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) ) + + enddo + !$OMP END DO + !$OMP END PARALLEL + + mo_x_v_ki_bi_ortho_erf_rk_cst_mu = 0.5d0 * mo_x_v_ki_bi_ortho_erf_rk_cst_mu + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int2_grad1_u12_bimo(:,k,i,ipoint) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \chi_k(r2) \phi_i(r2) + ! + END_DOC + + implicit none + integer :: ipoint + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint) & + !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao,int2_grad1_u12_bimo) + !$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) ) + + enddo + !$OMP END DO + !$OMP END PARALLEL + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, 3, mo_num, mo_num)] + + implicit none + integer :: i, j, ipoint + + do i = 1, mo_num + do j = 1, mo_num + do ipoint = 1, n_points_final_grid + mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,1,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,1,ipoint) + mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,2,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,2,ipoint) + mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,3,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,3,ipoint) enddo enddo enddo @@ -83,14 +156,15 @@ END_PROVIDER ! --- - BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, 3, mo_num, mo_num)] + BEGIN_DOC - ! x_W_ki_bi_ortho_erf_rk(ip,m,k,i) = \int dr chi_k(r) (1 - erf(mu |r-R_ip|)) (x(m)-X(m)_ip) phi_i(r) ON THE BI-ORTHO MO BASIS -! -! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z, -! -! R_ip = the "ip"-th point of the DFT Grid + ! + ! x_W_ki_bi_ortho_erf_rk(ip,m,k,i) = \int dr chi_k(r) \frac{(1 - erf(mu |r-R_ip|))}{2|r-R_ip|} (x(m)-R_ip(m)) phi_i(r) ON THE BI-ORTHO MO BASIS + ! + ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z, + ! + ! R_ip = the "ip"-th point of the DFT Grid END_DOC implicit none @@ -100,7 +174,7 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, double precision :: xyz double precision :: wall0, wall1 - print*,'providing x_W_ki_bi_ortho_erf_rk ...' + print*, ' providing x_W_ki_bi_ortho_erf_rk ...' call wall_time(wall0) !$OMP PARALLEL & @@ -126,7 +200,7 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, ! FREE mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp call wall_time(wall1) - print*,'time to provide x_W_ki_bi_ortho_erf_rk = ',wall1 - wall0 + print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0 END_PROVIDER diff --git a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f index 1fe27ab1..12361ace 100644 --- a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f +++ b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f @@ -53,26 +53,55 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n END_PROVIDER -subroutine give_integrals_3_body_bi_ort(n,l,k,m,j,i,integral) - implicit none - double precision, intent(out) :: integral - integer, intent(in) :: n,l,k,m,j,i - double precision :: weight - BEGIN_DOC -! with a BI ORTHONORMAL ORBITALS - END_DOC - integer :: ipoint,mm - integral = 0.d0 - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & - * x_W_ki_bi_ortho_erf_rk(ipoint,mm,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,mm,l,j) - integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & - * x_W_ki_bi_ortho_erf_rk(ipoint,mm,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,mm,k,i) - integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & - * x_W_ki_bi_ortho_erf_rk(ipoint,mm,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,mm,k,i) - enddo - enddo -end +! --- + +subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) + + BEGIN_DOC + ! + ! < n l k | -L | m j i > with a BI-ORTHONORMAL ORBITALS + ! + END_DOC + + implicit none + integer, intent(in) :: n,l,k,m,j,i + double precision, intent(out) :: integral + integer :: ipoint + double precision :: weight + + integral = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + +! integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & +! * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,1,l,j) & +! + x_W_ki_bi_ortho_erf_rk(ipoint,2,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,2,l,j) & +! + x_W_ki_bi_ortho_erf_rk(ipoint,3,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,3,l,j) ) +! integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & +! * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,1,k,i) & +! + x_W_ki_bi_ortho_erf_rk(ipoint,2,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,2,k,i) & +! + x_W_ki_bi_ortho_erf_rk(ipoint,3,n,m) * x_W_ki_bi_ortho_erf_rk(ipoint,3,k,i) ) +! integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & +! * ( x_W_ki_bi_ortho_erf_rk(ipoint,1,l,j) * x_W_ki_bi_ortho_erf_rk(ipoint,1,k,i) & +! + 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) ) + + enddo + +end subroutine give_integrals_3_body_bi_ort + +! --- diff --git a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f index e59b5f7a..7b99cc91 100644 --- a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f +++ b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f @@ -34,7 +34,7 @@ program debug_integ_jmu_modif !call test_int2_u2_j1b2() !call test_int2_grad1u2_grad2u2_j1b2() - !call test_grad_1_u_ij_mu() + !call test_int2_grad1_u12_ao() !call test_gradu_squared_u_ij_mu() end @@ -287,16 +287,16 @@ end subroutine test_int2_grad1u2_grad2u2_j1b2 ! --- -subroutine test_grad_1_u_ij_mu() +subroutine test_int2_grad1_u12_ao() implicit none integer :: i, j, ipoint double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz double precision :: integ(3) - print*, ' test_grad_1_u_ij_mu ...' + print*, ' test_int2_grad1_u12_ao ...' - PROVIDE grad_1_u_ij_mu + PROVIDE int2_grad1_u12_ao eps_ij = 1d-6 acc_tot = 0.d0 @@ -305,13 +305,13 @@ subroutine test_grad_1_u_ij_mu() do j = 1, ao_num do i = 1, ao_num - call num_grad_1_u_ij_mu(i, j, ipoint, integ) + call num_int2_grad1_u12_ao(i, j, ipoint, integ) - i_exc = grad_1_u_ij_mu(i,j,ipoint,1) + i_exc = int2_grad1_u12_ao(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 grad_1_u_ij_mu on', i, j, ipoint + print *, ' problem in x part of int2_grad1_u12_ao on', i, j, ipoint print *, ' analyt integ = ', i_exc print *, ' numeri integ = ', i_num print *, ' diff = ', acc_ij @@ -319,11 +319,11 @@ subroutine test_grad_1_u_ij_mu() acc_tot += acc_ij normalz += dabs(i_num) - i_exc = grad_1_u_ij_mu(i,j,ipoint,2) + i_exc = int2_grad1_u12_ao(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 grad_1_u_ij_mu on', i, j, ipoint + print *, ' problem in y part of int2_grad1_u12_ao on', i, j, ipoint print *, ' analyt integ = ', i_exc print *, ' numeri integ = ', i_num print *, ' diff = ', acc_ij @@ -331,11 +331,11 @@ subroutine test_grad_1_u_ij_mu() acc_tot += acc_ij normalz += dabs(i_num) - i_exc = grad_1_u_ij_mu(i,j,ipoint,3) + i_exc = int2_grad1_u12_ao(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 grad_1_u_ij_mu on', i, j, ipoint + print *, ' problem in z part of int2_grad1_u12_ao on', i, j, ipoint print *, ' analyt integ = ', i_exc print *, ' numeri integ = ', i_num print *, ' diff = ', acc_ij @@ -352,7 +352,7 @@ subroutine test_grad_1_u_ij_mu() print*, ' normalz = ', normalz return -end subroutine test_grad_1_u_ij_mu +end subroutine test_int2_grad1_u12_ao ! --- diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index b9c98ea9..bf37c551 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -4,20 +4,34 @@ ! TODO : strong optmization : write the loops in a different way ! : for each couple of AO, the gaussian product are done once for all -BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num,n_points_final_grid)] +BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_points_final_grid) ] BEGIN_DOC - ! - ! -1/2 [ (grad_1 u)^2 + (grad_2 u^2)] = - 1/4 * (1 - erf(mu*r12))^2 ! + ! if J(r1,r2) = u12: + ! + ! gradu_squared_u_ij_mu = -0.50 x \int r2 [ (grad_1 u12)^2 + (grad_2 u12^2)] \phi_i(2) \phi_j(2) + ! = -0.25 x \int r2 (1 - erf(mu*r12))^2 \phi_i(2) \phi_j(2) ! and ! (1 - erf(mu*r12))^2 = \sum_i coef_gauss_1_erf_x_2(i) * exp(-expo_gauss_1_erf_x_2(i) * r12^2) ! + ! if J(r1,r2) = u12 x v1 x v2 + ! + ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2 v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + u12^2 v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ] + ! = -0.25 x v1^2 \int r2 \phi_i(2) \phi_j(2) [1 - erf(mu r12)]^2 v2^2 + ! + -0.50 x (grad_1 v1)^2 \int r2 \phi_i(2) \phi_j(2) u12^2 v2^2 + ! + -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 + ! + ! END_DOC implicit none integer :: ipoint, i, j, m, igauss - double precision :: r(3), delta, coef, tmp + double precision :: r(3), delta, coef + double precision :: tmp_v, tmp_x, tmp_y, tmp_z, tmp1, tmp2, tmp3, tmp4, tmp5 double precision :: time0, time1 double precision, external :: overlap_gauss_r12_ao @@ -27,13 +41,28 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num,n_poin PROVIDE j1b_type j1b_pen if(j1b_type .eq. 3) then - ! v1_1b^2 \int d2 \phi_i(2) \phi_j(2) \frac{-[1 - \erf(\mu r12)]^2}{4} v2_1b^2 do ipoint = 1, n_points_final_grid - tmp = v_1b(ipoint) * v_1b(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) + tmp3 = tmp_v * tmp_x + tmp4 = tmp_v * tmp_y + tmp5 = tmp_v * tmp_z + do j = 1, ao_num do i = 1, ao_num - gradu_squared_u_ij_mu(j,i,ipoint) += tmp * int2_grad1u2_grad2u2_j1b2(i,j,ipoint) + + 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) enddo enddo enddo @@ -74,22 +103,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 + double precision :: contrib, weight1, ao_k_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) + weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num - contrib = weight1 * 0.5d0 * (aos_in_r_array_transp(ipoint,k) * aos_in_r_array_transp(ipoint,i)) - ! \int dr1 phi_k(r1) phi_i(r1) . \int dr2 |\grad_1 u(r1,r2)|^2 \phi_l(r2) \phi_j(r2) - ac_mat(k,i,l,j) += gradu_squared_u_ij_mu(l,j,ipoint) * contrib + do i = 1, ao_num + ao_i_r = aos_in_r_array_transp(ipoint,i) + + do k = 1, ao_num + ao_k_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 enddo enddo enddo diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index 9832e81d..db659520 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -1,15 +1,28 @@ ! --- -BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)] +BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (3, ao_num, ao_num, n_points_final_grid)] BEGIN_DOC ! - ! grad_1_u_ij_mu(i,j,ipoint) = \int dr2 [-1 * \grad_r1 u(r1,r2)] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2) - ! = \int dr2 [(r1 - r2) (erf(mu * r12)-1)/2 r_12] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2) + ! int2_grad1_u12_ao(:,i,j,ipoint) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) ! ! where r1 = r(ipoint) ! + ! if J(r1,r2) = u12: + ! + ! int2_grad1_u12_ao(:,i,j,ipoint) = 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r2) \phi_j(r2) + ! = 0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ] + ! + ! 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) ] + ! = 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) + ! + ! END_DOC implicit none @@ -25,10 +38,10 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin y = final_grid_points(2,ipoint) z = final_grid_points(3,ipoint) - tmp0 = v_1b (ipoint) - tmp_x = v_1b_grad(1,ipoint) - tmp_y = v_1b_grad(2,ipoint) - tmp_z = v_1b_grad(3,ipoint) + tmp0 = 0.5d0 * v_1b(ipoint) + tmp_x = v_1b_grad(1,ipoint) + tmp_y = v_1b_grad(2,ipoint) + tmp_z = v_1b_grad(3,ipoint) do j = 1, ao_num do i = 1, ao_num @@ -36,9 +49,9 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) - grad_1_u_ij_mu(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) + tmp_x * tmp2 - grad_1_u_ij_mu(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) + tmp_y * tmp2 - grad_1_u_ij_mu(i,j,ipoint,3) = 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_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 enddo enddo enddo @@ -51,61 +64,87 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin z = final_grid_points(3,ipoint) do j = 1, ao_num do i = 1, ao_num - grad_1_u_ij_mu(i,j,ipoint,1) = v_ij_erf_rk_cst_mu(i,j,ipoint) * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1) - grad_1_u_ij_mu(i,j,ipoint,2) = v_ij_erf_rk_cst_mu(i,j,ipoint) * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2) - grad_1_u_ij_mu(i,j,ipoint,3) = v_ij_erf_rk_cst_mu(i,j,ipoint) * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3) + 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) enddo enddo enddo - endif + int2_grad1_u12_ao *= 0.5d0 - grad_1_u_ij_mu *= 0.5d0 + endif END_PROVIDER ! --- BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)] - implicit none - BEGIN_DOC - ! tc_grad_and_lapl_ao(k,i,l,j) = - ! - ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) - ! - ! This is obtained by integration by parts. - END_DOC - integer :: ipoint,i,j,k,l,m - double precision :: contrib,weight1 - double precision, allocatable :: ac_mat(:,:,:,:) - allocate(ac_mat(ao_num, ao_num, ao_num, ao_num)) - ac_mat = 0.d0 - do m = 1, 3 - do ipoint = 1, n_points_final_grid - weight1 = final_weight_at_r_vector(ipoint) - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num - contrib = weight1 *0.5D0* (aos_in_r_array_transp(ipoint,k) * aos_grad_in_r_array_transp_bis(ipoint,i,m) & - -aos_in_r_array_transp(ipoint,i) * aos_grad_in_r_array_transp_bis(ipoint,k,m) ) - ! \int dr1 phi_k(r1) \grad_r1 phi_i(r1) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) - ac_mat(k,i,l,j) += grad_1_u_ij_mu(l,j,ipoint,m) * contrib - enddo - 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_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + BEGIN_DOC + ! + ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij > + ! + ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) + ! + ! This is obtained by integration by parts. + ! + END_DOC + + 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, 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 = 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) + + 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) + + 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 = weight1 * ( contrib_x + contrib_y + contrib_z ) + + ac_mat(k,i,l,j) += contrib + enddo + enddo + 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_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + enddo + enddo + enddo + enddo + + deallocate(ac_mat) END_PROVIDER + +! --- + diff --git a/src/non_h_ints_mu/numerical_integ.irp.f b/src/non_h_ints_mu/numerical_integ.irp.f index 17b666aa..dae68649 100644 --- a/src/non_h_ints_mu/numerical_integ.irp.f +++ b/src/non_h_ints_mu/numerical_integ.irp.f @@ -1,51 +1,6 @@ ! --- -! -! \int dr2 [-1 * \grad_r1 u(r1,r2)] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2) -! - -!BEGIN_PROVIDER [ double precision, num_grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)] -! -! implicit none -! -! integer :: i, j, ipoint, jpoint -! double precision :: tmp, r1(3), r2(3), grad(3) -! -! double precision, external :: ao_value -! double precision, external :: j12_nucl -! -! num_grad_1_u_ij_mu = 0.d0 -! -! do j = 1, ao_num -! do i = 1, ao_num -! -! 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) -! -! 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 = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint) -! -! call grad1_j12_mu_exc(r1, r2, grad) -! -! num_grad_1_u_ij_mu(i,j,ipoint,1) += tmp * (-1.d0 * grad(1)) -! num_grad_1_u_ij_mu(i,j,ipoint,2) += tmp * (-1.d0 * grad(2)) -! num_grad_1_u_ij_mu(i,j,ipoint,3) += tmp * (-1.d0 * grad(3)) -! enddo -! -! enddo -! enddo -! enddo -! -!END_PROVIDER - -! --- - double precision function num_v_ij_u_cst_mu_j1b(i, j, ipoint) BEGIN_DOC @@ -289,7 +244,7 @@ end subroutine num_x_v_ij_erf_rk_cst_mu_j1b ! --- -subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ) +subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ) implicit none @@ -328,7 +283,7 @@ subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ) integ(3) = tmp_z return -end subroutine num_grad_1_u_ij_mu +end subroutine num_int2_grad1_u12_ao ! ---