mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-23 12:55:37 +01:00
num integ grad & grad squared
This commit is contained in:
AbdAmmar
parent
38f8b96d41
commit
1f56b5d0f4
@ -33,6 +33,10 @@ doc: Number of angular grid points given from input. Warning, this number cannot
|
|||||||
interface: ezfio,provider,ocaml
|
interface: ezfio,provider,ocaml
|
||||||
default: 1202
|
default: 1202
|
||||||
|
|
||||||
|
[n_points_extra_final_grid]
|
||||||
|
type: integer
|
||||||
|
doc: Total number of extra_grid points
|
||||||
|
interface: ezfio
|
||||||
|
|
||||||
[extra_grid_type_sgn]
|
[extra_grid_type_sgn]
|
||||||
type: integer
|
type: integer
|
||||||
|
|||||||
@ -14,7 +14,7 @@
|
|||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
if(.not.my_extra_grid_becke)then
|
if(.not. my_extra_grid_becke) then
|
||||||
select case (extra_grid_type_sgn)
|
select case (extra_grid_type_sgn)
|
||||||
case(0)
|
case(0)
|
||||||
n_points_extra_radial_grid = 23
|
n_points_extra_radial_grid = 23
|
||||||
@ -33,7 +33,7 @@
|
|||||||
stop
|
stop
|
||||||
end select
|
end select
|
||||||
else
|
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
|
n_points_extra_integration_angular = my_n_pt_a_extra_grid
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
|||||||
@ -23,29 +23,33 @@ BEGIN_PROVIDER [integer, n_points_extra_final_grid]
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
print*,'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)
|
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)
|
call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, final_grid_points_extra, (3,n_points_extra_final_grid)]
|
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 [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, (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_PROVIDER [integer, index_final_points_extra_reverse, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num)]
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! final_grid_points_extra(1:3,j) = (/ x, y, z /) of the jth grid point
|
! 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
|
! 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(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
|
! 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
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
integer :: i,j,k,l,i_count
|
integer :: i,j,k,l,i_count
|
||||||
double precision :: r(3)
|
double precision :: r(3)
|
||||||
|
|
||||||
i_count = 0
|
i_count = 0
|
||||||
do j = 1, nucl_num
|
do j = 1, nucl_num
|
||||||
do i = 1, n_points_extra_radial_grid -1
|
do i = 1, n_points_extra_radial_grid -1
|
||||||
@ -67,3 +71,5 @@ END_PROVIDER
|
|||||||
enddo
|
enddo
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@ -14,7 +14,7 @@
|
|||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
if(.not.my_grid_becke)then
|
if(.not. my_grid_becke) then
|
||||||
select case (grid_type_sgn)
|
select case (grid_type_sgn)
|
||||||
case(0)
|
case(0)
|
||||||
n_points_radial_grid = 23
|
n_points_radial_grid = 23
|
||||||
|
|||||||
@ -13,17 +13,27 @@ program debug_fit
|
|||||||
|
|
||||||
PROVIDE mu_erf j1b_pen
|
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_j1b_nucl()
|
||||||
!call test_grad_j1b_nucl()
|
!call test_grad_j1b_nucl()
|
||||||
!call test_lapl_j1b_nucl()
|
!call test_lapl_j1b_nucl()
|
||||||
|
|
||||||
!call test_list_b2()
|
!call test_list_b2()
|
||||||
call test_list_b3()
|
!call test_list_b3()
|
||||||
|
|
||||||
!call test_fit_u()
|
!call test_fit_u()
|
||||||
!call test_fit_u2()
|
!call test_fit_u2()
|
||||||
!call test_fit_ugradu()
|
!call test_fit_ugradu()
|
||||||
|
|
||||||
|
call test_grad1_u12_withsq_num()
|
||||||
|
|
||||||
end
|
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
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
File diff suppressed because it is too large
Load Diff
700
src/non_h_ints_mu/jast_deriv_utils.irp.f
Normal file
700
src/non_h_ints_mu/jast_deriv_utils.irp.f
Normal file
@ -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
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
332
src/non_h_ints_mu/jast_deriv_utils_vect.irp.f
Normal file
332
src/non_h_ints_mu/jast_deriv_utils_vect.irp.f
Normal file
@ -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
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
@ -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
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
244
src/non_h_ints_mu/tc_integ_an.irp.f
Normal file
244
src/non_h_ints_mu/tc_integ_an.irp.f
Normal file
@ -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
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
190
src/non_h_ints_mu/tc_integ_num.irp.f
Normal file
190
src/non_h_ints_mu/tc_integ_num.irp.f
Normal file
@ -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
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
@ -618,13 +618,14 @@ subroutine test_int2_grad1_u12_square_ao()
|
|||||||
|
|
||||||
I_old = int2_grad1_u12_square_ao_num_1shot(j,i,ipoint)
|
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 (j,i,ipoint)
|
||||||
|
!I_new = int2_grad1_u12_square_ao_num (j,i,ipoint)
|
||||||
|
|
||||||
diff = dabs(I_new-I_old)
|
diff = dabs(I_new-I_old)
|
||||||
if(diff .gt. thr) then
|
if(diff .gt. thr) then
|
||||||
print *, ' problem on:', j, i, ipoint
|
print *, ' problem on:', j, i, ipoint
|
||||||
print *, ' old value :', I_old
|
print *, ' old value :', I_old
|
||||||
print *, ' new value :', I_new
|
print *, ' new value :', I_new
|
||||||
stop
|
!stop
|
||||||
endif
|
endif
|
||||||
|
|
||||||
accu += diff
|
accu += diff
|
||||||
@ -660,13 +661,14 @@ subroutine test_int2_grad1_u12_ao()
|
|||||||
do m = 1, 3
|
do m = 1, 3
|
||||||
I_old = int2_grad1_u12_ao_num_1shot(j,i,ipoint,m)
|
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 (j,i,ipoint,m)
|
||||||
|
!I_new = int2_grad1_u12_ao_num (j,i,ipoint,m)
|
||||||
|
|
||||||
diff = dabs(I_new-I_old)
|
diff = dabs(I_new-I_old)
|
||||||
if(diff .gt. thr) then
|
if(diff .gt. thr) then
|
||||||
print *, ' problem on:', j, i, ipoint, m
|
print *, ' problem on:', j, i, ipoint, m
|
||||||
print *, ' old value :', I_old
|
print *, ' old value :', I_old
|
||||||
print *, ' new value :', I_new
|
print *, ' new value :', I_new
|
||||||
stop
|
!stop
|
||||||
endif
|
endif
|
||||||
|
|
||||||
accu += diff
|
accu += diff
|
||||||
|
|||||||
@ -274,4 +274,15 @@ doc: size of radial grid over r1
|
|||||||
interface: ezfio,provider,ocaml
|
interface: ezfio,provider,ocaml
|
||||||
default: 30
|
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
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user