9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-10-05 07:35:59 +02:00
qp2/plugins/local/non_h_ints_mu/jast_deriv.irp.f
2024-01-16 19:07:20 +01:00

236 lines
7.4 KiB
Fortran

! ---
BEGIN_PROVIDER [double precision, grad1_u12_num, (n_points_extra_final_grid, n_points_final_grid, 3)]
&BEGIN_PROVIDER [double precision, grad1_u12_squared_num, (n_points_extra_final_grid, n_points_final_grid)]
BEGIN_DOC
!
!
! grad_1 u(r1,r2)
! numerical integration over r1 & r2
!
END_DOC
implicit none
integer :: ipoint, jpoint
double precision :: r1(3), r2(3)
double precision :: v_r1, v_r2, u2b_r12
double precision :: grad1_v(3), grad1_u2b(3)
double precision :: dx, dy, dz
double precision :: time0, time1
double precision, external :: j12_mu, env_nucl
PROVIDE env_type
PROVIDE final_grid_points_extra
print*, ' providing grad1_u12_num & grad1_u12_squared_num ...'
call wall_time(time0)
grad1_u12_num = 0.d0
grad1_u12_squared_num = 0.d0
if( ((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) .or. &
(j2e_type .eq. "Mur") ) then
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, jpoint, r1, r2, v_r1, v_r2, u2b_r12, grad1_v, grad1_u2b, dx, dy, dz) &
!$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points, &
!$OMP final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num)
!$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid ! r1
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
do jpoint = 1, n_points_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(r2, r1, grad1_u2b)
dx = grad1_u2b(1)
dy = grad1_u2b(2)
dz = grad1_u2b(3)
grad1_u12_num(jpoint,ipoint,1) = dx
grad1_u12_num(jpoint,ipoint,2) = dy
grad1_u12_num(jpoint,ipoint,3) = dz
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
elseif((j2e_type .eq. "Mu") .and. (env_type .ne. "None")) then
PROVIDE final_grid_points
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, jpoint, r1, r2, v_r1, v_r2, u2b_r12, grad1_v, grad1_u2b, dx, dy, dz) &
!$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points, &
!$OMP final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num)
!$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid ! r1
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
v_r1 = env_nucl(r1)
call grad1_env_nucl(r1, grad1_v)
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)
v_r2 = env_nucl(r2)
u2b_r12 = j12_mu(r1, r2)
call grad1_j12_mu(r2, r1, grad1_u2b)
dx = (grad1_u2b(1) * v_r1 + u2b_r12 * grad1_v(1)) * v_r2
dy = (grad1_u2b(2) * v_r1 + u2b_r12 * grad1_v(2)) * v_r2
dz = (grad1_u2b(3) * v_r1 + u2b_r12 * grad1_v(3)) * v_r2
grad1_u12_num(jpoint,ipoint,1) = dx
grad1_u12_num(jpoint,ipoint,2) = dy
grad1_u12_num(jpoint,ipoint,3) = dz
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
elseif(j2e_type .eq. "Qmckl") then
double precision :: f
f = 1.d0 / dble(elec_num - 1)
integer*8 :: n_points, n_points_max, k
integer :: ipoint_block, ipoint_end
n_points_max = n_points_extra_final_grid * n_points_final_grid
n_points = 100_8*n_points_extra_final_grid
double precision, allocatable :: rij(:,:,:)
allocate( rij(3, 2, n_points) )
use qmckl
integer(qmckl_exit_code) :: rc
double precision, allocatable :: gl(:,:,:)
allocate( gl(2,4,n_points) )
do ipoint_block = 1, n_points_final_grid, 100 ! r1
ipoint_end = min(n_points_final_grid, ipoint_block+99)
k=0
do ipoint = ipoint_block, ipoint_end
do jpoint = 1, n_points_extra_final_grid ! r2
k=k+1
rij(1:3, 1, k) = final_grid_points (1:3, ipoint)
rij(1:3, 2, k) = final_grid_points_extra(1:3, jpoint)
end do
enddo
rc = qmckl_set_electron_coord(qmckl_ctx_jastrow, 'N', n_points, rij, n_points*6_8)
if (rc /= QMCKL_SUCCESS) then
print *, irp_here, 'qmckl error in set_electron_coord'
rc = qmckl_check(qmckl_ctx_jastrow, rc)
stop -1
endif
! ---
! e-e term
rc = qmckl_get_jastrow_champ_factor_ee_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
if (rc /= QMCKL_SUCCESS) then
print *, irp_here, ' qmckl error in fact_ee_gl'
rc = qmckl_check(qmckl_ctx_jastrow, rc)
stop -1
endif
k=0
do ipoint = ipoint_block, ipoint_end
do jpoint = 1, n_points_extra_final_grid ! r2
k=k+1
grad1_u12_num(jpoint,ipoint,1) = gl(1,1,k)
grad1_u12_num(jpoint,ipoint,2) = gl(1,2,k)
grad1_u12_num(jpoint,ipoint,3) = gl(1,3,k)
enddo
enddo
! ---
! e-e-n term
! rc = qmckl_get_jastrow_champ_factor_een_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
! if (rc /= QMCKL_SUCCESS) then
! print *, irp_here, 'qmckl error in fact_een_gl'
! rc = qmckl_check(qmckl_ctx_jastrow, rc)
! stop -1
! endif
!
! k=0
! do ipoint = 1, n_points_final_grid ! r1
! do jpoint = 1, n_points_extra_final_grid ! r2
! k=k+1
! grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + gl(1,1,k)
! grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + gl(1,2,k)
! grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + gl(1,3,k)
! enddo
! enddo
! ---
! e-n term
rc = qmckl_get_jastrow_champ_factor_en_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
if (rc /= QMCKL_SUCCESS) then
print *, irp_here, 'qmckl error in fact_en_gl'
rc = qmckl_check(qmckl_ctx_jastrow, rc)
stop -1
endif
k=0
do ipoint = ipoint_block, ipoint_end ! r1
do jpoint = 1, n_points_extra_final_grid ! r2
k = k+1
grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + f * gl(1,1,k)
grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + f * gl(1,2,k)
grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + f * gl(1,3,k)
dx = grad1_u12_num(jpoint,ipoint,1)
dy = grad1_u12_num(jpoint,ipoint,2)
dz = grad1_u12_num(jpoint,ipoint,3)
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
enddo
enddo
enddo !ipoint_block
deallocate(gl, rij)
else
print *, ' Error in grad1_u12_num & grad1_u12_squared_num: Unknown Jastrow'
stop
endif ! j2e_type
call wall_time(time1)
print*, ' Wall time for grad1_u12_num & grad1_u12_squared_num (min) = ', (time1-time0)/60.d0
END_PROVIDER
! ---