From 1f56b5d0f4eb8f6f68ef83ef6744a26807514807 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Tue, 5 Sep 2023 11:52:08 +0200 Subject: [PATCH] num integ grad & grad squared --- src/becke_numerical_grid/EZFIO.cfg | 4 + src/becke_numerical_grid/extra_grid.irp.f | 4 +- .../extra_grid_vector.irp.f | 34 +- src/becke_numerical_grid/grid_becke.irp.f | 2 +- src/non_h_ints_mu/debug_fit.irp.f | 77 +- src/non_h_ints_mu/jast_deriv.irp.f | 1165 ----------------- src/non_h_ints_mu/jast_deriv_utils.irp.f | 700 ++++++++++ src/non_h_ints_mu/jast_deriv_utils_vect.irp.f | 332 +++++ src/non_h_ints_mu/tc_integ.irp.f | 530 -------- src/non_h_ints_mu/tc_integ_an.irp.f | 244 ++++ src/non_h_ints_mu/tc_integ_num.irp.f | 190 +++ src/non_h_ints_mu/test_non_h_ints.irp.f | 6 +- src/tc_keywords/EZFIO.cfg | 11 + 13 files changed, 1584 insertions(+), 1715 deletions(-) create mode 100644 src/non_h_ints_mu/jast_deriv_utils.irp.f create mode 100644 src/non_h_ints_mu/jast_deriv_utils_vect.irp.f delete mode 100644 src/non_h_ints_mu/tc_integ.irp.f create mode 100644 src/non_h_ints_mu/tc_integ_an.irp.f create mode 100644 src/non_h_ints_mu/tc_integ_num.irp.f diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg index 7861f074..e660fd6d 100644 --- a/src/becke_numerical_grid/EZFIO.cfg +++ b/src/becke_numerical_grid/EZFIO.cfg @@ -33,6 +33,10 @@ doc: Number of angular grid points given from input. Warning, this number cannot interface: ezfio,provider,ocaml default: 1202 +[n_points_extra_final_grid] +type: integer +doc: Total number of extra_grid points +interface: ezfio [extra_grid_type_sgn] type: integer diff --git a/src/becke_numerical_grid/extra_grid.irp.f b/src/becke_numerical_grid/extra_grid.irp.f index 9bd24f22..7df4dd6d 100644 --- a/src/becke_numerical_grid/extra_grid.irp.f +++ b/src/becke_numerical_grid/extra_grid.irp.f @@ -14,7 +14,7 @@ implicit none - if(.not.my_extra_grid_becke)then + if(.not. my_extra_grid_becke) then select case (extra_grid_type_sgn) case(0) n_points_extra_radial_grid = 23 @@ -33,7 +33,7 @@ stop end select else - n_points_extra_radial_grid = my_n_pt_r_extra_grid + n_points_extra_radial_grid = my_n_pt_r_extra_grid n_points_extra_integration_angular = my_n_pt_a_extra_grid endif diff --git a/src/becke_numerical_grid/extra_grid_vector.irp.f b/src/becke_numerical_grid/extra_grid_vector.irp.f index e4fc03b5..ae167282 100644 --- a/src/becke_numerical_grid/extra_grid_vector.irp.f +++ b/src/becke_numerical_grid/extra_grid_vector.irp.f @@ -23,29 +23,33 @@ BEGIN_PROVIDER [integer, n_points_extra_final_grid] enddo enddo - print*,'n_points_extra_final_grid = ',n_points_extra_final_grid - print*,'n max point = ',n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1) -! call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid) + print*, ' n_points_extra_final_grid = ', n_points_extra_final_grid + print*, ' n max point = ', n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1) + call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid) + END_PROVIDER ! --- BEGIN_PROVIDER [double precision, final_grid_points_extra, (3,n_points_extra_final_grid)] -&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_extra, (n_points_extra_final_grid) ] -&BEGIN_PROVIDER [integer, index_final_points_extra, (3,n_points_extra_final_grid) ] -&BEGIN_PROVIDER [integer, index_final_points_extra_reverse, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ] - implicit none +&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_extra, (n_points_extra_final_grid)] +&BEGIN_PROVIDER [integer, index_final_points_extra, (3,n_points_extra_final_grid)] +&BEGIN_PROVIDER [integer, index_final_points_extra_reverse, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num)] + BEGIN_DOC -! final_grid_points_extra(1:3,j) = (/ x, y, z /) of the jth grid point -! -! final_weight_at_r_vector_extra(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions -! -! index_final_points_extra(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point -! -! index_final_points_extra_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices + ! final_grid_points_extra(1:3,j) = (/ x, y, z /) of the jth grid point + ! + ! final_weight_at_r_vector_extra(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions + ! + ! index_final_points_extra(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point + ! + ! index_final_points_extra_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices END_DOC + + implicit none integer :: i,j,k,l,i_count double precision :: r(3) + i_count = 0 do j = 1, nucl_num do i = 1, n_points_extra_radial_grid -1 @@ -67,3 +71,5 @@ END_PROVIDER enddo END_PROVIDER + + diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f index 21b9f98d..91fdc563 100644 --- a/src/becke_numerical_grid/grid_becke.irp.f +++ b/src/becke_numerical_grid/grid_becke.irp.f @@ -14,7 +14,7 @@ implicit none - if(.not.my_grid_becke)then + if(.not. my_grid_becke) then select case (grid_type_sgn) case(0) n_points_radial_grid = 23 diff --git a/src/non_h_ints_mu/debug_fit.irp.f b/src/non_h_ints_mu/debug_fit.irp.f index 05d2db68..d3152836 100644 --- a/src/non_h_ints_mu/debug_fit.irp.f +++ b/src/non_h_ints_mu/debug_fit.irp.f @@ -13,17 +13,27 @@ program debug_fit PROVIDE mu_erf j1b_pen + if(j1b_type .ge. 100) then + my_extra_grid_becke = .True. + PROVIDE tc_grid2_a tc_grid2_r + my_n_pt_r_extra_grid = tc_grid2_r + my_n_pt_a_extra_grid = tc_grid2_a + touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid + endif + !call test_j1b_nucl() !call test_grad_j1b_nucl() !call test_lapl_j1b_nucl() !call test_list_b2() - call test_list_b3() + !call test_list_b3() !call test_fit_u() !call test_fit_u2() !call test_fit_ugradu() + call test_grad1_u12_withsq_num() + end ! --- @@ -643,4 +653,69 @@ end subroutine test_fit_u2 ! --- +subroutine test_grad1_u12_withsq_num() + + implicit none + integer :: ipoint, jpoint, m + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision, allocatable :: tmp_grad1_u12_squared(:,:), tmp_grad1_u12(:,:,:) + + print*, ' test_grad1_u12_withsq_num ...' + + PROVIDE grad1_u12_num grad1_u12_squared_num + + allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_points_final_grid)) + allocate(tmp_grad1_u12(n_points_extra_final_grid,n_points_final_grid,3)) + + eps_ij = 1d-7 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + + call get_grad1_u12_withsq_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12(1,ipoint,1) & + , tmp_grad1_u12(1,ipoint,2) & + , tmp_grad1_u12(1,ipoint,3) & + , tmp_grad1_u12_squared(1,ipoint)) + do jpoint = 1, n_points_extra_final_grid + + i_exc = grad1_u12_squared_num(jpoint,ipoint) + i_num = tmp_grad1_u12_squared(jpoint,ipoint) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in grad1_u12_squared_num on', ipoint, jpoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + stop + endif + acc_tot += acc_ij + normalz += dabs(i_num) + + do m = 1, 3 + i_exc = grad1_u12_num(jpoint,ipoint,m) + i_num = tmp_grad1_u12(jpoint,ipoint,m) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in grad1_u12_num on', ipoint, jpoint, m + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + stop + endif + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + enddo + enddo + + !print*, ' acc_tot = ', acc_tot + !print*, ' normalz = ', normalz + print*, ' accuracy (%) = ', 100.d0 * acc_tot / normalz + + return +end subroutine test_grad1_u12_withsq_num + +! --- + diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index 082412f9..ee01886c 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -123,1168 +123,3 @@ END_PROVIDER ! --- -subroutine get_grad1_u12_r1_seq(r1, n_grid2, resx, resy, resz) - - BEGIN_DOC - ! - ! grad_1 u(r1,r2) - ! - ! this will be integrated numerically over r2: - ! we use grid for r1 and extra_grid for r2 - ! - ! for 99 < j1b_type < 199 - ! - ! u(r1,r2) = j12_mu(r12) x v(r1) x v(r2) - ! grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2) - ! - END_DOC - - implicit none - integer, intent(in) :: n_grid2 - double precision, intent(in) :: r1(3) - double precision, intent(out) :: resx(n_grid2), resy(n_grid2), resz(n_grid2) - - double precision :: v1b_r1 - double precision :: grad1_v1b(3) - double precision, allocatable :: v1b_r2(:) - double precision, allocatable :: u2b_r12(:) - double precision, allocatable :: gradx1_u2b(:), grady1_u2b(:), gradz1_u2b(:) - double precision, external :: j1b_nucl - - PROVIDE j1b_type - PROVIDE final_grid_points_extra - - if( (j1b_type .eq. 100) .or. & - (j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then - - call grad1_j12_mu_r1_seq(r1, n_grid2, resx, resy, resz) - - elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then - - allocate(v1b_r2(n_grid2)) - allocate(u2b_r12(n_grid2)) - allocate(gradx1_u2b(n_grid2)) - allocate(grady1_u2b(n_grid2)) - allocate(gradz1_u2b(n_grid2)) - - v1b_r1 = j1b_nucl(r1) - call grad1_j1b_nucl(r1, grad1_v1b) - - call j1b_nucl_r1_seq(n_grid2, v1b_r2) - call j12_mu_r1_seq(r1, n_grid2, u2b_r12) - call grad1_j12_mu_r1_seq(r1, n_grid2, gradx1_u2b, grady1_u2b, gradz1_u2b) - - resx(:) = (gradx1_u2b(:) * v1b_r1 + u2b_r12(:) * grad1_v1b(1)) * v1b_r2(:) - resy(:) = (grady1_u2b(:) * v1b_r1 + u2b_r12(:) * grad1_v1b(2)) * v1b_r2(:) - resz(:) = (gradz1_u2b(:) * v1b_r1 + u2b_r12(:) * grad1_v1b(3)) * v1b_r2(:) - - deallocate(v1b_r2, u2b_r12, gradx1_u2b, grady1_u2b, gradz1_u2b) - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - return -end subroutine get_grad1_u12_r1_seq - -! --- - -subroutine get_grad1_u12_squared_r1_seq(r1, n_grid2, res) - - BEGIN_DOC - ! - ! grad_1 u(r1,r2) - ! - ! this will be integrated numerically over r2: - ! we use grid for r1 and extra_grid for r2 - ! - ! for 99 < j1b_type < 199 - ! - ! u(r1,r2) = j12_mu(r12) x v(r1) x v(r2) - ! grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2) - ! - END_DOC - - implicit none - integer, intent(in) :: n_grid2 - double precision, intent(in) :: r1(3) - double precision, intent(out) :: res(n_grid2) - - - integer :: jpoint - double precision :: r2(3) - double precision :: v1b_r1, v1b_r2, u2b_r12 - double precision :: grad1_v1b(3), grad1_u2b(3) - double precision :: dx, dy, dz - double precision, external :: j12_mu, j1b_nucl - - PROVIDE j1b_type - PROVIDE final_grid_points_extra - - if( (j1b_type .eq. 100) .or. & - (j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then - - do jpoint = 1, n_points_extra_final_grid ! r2 - - r2(1) = final_grid_points_extra(1,jpoint) - r2(2) = final_grid_points_extra(2,jpoint) - r2(3) = final_grid_points_extra(3,jpoint) - - call grad1_j12_mu(r1, r2, grad1_u2b) - - dx = grad1_u2b(1) - dy = grad1_u2b(2) - dz = grad1_u2b(3) - - res(jpoint) = dx*dx + dy*dy + dz*dz - enddo - - elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then - - v1b_r1 = j1b_nucl(r1) - call grad1_j1b_nucl(r1, grad1_v1b) - - do jpoint = 1, n_points_extra_final_grid ! r2 - - r2(1) = final_grid_points_extra(1,jpoint) - r2(2) = final_grid_points_extra(2,jpoint) - r2(3) = final_grid_points_extra(3,jpoint) - - v1b_r2 = j1b_nucl(r2) - u2b_r12 = j12_mu(r1, r2) - call grad1_j12_mu(r1, r2, grad1_u2b) - - dx = (grad1_u2b(1) * v1b_r1 + u2b_r12 * grad1_v1b(1)) * v1b_r2 - dy = (grad1_u2b(2) * v1b_r1 + u2b_r12 * grad1_v1b(2)) * v1b_r2 - dz = (grad1_u2b(3) * v1b_r1 + u2b_r12 * grad1_v1b(3)) * v1b_r2 - - res(jpoint) = dx*dx + dy*dy + dz*dz - enddo - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - return -end subroutine get_grad1_u12_squared_r1_seq - -! --- - -subroutine get_grad1_u12_withsq_r1_seq(r1, n_grid2, resx, resy, resz, res) - - BEGIN_DOC - ! - ! grad_1 u(r1,r2) - ! - ! this will be integrated numerically over r2: - ! we use grid for r1 and extra_grid for r2 - ! - ! for 99 < j1b_type < 199 - ! - ! u(r1,r2) = j12_mu(r12) x v(r1) x v(r2) - ! grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2) - ! - END_DOC - - implicit none - integer, intent(in) :: n_grid2 - double precision, intent(in) :: r1(3) - double precision, intent(out) :: resx(n_grid2), resy(n_grid2), resz(n_grid2), res(n_grid2) - - integer :: jpoint - double precision :: v1b_r1 - double precision :: grad1_v1b(3) - double precision, allocatable :: v1b_r2(:) - double precision, allocatable :: u2b_r12(:) - double precision, allocatable :: gradx1_u2b(:), grady1_u2b(:), gradz1_u2b(:) - double precision, external :: j1b_nucl - - PROVIDE j1b_type - PROVIDE final_grid_points_extra - - if( (j1b_type .eq. 100) .or. & - (j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then - - call grad1_j12_mu_r1_seq(r1, n_grid2, resx, resy, resz) - do jpoint = 1, n_points_extra_final_grid - res(jpoint) = resx(jpoint) * resx(jpoint) & - + resy(jpoint) * resy(jpoint) & - + resz(jpoint) * resz(jpoint) - enddo - - elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then - - allocate(v1b_r2(n_grid2)) - allocate(u2b_r12(n_grid2)) - allocate(gradx1_u2b(n_grid2)) - allocate(grady1_u2b(n_grid2)) - allocate(gradz1_u2b(n_grid2)) - - v1b_r1 = j1b_nucl(r1) - call grad1_j1b_nucl(r1, grad1_v1b) - - call j1b_nucl_r1_seq(n_grid2, v1b_r2) - call j12_mu_r1_seq(r1, n_grid2, u2b_r12) - call grad1_j12_mu_r1_seq(r1, n_grid2, gradx1_u2b, grady1_u2b, gradz1_u2b) - - do jpoint = 1, n_points_extra_final_grid - resx(jpoint) = (gradx1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(1)) * v1b_r2(jpoint) - resy(jpoint) = (grady1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(2)) * v1b_r2(jpoint) - resz(jpoint) = (gradz1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(3)) * v1b_r2(jpoint) - res (jpoint) = resx(jpoint) * resx(jpoint) & - + resy(jpoint) * resy(jpoint) & - + resz(jpoint) * resz(jpoint) - enddo - - deallocate(v1b_r2, u2b_r12, gradx1_u2b, grady1_u2b, gradz1_u2b) - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - return -end subroutine get_grad1_u12_withsq_r1_seq - -! --- - -double precision function j12_mu(r1, r2) - - include 'constants.include.F' - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: mu_tmp, r12 - - if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then - - r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & - + (r1(2) - r2(2)) * (r1(2) - r2(2)) & - + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) - mu_tmp = mu_erf * r12 - - j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu' - stop - - endif - - return -end function j12_mu - -! --- - -subroutine grad1_j12_mu(r1, r2, grad) - - BEGIN_DOC - ! - ! gradient of j(mu(r1,r2),r12) form of jastrow. - ! - ! if mu(r1,r2) = cst ---> j1b_type < 200 and - ! - ! d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2) - ! - ! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and - ! - ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) - ! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2) - ! - END_DOC - - include 'constants.include.F' - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision, intent(out) :: grad(3) - double precision :: dx, dy, dz, r12, tmp - - grad = 0.d0 - - if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then - - dx = r1(1) - r2(1) - dy = r1(2) - r2(2) - dz = r1(3) - r2(3) - - r12 = dsqrt(dx * dx + dy * dy + dz * dz) - if(r12 .lt. 1d-10) return - - tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 - - grad(1) = tmp * dx - grad(2) = tmp * dy - grad(3) = tmp * dz - - elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then - - double precision :: mu_val, mu_tmp, mu_der(3) - - dx = r1(1) - r2(1) - dy = r1(2) - r2(2) - dz = r1(3) - r2(3) - r12 = dsqrt(dx * dx + dy * dy + dz * dz) - - call mu_r_val_and_grad(r1, r2, mu_val, mu_der) - mu_tmp = mu_val * r12 - tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) - grad(1) = tmp * mu_der(1) - grad(2) = tmp * mu_der(2) - grad(3) = tmp * mu_der(3) - - if(r12 .lt. 1d-10) return - tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 - grad(1) = grad(1) + tmp * dx - grad(2) = grad(2) + tmp * dy - grad(3) = grad(3) + tmp * dz - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - return -end subroutine grad1_j12_mu - -! --- - -double precision function j1b_nucl(r) - - implicit none - double precision, intent(in) :: r(3) - integer :: i - double precision :: a, d, e, x, y, z - - if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then - - j1b_nucl = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - j1b_nucl = j1b_nucl - dexp(-a*dsqrt(d)) - enddo - - elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then - - j1b_nucl = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - e = 1.d0 - dexp(-a*d) - j1b_nucl = j1b_nucl * e - enddo - - elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then - - j1b_nucl = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - j1b_nucl = j1b_nucl - j1b_pen_coef(i) * dexp(-a*d) - enddo - - elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then - - j1b_nucl = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - x = r(1) - nucl_coord(i,1) - y = r(2) - nucl_coord(i,2) - z = r(3) - nucl_coord(i,3) - d = x*x + y*y + z*z - j1b_nucl = j1b_nucl - dexp(-a*d*d) - enddo - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl' - stop - - endif - - return -end function j1b_nucl - -! --- - -double precision function j1b_nucl_square(r) - - implicit none - double precision, intent(in) :: r(3) - integer :: i - double precision :: a, d, e, x, y, z - - if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then - - j1b_nucl_square = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - j1b_nucl_square = j1b_nucl_square - dexp(-a*dsqrt(d)) - enddo - j1b_nucl_square = j1b_nucl_square * j1b_nucl_square - - elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then - - j1b_nucl_square = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - e = 1.d0 - dexp(-a*d) - j1b_nucl_square = j1b_nucl_square * e - enddo - j1b_nucl_square = j1b_nucl_square * j1b_nucl_square - - elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then - - j1b_nucl_square = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - j1b_nucl_square = j1b_nucl_square - j1b_pen_coef(i) * dexp(-a*d) - enddo - j1b_nucl_square = j1b_nucl_square * j1b_nucl_square - - elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then - - j1b_nucl_square = 1.d0 - do i = 1, nucl_num - a = j1b_pen(i) - x = r(1) - nucl_coord(i,1) - y = r(2) - nucl_coord(i,2) - z = r(3) - nucl_coord(i,3) - d = x*x + y*y + z*z - j1b_nucl_square = j1b_nucl_square - dexp(-a*d*d) - enddo - j1b_nucl_square = j1b_nucl_square * j1b_nucl_square - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_square' - stop - - endif - - return -end function j1b_nucl_square - -! --- - -subroutine grad1_j1b_nucl(r, grad) - - implicit none - double precision, intent(in) :: r(3) - double precision, intent(out) :: grad(3) - integer :: ipoint, i, j, phase - double precision :: x, y, z, dx, dy, dz - double precision :: a, d, e - double precision :: fact_x, fact_y, fact_z - double precision :: ax_der, ay_der, az_der, a_expo - - if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then - - fact_x = 0.d0 - fact_y = 0.d0 - fact_z = 0.d0 - do i = 1, nucl_num - a = j1b_pen(i) - x = r(1) - nucl_coord(i,1) - y = r(2) - nucl_coord(i,2) - z = r(3) - nucl_coord(i,3) - d = dsqrt(x*x + y*y + z*z) - e = a * dexp(-a*d) / d - - fact_x += e * x - fact_y += e * y - fact_z += e * z - enddo - - grad(1) = fact_x - grad(2) = fact_y - grad(3) = fact_z - - elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then - - x = r(1) - y = r(2) - z = r(3) - - fact_x = 0.d0 - fact_y = 0.d0 - fact_z = 0.d0 - do i = 1, List_all_comb_b2_size - - phase = 0 - a_expo = 0.d0 - ax_der = 0.d0 - ay_der = 0.d0 - az_der = 0.d0 - do j = 1, nucl_num - a = dble(List_all_comb_b2(j,i)) * j1b_pen(j) - dx = x - nucl_coord(j,1) - dy = y - nucl_coord(j,2) - dz = z - nucl_coord(j,3) - - phase += List_all_comb_b2(j,i) - a_expo += a * (dx*dx + dy*dy + dz*dz) - ax_der += a * dx - ay_der += a * dy - az_der += a * dz - enddo - e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo) - - fact_x += e * ax_der - fact_y += e * ay_der - fact_z += e * az_der - enddo - - grad(1) = fact_x - grad(2) = fact_y - grad(3) = fact_z - - elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then - - fact_x = 0.d0 - fact_y = 0.d0 - fact_z = 0.d0 - do i = 1, nucl_num - a = j1b_pen(i) - x = r(1) - nucl_coord(i,1) - y = r(2) - nucl_coord(i,2) - z = r(3) - nucl_coord(i,3) - d = x*x + y*y + z*z - e = a * j1b_pen_coef(i) * dexp(-a*d) - - fact_x += e * x - fact_y += e * y - fact_z += e * z - enddo - - grad(1) = 2.d0 * fact_x - grad(2) = 2.d0 * fact_y - grad(3) = 2.d0 * fact_z - - elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then - - fact_x = 0.d0 - fact_y = 0.d0 - fact_z = 0.d0 - do i = 1, nucl_num - a = j1b_pen(i) - x = r(1) - nucl_coord(i,1) - y = r(2) - nucl_coord(i,2) - z = r(3) - nucl_coord(i,3) - d = x*x + y*y + z*z - e = a * d * dexp(-a*d*d) - - fact_x += e * x - fact_y += e * y - fact_z += e * z - enddo - - grad(1) = 4.d0 * fact_x - grad(2) = 4.d0 * fact_y - grad(3) = 4.d0 * fact_z - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented for grad1_j1b_nucl' - stop - - endif - - return -end subroutine grad1_j1b_nucl - -! --- - -subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision, intent(out) :: mu_val, mu_der(3) - double precision :: r(3) - double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1) - double precision :: dm_tot, tmp1, tmp2, tmp3 - double precision :: rho1, grad_rho1(3),rho2,rho_tot,inv_rho_tot - double precision :: f_rho1, f_rho2, d_drho_f_rho1 - double precision :: d_dx1_f_rho1(3),d_dx_rho_f_rho(3),nume - - if(j1b_type .eq. 200) then - - ! - ! r = 0.5 (r1 + r2) - ! - ! mu[rho(r)] = alpha sqrt(rho) + mu0 exp(-rho) - ! - ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx - ! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx) - ! - - PROVIDE mu_r_ct mu_erf - - r(1) = 0.5d0 * (r1(1) + r2(1)) - r(2) = 0.5d0 * (r1(2) + r2(2)) - r(3) = 0.5d0 * (r1(3) + r2(3)) - - call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b) - - dm_tot = dm_a(1) + dm_b(1) - tmp1 = dsqrt(dm_tot) - tmp2 = mu_erf * dexp(-dm_tot) - - mu_val = mu_r_ct * tmp1 + tmp2 - - mu_der = 0.d0 - if(dm_tot .lt. 1d-7) return - - tmp3 = 0.25d0 * mu_r_ct / tmp1 - 0.5d0 * tmp2 - mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1)) - mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1)) - mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1)) - - elseif(j1b_type .eq. 201) then - - ! - ! r = 0.5 (r1 + r2) - ! - ! mu[rho(r)] = alpha rho + mu0 exp(-rho) - ! - ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx - ! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx) - ! - - PROVIDE mu_r_ct mu_erf - - r(1) = 0.5d0 * (r1(1) + r2(1)) - r(2) = 0.5d0 * (r1(2) + r2(2)) - r(3) = 0.5d0 * (r1(3) + r2(3)) - - call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b) - - dm_tot = dm_a(1) + dm_b(1) - tmp2 = mu_erf * dexp(-dm_tot) - - mu_val = mu_r_ct * dm_tot + tmp2 - - tmp3 = 0.5d0 * (mu_r_ct - tmp2) - mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1)) - mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1)) - mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1)) - - elseif(j1b_type .eq. 202) then - - ! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO - ! - ! RHO = rho(r1) + rho(r2) - ! - ! f[rho] = alpha rho^beta + mu0 exp(-rho) - ! - ! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)]) - ! - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] } - ! - ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) - mu0 exp(-rho(r1))] (d rho(r1) / dx1) - ! - ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1) - - !!!!!!!!! rho1,rho2,rho1+rho2 - call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) - rho_tot = rho1 + rho2 - if(rho_tot.lt.1.d-10)rho_tot = 1.d-10 - inv_rho_tot = 1.d0/rho_tot - ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf * exp(-rho) - call get_all_f_rho(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2) - d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3) - d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3) - nume = rho1 * f_rho1 + rho2 * f_rho2 - mu_val = nume * inv_rho_tot - mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume) - elseif(j1b_type .eq. 203) then - - ! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO - ! - ! RHO = rho(r1) + rho(r2) - ! - ! f[rho] = alpha rho^beta + mu0 - ! - ! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)]) - ! - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] } - ! - ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1) - ! - ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1) - - !!!!!!!!! rho1,rho2,rho1+rho2 - call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) - rho_tot = rho1 + rho2 - if(rho_tot.lt.1.d-10)rho_tot = 1.d-10 - inv_rho_tot = 1.d0/rho_tot - ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf - call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2) - d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3) - d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3) - nume = rho1 * f_rho1 + rho2 * f_rho2 - mu_val = nume * inv_rho_tot - mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume) - elseif(j1b_type .eq. 204) then - - ! mu(r1,r2) = 1/2 * (f[rho(r1)] + f[rho(r2)]} - ! - ! f[rho] = alpha rho^beta + mu0 - ! - ! d/dx1 mu(r1,r2) = 1/2 * d/dx1 (rho(r1) f[rho(r1)]) - ! - ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1) - ! - ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1) - - !!!!!!!!! rho1,rho2,rho1+rho2 - call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) - rho_tot = rho1 + rho2 - if(rho_tot.lt.1.d-10)rho_tot = 1.d-10 - inv_rho_tot = 1.d0/rho_tot - ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf - call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2) - d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3) - d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3) - mu_val = 0.5d0 * ( f_rho1 + f_rho2) - mu_der(1:3) = d_dx_rho_f_rho(1:3) - else - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - return -end subroutine mu_r_val_and_grad - -! --- - -subroutine grad1_j1b_nucl_square_num(r1, grad) - - implicit none - double precision, intent(in) :: r1(3) - double precision, intent(out) :: grad(3) - double precision :: r(3), eps, tmp_eps, vp, vm - double precision, external :: j1b_nucl_square - - eps = 1d-5 - tmp_eps = 0.5d0 / eps - - r(1:3) = r1(1:3) - - r(1) = r(1) + eps - vp = j1b_nucl_square(r) - r(1) = r(1) - 2.d0 * eps - vm = j1b_nucl_square(r) - r(1) = r(1) + eps - grad(1) = tmp_eps * (vp - vm) - - r(2) = r(2) + eps - vp = j1b_nucl_square(r) - r(2) = r(2) - 2.d0 * eps - vm = j1b_nucl_square(r) - r(2) = r(2) + eps - grad(2) = tmp_eps * (vp - vm) - - r(3) = r(3) + eps - vp = j1b_nucl_square(r) - r(3) = r(3) - 2.d0 * eps - vm = j1b_nucl_square(r) - r(3) = r(3) + eps - grad(3) = tmp_eps * (vp - vm) - - return -end subroutine grad1_j1b_nucl_square_num - -! --- - -subroutine grad1_j12_mu_square_num(r1, r2, grad) - - include 'constants.include.F' - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision, intent(out) :: grad(3) - double precision :: r(3) - double precision :: eps, tmp_eps, vp, vm - double precision, external :: j12_mu_square - - eps = 1d-5 - tmp_eps = 0.5d0 / eps - - r(1:3) = r1(1:3) - - r(1) = r(1) + eps - vp = j12_mu_square(r, r2) - r(1) = r(1) - 2.d0 * eps - vm = j12_mu_square(r, r2) - r(1) = r(1) + eps - grad(1) = tmp_eps * (vp - vm) - - r(2) = r(2) + eps - vp = j12_mu_square(r, r2) - r(2) = r(2) - 2.d0 * eps - vm = j12_mu_square(r, r2) - r(2) = r(2) + eps - grad(2) = tmp_eps * (vp - vm) - - r(3) = r(3) + eps - vp = j12_mu_square(r, r2) - r(3) = r(3) - 2.d0 * eps - vm = j12_mu_square(r, r2) - r(3) = r(3) + eps - grad(3) = tmp_eps * (vp - vm) - - return -end subroutine grad1_j12_mu_square_num - -! --- - -double precision function j12_mu_square(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision, external :: j12_mu - - j12_mu_square = j12_mu(r1, r2) * j12_mu(r1, r2) - - return -end function j12_mu_square - -! --- - -subroutine f_mu_and_deriv_mu(rho,alpha,mu0,beta,f_mu,d_drho_f_mu) - implicit none - BEGIN_DOC -! function giving mu as a function of rho -! -! f_mu = alpha * rho**beta + mu0 * exp(-rho) -! -! and its derivative with respect to rho d_drho_f_mu - END_DOC - double precision, intent(in) :: rho,alpha,mu0,beta - double precision, intent(out) :: f_mu,d_drho_f_mu - f_mu = alpha * (rho)**beta + mu0 * dexp(-rho) - d_drho_f_mu = alpha * beta * rho**(beta-1.d0) - mu0 * dexp(-rho) - -end - - -subroutine get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) - implicit none - BEGIN_DOC -! returns the density in r1,r2 and grad_rho at r1 - END_DOC - double precision, intent(in) :: r1(3),r2(3) - double precision, intent(out):: grad_rho1(3),rho1,rho2 - double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1) - call density_and_grad_alpha_beta(r1, dm_a, dm_b, grad_dm_a, grad_dm_b) - rho1 = dm_a(1) + dm_b(1) - grad_rho1(1:3) = grad_dm_a(1:3,1) + grad_dm_b(1:3,1) - call density_and_grad_alpha_beta(r2, dm_a, dm_b, grad_dm_a, grad_dm_b) - rho2 = dm_a(1) + dm_b(1) -end - -subroutine get_all_f_rho(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2) - implicit none - BEGIN_DOC -! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1)) - END_DOC - double precision, intent(in) :: rho1,rho2,alpha,mu0,beta - double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2 - double precision :: tmp - call f_mu_and_deriv_mu(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1) - call f_mu_and_deriv_mu(rho2,alpha,mu0,beta,f_rho2,tmp) -end - - -subroutine get_all_f_rho_simple(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2) - implicit none - BEGIN_DOC -! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1)) - END_DOC - double precision, intent(in) :: rho1,rho2,alpha,mu0,beta - double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2 - double precision :: tmp - call f_mu_and_deriv_mu_simple(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1) - call f_mu_and_deriv_mu_simple(rho2,alpha,mu0,beta,f_rho2,tmp) -end - -subroutine f_mu_and_deriv_mu_simple(rho,alpha,mu0,beta,f_mu,d_drho_f_mu) - implicit none - BEGIN_DOC -! function giving mu as a function of rho -! -! f_mu = alpha * rho**beta + mu0 -! -! and its derivative with respect to rho d_drho_f_mu - END_DOC - double precision, intent(in) :: rho,alpha,mu0,beta - double precision, intent(out) :: f_mu,d_drho_f_mu - f_mu = alpha * (rho)**beta + mu0 - d_drho_f_mu = alpha * beta * rho**(beta-1.d0) - -end - -! --- - -subroutine grad1_j12_mu_r1_seq(r1, n_grid2, gradx, grady, gradz) - - BEGIN_DOC - ! - ! gradient of j(mu(r1,r2),r12) form of jastrow. - ! - ! if mu(r1,r2) = cst ---> j1b_type < 200 and - ! - ! d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2) - ! - ! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and - ! - ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) - ! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2) - ! - END_DOC - - include 'constants.include.F' - - implicit none - integer , intent(in) :: n_grid2 - double precision, intent(in) :: r1(3) - double precision, intent(out) :: gradx(n_grid2) - double precision, intent(out) :: grady(n_grid2) - double precision, intent(out) :: gradz(n_grid2) - - integer :: jpoint - double precision :: r2(3) - double precision :: dx, dy, dz, r12, tmp - - if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then - - do jpoint = 1, n_points_extra_final_grid ! r2 - - r2(1) = final_grid_points_extra(1,jpoint) - r2(2) = final_grid_points_extra(2,jpoint) - r2(3) = final_grid_points_extra(3,jpoint) - - dx = r1(1) - r2(1) - dy = r1(2) - r2(2) - dz = r1(3) - r2(3) - - r12 = dsqrt(dx * dx + dy * dy + dz * dz) - if(r12 .lt. 1d-10) return - - tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 - - gradx(jpoint) = tmp * dx - grady(jpoint) = tmp * dy - gradz(jpoint) = tmp * dz - enddo - - elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then - - double precision :: mu_val, mu_tmp, mu_der(3) - - do jpoint = 1, n_points_extra_final_grid ! r2 - - r2(1) = final_grid_points_extra(1,jpoint) - r2(2) = final_grid_points_extra(2,jpoint) - r2(3) = final_grid_points_extra(3,jpoint) - - dx = r1(1) - r2(1) - dy = r1(2) - r2(2) - dz = r1(3) - r2(3) - r12 = dsqrt(dx * dx + dy * dy + dz * dz) - - call mu_r_val_and_grad(r1, r2, mu_val, mu_der) - mu_tmp = mu_val * r12 - tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) - gradx(jpoint) = tmp * mu_der(1) - grady(jpoint) = tmp * mu_der(2) - gradz(jpoint) = tmp * mu_der(3) - - if(r12 .lt. 1d-10) return - tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 - - gradx(jpoint) = gradx(jpoint) + tmp * dx - grady(jpoint) = grady(jpoint) + tmp * dy - gradz(jpoint) = gradz(jpoint) + tmp * dz - enddo - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - return -end subroutine grad1_j12_mu_r1_seq - -! --- - -subroutine j12_mu_r1_seq(r1, n_grid2, res) - - include 'constants.include.F' - - implicit none - integer, intent(in) :: n_grid2 - double precision, intent(in) :: r1(3) - double precision, intent(out) :: res(n_grid2) - - integer :: jpoint - double precision :: r2(3) - double precision :: mu_tmp, r12 - - PROVIDE final_grid_points_extra - - if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then - - do jpoint = 1, n_points_extra_final_grid ! r2 - - r2(1) = final_grid_points_extra(1,jpoint) - r2(2) = final_grid_points_extra(2,jpoint) - r2(3) = final_grid_points_extra(3,jpoint) - - r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & - + (r1(2) - r2(2)) * (r1(2) - r2(2)) & - + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) - mu_tmp = mu_erf * r12 - - res(jpoint) = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf - enddo - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu_r1_seq' - stop - - endif - - return -end subroutine j12_mu_r1_seq - -! --- - -subroutine j1b_nucl_r1_seq(n_grid2, res) - - implicit none - integer, intent(in) :: n_grid2 - double precision, intent(out) :: res(n_grid2) - - double precision :: r(3) - integer :: i, jpoint - double precision :: a, d, e, x, y, z - - if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then - - res = 1.d0 - - do jpoint = 1, n_points_extra_final_grid ! r2 - r(1) = final_grid_points_extra(1,jpoint) - r(2) = final_grid_points_extra(2,jpoint) - r(3) = final_grid_points_extra(3,jpoint) - - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - - res(jpoint) -= dexp(-a*dsqrt(d)) - enddo - enddo - - elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then - - res = 1.d0 - - do jpoint = 1, n_points_extra_final_grid ! r2 - r(1) = final_grid_points_extra(1,jpoint) - r(2) = final_grid_points_extra(2,jpoint) - r(3) = final_grid_points_extra(3,jpoint) - - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - e = 1.d0 - dexp(-a*d) - - res(jpoint) *= e - enddo - enddo - - elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then - - res = 1.d0 - - do jpoint = 1, n_points_extra_final_grid ! r2 - r(1) = final_grid_points_extra(1,jpoint) - r(2) = final_grid_points_extra(2,jpoint) - r(3) = final_grid_points_extra(3,jpoint) - - do i = 1, nucl_num - a = j1b_pen(i) - d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & - + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & - + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) - res(jpoint) -= j1b_pen_coef(i) * dexp(-a*d) - enddo - enddo - - elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then - - res = 1.d0 - - do jpoint = 1, n_points_extra_final_grid ! r2 - r(1) = final_grid_points_extra(1,jpoint) - r(2) = final_grid_points_extra(2,jpoint) - r(3) = final_grid_points_extra(3,jpoint) - - do i = 1, nucl_num - a = j1b_pen(i) - x = r(1) - nucl_coord(i,1) - y = r(2) - nucl_coord(i,2) - z = r(3) - nucl_coord(i,3) - d = x*x + y*y + z*z - res(jpoint) -= dexp(-a*d*d) - enddo - enddo - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_r1_seq' - stop - - endif - - return -end subroutine j1b_nucl_r1_seq - -! --- - diff --git a/src/non_h_ints_mu/jast_deriv_utils.irp.f b/src/non_h_ints_mu/jast_deriv_utils.irp.f new file mode 100644 index 00000000..bcbe16af --- /dev/null +++ b/src/non_h_ints_mu/jast_deriv_utils.irp.f @@ -0,0 +1,700 @@ + +! --- + +double precision function j12_mu(r1, r2) + + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: mu_tmp, r12 + + if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then + + r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & + + (r1(2) - r2(2)) * (r1(2) - r2(2)) & + + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) + mu_tmp = mu_erf * r12 + + j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu' + stop + + endif + + return +end function j12_mu + +! --- + +subroutine grad1_j12_mu(r1, r2, grad) + + BEGIN_DOC + ! + ! gradient of j(mu(r1,r2),r12) form of jastrow. + ! + ! if mu(r1,r2) = cst ---> j1b_type < 200 and + ! + ! d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2) + ! + ! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and + ! + ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) + ! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2) + ! + END_DOC + + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out) :: grad(3) + double precision :: dx, dy, dz, r12, tmp + + grad = 0.d0 + + if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + if(r12 .lt. 1d-10) return + + tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 + + grad(1) = tmp * dx + grad(2) = tmp * dy + grad(3) = tmp * dz + + elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then + + double precision :: mu_val, mu_tmp, mu_der(3) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + + call mu_r_val_and_grad(r1, r2, mu_val, mu_der) + mu_tmp = mu_val * r12 + tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) + grad(1) = tmp * mu_der(1) + grad(2) = tmp * mu_der(2) + grad(3) = tmp * mu_der(3) + + if(r12 .lt. 1d-10) return + tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 + grad(1) = grad(1) + tmp * dx + grad(2) = grad(2) + tmp * dy + grad(3) = grad(3) + tmp * dz + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + + return +end subroutine grad1_j12_mu + +! --- + +double precision function j1b_nucl(r) + + implicit none + double precision, intent(in) :: r(3) + integer :: i + double precision :: a, d, e, x, y, z + + if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then + + j1b_nucl = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + j1b_nucl = j1b_nucl - dexp(-a*dsqrt(d)) + enddo + + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then + + j1b_nucl = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + e = 1.d0 - dexp(-a*d) + j1b_nucl = j1b_nucl * e + enddo + + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then + + j1b_nucl = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + j1b_nucl = j1b_nucl - j1b_pen_coef(i) * dexp(-a*d) + enddo + + elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then + + j1b_nucl = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = x*x + y*y + z*z + j1b_nucl = j1b_nucl - dexp(-a*d*d) + enddo + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl' + stop + + endif + + return +end function j1b_nucl + +! --- + +double precision function j1b_nucl_square(r) + + implicit none + double precision, intent(in) :: r(3) + integer :: i + double precision :: a, d, e, x, y, z + + if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + j1b_nucl_square = j1b_nucl_square - dexp(-a*dsqrt(d)) + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + e = 1.d0 - dexp(-a*d) + j1b_nucl_square = j1b_nucl_square * e + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + j1b_nucl_square = j1b_nucl_square - j1b_pen_coef(i) * dexp(-a*d) + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then + + j1b_nucl_square = 1.d0 + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = x*x + y*y + z*z + j1b_nucl_square = j1b_nucl_square - dexp(-a*d*d) + enddo + j1b_nucl_square = j1b_nucl_square * j1b_nucl_square + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_square' + stop + + endif + + return +end function j1b_nucl_square + +! --- + +subroutine grad1_j1b_nucl(r, grad) + + implicit none + double precision, intent(in) :: r(3) + double precision, intent(out) :: grad(3) + integer :: ipoint, i, j, phase + double precision :: x, y, z, dx, dy, dz + double precision :: a, d, e + double precision :: fact_x, fact_y, fact_z + double precision :: ax_der, ay_der, az_der, a_expo + + if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = dsqrt(x*x + y*y + z*z) + e = a * dexp(-a*d) / d + + fact_x += e * x + fact_y += e * y + fact_z += e * z + enddo + + grad(1) = fact_x + grad(2) = fact_y + grad(3) = fact_z + + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then + + x = r(1) + y = r(2) + z = r(3) + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + do i = 1, List_all_comb_b2_size + + phase = 0 + a_expo = 0.d0 + ax_der = 0.d0 + ay_der = 0.d0 + az_der = 0.d0 + do j = 1, nucl_num + a = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + dx = x - nucl_coord(j,1) + dy = y - nucl_coord(j,2) + dz = z - nucl_coord(j,3) + + phase += List_all_comb_b2(j,i) + a_expo += a * (dx*dx + dy*dy + dz*dz) + ax_der += a * dx + ay_der += a * dy + az_der += a * dz + enddo + e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo) + + fact_x += e * ax_der + fact_y += e * ay_der + fact_z += e * az_der + enddo + + grad(1) = fact_x + grad(2) = fact_y + grad(3) = fact_z + + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = x*x + y*y + z*z + e = a * j1b_pen_coef(i) * dexp(-a*d) + + fact_x += e * x + fact_y += e * y + fact_z += e * z + enddo + + grad(1) = 2.d0 * fact_x + grad(2) = 2.d0 * fact_y + grad(3) = 2.d0 * fact_z + + elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = x*x + y*y + z*z + e = a * d * dexp(-a*d*d) + + fact_x += e * x + fact_y += e * y + fact_z += e * z + enddo + + grad(1) = 4.d0 * fact_x + grad(2) = 4.d0 * fact_y + grad(3) = 4.d0 * fact_z + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented for grad1_j1b_nucl' + stop + + endif + + return +end subroutine grad1_j1b_nucl + +! --- + +subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out) :: mu_val, mu_der(3) + double precision :: r(3) + double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1) + double precision :: dm_tot, tmp1, tmp2, tmp3 + double precision :: rho1, grad_rho1(3),rho2,rho_tot,inv_rho_tot + double precision :: f_rho1, f_rho2, d_drho_f_rho1 + double precision :: d_dx1_f_rho1(3),d_dx_rho_f_rho(3),nume + + if(j1b_type .eq. 200) then + + ! + ! r = 0.5 (r1 + r2) + ! + ! mu[rho(r)] = alpha sqrt(rho) + mu0 exp(-rho) + ! + ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx + ! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx) + ! + + PROVIDE mu_r_ct mu_erf + + r(1) = 0.5d0 * (r1(1) + r2(1)) + r(2) = 0.5d0 * (r1(2) + r2(2)) + r(3) = 0.5d0 * (r1(3) + r2(3)) + + call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b) + + dm_tot = dm_a(1) + dm_b(1) + tmp1 = dsqrt(dm_tot) + tmp2 = mu_erf * dexp(-dm_tot) + + mu_val = mu_r_ct * tmp1 + tmp2 + + mu_der = 0.d0 + if(dm_tot .lt. 1d-7) return + + tmp3 = 0.25d0 * mu_r_ct / tmp1 - 0.5d0 * tmp2 + mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1)) + mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1)) + mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1)) + + elseif(j1b_type .eq. 201) then + + ! + ! r = 0.5 (r1 + r2) + ! + ! mu[rho(r)] = alpha rho + mu0 exp(-rho) + ! + ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx + ! d mu[rho(r)] / dx = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx) + ! + + PROVIDE mu_r_ct mu_erf + + r(1) = 0.5d0 * (r1(1) + r2(1)) + r(2) = 0.5d0 * (r1(2) + r2(2)) + r(3) = 0.5d0 * (r1(3) + r2(3)) + + call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b) + + dm_tot = dm_a(1) + dm_b(1) + tmp2 = mu_erf * dexp(-dm_tot) + + mu_val = mu_r_ct * dm_tot + tmp2 + + tmp3 = 0.5d0 * (mu_r_ct - tmp2) + mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1)) + mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1)) + mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1)) + + elseif(j1b_type .eq. 202) then + + ! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO + ! + ! RHO = rho(r1) + rho(r2) + ! + ! f[rho] = alpha rho^beta + mu0 exp(-rho) + ! + ! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)]) + ! - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] } + ! + ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) - mu0 exp(-rho(r1))] (d rho(r1) / dx1) + ! + ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1) + + !!!!!!!!! rho1,rho2,rho1+rho2 + call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) + rho_tot = rho1 + rho2 + if(rho_tot.lt.1.d-10)rho_tot = 1.d-10 + inv_rho_tot = 1.d0/rho_tot + ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf * exp(-rho) + call get_all_f_rho(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2) + d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3) + d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3) + nume = rho1 * f_rho1 + rho2 * f_rho2 + mu_val = nume * inv_rho_tot + mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume) + elseif(j1b_type .eq. 203) then + + ! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO + ! + ! RHO = rho(r1) + rho(r2) + ! + ! f[rho] = alpha rho^beta + mu0 + ! + ! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)]) + ! - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] } + ! + ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1) + ! + ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1) + + !!!!!!!!! rho1,rho2,rho1+rho2 + call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) + rho_tot = rho1 + rho2 + if(rho_tot.lt.1.d-10)rho_tot = 1.d-10 + inv_rho_tot = 1.d0/rho_tot + ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf + call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2) + d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3) + d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3) + nume = rho1 * f_rho1 + rho2 * f_rho2 + mu_val = nume * inv_rho_tot + mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume) + elseif(j1b_type .eq. 204) then + + ! mu(r1,r2) = 1/2 * (f[rho(r1)] + f[rho(r2)]} + ! + ! f[rho] = alpha rho^beta + mu0 + ! + ! d/dx1 mu(r1,r2) = 1/2 * d/dx1 (rho(r1) f[rho(r1)]) + ! + ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1) + ! + ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1) + + !!!!!!!!! rho1,rho2,rho1+rho2 + call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) + rho_tot = rho1 + rho2 + if(rho_tot.lt.1.d-10)rho_tot = 1.d-10 + inv_rho_tot = 1.d0/rho_tot + ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf + call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2) + d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3) + d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3) + mu_val = 0.5d0 * ( f_rho1 + f_rho2) + mu_der(1:3) = d_dx_rho_f_rho(1:3) + else + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + + return +end subroutine mu_r_val_and_grad + +! --- + +subroutine grad1_j1b_nucl_square_num(r1, grad) + + implicit none + double precision, intent(in) :: r1(3) + double precision, intent(out) :: grad(3) + double precision :: r(3), eps, tmp_eps, vp, vm + double precision, external :: j1b_nucl_square + + eps = 1d-5 + tmp_eps = 0.5d0 / eps + + r(1:3) = r1(1:3) + + r(1) = r(1) + eps + vp = j1b_nucl_square(r) + r(1) = r(1) - 2.d0 * eps + vm = j1b_nucl_square(r) + r(1) = r(1) + eps + grad(1) = tmp_eps * (vp - vm) + + r(2) = r(2) + eps + vp = j1b_nucl_square(r) + r(2) = r(2) - 2.d0 * eps + vm = j1b_nucl_square(r) + r(2) = r(2) + eps + grad(2) = tmp_eps * (vp - vm) + + r(3) = r(3) + eps + vp = j1b_nucl_square(r) + r(3) = r(3) - 2.d0 * eps + vm = j1b_nucl_square(r) + r(3) = r(3) + eps + grad(3) = tmp_eps * (vp - vm) + + return +end subroutine grad1_j1b_nucl_square_num + +! --- + +subroutine grad1_j12_mu_square_num(r1, r2, grad) + + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out) :: grad(3) + double precision :: r(3) + double precision :: eps, tmp_eps, vp, vm + double precision, external :: j12_mu_square + + eps = 1d-5 + tmp_eps = 0.5d0 / eps + + r(1:3) = r1(1:3) + + r(1) = r(1) + eps + vp = j12_mu_square(r, r2) + r(1) = r(1) - 2.d0 * eps + vm = j12_mu_square(r, r2) + r(1) = r(1) + eps + grad(1) = tmp_eps * (vp - vm) + + r(2) = r(2) + eps + vp = j12_mu_square(r, r2) + r(2) = r(2) - 2.d0 * eps + vm = j12_mu_square(r, r2) + r(2) = r(2) + eps + grad(2) = tmp_eps * (vp - vm) + + r(3) = r(3) + eps + vp = j12_mu_square(r, r2) + r(3) = r(3) - 2.d0 * eps + vm = j12_mu_square(r, r2) + r(3) = r(3) + eps + grad(3) = tmp_eps * (vp - vm) + + return +end subroutine grad1_j12_mu_square_num + +! --- + +double precision function j12_mu_square(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, external :: j12_mu + + j12_mu_square = j12_mu(r1, r2) * j12_mu(r1, r2) + + return +end function j12_mu_square + +! --- + +subroutine f_mu_and_deriv_mu(rho,alpha,mu0,beta,f_mu,d_drho_f_mu) + implicit none + BEGIN_DOC +! function giving mu as a function of rho +! +! f_mu = alpha * rho**beta + mu0 * exp(-rho) +! +! and its derivative with respect to rho d_drho_f_mu + END_DOC + double precision, intent(in) :: rho,alpha,mu0,beta + double precision, intent(out) :: f_mu,d_drho_f_mu + f_mu = alpha * (rho)**beta + mu0 * dexp(-rho) + d_drho_f_mu = alpha * beta * rho**(beta-1.d0) - mu0 * dexp(-rho) + +end + + +subroutine get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1) + implicit none + BEGIN_DOC +! returns the density in r1,r2 and grad_rho at r1 + END_DOC + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: grad_rho1(3),rho1,rho2 + double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1) + call density_and_grad_alpha_beta(r1, dm_a, dm_b, grad_dm_a, grad_dm_b) + rho1 = dm_a(1) + dm_b(1) + grad_rho1(1:3) = grad_dm_a(1:3,1) + grad_dm_b(1:3,1) + call density_and_grad_alpha_beta(r2, dm_a, dm_b, grad_dm_a, grad_dm_b) + rho2 = dm_a(1) + dm_b(1) +end + +subroutine get_all_f_rho(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2) + implicit none + BEGIN_DOC +! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1)) + END_DOC + double precision, intent(in) :: rho1,rho2,alpha,mu0,beta + double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2 + double precision :: tmp + call f_mu_and_deriv_mu(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1) + call f_mu_and_deriv_mu(rho2,alpha,mu0,beta,f_rho2,tmp) +end + + +subroutine get_all_f_rho_simple(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2) + implicit none + BEGIN_DOC +! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1)) + END_DOC + double precision, intent(in) :: rho1,rho2,alpha,mu0,beta + double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2 + double precision :: tmp + call f_mu_and_deriv_mu_simple(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1) + call f_mu_and_deriv_mu_simple(rho2,alpha,mu0,beta,f_rho2,tmp) +end + +subroutine f_mu_and_deriv_mu_simple(rho,alpha,mu0,beta,f_mu,d_drho_f_mu) + implicit none + BEGIN_DOC +! function giving mu as a function of rho +! +! f_mu = alpha * rho**beta + mu0 +! +! and its derivative with respect to rho d_drho_f_mu + END_DOC + double precision, intent(in) :: rho,alpha,mu0,beta + double precision, intent(out) :: f_mu,d_drho_f_mu + f_mu = alpha * (rho)**beta + mu0 + d_drho_f_mu = alpha * beta * rho**(beta-1.d0) + +end + +! --- + diff --git a/src/non_h_ints_mu/jast_deriv_utils_vect.irp.f b/src/non_h_ints_mu/jast_deriv_utils_vect.irp.f new file mode 100644 index 00000000..f9512827 --- /dev/null +++ b/src/non_h_ints_mu/jast_deriv_utils_vect.irp.f @@ -0,0 +1,332 @@ + +! --- + +subroutine get_grad1_u12_withsq_r1_seq(r1, n_grid2, resx, resy, resz, res) + + BEGIN_DOC + ! + ! grad_1 u(r1,r2) + ! + ! this will be integrated numerically over r2: + ! we use grid for r1 and extra_grid for r2 + ! + ! for 99 < j1b_type < 199 + ! + ! u(r1,r2) = j12_mu(r12) x v(r1) x v(r2) + ! grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2) + ! + END_DOC + + implicit none + integer, intent(in) :: n_grid2 + double precision, intent(in) :: r1(3) + double precision, intent(out) :: resx(n_grid2), resy(n_grid2), resz(n_grid2), res(n_grid2) + + integer :: jpoint + double precision :: v1b_r1 + double precision :: grad1_v1b(3) + double precision, allocatable :: v1b_r2(:) + double precision, allocatable :: u2b_r12(:) + double precision, allocatable :: gradx1_u2b(:), grady1_u2b(:), gradz1_u2b(:) + double precision, external :: j1b_nucl + + PROVIDE j1b_type + PROVIDE final_grid_points_extra + + if( (j1b_type .eq. 100) .or. & + (j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then + + call grad1_j12_mu_r1_seq(r1, n_grid2, resx, resy, resz) + do jpoint = 1, n_points_extra_final_grid + res(jpoint) = resx(jpoint) * resx(jpoint) & + + resy(jpoint) * resy(jpoint) & + + resz(jpoint) * resz(jpoint) + enddo + + elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then + + allocate(v1b_r2(n_grid2)) + allocate(u2b_r12(n_grid2)) + allocate(gradx1_u2b(n_grid2)) + allocate(grady1_u2b(n_grid2)) + allocate(gradz1_u2b(n_grid2)) + + v1b_r1 = j1b_nucl(r1) + call grad1_j1b_nucl(r1, grad1_v1b) + + call j1b_nucl_r1_seq(n_grid2, v1b_r2) + call j12_mu_r1_seq(r1, n_grid2, u2b_r12) + call grad1_j12_mu_r1_seq(r1, n_grid2, gradx1_u2b, grady1_u2b, gradz1_u2b) + + do jpoint = 1, n_points_extra_final_grid + resx(jpoint) = (gradx1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(1)) * v1b_r2(jpoint) + resy(jpoint) = (grady1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(2)) * v1b_r2(jpoint) + resz(jpoint) = (gradz1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(3)) * v1b_r2(jpoint) + res (jpoint) = resx(jpoint) * resx(jpoint) & + + resy(jpoint) * resy(jpoint) & + + resz(jpoint) * resz(jpoint) + enddo + + deallocate(v1b_r2, u2b_r12, gradx1_u2b, grady1_u2b, gradz1_u2b) + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + + return +end subroutine get_grad1_u12_withsq_r1_seq + +! --- + +subroutine grad1_j12_mu_r1_seq(r1, n_grid2, gradx, grady, gradz) + + BEGIN_DOC + ! + ! gradient of j(mu(r1,r2),r12) form of jastrow. + ! + ! if mu(r1,r2) = cst ---> j1b_type < 200 and + ! + ! d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2) + ! + ! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and + ! + ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) + ! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2) + ! + END_DOC + + include 'constants.include.F' + + implicit none + integer , intent(in) :: n_grid2 + double precision, intent(in) :: r1(3) + double precision, intent(out) :: gradx(n_grid2) + double precision, intent(out) :: grady(n_grid2) + double precision, intent(out) :: gradz(n_grid2) + + integer :: jpoint + double precision :: r2(3) + double precision :: dx, dy, dz, r12, tmp + + if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + if(r12 .lt. 1d-10) then + gradx(jpoint) = 0.d0 + grady(jpoint) = 0.d0 + gradz(jpoint) = 0.d0 + cycle + endif + + tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 + + gradx(jpoint) = tmp * dx + grady(jpoint) = tmp * dy + gradz(jpoint) = tmp * dz + enddo + + elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then + + double precision :: mu_val, mu_tmp, mu_der(3) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + + call mu_r_val_and_grad(r1, r2, mu_val, mu_der) + mu_tmp = mu_val * r12 + tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) + gradx(jpoint) = tmp * mu_der(1) + grady(jpoint) = tmp * mu_der(2) + gradz(jpoint) = tmp * mu_der(3) + + if(r12 .lt. 1d-10) then + gradx(jpoint) = 0.d0 + grady(jpoint) = 0.d0 + gradz(jpoint) = 0.d0 + cycle + endif + + tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 + + gradx(jpoint) = gradx(jpoint) + tmp * dx + grady(jpoint) = grady(jpoint) + tmp * dy + gradz(jpoint) = gradz(jpoint) + tmp * dz + enddo + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + + return +end subroutine grad1_j12_mu_r1_seq + +! --- + +subroutine j12_mu_r1_seq(r1, n_grid2, res) + + include 'constants.include.F' + + implicit none + integer, intent(in) :: n_grid2 + double precision, intent(in) :: r1(3) + double precision, intent(out) :: res(n_grid2) + + integer :: jpoint + double precision :: r2(3) + double precision :: mu_tmp, r12 + + PROVIDE final_grid_points_extra + + if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & + + (r1(2) - r2(2)) * (r1(2) - r2(2)) & + + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) + mu_tmp = mu_erf * r12 + + res(jpoint) = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf + enddo + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu_r1_seq' + stop + + endif + + return +end subroutine j12_mu_r1_seq + +! --- + +subroutine j1b_nucl_r1_seq(n_grid2, res) + + ! TODO + ! change loops order + + implicit none + integer, intent(in) :: n_grid2 + double precision, intent(out) :: res(n_grid2) + + double precision :: r(3) + integer :: i, jpoint + double precision :: a, d, e, x, y, z + + if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then + + res = 1.d0 + + do jpoint = 1, n_points_extra_final_grid ! r2 + r(1) = final_grid_points_extra(1,jpoint) + r(2) = final_grid_points_extra(2,jpoint) + r(3) = final_grid_points_extra(3,jpoint) + + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + + res(jpoint) -= dexp(-a*dsqrt(d)) + enddo + enddo + + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then + + res = 1.d0 + + do jpoint = 1, n_points_extra_final_grid ! r2 + r(1) = final_grid_points_extra(1,jpoint) + r(2) = final_grid_points_extra(2,jpoint) + r(3) = final_grid_points_extra(3,jpoint) + + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + e = 1.d0 - dexp(-a*d) + + res(jpoint) *= e + enddo + enddo + + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then + + res = 1.d0 + + do jpoint = 1, n_points_extra_final_grid ! r2 + r(1) = final_grid_points_extra(1,jpoint) + r(2) = final_grid_points_extra(2,jpoint) + r(3) = final_grid_points_extra(3,jpoint) + + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + res(jpoint) -= j1b_pen_coef(i) * dexp(-a*d) + enddo + enddo + + elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then + + res = 1.d0 + + do jpoint = 1, n_points_extra_final_grid ! r2 + r(1) = final_grid_points_extra(1,jpoint) + r(2) = final_grid_points_extra(2,jpoint) + r(3) = final_grid_points_extra(3,jpoint) + + do i = 1, nucl_num + a = j1b_pen(i) + x = r(1) - nucl_coord(i,1) + y = r(2) - nucl_coord(i,2) + z = r(3) - nucl_coord(i,3) + d = x*x + y*y + z*z + res(jpoint) -= dexp(-a*d*d) + enddo + enddo + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_r1_seq' + stop + + endif + + return +end subroutine j1b_nucl_r1_seq + +! --- + diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f deleted file mode 100644 index dd21a67f..00000000 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ /dev/null @@ -1,530 +0,0 @@ - -BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)] - - BEGIN_DOC - ! - ! TODO - ! combine with int2_grad1_u12_square_ao to avoid repeated calculation ? - ! - ! 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 (j1b_type .eq. 1) - ! - ! 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 (j1b_type .eq. 3) - ! - ! 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 - integer :: ipoint, i, j, m, jpoint - double precision :: time0, time1 - double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 - - print*, ' providing int2_grad1_u12_ao ...' - call wall_time(time0) - - PROVIDE j1b_type - - if(read_tc_integ) then - - open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read") - read(11) int2_grad1_u12_ao - - else - - if(j1b_type .eq. 0) then - - PROVIDE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu - - int2_grad1_u12_ao = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp1) & - !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points & - !$OMP , v_ij_erf_rk_cst_mu, x_v_ij_erf_rk_cst_mu, int2_grad1_u12_ao) - !$OMP DO SCHEDULE (static) - do ipoint = 1, n_points_final_grid - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) - do j = 1, ao_num - do i = 1, ao_num - tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) - int2_grad1_u12_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)) - int2_grad1_u12_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)) - int2_grad1_u12_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then - - PROVIDE v_1b_grad - PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b_an x_v_ij_erf_rk_cst_mu_j1b - - int2_grad1_u12_ao = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp0, tmp1, tmp2, tmp_x, tmp_y, tmp_z) & - !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points, v_1b, v_1b_grad & - !$OMP , v_ij_erf_rk_cst_mu_j1b, v_ij_u_cst_mu_j1b_an, x_v_ij_erf_rk_cst_mu_j1b, int2_grad1_u12_ao) - !$OMP DO SCHEDULE (static) - do ipoint = 1, n_points_final_grid - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) - 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 - tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) - tmp2 = v_ij_u_cst_mu_j1b_an(i,j,ipoint) - int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x - int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y - int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - FREE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b_an x_v_ij_erf_rk_cst_mu_j1b - - elseif(j1b_type .ge. 100) then - - PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra - - double precision, allocatable :: tmp(:,:,:) - allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) - tmp = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (j, i, jpoint) & - !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) - !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do i = 1, ao_num - do jpoint = 1, n_points_extra_final_grid - tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - integer :: n_blocks, n_rest, n_pass - integer :: i_blocks, i_rest, i_pass, ii - double precision, allocatable :: tmp_grad1_u12(:,:,:) - - ! n_points_final_grid = n_blocks * n_pass + n_rest - n_blocks = 8 - n_rest = int(mod(n_points_final_grid, n_blocks)) - n_pass = int((n_points_final_grid - n_rest) / n_blocks) - - if(n_pass .le. 1) then - print*, ' blocks are to large or grid is very small !' - stop - endif - - allocate(tmp_grad1_u12(n_points_extra_final_grid,n_blocks,3)) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i_pass, i_blocks, ipoint, ii, m, tmp_grad1_u12) & - !$OMP SHARED (n_pass, n_blocks, n_points_extra_final_grid, ao_num, & - !$OMP final_grid_points, tmp, int2_grad1_u12_ao) - !$OMP DO - do i_pass = 1, n_pass - ii = (i_pass-1)*n_blocks + 1 - - do i_blocks = 1, n_blocks - ipoint = ii - 1 + i_blocks ! r1 - call get_grad1_u12_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12(1,i_blocks,1) & - , tmp_grad1_u12(1,i_blocks,2) & - , tmp_grad1_u12(1,i_blocks,3) ) - enddo - - do m = 1, 3 - call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, 1.d0 & - , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid & - , 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num) - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - deallocate(tmp_grad1_u12) - - ! TODO - ! OPENMP - if(n_rest .ne. 0) then - - allocate(tmp_grad1_u12(n_points_extra_final_grid,n_rest,3)) - - ii = n_pass*n_blocks + 1 - do i_rest = 1, n_rest - ipoint = ii - 1 + i_rest ! r1 - call get_grad1_u12_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12(1,i_rest,1) & - , tmp_grad1_u12(1,i_rest,2) & - , tmp_grad1_u12(1,i_rest,3) ) - enddo - - do m = 1, 3 - call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, 1.d0 & - , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid & - , 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num) - enddo - - deallocate(tmp_grad1_u12) - endif - - deallocate(tmp) - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - endif - - if(write_tc_integ.and.mpi_master) then - open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write") - call ezfio_set_work_empty(.False.) - write(11) int2_grad1_u12_ao - close(11) - call ezfio_set_tc_keywords_io_tc_integ('Read') - endif - - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_ao =', time1-time0 - call print_memory_usage() - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_num_1shot, (ao_num, ao_num, n_points_final_grid, 3)] - - BEGIN_DOC - ! - ! int2_grad1_u12_ao_num_1shot(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) - ! - END_DOC - - implicit none - integer :: ipoint, i, j, m, jpoint - double precision :: time0, time1 - double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 - - print*, ' providing int2_grad1_u12_ao_num_1shot ...' - call wall_time(time0) - - PROVIDE j1b_type - - if(j1b_type .ge. 100) then - - PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra - PROVIDE grad1_u12_num - - double precision, allocatable :: tmp(:,:,:) - allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) - tmp = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (j, i, jpoint) & - !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) - !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do i = 1, ao_num - do jpoint = 1, n_points_extra_final_grid - tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - do m = 1, 3 - !call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -1.d0 & - ! this work also because of the symmetry in K(1,2) and sign compensation in L(1,2,3) - call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0 & - , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid & - , 0.d0, int2_grad1_u12_ao_num_1shot(1,1,1,m), ao_num*ao_num) - enddo - - deallocate(tmp) - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_ao_num_1shot =', time1-time0 - call print_memory_usage() - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_points_final_grid)] - - BEGIN_DOC - ! - ! int2_grad1_u12_square_ao = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2 - ! - END_DOC - - implicit none - integer :: ipoint, i, j, m, jpoint - double precision :: time0, time1 - double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 - - print*, ' providing int2_grad1_u12_square_ao ...' - call wall_time(time0) - - PROVIDE j1b_type - - if(j1b_type .eq. 0) then - - PROVIDE int2_grad1u2_grad2u2 - - int2_grad1_u12_square_ao = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, ipoint) & - !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, int2_grad1u2_grad2u2) - !$OMP DO SCHEDULE (static) - do ipoint = 1, n_points_final_grid - do j = 1, ao_num - do i = 1, ao_num - int2_grad1_u12_square_ao(i,j,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then - - if(use_ipp) then - - ! the term u12_grad1_u12_j1b_grad1_j1b is added directly for performance - PROVIDE u12sq_j1bsq grad12_j12 - - int2_grad1_u12_square_ao = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, ipoint) & - !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12) - !$OMP DO SCHEDULE (static) - do ipoint = 1, n_points_final_grid - do j = 1, ao_num - do i = 1, ao_num - int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - FREE u12sq_j1bsq grad12_j12 - - else - - PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 - - int2_grad1_u12_square_ao = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, ipoint) & - !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12, u12_grad1_u12_j1b_grad1_j1b) - !$OMP DO SCHEDULE (static) - do ipoint = 1, n_points_final_grid - do j = 1, ao_num - do i = 1, ao_num - int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - FREE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 - - endif - - elseif(j1b_type .ge. 100) then - - PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra - - double precision, allocatable :: tmp(:,:,:) - allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (j, i, jpoint) & - !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) - !$OMP DO COLLAPSE(2) - do j = 1, ao_num - do i = 1, ao_num - do jpoint = 1, n_points_extra_final_grid - tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - integer :: n_blocks, n_rest, n_pass - integer :: i_blocks, i_rest, i_pass, ii - double precision, allocatable :: tmp_grad1_u12_squared(:,:) - - ! n_points_final_grid = n_blocks * n_pass + n_rest - n_blocks = 16 - n_rest = int(mod(n_points_final_grid, n_blocks)) - n_pass = int((n_points_final_grid - n_rest) / n_blocks) - - if(n_pass .le. 1) then - print*, ' blocks are to large or grid is very small !' - stop - endif - - allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_blocks)) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i_pass, i_blocks, ipoint, ii, tmp_grad1_u12_squared) & - !$OMP SHARED (n_pass, n_blocks, n_points_extra_final_grid, ao_num, & - !$OMP final_grid_points, tmp, int2_grad1_u12_square_ao) - !$OMP DO - do i_pass = 1, n_pass - ii = (i_pass-1)*n_blocks + 1 - - do i_blocks = 1, n_blocks - ipoint = ii - 1 + i_blocks ! r1 - call get_grad1_u12_squared_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12_squared(1,i_blocks)) - enddo - - call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, -0.5d0 & - , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid & - , 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num) - enddo - !$OMP END DO - !$OMP END PARALLEL - - deallocate(tmp_grad1_u12_squared) - - ! TODO - ! OPENMP - if(n_rest .ne. 0) then - - allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_rest)) - - ii = n_pass*n_blocks + 1 - do i_rest = 1, n_rest - ipoint = ii - 1 + i_rest ! r1 - call get_grad1_u12_squared_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12_squared(1,i_rest)) - enddo - - call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, -0.5d0 & - , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid & - , 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num) - - deallocate(tmp_grad1_u12_squared) - endif - - deallocate(tmp) - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_square_ao =', time1-time0 - call print_memory_usage() - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao_num_1shot, (ao_num, ao_num, n_points_final_grid)] - - BEGIN_DOC - ! - ! int2_grad1_u12_square_ao_num_1shot = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2 - ! - END_DOC - - implicit none - integer :: i, j, jpoint - double precision :: time0, time1 - - print*, ' providing int2_grad1_u12_square_ao_num_1shot ...' - call wall_time(time0) - - PROVIDE j1b_type - - if(j1b_type .ge. 100) then - - PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra - PROVIDE grad1_u12_squared_num - - double precision, allocatable :: tmp(:,:,:) - allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (j, i, jpoint) & - !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) - !$OMP DO COLLAPSE(2) - do j = 1, ao_num - do i = 1, ao_num - do jpoint = 1, n_points_extra_final_grid - tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 & - , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid & - , 0.d0, int2_grad1_u12_square_ao_num_1shot(1,1,1), ao_num*ao_num) - - FREE grad1_u12_squared_num - deallocate(tmp) - - else - - print *, ' j1b_type = ', j1b_type, 'not implemented yet' - stop - - endif - - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_square_ao_num_1shot =', time1-time0 - call print_memory_usage() - -END_PROVIDER - -! --- - diff --git a/src/non_h_ints_mu/tc_integ_an.irp.f b/src/non_h_ints_mu/tc_integ_an.irp.f new file mode 100644 index 00000000..e2f181a0 --- /dev/null +++ b/src/non_h_ints_mu/tc_integ_an.irp.f @@ -0,0 +1,244 @@ + +BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)] + + BEGIN_DOC + ! + ! TODO + ! combine with int2_grad1_u12_square_ao to avoid repeated calculation ? + ! + ! 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 (j1b_type .eq. 1) + ! + ! 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 (j1b_type .eq. 3) + ! + ! 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 + integer :: ipoint, i, j, m, jpoint + double precision :: time0, time1 + double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 + + print*, ' providing int2_grad1_u12_ao ...' + call wall_time(time0) + + PROVIDE j1b_type + + if(read_tc_integ) then + + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read") + read(11) int2_grad1_u12_ao + + else + + if(j1b_type .eq. 0) then + + PROVIDE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu + + int2_grad1_u12_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp1) & + !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points & + !$OMP , v_ij_erf_rk_cst_mu, x_v_ij_erf_rk_cst_mu, int2_grad1_u12_ao) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) + int2_grad1_u12_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1)) + int2_grad1_u12_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2)) + int2_grad1_u12_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3)) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then + + PROVIDE v_1b_grad + PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b_an x_v_ij_erf_rk_cst_mu_j1b + + int2_grad1_u12_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp0, tmp1, tmp2, tmp_x, tmp_y, tmp_z) & + !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points, v_1b, v_1b_grad & + !$OMP , v_ij_erf_rk_cst_mu_j1b, v_ij_u_cst_mu_j1b_an, x_v_ij_erf_rk_cst_mu_j1b, int2_grad1_u12_ao) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + 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 + tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) + tmp2 = v_ij_u_cst_mu_j1b_an(i,j,ipoint) + int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x + int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y + int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + FREE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b_an x_v_ij_erf_rk_cst_mu_j1b + + elseif(j1b_type .ge. 100) then + + PROVIDE int2_grad1_u12_ao_num + int2_grad1_u12_ao = int2_grad1_u12_ao_num + FREE int2_grad1_u12_ao_num + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + endif + + if(write_tc_integ.and.mpi_master) then + open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write") + call ezfio_set_work_empty(.False.) + write(11) int2_grad1_u12_ao + close(11) + call ezfio_set_tc_keywords_io_tc_integ('Read') + endif + + call wall_time(time1) + print*, ' wall time for int2_grad1_u12_ao =', time1-time0 + call print_memory_usage() + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int2_grad1_u12_square_ao = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2 + ! + END_DOC + + implicit none + integer :: ipoint, i, j, m, jpoint + double precision :: time0, time1 + double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 + + print*, ' providing int2_grad1_u12_square_ao ...' + call wall_time(time0) + + PROVIDE j1b_type + + if(j1b_type .eq. 0) then + + PROVIDE int2_grad1u2_grad2u2 + + int2_grad1_u12_square_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, ipoint) & + !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, int2_grad1u2_grad2u2) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + int2_grad1_u12_square_ao(i,j,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then + + if(use_ipp) then + + ! the term u12_grad1_u12_j1b_grad1_j1b is added directly for performance + PROVIDE u12sq_j1bsq grad12_j12 + + int2_grad1_u12_square_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, ipoint) & + !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + FREE u12sq_j1bsq grad12_j12 + + else + + PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 + + int2_grad1_u12_square_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, ipoint) & + !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12, u12_grad1_u12_j1b_grad1_j1b) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + FREE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 + + endif + + elseif(j1b_type .ge. 100) then + + PROVIDE int2_grad1_u12_square_ao_num + int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num + FREE int2_grad1_u12_square_ao_num + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + + call wall_time(time1) + print*, ' wall time for int2_grad1_u12_square_ao =', time1-time0 + call print_memory_usage() + +END_PROVIDER + +! --- + diff --git a/src/non_h_ints_mu/tc_integ_num.irp.f b/src/non_h_ints_mu/tc_integ_num.irp.f new file mode 100644 index 00000000..707d484c --- /dev/null +++ b/src/non_h_ints_mu/tc_integ_num.irp.f @@ -0,0 +1,190 @@ + + BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_num , (ao_num,ao_num,n_points_final_grid,3)] +&BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao_num, (ao_num,ao_num,n_points_final_grid) ] + + BEGIN_DOC + ! + ! int2_grad1_u12_ao_num(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) + ! + ! int2_grad1_u12_square_ao_num = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2 + ! + END_DOC + + implicit none + integer :: ipoint, i, j, m, jpoint + integer :: n_blocks, n_rest, n_pass + integer :: i_blocks, i_rest, i_pass, ii + double precision :: time0, time1 + double precision, allocatable :: tmp(:,:,:) + double precision, allocatable :: tmp_grad1_u12(:,:,:), tmp_grad1_u12_squared(:,:) + + ! TODO + ! tmp_grad1_u12_squared get be obtained from tmp_grad1_u12 + + print*, ' providing int2_grad1_u12_ao_num & int2_grad1_u12_square_ao_num ...' + call wall_time(time0) + + PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra + + allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, i, jpoint) & + !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) + !$OMP DO SCHEDULE (static) + do j = 1, ao_num + do i = 1, ao_num + do jpoint = 1, n_points_extra_final_grid + tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + ! n_points_final_grid = n_blocks * n_pass + n_rest + n_blocks = 4 + n_rest = int(mod(n_points_final_grid, n_blocks)) + n_pass = int((n_points_final_grid - n_rest) / n_blocks) + + if(n_pass .le. 1) then + print*, ' blocks are to large or grid is very small !' + stop + endif + + allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_blocks)) + allocate(tmp_grad1_u12(n_points_extra_final_grid,n_blocks,3)) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i_pass, i_blocks, ipoint, ii, m, tmp_grad1_u12, & + !$OMP tmp_grad1_u12_squared) & + !$OMP SHARED (n_pass, n_blocks, n_points_extra_final_grid, ao_num, & + !$OMP final_grid_points, tmp, int2_grad1_u12_ao_num, & + !$OMP int2_grad1_u12_square_ao_num) + !$OMP DO + do i_pass = 1, n_pass + ii = (i_pass-1)*n_blocks + 1 + + do i_blocks = 1, n_blocks + ipoint = ii - 1 + i_blocks ! r1 + call get_grad1_u12_withsq_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12(1,i_blocks,1) & + , tmp_grad1_u12(1,i_blocks,2) & + , tmp_grad1_u12(1,i_blocks,3) & + , tmp_grad1_u12_squared(1,i_blocks)) + enddo + + do m = 1, 3 + call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, 1.d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_ao_num(1,1,ii,m), ao_num*ao_num) + enddo + call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, -0.5d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_square_ao_num(1,1,ii), ao_num*ao_num) + enddo + !$OMP END DO + !$OMP END PARALLEL + + deallocate(tmp_grad1_u12, tmp_grad1_u12_squared) + + ! TODO + ! OPENMP + if(n_rest .ne. 0) then + + allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_rest)) + allocate(tmp_grad1_u12(n_points_extra_final_grid,n_rest,3)) + + ii = n_pass*n_blocks + 1 + do i_rest = 1, n_rest + ipoint = ii - 1 + i_rest ! r1 + call get_grad1_u12_withsq_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12(1,i_rest,1) & + , tmp_grad1_u12(1,i_rest,2) & + , tmp_grad1_u12(1,i_rest,3) & + , tmp_grad1_u12_squared(1,i_rest)) + enddo + + do m = 1, 3 + call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, 1.d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_ao_num(1,1,ii,m), ao_num*ao_num) + enddo + call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, -0.5d0 & + , tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_square_ao_num(1,1,ii), ao_num*ao_num) + + deallocate(tmp_grad1_u12, tmp_grad1_u12_squared) + endif + + deallocate(tmp) + + call wall_time(time1) + print*, ' wall time for int2_grad1_u12_ao_num & int2_grad1_u12_square_ao_num =', time1-time0 + call print_memory_usage() + +END_PROVIDER + +! --- + + BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_num_1shot , (ao_num,ao_num,n_points_final_grid,3)] +&BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao_num_1shot, (ao_num,ao_num,n_points_final_grid) ] + + BEGIN_DOC + ! + ! int2_grad1_u12_ao_num_1shot(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) + ! + ! int2_grad1_u12_square_ao_num_1shot = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2 + ! + END_DOC + + implicit none + integer :: ipoint, i, j, m, jpoint + double precision :: time0, time1 + double precision, allocatable :: tmp(:,:,:) + + print*, ' providing int2_grad1_u12_ao_num_1shot & int2_grad1_u12_square_ao_num_1shot ...' + call wall_time(time0) + + PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra + PROVIDE grad1_u12_num grad1_u12_squared_num + + allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, i, jpoint) & + !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) + !$OMP DO SCHEDULE (static) + do j = 1, ao_num + do i = 1, ao_num + do jpoint = 1, n_points_extra_final_grid + tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do m = 1, 3 + !call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -1.d0 & + ! this work also because of the symmetry in K(1,2) and sign compensation in L(1,2,3) + call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0 & + , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_ao_num_1shot(1,1,1,m), ao_num*ao_num) + enddo + FREE grad1_u12_num + + call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 & + , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_square_ao_num_1shot(1,1,1), ao_num*ao_num) + FREE grad1_u12_squared_num + + deallocate(tmp) + + call wall_time(time1) + print*, ' wall time for int2_grad1_u12_ao_num_1shot & int2_grad1_u12_square_ao_num_1shot =', time1-time0 + call print_memory_usage() + +END_PROVIDER + +! --- + diff --git a/src/non_h_ints_mu/test_non_h_ints.irp.f b/src/non_h_ints_mu/test_non_h_ints.irp.f index 7d84e73b..84674fa0 100644 --- a/src/non_h_ints_mu/test_non_h_ints.irp.f +++ b/src/non_h_ints_mu/test_non_h_ints.irp.f @@ -618,13 +618,14 @@ subroutine test_int2_grad1_u12_square_ao() I_old = int2_grad1_u12_square_ao_num_1shot(j,i,ipoint) I_new = int2_grad1_u12_square_ao (j,i,ipoint) + !I_new = int2_grad1_u12_square_ao_num (j,i,ipoint) diff = dabs(I_new-I_old) if(diff .gt. thr) then print *, ' problem on:', j, i, ipoint print *, ' old value :', I_old print *, ' new value :', I_new - stop + !stop endif accu += diff @@ -660,13 +661,14 @@ subroutine test_int2_grad1_u12_ao() do m = 1, 3 I_old = int2_grad1_u12_ao_num_1shot(j,i,ipoint,m) I_new = int2_grad1_u12_ao (j,i,ipoint,m) + !I_new = int2_grad1_u12_ao_num (j,i,ipoint,m) diff = dabs(I_new-I_old) if(diff .gt. thr) then print *, ' problem on:', j, i, ipoint, m print *, ' old value :', I_old print *, ' new value :', I_new - stop + !stop endif accu += diff diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 29c238e1..2d244586 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -274,4 +274,15 @@ doc: size of radial grid over r1 interface: ezfio,provider,ocaml default: 30 +[tc_grid2_a] +type: integer +doc: size of angular grid over r2 +interface: ezfio,provider,ocaml +default: 50 + +[tc_grid2_r] +type: integer +doc: size of radial grid over r2 +interface: ezfio,provider,ocaml +default: 30