mirror of
https://github.com/LCPQ/quantum_package
synced 2024-10-19 22:41:48 +02:00
added dft_utils_one_body
This commit is contained in:
parent
bc39df1a01
commit
430cbc91e2
6
src/dft_utils_on_grid/NEED
Normal file
6
src/dft_utils_on_grid/NEED
Normal file
@ -0,0 +1,6 @@
|
||||
dft_keywords
|
||||
determinants
|
||||
ao_basis
|
||||
mo_basis
|
||||
becke_numerical_grid
|
||||
density_for_dft
|
77
src/dft_utils_on_grid/ao_on_grid.irp.f
Normal file
77
src/dft_utils_on_grid/ao_on_grid.irp.f
Normal file
@ -0,0 +1,77 @@
|
||||
BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)]
|
||||
&BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! aos_in_r_array(i,j) = value of the ith ao on the jth grid point
|
||||
!
|
||||
! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point
|
||||
END_DOC
|
||||
integer :: i,j
|
||||
double precision :: aos_array(ao_num), r(3)
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
call give_all_aos_at_r(r,aos_array)
|
||||
do j = 1, ao_num
|
||||
aos_in_r_array(j,i) = aos_array(j)
|
||||
aos_in_r_array_transp(i,j) = aos_array(j)
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
|
||||
&BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp, (n_points_final_grid,ao_num,3)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! aos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith ao on the jth grid point
|
||||
!
|
||||
! aos_grad_in_r_array_transp(i,j,k) = value of the kth component of the gradient of jth ao on the ith grid point
|
||||
!
|
||||
! k = 1 : x, k= 2, y, k 3, z
|
||||
END_DOC
|
||||
integer :: i,j,m
|
||||
double precision :: aos_array(ao_num), r(3)
|
||||
double precision :: aos_grad_array(3,ao_num)
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
|
||||
do m = 1, 3
|
||||
do j = 1, ao_num
|
||||
aos_grad_in_r_array(j,i,m) = aos_grad_array(m,j)
|
||||
aos_grad_in_r_array_transp(i,j,m) = aos_grad_array(m,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)]
|
||||
&BEGIN_PROVIDER[double precision, aos_lapl_in_r_array_transp, (n_points_final_grid,ao_num,3)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! aos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith ao on the jth grid point
|
||||
!
|
||||
! aos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point
|
||||
!
|
||||
! k = 1 : x, k= 2, y, k 3, z
|
||||
END_DOC
|
||||
integer :: i,j,m
|
||||
double precision :: aos_array(ao_num), r(3)
|
||||
double precision :: aos_grad_array(ao_num,3)
|
||||
double precision :: aos_lapl_array(ao_num,3)
|
||||
do m = 1, 3
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
|
||||
do j = 1, ao_num
|
||||
aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m)
|
||||
aos_lapl_in_r_array_transp(i,j,m) = aos_lapl_array(j,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
297
src/dft_utils_on_grid/density_matrices_in_real_space.irp.f
Normal file
297
src/dft_utils_on_grid/density_matrices_in_real_space.irp.f
Normal file
@ -0,0 +1,297 @@
|
||||
subroutine dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! input: r(1) ==> r(1) = x, r(2) = y, r(3) = z
|
||||
! output : dm_a = alpha density evaluated at r(3)
|
||||
! output : dm_b = beta density evaluated at r(3)
|
||||
END_DOC
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
|
||||
integer :: istate
|
||||
double precision :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v
|
||||
call give_all_aos_at_r(r,aos_array)
|
||||
do istate = 1, N_states
|
||||
aos_array_bis = aos_array
|
||||
! alpha density
|
||||
call dgemv('N',ao_num,ao_num,1.d0,one_body_dm_alpha_ao_for_dft(1,1,istate),ao_num,aos_array,1,0.d0,aos_array_bis,1)
|
||||
dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
|
||||
! beta density
|
||||
aos_array_bis = aos_array
|
||||
call dgemv('N',ao_num,ao_num,1.d0,one_body_dm_beta_ao_for_dft(1,1,istate),ao_num,aos_array,1,0.d0,aos_array_bis,1)
|
||||
dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
subroutine dm_dft_alpha_beta_and_all_aos_at_r(r,dm_a,dm_b,aos_array)
|
||||
BEGIN_DOC
|
||||
! input: r(1) ==> r(1) = x, r(2) = y, r(3) = z
|
||||
! output : dm_a = alpha density evaluated at r
|
||||
! output : dm_b = beta density evaluated at r
|
||||
! output : aos_array(i) = ao(i) evaluated at r
|
||||
END_DOC
|
||||
implicit none
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
|
||||
double precision, intent(out) :: aos_array(ao_num)
|
||||
integer :: istate
|
||||
double precision :: aos_array_bis(ao_num),u_dot_v
|
||||
call give_all_aos_at_r(r,aos_array)
|
||||
do istate = 1, N_states
|
||||
aos_array_bis = aos_array
|
||||
! alpha density
|
||||
call dsymv('U',ao_num,1.d0,one_body_dm_alpha_ao_for_dft(1,1,istate),size(one_body_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
|
||||
dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
|
||||
! beta density
|
||||
aos_array_bis = aos_array
|
||||
call dsymv('U',ao_num,1.d0,one_body_dm_beta_ao_for_dft(1,1,istate),size(one_body_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
|
||||
dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, grad_dm_a, grad_dm_b, aos_array, grad_aos_array)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
|
||||
! output : dm_a = alpha density evaluated at r
|
||||
! : dm_b = beta density evaluated at r
|
||||
! : aos_array(i) = ao(i) evaluated at r
|
||||
! : grad_dm_a(1) = X gradient of the alpha density evaluated in r
|
||||
! : grad_dm_a(1) = X gradient of the beta density evaluated in r
|
||||
! : grad_aos_array(1) = X gradient of the aos(i) evaluated at r
|
||||
END_DOC
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
|
||||
double precision, intent(out) :: grad_dm_a(3,N_states),grad_dm_b(3,N_states)
|
||||
double precision, intent(out) :: grad_aos_array(3,ao_num)
|
||||
integer :: i,j,istate
|
||||
double precision :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v
|
||||
double precision :: aos_grad_array(ao_num,3), aos_grad_array_bis(ao_num,3)
|
||||
|
||||
call give_all_aos_and_grad_at_r(r,aos_array,grad_aos_array)
|
||||
do i = 1, ao_num
|
||||
do j = 1, 3
|
||||
aos_grad_array(i,j) = grad_aos_array(j,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do istate = 1, N_states
|
||||
! alpha density
|
||||
! aos_array_bis = \rho_ao * aos_array
|
||||
call dsymv('U',ao_num,1.d0,one_body_dm_alpha_ao_for_dft(1,1,istate),size(one_body_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
|
||||
dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
|
||||
|
||||
! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
|
||||
grad_dm_a(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
|
||||
grad_dm_a(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
|
||||
grad_dm_a(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
|
||||
grad_dm_a *= 2.d0
|
||||
! aos_grad_array_bis = \rho_ao * aos_grad_array
|
||||
|
||||
! beta density
|
||||
call dsymv('U',ao_num,1.d0,one_body_dm_beta_ao_for_dft(1,1,istate),size(one_body_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
|
||||
dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
|
||||
|
||||
! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
|
||||
grad_dm_b(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
|
||||
grad_dm_b(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
|
||||
grad_dm_b(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
|
||||
grad_dm_b *= 2.d0
|
||||
! aos_grad_array_bis = \rho_ao * aos_grad_array
|
||||
enddo
|
||||
end
|
||||
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
|
||||
implicit none
|
||||
integer :: i,j,k,l,m,istate
|
||||
double precision :: contrib
|
||||
double precision :: r(3)
|
||||
double precision :: aos_array(ao_num),mos_array(mo_tot_num)
|
||||
do j = 1, nucl_num
|
||||
do k = 1, n_points_radial_grid -1
|
||||
do l = 1, n_points_integration_angular
|
||||
do istate = 1, N_States
|
||||
one_body_dm_mo_alpha_at_grid_points(l,k,j,istate) = 0.d0
|
||||
one_body_dm_mo_beta_at_grid_points(l,k,j,istate) = 0.d0
|
||||
enddo
|
||||
r(1) = grid_points_per_atom(1,l,k,j)
|
||||
r(2) = grid_points_per_atom(2,l,k,j)
|
||||
r(3) = grid_points_per_atom(3,l,k,j)
|
||||
|
||||
double precision :: dm_a,dm_b
|
||||
call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
|
||||
one_body_dm_mo_alpha_at_grid_points(l,k,j,1) = dm_a
|
||||
one_body_dm_mo_beta_at_grid_points(l,k,j,1) = dm_b
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_alpha_at_r, (n_points_final_grid,N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_beta_at_r, (n_points_final_grid,N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! one_body_dm_alpha_at_r(i,istate) = n_alpha(r_i,istate)
|
||||
! one_body_dm_beta_at_r(i,istate) = n_beta(r_i,istate)
|
||||
! where r_i is the ith point of the grid and istate is the state number
|
||||
END_DOC
|
||||
integer :: i,istate
|
||||
double precision :: r(3)
|
||||
double precision, allocatable :: dm_a(:),dm_b(:)
|
||||
allocate(dm_a(N_states),dm_b(N_states))
|
||||
do istate = 1, N_states
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
|
||||
one_body_dm_alpha_at_r(i,istate) = dm_a(istate)
|
||||
one_body_dm_beta_at_r(i,istate) = dm_b(istate)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_alpha_and_grad_at_r, (4,n_points_final_grid,N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_beta_and_grad_at_r, (4,n_points_final_grid,N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_grad_2_dm_alpha_at_r, (n_points_final_grid,N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_grad_2_dm_beta_at_r, (n_points_final_grid,N_states) ]
|
||||
BEGIN_DOC
|
||||
! one_body_dm_alpha_and_grad_at_r(1,i,i_state) = d\dx n_alpha(r_i,istate)
|
||||
! one_body_dm_alpha_and_grad_at_r(2,i,i_state) = d\dy n_alpha(r_i,istate)
|
||||
! one_body_dm_alpha_and_grad_at_r(3,i,i_state) = d\dz n_alpha(r_i,istate)
|
||||
! one_body_dm_alpha_and_grad_at_r(4,i,i_state) = n_alpha(r_i,istate)
|
||||
! one_body_grad_2_dm_alpha_at_r(i,istate) = d\dx n_alpha(r_i,istate)^2 + d\dy n_alpha(r_i,istate)^2 + d\dz n_alpha(r_i,istate)^2
|
||||
! where r_i is the ith point of the grid and istate is the state number
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j,k,l,m,istate
|
||||
double precision :: contrib
|
||||
double precision :: r(3)
|
||||
double precision, allocatable :: aos_array(:),grad_aos_array(:,:)
|
||||
double precision, allocatable :: dm_a(:),dm_b(:), dm_a_grad(:,:), dm_b_grad(:,:)
|
||||
allocate(dm_a(N_states),dm_b(N_states), dm_a_grad(3,N_states), dm_b_grad(3,N_states))
|
||||
allocate(aos_array(ao_num),grad_aos_array(3,ao_num))
|
||||
do istate = 1, N_states
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
!!!! Works also with the ao basis
|
||||
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, dm_a_grad, dm_b_grad, aos_array, grad_aos_array)
|
||||
one_body_dm_alpha_and_grad_at_r(1,i,istate) = dm_a_grad(1,istate)
|
||||
one_body_dm_alpha_and_grad_at_r(2,i,istate) = dm_a_grad(2,istate)
|
||||
one_body_dm_alpha_and_grad_at_r(3,i,istate) = dm_a_grad(3,istate)
|
||||
one_body_dm_alpha_and_grad_at_r(4,i,istate) = dm_a(istate)
|
||||
one_body_grad_2_dm_alpha_at_r(i,istate) = dm_a_grad(1,istate) * dm_a_grad(1,istate) + dm_a_grad(2,istate) * dm_a_grad(2,istate) + dm_a_grad(3,istate) * dm_a_grad(3,istate)
|
||||
|
||||
one_body_dm_beta_and_grad_at_r(1,i,istate) = dm_b_grad(1,istate)
|
||||
one_body_dm_beta_and_grad_at_r(2,i,istate) = dm_b_grad(2,istate)
|
||||
one_body_dm_beta_and_grad_at_r(3,i,istate) = dm_b_grad(3,istate)
|
||||
one_body_dm_beta_and_grad_at_r(4,i,istate) = dm_b(istate)
|
||||
one_body_grad_2_dm_beta_at_r(i,istate) = dm_b_grad(1,istate) * dm_b_grad(1,istate) + dm_b_grad(2,istate) * dm_b_grad(2,istate) + dm_b_grad(3,istate) * dm_b_grad(3,istate)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_and_grad_at_grid_points, (4,n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_and_grad_at_grid_points, (4,n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
|
||||
BEGIN_DOC
|
||||
! one_body_dm_mo_alpha_and_grad_at_grid_points(1,.....,i_state) = d\dx \rho_{\alpha}^{\istate}(r)
|
||||
! one_body_dm_mo_alpha_and_grad_at_grid_points(2,.....,i_state) = d\dy \rho_{\alpha}^{\istate}(r)
|
||||
! one_body_dm_mo_alpha_and_grad_at_grid_points(3,.....,i_state) = d\dz \rho_{\alpha}^{\istate}(r)
|
||||
! one_body_dm_mo_alpha_and_grad_at_grid_points(4,.....,i_state) = \rho_{\alpha}^{\istate}(r)
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j,k,l,m,istate
|
||||
double precision :: contrib
|
||||
double precision :: r(3)
|
||||
double precision :: aos_array(ao_num),grad_aos_array(3,ao_num)
|
||||
do istate = 1, N_States
|
||||
do j = 1, nucl_num
|
||||
do k = 1, n_points_radial_grid -1
|
||||
do l = 1, n_points_integration_angular
|
||||
do m = 1, 4
|
||||
one_body_dm_mo_alpha_and_grad_at_grid_points(m,l,k,j,istate) = 0.d0
|
||||
one_body_dm_mo_beta_and_grad_at_grid_points(m,l,k,j,istate) = 0.d0
|
||||
enddo
|
||||
r(1) = grid_points_per_atom(1,l,k,j)
|
||||
r(2) = grid_points_per_atom(2,l,k,j)
|
||||
r(3) = grid_points_per_atom(3,l,k,j)
|
||||
|
||||
!!!!! Works also with the ao basis
|
||||
double precision :: dm_a(N_states),dm_b(N_states), dm_a_grad(3,N_states), dm_b_grad(3,N_states)
|
||||
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, dm_a_grad, dm_b_grad, aos_array, grad_aos_array)
|
||||
one_body_dm_mo_alpha_and_grad_at_grid_points(1,l,k,j,istate) = dm_a_grad(1,istate)
|
||||
one_body_dm_mo_alpha_and_grad_at_grid_points(2,l,k,j,istate) = dm_a_grad(2,istate)
|
||||
one_body_dm_mo_alpha_and_grad_at_grid_points(3,l,k,j,istate) = dm_a_grad(3,istate)
|
||||
one_body_dm_mo_alpha_and_grad_at_grid_points(4,l,k,j,istate) = dm_a(istate)
|
||||
|
||||
|
||||
one_body_dm_mo_beta_and_grad_at_grid_points(1,l,k,j,istate) = dm_b_grad(1,istate)
|
||||
one_body_dm_mo_beta_and_grad_at_grid_points(2,l,k,j,istate) = dm_b_grad(2,istate)
|
||||
one_body_dm_mo_beta_and_grad_at_grid_points(3,l,k,j,istate) = dm_b_grad(3,istate)
|
||||
one_body_dm_mo_beta_and_grad_at_grid_points(4,l,k,j,istate) = dm_b(istate)
|
||||
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, integral_density_alpha_knowles_becke_per_atom, (nucl_num)]
|
||||
&BEGIN_PROVIDER [ double precision, integral_density_beta_knowles_becke_per_atom, (nucl_num)]
|
||||
implicit none
|
||||
double precision :: accu
|
||||
integer :: i,j,k,l,istate
|
||||
double precision :: x
|
||||
double precision :: integrand(n_points_integration_angular), weights(n_points_integration_angular)
|
||||
double precision :: f_average_angular_alpha,f_average_angular_beta
|
||||
double precision :: derivative_knowles_function,knowles_function
|
||||
|
||||
! Run over all nuclei in order to perform the Voronoi partition
|
||||
! according ot equation (6) of the paper of Becke (JCP, (88), 1988)
|
||||
! Here the m index is referred to the w_m(r) weight functions of equation (22)
|
||||
! Run over all points of integrations : there are
|
||||
! n_points_radial_grid (i) * n_points_integration_angular (k)
|
||||
do j = 1, nucl_num
|
||||
integral_density_alpha_knowles_becke_per_atom(j) = 0.d0
|
||||
integral_density_beta_knowles_becke_per_atom(j) = 0.d0
|
||||
do i = 1, n_points_radial_grid-1
|
||||
! Angular integration over the solid angle Omega for a FIXED angular coordinate "r"
|
||||
f_average_angular_alpha = 0.d0
|
||||
f_average_angular_beta = 0.d0
|
||||
do istate = 1, N_states
|
||||
do k = 1, n_points_integration_angular
|
||||
f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j,istate) * weight_functions_at_grid_points(k,i,j)
|
||||
f_average_angular_beta += weights_angular_points(k) * one_body_dm_mo_beta_at_grid_points(k,i,j,istate) * weight_functions_at_grid_points(k,i,j)
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
|
||||
double precision :: contrib_integration
|
||||
! print*,m_knowles
|
||||
contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x) &
|
||||
*knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
|
||||
integral_density_alpha_knowles_becke_per_atom(j) += contrib_integration *f_average_angular_alpha
|
||||
integral_density_beta_knowles_becke_per_atom(j) += contrib_integration *f_average_angular_beta
|
||||
enddo
|
||||
integral_density_alpha_knowles_becke_per_atom(j) *= dr_radial_integral
|
||||
integral_density_beta_knowles_becke_per_atom(j) *= dr_radial_integral
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
56
src/dft_utils_on_grid/mo_on_grid.irp.f
Normal file
56
src/dft_utils_on_grid/mo_on_grid.irp.f
Normal file
@ -0,0 +1,56 @@
|
||||
BEGIN_PROVIDER[double precision, mos_in_r_array, (mo_tot_num,n_points_final_grid)]
|
||||
&BEGIN_PROVIDER[double precision, mos_in_r_array_transp,(n_points_final_grid,mo_tot_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! mos_in_r_array(i,j) = value of the ith mo on the jth grid point
|
||||
!
|
||||
! mos_in_r_array_transp(i,j) = value of the jth mo on the ith grid point
|
||||
END_DOC
|
||||
integer :: i,j
|
||||
double precision :: mos_array(mo_tot_num), r(3)
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
call give_all_mos_at_r(r,mos_array)
|
||||
do j = 1, mo_tot_num
|
||||
mos_in_r_array(j,i) = mos_array(j)
|
||||
mos_in_r_array_transp(i,j) = mos_array(j)
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER[double precision, mos_grad_in_r_array,(mo_tot_num,n_points_final_grid,3)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! mos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith mo on the jth grid point
|
||||
!
|
||||
! mos_grad_in_r_array_transp(i,j,k) = value of the kth component of the gradient of jth mo on the ith grid point
|
||||
!
|
||||
! k = 1 : x, k= 2, y, k 3, z
|
||||
END_DOC
|
||||
integer :: m
|
||||
mos_grad_in_r_array = 0.d0
|
||||
do m=1,3
|
||||
call dgemm('N','N',mo_tot_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_tot_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_grad_in_r_array(1,1,m),mo_tot_num)
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER[double precision, mos_lapl_in_r_array,(mo_tot_num,n_points_final_grid,3)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! mos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith mo on the jth grid point
|
||||
!
|
||||
! mos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth mo on the ith grid point
|
||||
!
|
||||
! k = 1 : x, k= 2, y, k 3, z
|
||||
END_DOC
|
||||
integer :: m
|
||||
mos_lapl_in_r_array = 0.d0
|
||||
do m=1,3
|
||||
call dgemm('N','N',mo_tot_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_tot_num,aos_lapl_in_r_array(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_tot_num)
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
4
src/dft_utils_one_body/NEED
Normal file
4
src/dft_utils_one_body/NEED
Normal file
@ -0,0 +1,4 @@
|
||||
density_for_dft
|
||||
dft_utils_on_grid
|
||||
integrals_bielec
|
||||
integrals_bielec_erf
|
12
src/dft_utils_one_body/README.rst
Normal file
12
src/dft_utils_one_body/README.rst
Normal file
@ -0,0 +1,12 @@
|
||||
===========
|
||||
RSDFT_Utils
|
||||
===========
|
||||
|
||||
Needed Modules
|
||||
==============
|
||||
.. Do not edit this section It was auto-generated
|
||||
.. by the `update_README.py` script.
|
||||
Documentation
|
||||
=============
|
||||
.. Do not edit this section It was auto-generated
|
||||
.. by the `update_README.py` script.
|
125
src/dft_utils_one_body/e_c_e_x_and_potentials_mo.irp.f
Normal file
125
src/dft_utils_one_body/e_c_e_x_and_potentials_mo.irp.f
Normal file
@ -0,0 +1,125 @@
|
||||
|
||||
BEGIN_PROVIDER [double precision, potential_x_alpha_ao,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_x_beta_ao,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_alpha_ao,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_beta_ao,(ao_num,ao_num,N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! alpha/beta exchange/correlation potentials on the AO basis
|
||||
END_DOC
|
||||
|
||||
if(trim(exchange_functional)=="short_range_LDA")then
|
||||
potential_x_alpha_ao = potential_x_alpha_ao_LDA
|
||||
potential_x_beta_ao = potential_x_beta_ao_LDA
|
||||
else if(exchange_functional.EQ."short_range_PBE")then
|
||||
potential_x_alpha_ao = potential_x_alpha_ao_PBE
|
||||
potential_x_beta_ao = potential_x_beta_ao_PBE
|
||||
else if(exchange_functional.EQ."None")then
|
||||
potential_x_alpha_ao = 0.d0
|
||||
potential_x_beta_ao = 0.d0
|
||||
else
|
||||
print*, 'Exchange functional required does not exist ...'
|
||||
print*,'exchange_functional',exchange_functional
|
||||
stop
|
||||
endif
|
||||
|
||||
if(trim(correlation_functional)=="short_range_LDA")then
|
||||
potential_c_alpha_ao = potential_c_alpha_ao_LDA
|
||||
potential_c_beta_ao = potential_c_beta_ao_LDA
|
||||
else if(correlation_functional.EQ."short_range_PBE")then
|
||||
potential_c_alpha_ao = potential_c_alpha_ao_PBE
|
||||
potential_c_beta_ao = potential_c_beta_ao_PBE
|
||||
else if(correlation_functional.EQ."None")then
|
||||
potential_c_alpha_ao = 0.d0
|
||||
potential_c_beta_ao = 0.d0
|
||||
else
|
||||
print*, 'Correlation functional required does not ecist ...'
|
||||
print*,'correlation_functional',correlation_functional
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, potential_x_alpha_mo,(mo_tot_num,mo_tot_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_x_beta_mo,(mo_tot_num,mo_tot_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_alpha_mo,(mo_tot_num,mo_tot_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_beta_mo,(mo_tot_num,mo_tot_num,N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! alpha/beta exchange/correlation potentials on the MO basis
|
||||
END_DOC
|
||||
integer :: istate
|
||||
do istate = 1, N_states
|
||||
call ao_to_mo( &
|
||||
potential_x_alpha_ao(1,1,istate), &
|
||||
size(potential_x_alpha_ao,1), &
|
||||
potential_x_alpha_mo(1,1,istate), &
|
||||
size(potential_x_alpha_mo,1) &
|
||||
)
|
||||
|
||||
call ao_to_mo( &
|
||||
potential_x_beta_ao(1,1,istate), &
|
||||
size(potential_x_beta_ao,1), &
|
||||
potential_x_beta_mo(1,1,istate), &
|
||||
size(potential_x_beta_mo,1) &
|
||||
)
|
||||
|
||||
|
||||
call ao_to_mo( &
|
||||
potential_c_alpha_ao(1,1,istate), &
|
||||
size(potential_c_alpha_ao,1), &
|
||||
potential_c_alpha_mo(1,1,istate), &
|
||||
size(potential_c_alpha_mo,1) &
|
||||
)
|
||||
|
||||
call ao_to_mo( &
|
||||
potential_c_beta_ao(1,1,istate), &
|
||||
size(potential_c_beta_ao,1), &
|
||||
potential_c_beta_mo(1,1,istate), &
|
||||
size(potential_c_beta_mo,1) &
|
||||
)
|
||||
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, energy_x, (N_states)]
|
||||
&BEGIN_PROVIDER [double precision, energy_c, (N_states)]
|
||||
implicit none
|
||||
if(trim(exchange_functional)=="short_range_LDA")then
|
||||
energy_x = energy_x_LDA
|
||||
energy_x = energy_x_LDA
|
||||
else if(exchange_functional.EQ."short_range_PBE")then
|
||||
energy_x = energy_x_PBE
|
||||
energy_x = energy_x_PBE
|
||||
else if(exchange_functional.EQ."None")then
|
||||
energy_x = 0.d0
|
||||
energy_x = 0.d0
|
||||
else
|
||||
print*, 'Exchange functional required does not exist ...'
|
||||
print*,'exchange_functional',exchange_functional
|
||||
stop
|
||||
endif
|
||||
|
||||
if(trim(correlation_functional)=="short_range_LDA")then
|
||||
energy_c = energy_c_LDA
|
||||
energy_c = energy_c_LDA
|
||||
else if(correlation_functional.EQ."short_range_PBE")then
|
||||
energy_c = energy_c_PBE
|
||||
energy_c = energy_c_PBE
|
||||
else if(correlation_functional.EQ."None")then
|
||||
energy_c = 0.d0
|
||||
energy_c = 0.d0
|
||||
else
|
||||
print*, 'Correlation functional required does not ecist ...'
|
||||
print*,'correlation_functional',correlation_functional
|
||||
stop
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
1007
src/dft_utils_one_body/exc_sr_lda.irp.f
Normal file
1007
src/dft_utils_one_body/exc_sr_lda.irp.f
Normal file
File diff suppressed because it is too large
Load Diff
477
src/dft_utils_one_body/exc_sr_pbe.irp.f
Normal file
477
src/dft_utils_one_body/exc_sr_pbe.irp.f
Normal file
@ -0,0 +1,477 @@
|
||||
subroutine ec_pbe_sr(mu,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
|
||||
!************************************************************************
|
||||
! Short-range PBE correlation energy functional for erf interaction
|
||||
!
|
||||
!
|
||||
!************************************************************************
|
||||
include 'constants.include.F'
|
||||
implicit none
|
||||
! input
|
||||
double precision, intent(in) :: rhoc,rhoo,mu
|
||||
double precision, intent(in) :: sigmacc,sigmaco,sigmaoo
|
||||
! output
|
||||
double precision, intent(out) :: ec
|
||||
double precision, intent(out) :: vrhoc,vrhoo
|
||||
double precision, intent(out) :: vsigmacc,vsigmaco,vsigmaoo
|
||||
! local
|
||||
double precision tol
|
||||
parameter(tol=1d-12)
|
||||
|
||||
character(len=30) namedummy
|
||||
|
||||
double precision eccerflda
|
||||
double precision vrhoccerflda
|
||||
double precision vrhoocerflda
|
||||
|
||||
double precision ecclda
|
||||
double precision vrhocclda
|
||||
double precision vrhooclda
|
||||
|
||||
integer i,igrad
|
||||
double precision rho,drho2,rhoa,rhob
|
||||
double precision ecerflda,decerfldadrho
|
||||
double precision eclda,decldadrho
|
||||
double precision ecerfpbe,decerfpbedrho,decerfpbedrhoo
|
||||
double precision decerfpbeddrho2
|
||||
double precision arglog,arglogs,arglogss,alpha,beta,betas,gamma
|
||||
double precision Aa,Ab,Ac,Aas,tq,tqs,tqss,decerfpur,decpur
|
||||
double precision t1,t2,t3,t4,t5,t6,t7,t8,t9,t10
|
||||
double precision t11,t12,t13,t14,t15,t16,t17,t18,t19
|
||||
double precision zeta,phi,phi2,phi3,phi4,phis,arglogsc
|
||||
double precision dlogarglog
|
||||
double precision, parameter :: f13=0.333333333333333d0
|
||||
|
||||
|
||||
! Parameter of the modified interaction
|
||||
|
||||
ec = 0.d0
|
||||
vrhoc = 0.d0
|
||||
vrhoo = 0.d0
|
||||
vsigmacc = 0.d0
|
||||
vsigmaco = 0.d0
|
||||
vsigmaoo = 0.d0
|
||||
|
||||
! First-type gradient functional
|
||||
igrad=1
|
||||
|
||||
alpha=2.78d0
|
||||
gamma=3.1091d-2
|
||||
|
||||
|
||||
! test on density
|
||||
if (dabs(rhoc).lt.tol) return
|
||||
double precision :: vc_a,vc_b
|
||||
! Spin polarisation
|
||||
rhoa=max((rhoc+rhoo)*.5d0,1.0d-15)
|
||||
rhob=max((rhoc-rhoo)*.5d0,1.0d-15)
|
||||
|
||||
call ec_lda_sr(mu,rhoa,rhob,eccerflda,vc_a,vc_b)
|
||||
ecerflda = eccerflda
|
||||
vrhoccerflda = 0.5d0 * (vc_a + vc_b)
|
||||
vrhoocerflda = 0.5d0 * (vc_a - vc_b)
|
||||
|
||||
! Density
|
||||
rho = rhoc
|
||||
|
||||
! Square of density gradient
|
||||
drho2 = sigmacc
|
||||
|
||||
zeta = (rhoa-rhob)/(rhoa+rhob)
|
||||
|
||||
! LDA energy density
|
||||
double precision :: vc_a_lda,vc_b_lda
|
||||
call ec_lda(rhoa,rhob,ecclda,vc_a_lda,vc_b_lda)
|
||||
eclda = ecclda
|
||||
!decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
|
||||
!decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
|
||||
|
||||
if ((ecerflda/eclda).le.0d0) then
|
||||
beta=0d0
|
||||
else
|
||||
beta=6.6725d-2*(ecerflda/eclda)**alpha
|
||||
endif
|
||||
phi=((1d0+zeta)**(2d0/3d0)+(1d0-zeta)**(2d0/3d0))/2d0
|
||||
phi2=phi*phi
|
||||
phi3=phi2*phi
|
||||
phi4=phi3*phi
|
||||
tq=drho2*6.346820607d-2*rho**(-7d0/3d0)/phi2
|
||||
Ab=dexp(-ecerflda/(rho*gamma*phi3))-1d0
|
||||
if (dabs(Ab).le.dabs(beta*tol)) then
|
||||
ecerfpbe=ecerflda
|
||||
else
|
||||
Aa=beta/(gamma*Ab)
|
||||
Ac=1d0+Aa*tq+Aa**2*tq**2
|
||||
if (Aa.lt.tol) Aa=tol
|
||||
arglog=1d0+beta*(1d0-1d0/Ac)/(gamma*Aa)
|
||||
ecerfpbe=ecerflda+rho*phi3*gamma*dlog(arglog)
|
||||
end if
|
||||
|
||||
ec = ecerfpbe
|
||||
|
||||
! if(ldebug) write(*,*)"ecerfpbe=",ecerfpbe
|
||||
|
||||
! Derive
|
||||
|
||||
|
||||
! LDA energy density derivative
|
||||
decerfldadrho = vrhoccerflda
|
||||
decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
|
||||
|
||||
decerfpur=(decerfldadrho-ecerflda/rho)/rho
|
||||
decpur=(decldadrho-eclda/rho)/rho
|
||||
betas=alpha*beta*(decerfpur*rho/ecerflda-decpur*rho/eclda)
|
||||
phis=((rhoa - rhob)*((rhoa/(rhoa + rhob))**f13 - (rhob/(rhoa + rhob))**f13))/(3d0*2d0**f13*(rhoa/(rhoa + rhob))**f13*(rhob/(rhoa + rhob))**f13*(rhoa + rhob)**2)
|
||||
if (dabs(Ab).le.dabs(beta*tol)) then
|
||||
decerfpbedrho=decerfldadrho
|
||||
else
|
||||
Aas=betas/(gamma*Ab)+Aa*(1d0+1d0/Ab)*(decerfpur/phi3-3d0*phis*ecerflda/(rho*phi4))/gamma
|
||||
tqs=-7d0*tq/(3d0*rho)-2d0*tq*phis/phi
|
||||
arglogs=betas*tq*(1d0+Aa*tq)/(Ac*gamma)+beta*tqs*(1d0+Aa*tq)/(Ac*gamma)-beta*tq*Aa*tq*(Aas*tq+Aa*tqs)*(2d0+Aa*tq)/(Ac**2*gamma)
|
||||
dlogarglog=dlog(arglog)
|
||||
decerfpbedrho=decerfldadrho+gamma*(phi3*dlogarglog+3d0*rho*phis*phi2*dlogarglog+rho*phi3*arglogs/arglog)
|
||||
end if
|
||||
|
||||
if (dabs(Ab).le.dabs(beta*tol)) then
|
||||
decerfpbeddrho2=0.0d0
|
||||
else
|
||||
arglogsc=Ab*(Aa+2d0*Aa*Aa*tq)/(Ac*Ac)
|
||||
tqss=6.346820607d-2*rho**(-7d0/3d0)/phi2
|
||||
arglogss=tqss*arglogsc
|
||||
decerfpbeddrho2=rho*gamma*phi3*arglogss/arglog
|
||||
end if
|
||||
|
||||
! LDA energy density derivative
|
||||
decerfldadrho = vrhoocerflda
|
||||
decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
|
||||
|
||||
decerfpur=decerfldadrho/rho
|
||||
decpur=decldadrho/rho
|
||||
betas=alpha*beta*(decerfpur*rho/ecerflda-decpur*rho/eclda)
|
||||
phis=(rhob*(rhoa/(rhoa + rhob))**(2d0*f13)-rhoa*(rhob/(rhoa + rhob))**(2d0*f13))/(3d0*2d0**f13*rhoa*rhob)
|
||||
|
||||
if (dabs(Ab).le.dabs(beta*tol)) then
|
||||
decerfpbedrhoo=decerfldadrho
|
||||
else
|
||||
Aas=betas/(gamma*Ab)+Aa*(1d0+1d0/Ab)*(decerfpur/phi3-3d0*phis*ecerflda/(rho*phi4))/gamma
|
||||
tqs=-2d0*tq*phis/phi
|
||||
arglogs=betas*tq*(1d0+Aa*tq)/(Ac*gamma)+beta*tqs*(1d0+Aa*tq)/(Ac*gamma)-beta*tq*Aa*tq*(Aas*tq+Aa*tqs)*(2d0+Aa*tq)/(Ac**2*gamma)
|
||||
decerfpbedrhoo=decerfldadrho+gamma*(3d0*rho*phis*phi2*dlog(arglog)+rho*phi3*arglogs/arglog)
|
||||
end if
|
||||
|
||||
! derivatives
|
||||
vrhoc = vrhoc + decerfpbedrho
|
||||
vrhoo = vrhoo + decerfpbedrhoo
|
||||
vsigmacc = vsigmacc + decerfpbeddrho2
|
||||
|
||||
|
||||
end
|
||||
|
||||
subroutine ex_pbe_sr(mu,rho_a,rho_b,grd_rho_a_2,grd_rho_b_2,grd_rho_a_b,ex,vx_rho_a,vx_rho_b,vx_grd_rho_a_2,vx_grd_rho_b_2,vx_grd_rho_a_b)
|
||||
BEGIN_DOC
|
||||
!rho_a = density alpha
|
||||
!rho_b = density beta
|
||||
!grd_rho_a_2 = (gradient rho_a)^2
|
||||
!grd_rho_b_2 = (gradient rho_b)^2
|
||||
!grd_rho_a_b = (gradient rho_a).(gradient rho_b)
|
||||
!ex = exchange energy density at point r
|
||||
!vx_rho_a = d ex / d rho_a
|
||||
!vx_rho_b = d ex / d rho_b
|
||||
!vx_grd_rho_a_2 = d ex / d grd_rho_a_2
|
||||
!vx_grd_rho_b_2 = d ex / d grd_rho_b_2
|
||||
!vx_grd_rho_a_b = d ex / d grd_rho_a_b
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
! input
|
||||
double precision, intent(in) :: mu,rho_a, rho_b
|
||||
double precision, intent(in) :: grd_rho_a_2, grd_rho_b_2, grd_rho_a_b
|
||||
|
||||
! output
|
||||
double precision, intent(out) :: ex
|
||||
double precision, intent(out) :: vx_rho_a, vx_rho_b
|
||||
double precision, intent(out) :: vx_grd_rho_a_2, vx_grd_rho_b_2, vx_grd_rho_a_b
|
||||
|
||||
! function
|
||||
double precision berf
|
||||
double precision dberfda
|
||||
|
||||
! local
|
||||
double precision, parameter :: tol=1d-12
|
||||
double precision, parameter :: f13=0.333333333333333d0
|
||||
|
||||
double precision exerflda,vxerflda_a,vxerflda_b
|
||||
double precision dexerfldadrho
|
||||
double precision exerfpbe_a, exerfpbe_b
|
||||
double precision dexerfpbedrho_a, dexerfpbedrho_b
|
||||
double precision dexerfpbeddrho2_a, dexerfpbeddrho2_b
|
||||
|
||||
double precision rho,drho2
|
||||
double precision rho_a_2, rho_b_2
|
||||
double precision t1,t2,t3,t4
|
||||
double precision kappa,sq,sqs,sqss,fx,fxs,ksig
|
||||
|
||||
! Parameter of the modified interaction
|
||||
|
||||
! initialization
|
||||
ex=0.d0
|
||||
vx_rho_a=0.d0
|
||||
vx_rho_b=0.d0
|
||||
vx_grd_rho_a_2=0.d0
|
||||
vx_grd_rho_b_2=0.d0
|
||||
vx_grd_rho_a_b=0.d0
|
||||
|
||||
|
||||
! spin scaling relation Ex[rho_a,rho_b] = (1/2) (Ex[2rho_a,2rho_a] + Ex[2rho_b,2rho_b])
|
||||
|
||||
! two times spin alpha density
|
||||
rho = max(rho_a,tol)*2.d0
|
||||
|
||||
! test on density
|
||||
if (rho >= tol) then
|
||||
|
||||
! call srLDA Ex[2*rho_a,2*rho_a]
|
||||
call ex_lda_sr(mu,rho_a,rho_a,exerflda,vxerflda_a,vxerflda_b)
|
||||
dexerfldadrho = (vxerflda_a + vxerflda_b)*0.5d0
|
||||
|
||||
! square of two times spin alpha density gradient
|
||||
drho2=max(grd_rho_a_2,0d0)*4.0d0
|
||||
|
||||
kappa=0.804d0
|
||||
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
|
||||
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
|
||||
exerfpbe_a=exerflda*fx
|
||||
|
||||
! Derivatives
|
||||
sqs=-8d0*sq/(3d0*rho)
|
||||
fxs=kappa**2*(-1.616204596739954813d-1*mu*rho**(-4d0*f13)/3d0*dberfda(1.616204596739954813d-1*mu*rho**(-f13))*sq+berf(1.616204596739954813d-1*mu*rho**(-f13))*sqs)/(kappa+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq)**2
|
||||
dexerfpbedrho_a=dexerfldadrho*fx+exerflda*fxs
|
||||
sqss=2.6121172985233599567768d-2*rho**(-8d0/3d0)
|
||||
dexerfpbeddrho2_a=exerflda*berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sqss*kappa**2/(kappa+berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sq)**2
|
||||
|
||||
endif
|
||||
|
||||
|
||||
! two times spin beta density
|
||||
rho = max(rho_b,tol)*2.d0
|
||||
|
||||
! test on density
|
||||
if (rho >= tol) then
|
||||
|
||||
! call srLDA Ex[2*rho_b,2*rho_b]
|
||||
call ex_lda_sr(mu,rho_b,rho_b,exerflda,vxerflda_a,vxerflda_b)
|
||||
dexerfldadrho = (vxerflda_a + vxerflda_b)*0.5d0
|
||||
|
||||
! square of two times spin beta density gradient
|
||||
drho2=max(grd_rho_b_2,0d0)*4.0d0
|
||||
|
||||
kappa=0.804d0
|
||||
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
|
||||
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
|
||||
exerfpbe_b=exerflda*fx
|
||||
|
||||
! Derivatives
|
||||
sqs=-8d0*sq/(3d0*rho)
|
||||
fxs=kappa**2*(-1.616204596739954813d-1*mu*rho**(-4d0*f13)/3d0*dberfda(1.616204596739954813d-1*mu*rho**(-f13))*sq+berf(1.616204596739954813d-1*mu*rho**(-f13))*sqs)/(kappa+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq)**2
|
||||
dexerfpbedrho_b=dexerfldadrho*fx+exerflda*fxs
|
||||
sqss=2.6121172985233599567768d-2*rho**(-8d0/3d0)
|
||||
dexerfpbeddrho2_b=exerflda*berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sqss*kappa**2/(kappa+berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sq)**2
|
||||
|
||||
endif
|
||||
|
||||
|
||||
ex = (exerfpbe_a+exerfpbe_b)*0.5d0
|
||||
vx_rho_a = dexerfpbedrho_a
|
||||
vx_rho_b = dexerfpbedrho_a
|
||||
vx_grd_rho_a_2 = 2.d0*dexerfpbeddrho2_a
|
||||
vx_grd_rho_b_2 = 2.d0*dexerfpbeddrho2_b
|
||||
vx_grd_rho_a_b = 0.d0
|
||||
|
||||
end
|
||||
|
||||
subroutine ex_pbe_sr_only(mu,rho_a,rho_b,grd_rho_a_2,grd_rho_b_2,grd_rho_a_b,ex)
|
||||
BEGIN_DOC
|
||||
!rho_a = density alpha
|
||||
!rho_b = density beta
|
||||
!grd_rho_a_2 = (gradient rho_a)^2
|
||||
!grd_rho_b_2 = (gradient rho_b)^2
|
||||
!grd_rho_a_b = (gradient rho_a).(gradient rho_b)
|
||||
!ex = exchange energy density at point r
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
! input
|
||||
double precision, intent(in) :: mu,rho_a, rho_b
|
||||
double precision, intent(in) :: grd_rho_a_2, grd_rho_b_2, grd_rho_a_b
|
||||
|
||||
! output
|
||||
double precision, intent(out) :: ex
|
||||
|
||||
! function
|
||||
double precision berf
|
||||
|
||||
! local
|
||||
double precision, parameter :: tol=1d-12
|
||||
double precision, parameter :: f13=0.333333333333333d0
|
||||
|
||||
double precision exerflda,vxerflda_a,vxerflda_b
|
||||
double precision exerfpbe_a, exerfpbe_b
|
||||
|
||||
double precision rho,drho2
|
||||
double precision kappa,sq,fx
|
||||
|
||||
|
||||
! initialization
|
||||
ex=0.d0
|
||||
|
||||
|
||||
! spin scaling relation Ex[rho_a,rho_b] = (1/2) (Ex[2rho_a,2rho_a] + Ex[2rho_b,2rho_b])
|
||||
|
||||
! two times spin alpha density
|
||||
rho = max(rho_a,tol)*2.d0
|
||||
|
||||
! test on density
|
||||
if (rho >= tol) then
|
||||
|
||||
! call srLDA Ex[2*rho_a,2*rho_a]
|
||||
call ex_lda_sr(mu,rho_a,rho_a,exerflda,vxerflda_a,vxerflda_b)
|
||||
|
||||
! square of two times spin alpha density gradient
|
||||
drho2=max(grd_rho_a_2,0d0)*4.0d0
|
||||
|
||||
kappa=0.804d0
|
||||
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
|
||||
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
|
||||
exerfpbe_a=exerflda*fx
|
||||
|
||||
endif
|
||||
|
||||
|
||||
! two times spin beta density
|
||||
rho = max(rho_b,tol)*2.d0
|
||||
|
||||
! test on density
|
||||
if (rho >= tol) then
|
||||
|
||||
! call srLDA Ex[2*rho_b,2*rho_b]
|
||||
call ex_lda_sr(mu,rho_b,rho_b,exerflda,vxerflda_a,vxerflda_b)
|
||||
|
||||
! square of two times spin beta density gradient
|
||||
drho2=max(grd_rho_b_2,0d0)*4.0d0
|
||||
|
||||
kappa=0.804d0
|
||||
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
|
||||
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
|
||||
exerfpbe_b=exerflda*fx
|
||||
|
||||
endif
|
||||
|
||||
ex = (exerfpbe_a+exerfpbe_b)*0.5d0
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine ec_pbe_only(mu,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec)
|
||||
!************************************************************************
|
||||
! Short-range PBE correlation energy functional for erf interaction
|
||||
!
|
||||
!
|
||||
!************************************************************************
|
||||
include 'constants.include.F'
|
||||
implicit none
|
||||
! input
|
||||
double precision, intent(in) :: rhoc,rhoo,mu
|
||||
double precision, intent(in) :: sigmacc,sigmaco,sigmaoo
|
||||
! output
|
||||
double precision, intent(out) :: ec
|
||||
! local
|
||||
double precision tol
|
||||
parameter(tol=1d-12)
|
||||
|
||||
character(len=30) namedummy
|
||||
|
||||
double precision eccerflda
|
||||
double precision vrhoccerflda
|
||||
double precision vrhoocerflda
|
||||
|
||||
double precision ecclda
|
||||
double precision vrhocclda
|
||||
double precision vrhooclda
|
||||
|
||||
integer i,igrad
|
||||
double precision rho,drho2,rhoa,rhob
|
||||
double precision ecerflda
|
||||
double precision eclda,decldadrho
|
||||
double precision ecerfpbe
|
||||
double precision arglog,alpha,beta,gamma
|
||||
double precision Aa,Ab,Ac,tq
|
||||
double precision zeta,phi,phi2,phi3,phi4
|
||||
double precision, parameter :: f13=0.333333333333333d0
|
||||
|
||||
|
||||
! Parameter of the modified interaction
|
||||
|
||||
ec = 0.d0
|
||||
|
||||
! First-type gradient functional
|
||||
igrad=1
|
||||
|
||||
alpha=2.78d0
|
||||
gamma=3.1091d-2
|
||||
|
||||
! test on density
|
||||
if (dabs(rhoc).lt.tol) return
|
||||
double precision :: vc_a,vc_b
|
||||
! Spin polarisation
|
||||
rhoa=max((rhoc+rhoo)*.5d0,1.0d-15)
|
||||
rhob=max((rhoc-rhoo)*.5d0,1.0d-15)
|
||||
|
||||
call ec_lda_sr(mu,rhoa,rhob,eccerflda,vc_a,vc_b)
|
||||
ecerflda = eccerflda
|
||||
vrhoccerflda = 0.5d0 * (vc_a + vc_b)
|
||||
vrhoocerflda = 0.5d0 * (vc_a - vc_b)
|
||||
|
||||
! Density
|
||||
rho = rhoc
|
||||
rho = max(rho,1.d-10)
|
||||
|
||||
! Square of density gradient
|
||||
drho2 = sigmacc
|
||||
|
||||
zeta = (rhoa-rhob)/(rhoa+rhob)
|
||||
zeta = max(zeta,1.d-10)
|
||||
|
||||
! LDA energy density
|
||||
double precision :: vc_a_lda,vc_b_lda
|
||||
call ec_lda(rhoa,rhob,ecclda,vc_a_lda,vc_b_lda)
|
||||
eclda = ecclda
|
||||
decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
|
||||
decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
|
||||
|
||||
if ((ecerflda/eclda).le.0d0) then
|
||||
beta=0d0
|
||||
else
|
||||
beta=6.6725d-2*(ecerflda/eclda)**alpha
|
||||
endif
|
||||
phi=((1d0+zeta)**(2d0/3d0)+(1d0-zeta)**(2d0/3d0))/2d0
|
||||
phi2=phi*phi
|
||||
phi3=phi2*phi
|
||||
phi4=phi3*phi
|
||||
tq=drho2*6.346820607d-2*rho**(-7d0/3d0)/phi2
|
||||
Ab=dexp(-ecerflda/(rho*gamma*phi3))-1d0
|
||||
if (dabs(Ab).le.dabs(beta*tol)) then
|
||||
ecerfpbe=ecerflda
|
||||
else
|
||||
Aa=beta/(gamma*Ab)
|
||||
Ac=1d0+Aa*tq+Aa**2*tq**2
|
||||
if (Aa.lt.tol) Aa=tol
|
||||
arglog=1d0+beta*(1d0-1d0/Ac)/(gamma*Aa)
|
||||
arglog=max(arglog,1.d-10)
|
||||
ecerfpbe=ecerflda+rho*phi3*gamma*dlog(arglog)
|
||||
end if
|
||||
|
||||
ec = ecerfpbe
|
||||
|
||||
end
|
29
src/dft_utils_one_body/one_body_psi_energy.irp.f
Normal file
29
src/dft_utils_one_body/one_body_psi_energy.irp.f
Normal file
@ -0,0 +1,29 @@
|
||||
BEGIN_PROVIDER [double precision, psi_dft_energy_kinetic, (N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, psi_dft_energy_nuclear_elec, (N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, psi_dft_energy_h_core, (N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! kinetic, electron-nuclear and total h_core energy computed with the density matrix one_body_dm_mo_beta_for_dft+one_body_dm_mo_alpha_for_dft
|
||||
END_DOC
|
||||
integer :: i,j,istate
|
||||
double precision :: accu
|
||||
psi_dft_energy_kinetic = 0.d0
|
||||
psi_dft_energy_nuclear_elec = 0.d0
|
||||
do istate = 1, N_states
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
psi_dft_energy_kinetic(istate) += ( one_body_dm_mo_alpha_for_dft(j,i,istate)+one_body_dm_mo_beta_for_dft(j,i,istate)) * mo_kinetic_integral(j,i)
|
||||
psi_dft_energy_nuclear_elec(istate) += ( one_body_dm_mo_alpha_for_dft(j,i,istate)+one_body_dm_mo_beta_for_dft(j,i,istate)) * mo_nucl_elec_integral(j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i = 1, N_states
|
||||
do j = 1, mo_tot_num
|
||||
accu += one_body_dm_mo_alpha_for_dft(j,j,i) + one_body_dm_mo_beta_for_dft(j,j,i)
|
||||
enddo
|
||||
accu = (elec_alpha_num + elec_beta_num ) / accu
|
||||
psi_energy_h_core(i) = psi_dft_energy_h_core(i) * accu
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
272
src/dft_utils_one_body/potentials_ao_on_grids.irp.f
Normal file
272
src/dft_utils_one_body/potentials_ao_on_grids.irp.f
Normal file
@ -0,0 +1,272 @@
|
||||
BEGIN_PROVIDER[double precision, aos_vc_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_vc_beta_LDA_w, (n_points_final_grid,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_vx_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_vx_beta_LDA_w, (n_points_final_grid,ao_num,N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! aos_vxc_alpha_LDA_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
|
||||
END_DOC
|
||||
integer :: istate,i,j
|
||||
double precision :: r(3)
|
||||
double precision :: mu,weight
|
||||
double precision :: e_c,vc_a,vc_b,e_x,vx_a,vx_b
|
||||
double precision, allocatable :: rhoa(:),rhob(:)
|
||||
allocate(rhoa(N_states), rhob(N_states))
|
||||
do istate = 1, N_states
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
weight=final_weight_functions_at_final_grid_points(i)
|
||||
rhoa(istate) = one_body_dm_alpha_at_r(i,istate)
|
||||
rhob(istate) = one_body_dm_beta_at_r(i,istate)
|
||||
call ec_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
|
||||
call ex_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
|
||||
do j =1, ao_num
|
||||
aos_vc_alpha_LDA_w(i,j,istate) = vc_a * aos_in_r_array(j,i)*weight
|
||||
aos_vc_beta_LDA_w(i,j,istate) = vc_b * aos_in_r_array(j,i)*weight
|
||||
aos_vx_alpha_LDA_w(i,j,istate) = vx_a * aos_in_r_array(j,i)*weight
|
||||
aos_vx_beta_LDA_w(i,j,istate) = vx_b * aos_in_r_array(j,i)*weight
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER[double precision, energy_x_LDA, (N_states) ]
|
||||
&BEGIN_PROVIDER[double precision, energy_c_LDA, (N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! exchange/correlation energy with the short range LDA functional
|
||||
END_DOC
|
||||
integer :: istate,i,j
|
||||
double precision :: r(3)
|
||||
double precision :: mu,weight
|
||||
double precision :: e_c,vc_a,vc_b,e_x,vx_a,vx_b
|
||||
double precision, allocatable :: rhoa(:),rhob(:)
|
||||
allocate(rhoa(N_states), rhob(N_states))
|
||||
energy_x_LDA = 0.d0
|
||||
energy_c_LDA = 0.d0
|
||||
do istate = 1, N_states
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
weight=final_weight_functions_at_final_grid_points(i)
|
||||
rhoa(istate) = one_body_dm_alpha_at_r(i,istate)
|
||||
rhob(istate) = one_body_dm_beta_at_r(i,istate)
|
||||
call ec_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
|
||||
call ex_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
|
||||
energy_x_LDA(istate) += weight * e_x
|
||||
energy_c_LDA(istate) += weight * e_c
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, potential_x_alpha_ao_LDA,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_x_beta_ao_LDA,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_LDA,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_beta_ao_LDA,(ao_num,ao_num,N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis
|
||||
END_DOC
|
||||
integer :: istate
|
||||
double precision :: wall_1,wall_2
|
||||
call wall_time(wall_1)
|
||||
do istate = 1, N_states
|
||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vc_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_c_alpha_ao_LDA(1,1,istate),ao_num)
|
||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vc_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_c_beta_ao_LDA(1,1,istate),ao_num)
|
||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vx_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_x_alpha_ao_LDA(1,1,istate),ao_num)
|
||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vx_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_x_beta_ao_LDA(1,1,istate),ao_num)
|
||||
enddo
|
||||
call wall_time(wall_2)
|
||||
print*,'time to provide potential_x/c_alpha/beta_ao_LDA = ',wall_2 - wall_1
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER[double precision, aos_vc_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_vc_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_vx_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_vx_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_dvc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_dvc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_dvx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, aos_dvx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, grad_aos_dvc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, grad_aos_dvc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, grad_aos_dvx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
&BEGIN_PROVIDER[double precision, grad_aos_dvx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! aos_vxc_alpha_PBE_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
|
||||
END_DOC
|
||||
integer :: istate,i,j,m
|
||||
double precision :: r(3)
|
||||
double precision :: mu,weight
|
||||
double precision, allocatable :: ex(:), ec(:)
|
||||
double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
|
||||
double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
|
||||
double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
|
||||
double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
|
||||
allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
|
||||
allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
|
||||
|
||||
|
||||
allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
|
||||
allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
|
||||
allocate(contrib_grad_xa(3,N_states),contrib_grad_xb(3,N_states),contrib_grad_ca(3,N_states),contrib_grad_cb(3,N_states))
|
||||
do istate = 1, N_states
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
weight=final_weight_functions_at_final_grid_points(i)
|
||||
rho_a(istate) = one_body_dm_alpha_and_grad_at_r(4,i,istate)
|
||||
rho_b(istate) = one_body_dm_beta_and_grad_at_r(4,i,istate)
|
||||
grad_rho_a(1:3,istate) = one_body_dm_alpha_and_grad_at_r(1:3,i,istate)
|
||||
grad_rho_b(1:3,istate) = one_body_dm_beta_and_grad_at_r(1:3,i,istate)
|
||||
grad_rho_a_2 = 0.d0
|
||||
grad_rho_b_2 = 0.d0
|
||||
grad_rho_a_b = 0.d0
|
||||
do m = 1, 3
|
||||
grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
|
||||
grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
|
||||
grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
|
||||
enddo
|
||||
|
||||
! inputs
|
||||
call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
|
||||
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
|
||||
vx_rho_a(istate) *= weight
|
||||
vc_rho_a(istate) *= weight
|
||||
vx_rho_b(istate) *= weight
|
||||
vc_rho_b(istate) *= weight
|
||||
do m= 1,3
|
||||
contrib_grad_ca(m,istate) = weight * (2.d0 * vc_grad_rho_a_2(istate) * grad_rho_a(m,istate) + vc_grad_rho_a_b(istate) * grad_rho_b(m,istate))
|
||||
contrib_grad_xa(m,istate) = weight * (2.d0 * vx_grad_rho_a_2(istate) * grad_rho_a(m,istate) + vx_grad_rho_a_b(istate) * grad_rho_b(m,istate))
|
||||
contrib_grad_cb(m,istate) = weight * (2.d0 * vc_grad_rho_b_2(istate) * grad_rho_b(m,istate) + vc_grad_rho_a_b(istate) * grad_rho_a(m,istate))
|
||||
contrib_grad_xb(m,istate) = weight * (2.d0 * vx_grad_rho_b_2(istate) * grad_rho_b(m,istate) + vx_grad_rho_a_b(istate) * grad_rho_a(m,istate))
|
||||
enddo
|
||||
do j = 1, ao_num
|
||||
aos_vc_alpha_PBE_w(j,i,istate) = vc_rho_a(istate) * aos_in_r_array(j,i)
|
||||
aos_vc_beta_PBE_w (j,i,istate) = vc_rho_b(istate) * aos_in_r_array(j,i)
|
||||
aos_vx_alpha_PBE_w(j,i,istate) = vx_rho_a(istate) * aos_in_r_array(j,i)
|
||||
aos_vx_beta_PBE_w (j,i,istate) = vx_rho_b(istate) * aos_in_r_array(j,i)
|
||||
do m = 1,3
|
||||
aos_dvc_alpha_PBE_w(j,i,m,istate) = contrib_grad_ca(m,istate) * aos_in_r_array(j,i)
|
||||
aos_dvc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_in_r_array(j,i)
|
||||
aos_dvx_alpha_PBE_w(j,i,m,istate) = contrib_grad_xa(m,istate) * aos_in_r_array(j,i)
|
||||
aos_dvx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_in_r_array(j,i)
|
||||
|
||||
grad_aos_dvc_alpha_PBE_w (j,i,m,istate) = contrib_grad_ca(m,istate) * aos_grad_in_r_array(j,i,m)
|
||||
grad_aos_dvc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_grad_in_r_array(j,i,m)
|
||||
grad_aos_dvx_alpha_PBE_w (j,i,m,istate) = contrib_grad_xa(m,istate) * aos_grad_in_r_array(j,i,m)
|
||||
grad_aos_dvx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_grad_in_r_array(j,i,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER[double precision, energy_x_PBE, (N_states) ]
|
||||
&BEGIN_PROVIDER[double precision, energy_c_PBE, (N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! exchange/correlation energy with the short range PBE functional
|
||||
END_DOC
|
||||
integer :: istate,i,j,m
|
||||
double precision :: r(3)
|
||||
double precision :: mu,weight
|
||||
double precision, allocatable :: ex(:), ec(:)
|
||||
double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
|
||||
double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
|
||||
double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
|
||||
double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
|
||||
allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
|
||||
allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
|
||||
|
||||
|
||||
allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
|
||||
allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
|
||||
energy_x_PBE = 0.d0
|
||||
energy_c_PBE = 0.d0
|
||||
do istate = 1, N_states
|
||||
do i = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,i)
|
||||
r(2) = final_grid_points(2,i)
|
||||
r(3) = final_grid_points(3,i)
|
||||
weight=final_weight_functions_at_final_grid_points(i)
|
||||
rho_a(istate) = one_body_dm_alpha_and_grad_at_r(4,i,istate)
|
||||
rho_b(istate) = one_body_dm_beta_and_grad_at_r(4,i,istate)
|
||||
grad_rho_a(1:3,istate) = one_body_dm_alpha_and_grad_at_r(1:3,i,istate)
|
||||
grad_rho_b(1:3,istate) = one_body_dm_beta_and_grad_at_r(1:3,i,istate)
|
||||
grad_rho_a_2 = 0.d0
|
||||
grad_rho_b_2 = 0.d0
|
||||
grad_rho_a_b = 0.d0
|
||||
do m = 1, 3
|
||||
grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
|
||||
grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
|
||||
grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
|
||||
enddo
|
||||
|
||||
! inputs
|
||||
call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
|
||||
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
|
||||
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
|
||||
energy_x_PBE += ex * weight
|
||||
energy_c_PBE += ec * weight
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, potential_x_alpha_ao_PBE,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_x_beta_ao_PBE,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_PBE,(ao_num,ao_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, potential_c_beta_ao_PBE,(ao_num,ao_num,N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! exchange/correlation alpha/beta potentials with the short range PBE functional on the AO basis
|
||||
END_DOC
|
||||
integer :: istate, m
|
||||
double precision :: wall_1,wall_2
|
||||
call wall_time(wall_1)
|
||||
potential_c_alpha_ao_PBE = 0.d0
|
||||
potential_x_alpha_ao_PBE = 0.d0
|
||||
potential_c_beta_ao_PBE = 0.d0
|
||||
potential_x_beta_ao_PBE = 0.d0
|
||||
do istate = 1, N_states
|
||||
! correlation alpha
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vc_alpha_PBE_w(1,1,istate),size(aos_vc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1))
|
||||
! correlation beta
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vc_beta_PBE_w(1,1,istate),size(aos_vc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1))
|
||||
! exchange alpha
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vx_alpha_PBE_w(1,1,istate),size(aos_vx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1))
|
||||
! exchange beta
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vx_beta_PBE_w(1,1,istate),size(aos_vx_beta_PBE_w,1), aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate), size(potential_x_beta_ao_PBE,1))
|
||||
do m= 1,3
|
||||
! correlation alpha
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvc_alpha_PBE_w(1,1,m,istate),size(aos_dvc_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1))
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvc_alpha_PBE_w(1,1,m,istate),size(grad_aos_dvc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1))
|
||||
! correlation beta
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvc_beta_PBE_w(1,1,m,istate),size(aos_dvc_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1))
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvc_beta_PBE_w(1,1,m,istate),size(grad_aos_dvc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1))
|
||||
! exchange alpha
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvx_alpha_PBE_w(1,1,m,istate),size(aos_dvx_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1))
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvx_alpha_PBE_w(1,1,m,istate),size(grad_aos_dvx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1))
|
||||
! exchange beta
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvx_beta_PBE_w(1,1,m,istate),size(aos_dvx_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate),size(potential_x_beta_ao_PBE,1))
|
||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvx_beta_PBE_w(1,1,m,istate),size(grad_aos_dvx_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate),size(potential_x_beta_ao_PBE,1))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall_2)
|
||||
|
||||
END_PROVIDER
|
@ -0,0 +1,74 @@
|
||||
subroutine rho_ab_to_rho_oc(rho_a,rho_b,rho_o,rho_c)
|
||||
implicit none
|
||||
double precision, intent(in) :: rho_a,rho_b
|
||||
double precision, intent(out) :: rho_o,rho_c
|
||||
rho_c=rho_a+rho_b
|
||||
rho_o=rho_a-rho_b
|
||||
end
|
||||
|
||||
subroutine rho_oc_to_rho_ab(rho_o,rho_c,rho_a,rho_b)
|
||||
implicit none
|
||||
double precision, intent(in) :: rho_o,rho_c
|
||||
double precision, intent(out) :: rho_a,rho_b
|
||||
rho_a= 0.5d0*(rho_c+rho_o)
|
||||
rho_b= 0.5d0*(rho_c-rho_o)
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,grad_rho_o_2,grad_rho_c_2,grad_rho_o_c)
|
||||
implicit none
|
||||
double precision, intent(in) :: grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
|
||||
double precision, intent(out) :: grad_rho_o_2,grad_rho_c_2,grad_rho_o_c
|
||||
grad_rho_c_2 = grad_rho_a_2 + grad_rho_b_2 + 2d0*grad_rho_a_b
|
||||
grad_rho_o_2 = grad_rho_a_2 + grad_rho_b_2 - 2d0*grad_rho_a_b
|
||||
grad_rho_o_c = grad_rho_a_2 - grad_rho_b_2
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine v_rho_ab_to_v_rho_oc(v_rho_a,v_rho_b,v_rho_o,v_rho_c)
|
||||
implicit none
|
||||
double precision, intent(in) :: v_rho_a,v_rho_b
|
||||
double precision, intent(out) :: v_rho_o,v_rho_c
|
||||
v_rho_c = 0.5d0*(v_rho_a + v_rho_b)
|
||||
v_rho_o = 0.5d0*(v_rho_a - v_rho_b)
|
||||
end
|
||||
|
||||
subroutine v_rho_oc_to_v_rho_ab(v_rho_o,v_rho_c,v_rho_a,v_rho_b)
|
||||
implicit none
|
||||
double precision, intent(in) :: v_rho_o,v_rho_c
|
||||
double precision, intent(out) :: v_rho_a,v_rho_b
|
||||
v_rho_a = v_rho_c + v_rho_o
|
||||
v_rho_b = v_rho_c - v_rho_o
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine v_grad_rho_oc_to_v_grad_rho_ab(v_grad_rho_o_2,v_grad_rho_c_2,v_grad_rho_o_c,v_grad_rho_a_2,v_grad_rho_b_2,v_grad_rho_a_b)
|
||||
implicit none
|
||||
double precision, intent(in) :: v_grad_rho_o_2,v_grad_rho_c_2,v_grad_rho_o_c
|
||||
double precision, intent(out) :: v_grad_rho_a_2,v_grad_rho_b_2,v_grad_rho_a_b
|
||||
v_grad_rho_a_2 = v_grad_rho_o_2 + v_grad_rho_c_2 + v_grad_rho_o_c
|
||||
v_grad_rho_b_2 = v_grad_rho_o_2 + v_grad_rho_c_2 - v_grad_rho_o_c
|
||||
v_grad_rho_a_b = -2d0 * v_grad_rho_o_2 + 2d0 * v_grad_rho_c_2
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
16
src/dft_utils_one_body/shifted_potential.irp.f
Normal file
16
src/dft_utils_one_body/shifted_potential.irp.f
Normal file
@ -0,0 +1,16 @@
|
||||
BEGIN_PROVIDER [double precision, shifting_constant, (N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! shifting_constant = (E_{Hxc} - <\Psi | V_{Hxc} | \Psi>) / N_elec
|
||||
! constant to add to the potential in order to obtain the variational energy as
|
||||
! the eigenvalue of the effective long-range Hamiltonian
|
||||
! (see original paper of Levy PRL 113, 113002 (2014), equation (17) )
|
||||
END_DOC
|
||||
integer :: istate
|
||||
do istate = 1, N_states
|
||||
shifting_constant(istate) = energy_x(istate) + energy_c(istate) + short_range_Hartree(istate) - Trace_v_Hxc(istate)
|
||||
enddo
|
||||
shifting_constant = shifting_constant / dble(elec_num)
|
||||
|
||||
|
||||
END_PROVIDER
|
118
src/dft_utils_one_body/short_range_coulomb.irp.f
Normal file
118
src/dft_utils_one_body/short_range_coulomb.irp.f
Normal file
@ -0,0 +1,118 @@
|
||||
BEGIN_PROVIDER [double precision, short_range_Hartree_operator, (mo_tot_num,mo_tot_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, short_range_Hartree, (N_states)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! short_range_Hartree_operator(i,j) = \int dr i(r)j(r) \int r' \rho(r') W_{ee}^{sr}
|
||||
! short_range_Hartree = 0.5 * \sum_{i,j} \rho_{ij} short_range_Hartree_operator(i,j)
|
||||
! = 0.5 * \int dr \int r' \rho(r) \rho(r') W_{ee}^{sr}
|
||||
END_DOC
|
||||
integer :: i,j,k,l,m,n,istate
|
||||
double precision :: get_mo_bielec_integral,get_mo_bielec_integral_erf
|
||||
double precision :: integral, integral_erf, contrib
|
||||
double precision :: integrals_array(mo_tot_num,mo_tot_num),integrals_erf_array(mo_tot_num,mo_tot_num)
|
||||
short_range_Hartree_operator = 0.d0
|
||||
short_range_Hartree = 0.d0
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
if(dabs(one_body_dm_average_mo_for_dft(j,i)).le.1.d-12)cycle
|
||||
call get_mo_bielec_integrals_i1j1(i,j,mo_tot_num,integrals_array,mo_integrals_map)
|
||||
call get_mo_bielec_integrals_erf_i1j1(i,j,mo_tot_num,integrals_erf_array,mo_integrals_erf_map)
|
||||
do istate = 1, N_states
|
||||
do k = 1, mo_tot_num
|
||||
do l = 1, mo_tot_num
|
||||
integral = integrals_array(l,k)
|
||||
integral_erf = integrals_erf_array(l,k)
|
||||
contrib = one_body_dm_mo_for_dft(i,j,istate) * (integral - integral_erf)
|
||||
short_range_Hartree_operator(l,k,istate) += contrib
|
||||
short_range_Hartree(istate) += contrib * one_body_dm_mo_for_dft(k,l,istate)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
short_range_Hartree = short_range_Hartree * 0.5d0
|
||||
print*, 'short_range_Hartree',short_range_Hartree
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, effective_one_e_potential, (mo_tot_num, mo_tot_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, effective_one_e_potential_without_kin, (mo_tot_num, mo_tot_num,N_states)]
|
||||
&BEGIN_PROVIDER [double precision, shifted_effective_one_e_potential_without_kin, (mo_tot_num, mo_tot_num,N_states)]
|
||||
implicit none
|
||||
integer :: i,j,istate
|
||||
effective_one_e_potential = 0.d0
|
||||
BEGIN_DOC
|
||||
! effective_one_e_potential(i,j) = <i| v_{H}^{sr} |j> + <i| h_{core} |j> + <i|v_{xc} |j>
|
||||
! Taking the expectation value does not provide any energy
|
||||
! but effective_one_e_potential(i,j) is the potential coupling DFT and WFT part to be used in any WFT calculation
|
||||
! shifted_effective_one_e_potential_without_kin = effective_one_e_potential_without_kin + shifting_constant on the diagonal
|
||||
END_DOC
|
||||
do istate = 1, N_states
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
effective_one_e_potential(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) &
|
||||
+ 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate) &
|
||||
+ potential_x_beta_mo(i,j,istate) + potential_c_beta_mo(i,j,istate) )
|
||||
effective_one_e_potential_without_kin(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_nucl_elec_integral(i,j) &
|
||||
+ 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate) &
|
||||
+ potential_x_beta_mo(i,j,istate) + potential_c_beta_mo(i,j,istate) )
|
||||
shifted_effective_one_e_potential_without_kin(j,i,istate) = effective_one_e_potential_without_kin(j,i,istate)
|
||||
enddo
|
||||
enddo
|
||||
do i = 1, mo_tot_num
|
||||
shifted_effective_one_e_potential_without_kin(i,i,istate) += shifting_constant(istate)
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, Fock_matrix_expectation_value]
|
||||
implicit none
|
||||
call get_average(effective_one_e_potential,one_body_dm_average_mo_for_dft,Fock_matrix_expectation_value)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, Trace_v_xc, (N_states)]
|
||||
&BEGIN_PROVIDER [double precision, Trace_v_H, (N_states)]
|
||||
&BEGIN_PROVIDER [double precision, Trace_v_Hxc, (N_states)]
|
||||
implicit none
|
||||
integer :: i,j,istate
|
||||
double precision :: dm
|
||||
BEGIN_DOC
|
||||
! Trace_v_xc = \sum_{i,j} (rho_{ij}_\alpha v^{xc}_{ij}^\alpha + rho_{ij}_\beta v^{xc}_{ij}^\beta)
|
||||
! Trace_v_Hxc = \sum_{i,j} v^{H}_{ij} (rho_{ij}_\alpha + rho_{ij}_\beta)
|
||||
! Trace_v_Hxc = \sum_{i,j} rho_{ij} v^{Hxc}_{ij}
|
||||
END_DOC
|
||||
do istate = 1, N_states
|
||||
Trace_v_xc(istate) = 0.d0
|
||||
Trace_v_H(istate) = 0.d0
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
Trace_v_xc(istate) += (potential_x_alpha_mo(j,i,istate) + potential_c_alpha_mo(j,i,istate)) * one_body_dm_mo_alpha_for_dft(j,i,istate)
|
||||
Trace_v_xc(istate) += (potential_x_beta_mo(j,i,istate) + potential_c_beta_mo(j,i,istate) ) * one_body_dm_mo_beta_for_dft(j,i,istate)
|
||||
dm = one_body_dm_mo_alpha_for_dft(j,i,istate) + one_body_dm_mo_beta_for_dft(j,i,istate)
|
||||
Trace_v_H(istate) += dm * short_range_Hartree_operator(j,i,istate)
|
||||
enddo
|
||||
enddo
|
||||
Trace_v_Hxc(istate) = Trace_v_xc(istate) + Trace_v_H(istate)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, DFT_one_e_energy_potential, (mo_tot_num, mo_tot_num,N_states)]
|
||||
implicit none
|
||||
integer :: i,j,istate
|
||||
BEGIN_DOC
|
||||
! one_e_energy_potential(i,j) = <i|h_{core}|j> + \int dr i(r)j(r) \int r' \rho(r') W_{ee}^{sr}
|
||||
! If one take the expectation value over Psi, one gets the total one body energy
|
||||
END_DOC
|
||||
do istate = 1, N_states
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
DFT_one_e_energy_potential(j,i,istate) = mo_nucl_elec_integral(j,i) + mo_kinetic_integral(j,i) + short_range_Hartree_operator(j,i,istate) * 0.5d0
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
51
src/dft_utils_one_body/utils.irp.f
Normal file
51
src/dft_utils_one_body/utils.irp.f
Normal file
@ -0,0 +1,51 @@
|
||||
|
||||
subroutine GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, &
|
||||
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &
|
||||
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
|
||||
implicit none
|
||||
double precision, intent(in) :: r(3),rho_a(N_states),rho_b(N_states),grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
|
||||
double precision, intent(out) :: ex(N_states),vx_rho_a(N_states),vx_rho_b(N_states),vx_grad_rho_a_2(N_states),vx_grad_rho_b_2(N_states),vx_grad_rho_a_b(N_states)
|
||||
double precision, intent(out) :: ec(N_states),vc_rho_a(N_states),vc_rho_b(N_states),vc_grad_rho_a_2(N_states),vc_grad_rho_b_2(N_states),vc_grad_rho_a_b(N_states)
|
||||
integer :: istate
|
||||
double precision :: r2(3),dr2(3), local_potential,r12,dx2,mu
|
||||
do istate = 1, N_states
|
||||
if(exchange_functional.EQ."short_range_PBE")then
|
||||
call ex_pbe_sr(mu_erf,rho_a(istate),rho_b(istate),grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),ex(istate),vx_rho_a(istate),vx_rho_b(istate),vx_grad_rho_a_2(istate),vx_grad_rho_b_2(istate),vx_grad_rho_a_b(istate))
|
||||
else if(exchange_functional.EQ."None")then
|
||||
ex = 0.d0
|
||||
vx_rho_a = 0.d0
|
||||
vx_rho_b = 0.d0
|
||||
vx_grad_rho_a_2 = 0.d0
|
||||
vx_grad_rho_a_b = 0.d0
|
||||
vx_grad_rho_b_2 = 0.d0
|
||||
else
|
||||
print*, 'Exchange functional required does not exist ...'
|
||||
print*,'exchange_functional',exchange_functional
|
||||
stop
|
||||
endif
|
||||
|
||||
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
|
||||
if(correlation_functional.EQ."short_range_PBE")then
|
||||
! convertion from (alpha,beta) formalism to (closed, open) formalism
|
||||
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
|
||||
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
|
||||
|
||||
call ec_pbe_sr(mu_erf,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec(istate),vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
|
||||
|
||||
call v_rho_oc_to_v_rho_ab(vrhoo,vrhoc,vc_rho_a(istate),vc_rho_b(istate))
|
||||
call v_grad_rho_oc_to_v_grad_rho_ab(vsigmaoo,vsigmacc,vsigmaco,vc_grad_rho_a_2(istate),vc_grad_rho_b_2(istate),vc_grad_rho_a_b(istate))
|
||||
else if(correlation_functional.EQ."None")then
|
||||
ec = 0.d0
|
||||
vc_rho_a = 0.d0
|
||||
vc_rho_b = 0.d0
|
||||
vc_grad_rho_a_2 = 0.d0
|
||||
vc_grad_rho_a_b = 0.d0
|
||||
vc_grad_rho_b_2 = 0.d0
|
||||
else
|
||||
print*, 'Correlation functional required does not exist ...'
|
||||
print*, 'correlation_functional',correlation_functional
|
||||
stop
|
||||
endif
|
||||
enddo
|
||||
end
|
||||
|
25
src/integrals_bielec_erf/EZFIO.cfg
Normal file
25
src/integrals_bielec_erf/EZFIO.cfg
Normal file
@ -0,0 +1,25 @@
|
||||
[disk_access_ao_integrals_erf]
|
||||
type: Disk_access
|
||||
doc: Read/Write AO integrals with the long range interaction from/to disk [ Write | Read | None ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: None
|
||||
|
||||
|
||||
[disk_access_mo_integrals_erf]
|
||||
type: Disk_access
|
||||
doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: None
|
||||
|
||||
[disk_access_mo_integrals_sr]
|
||||
type: Disk_access
|
||||
doc: Read/Write MO integrals with the short range interaction from/to disk [ Write | Read | None ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: None
|
||||
|
||||
[disk_access_ao_integrals_sr]
|
||||
type: Disk_access
|
||||
doc: Read/Write AO integrals with the short range interaction from/to disk [ Write | Read | None ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: None
|
||||
|
5
src/integrals_bielec_erf/NEED
Normal file
5
src/integrals_bielec_erf/NEED
Normal file
@ -0,0 +1,5 @@
|
||||
pseudo
|
||||
bitmask
|
||||
zmq
|
||||
integrals_bielec
|
||||
dft_keywords
|
649
src/integrals_bielec_erf/ao_bi_integrals_erf.irp.f
Normal file
649
src/integrals_bielec_erf/ao_bi_integrals_erf.irp.f
Normal file
@ -0,0 +1,649 @@
|
||||
double precision function ao_bielec_integral_erf(i,j,k,l)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! integral of the AO basis <ik|jl> or (ij|kl)
|
||||
! i(r1) j(r1) 1/r12 k(r2) l(r2)
|
||||
END_DOC
|
||||
|
||||
integer,intent(in) :: i,j,k,l
|
||||
integer :: p,q,r,s
|
||||
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
|
||||
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
|
||||
double precision :: integral
|
||||
include 'Utils/constants.include.F'
|
||||
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
|
||||
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
|
||||
integer :: iorder_p(3), iorder_q(3)
|
||||
double precision :: ao_bielec_integral_schwartz_accel_erf
|
||||
|
||||
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
|
||||
ao_bielec_integral_erf = ao_bielec_integral_schwartz_accel_erf(i,j,k,l)
|
||||
return
|
||||
endif
|
||||
|
||||
dim1 = n_pt_max_integrals
|
||||
|
||||
num_i = ao_nucl(i)
|
||||
num_j = ao_nucl(j)
|
||||
num_k = ao_nucl(k)
|
||||
num_l = ao_nucl(l)
|
||||
ao_bielec_integral_erf = 0.d0
|
||||
|
||||
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
I_center(p) = nucl_coord(num_i,p)
|
||||
J_center(p) = nucl_coord(num_j,p)
|
||||
K_center(p) = nucl_coord(num_k,p)
|
||||
L_center(p) = nucl_coord(num_l,p)
|
||||
enddo
|
||||
|
||||
double precision :: coef1, coef2, coef3, coef4
|
||||
double precision :: p_inv,q_inv
|
||||
double precision :: general_primitive_integral_erf
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
|
||||
I_power,J_power,I_center,J_center,dim1)
|
||||
p_inv = 1.d0/pp
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
|
||||
K_power,L_power,K_center,L_center,dim1)
|
||||
q_inv = 1.d0/qq
|
||||
integral = general_primitive_integral_erf(dim1, &
|
||||
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
||||
ao_bielec_integral_erf = ao_bielec_integral_erf + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
|
||||
else
|
||||
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
enddo
|
||||
double precision :: ERI_erf
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
integral = ERI_erf( &
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
|
||||
I_power(1),J_power(1),K_power(1),L_power(1), &
|
||||
I_power(2),J_power(2),K_power(2),L_power(2), &
|
||||
I_power(3),J_power(3),K_power(3),L_power(3))
|
||||
ao_bielec_integral_erf = ao_bielec_integral_erf + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
double precision function ao_bielec_integral_schwartz_accel_erf(i,j,k,l)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! integral of the AO basis <ik|jl> or (ij|kl)
|
||||
! i(r1) j(r1) 1/r12 k(r2) l(r2)
|
||||
END_DOC
|
||||
integer,intent(in) :: i,j,k,l
|
||||
integer :: p,q,r,s
|
||||
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
|
||||
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
|
||||
double precision :: integral
|
||||
include 'Utils/constants.include.F'
|
||||
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
|
||||
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
|
||||
integer :: iorder_p(3), iorder_q(3)
|
||||
double precision, allocatable :: schwartz_kl(:,:)
|
||||
double precision :: schwartz_ij
|
||||
|
||||
dim1 = n_pt_max_integrals
|
||||
|
||||
num_i = ao_nucl(i)
|
||||
num_j = ao_nucl(j)
|
||||
num_k = ao_nucl(k)
|
||||
num_l = ao_nucl(l)
|
||||
ao_bielec_integral_schwartz_accel_erf = 0.d0
|
||||
double precision :: thr
|
||||
thr = ao_integrals_threshold*ao_integrals_threshold
|
||||
|
||||
allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k)))
|
||||
|
||||
double precision :: coef3
|
||||
double precision :: coef2
|
||||
double precision :: p_inv,q_inv
|
||||
double precision :: coef1
|
||||
double precision :: coef4
|
||||
|
||||
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
I_center(p) = nucl_coord(num_i,p)
|
||||
J_center(p) = nucl_coord(num_j,p)
|
||||
K_center(p) = nucl_coord(num_k,p)
|
||||
L_center(p) = nucl_coord(num_l,p)
|
||||
enddo
|
||||
|
||||
schwartz_kl(0,0) = 0.d0
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k)
|
||||
schwartz_kl(0,r) = 0.d0
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef2 = coef1 * ao_coef_normalized_ordered_transp(s,l) * ao_coef_normalized_ordered_transp(s,l)
|
||||
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
|
||||
K_power,L_power,K_center,L_center,dim1)
|
||||
q_inv = 1.d0/qq
|
||||
schwartz_kl(s,r) = general_primitive_integral_erf(dim1, &
|
||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q, &
|
||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q) &
|
||||
* coef2
|
||||
schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
|
||||
enddo
|
||||
schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0))
|
||||
enddo
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
|
||||
I_power,J_power,I_center,J_center,dim1)
|
||||
p_inv = 1.d0/pp
|
||||
schwartz_ij = general_primitive_integral_erf(dim1, &
|
||||
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||
P_new,P_center,fact_p,pp,p_inv,iorder_p) * &
|
||||
coef2*coef2
|
||||
if (schwartz_kl(0,0)*schwartz_ij < thr) then
|
||||
cycle
|
||||
endif
|
||||
do r = 1, ao_prim_num(k)
|
||||
if (schwartz_kl(0,r)*schwartz_ij < thr) then
|
||||
cycle
|
||||
endif
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
if (schwartz_kl(s,r)*schwartz_ij < thr) then
|
||||
cycle
|
||||
endif
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
double precision :: general_primitive_integral_erf
|
||||
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
|
||||
K_power,L_power,K_center,L_center,dim1)
|
||||
q_inv = 1.d0/qq
|
||||
integral = general_primitive_integral_erf(dim1, &
|
||||
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
||||
ao_bielec_integral_schwartz_accel_erf = ao_bielec_integral_schwartz_accel_erf + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
|
||||
else
|
||||
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
enddo
|
||||
double precision :: ERI_erf
|
||||
|
||||
schwartz_kl(0,0) = 0.d0
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef1 = ao_coef_normalized_ordered_transp(r,k)*ao_coef_normalized_ordered_transp(r,k)
|
||||
schwartz_kl(0,r) = 0.d0
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(s,l)*ao_coef_normalized_ordered_transp(s,l)
|
||||
schwartz_kl(s,r) = ERI_erf( &
|
||||
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
|
||||
K_power(1),L_power(1),K_power(1),L_power(1), &
|
||||
K_power(2),L_power(2),K_power(2),L_power(2), &
|
||||
K_power(3),L_power(3),K_power(3),L_power(3)) * &
|
||||
coef2
|
||||
schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
|
||||
enddo
|
||||
schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0))
|
||||
enddo
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
schwartz_ij = ERI_erf( &
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),&
|
||||
I_power(1),J_power(1),I_power(1),J_power(1), &
|
||||
I_power(2),J_power(2),I_power(2),J_power(2), &
|
||||
I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2
|
||||
if (schwartz_kl(0,0)*schwartz_ij < thr) then
|
||||
cycle
|
||||
endif
|
||||
do r = 1, ao_prim_num(k)
|
||||
if (schwartz_kl(0,r)*schwartz_ij < thr) then
|
||||
cycle
|
||||
endif
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
if (schwartz_kl(s,r)*schwartz_ij < thr) then
|
||||
cycle
|
||||
endif
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
integral = ERI_erf( &
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
|
||||
I_power(1),J_power(1),K_power(1),L_power(1), &
|
||||
I_power(2),J_power(2),K_power(2),L_power(2), &
|
||||
I_power(3),J_power(3),K_power(3),L_power(3))
|
||||
ao_bielec_integral_schwartz_accel_erf = ao_bielec_integral_schwartz_accel_erf + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
|
||||
endif
|
||||
deallocate (schwartz_kl)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine compute_ao_bielec_integrals_erf(j,k,l,sze,buffer_value)
|
||||
implicit none
|
||||
use map_module
|
||||
|
||||
BEGIN_DOC
|
||||
! Compute AO 1/r12 integrals for all i and fixed j,k,l
|
||||
END_DOC
|
||||
|
||||
include 'Utils/constants.include.F'
|
||||
integer, intent(in) :: j,k,l,sze
|
||||
real(integral_kind), intent(out) :: buffer_value(sze)
|
||||
double precision :: ao_bielec_integral_erf
|
||||
|
||||
integer :: i
|
||||
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
buffer_value = 0._integral_kind
|
||||
return
|
||||
endif
|
||||
if (ao_bielec_integral_erf_schwartz(j,l) < thresh ) then
|
||||
buffer_value = 0._integral_kind
|
||||
return
|
||||
endif
|
||||
|
||||
do i = 1, ao_num
|
||||
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
|
||||
buffer_value(i) = 0._integral_kind
|
||||
cycle
|
||||
endif
|
||||
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh ) then
|
||||
buffer_value(i) = 0._integral_kind
|
||||
cycle
|
||||
endif
|
||||
!DIR$ FORCEINLINE
|
||||
buffer_value(i) = ao_bielec_integral_erf(i,k,j,l)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
double precision function general_primitive_integral_erf(dim, &
|
||||
P_new,P_center,fact_p,p,p_inv,iorder_p, &
|
||||
Q_new,Q_center,fact_q,q,q_inv,iorder_q)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes the integral <pq|rs> where p,q,r,s are Gaussian primitives
|
||||
END_DOC
|
||||
integer,intent(in) :: dim
|
||||
include 'Utils/constants.include.F'
|
||||
double precision, intent(in) :: P_new(0:max_dim,3),P_center(3),fact_p,p,p_inv
|
||||
double precision, intent(in) :: Q_new(0:max_dim,3),Q_center(3),fact_q,q,q_inv
|
||||
integer, intent(in) :: iorder_p(3)
|
||||
integer, intent(in) :: iorder_q(3)
|
||||
|
||||
double precision :: r_cut,gama_r_cut,rho,dist
|
||||
double precision :: dx(0:max_dim),Ix_pol(0:max_dim),dy(0:max_dim),Iy_pol(0:max_dim),dz(0:max_dim),Iz_pol(0:max_dim)
|
||||
integer :: n_Ix,n_Iy,n_Iz,nx,ny,nz
|
||||
double precision :: bla
|
||||
integer :: ix,iy,iz,jx,jy,jz,i
|
||||
double precision :: a,b,c,d,e,f,accu,pq,const
|
||||
double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2
|
||||
integer :: n_pt_tmp,n_pt_out, iorder
|
||||
double precision :: d1(0:max_dim),d_poly(0:max_dim),rint,d1_screened(0:max_dim)
|
||||
|
||||
general_primitive_integral_erf = 0.d0
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx,Ix_pol,dy,Iy_pol,dz,Iz_pol
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: d1, d_poly
|
||||
|
||||
! Gaussian Product
|
||||
! ----------------
|
||||
double precision :: p_plus_q
|
||||
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
|
||||
pq = p_inv*0.5d0*q_inv
|
||||
|
||||
pq_inv = 0.5d0/p_plus_q
|
||||
p10_1 = q*pq ! 1/(2p)
|
||||
p01_1 = p*pq ! 1/(2q)
|
||||
pq_inv_2 = pq_inv+pq_inv
|
||||
p10_2 = pq_inv_2 * p10_1*q !0.5d0*q/(pq + p*p)
|
||||
p01_2 = pq_inv_2 * p01_1*p !0.5d0*p/(q*q + pq)
|
||||
|
||||
|
||||
accu = 0.d0
|
||||
iorder = iorder_p(1)+iorder_q(1)+iorder_p(1)+iorder_q(1)
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do ix=0,iorder
|
||||
Ix_pol(ix) = 0.d0
|
||||
enddo
|
||||
n_Ix = 0
|
||||
do ix = 0, iorder_p(1)
|
||||
if (abs(P_new(ix,1)) < thresh) cycle
|
||||
a = P_new(ix,1)
|
||||
do jx = 0, iorder_q(1)
|
||||
d = a*Q_new(jx,1)
|
||||
if (abs(d) < thresh) cycle
|
||||
!DEC$ FORCEINLINE
|
||||
call give_polynom_mult_center_x(P_center(1),Q_center(1),ix,jx,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dx,nx)
|
||||
!DEC$ FORCEINLINE
|
||||
call add_poly_multiply(dx,nx,d,Ix_pol,n_Ix)
|
||||
enddo
|
||||
enddo
|
||||
if (n_Ix == -1) then
|
||||
return
|
||||
endif
|
||||
iorder = iorder_p(2)+iorder_q(2)+iorder_p(2)+iorder_q(2)
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do ix=0, iorder
|
||||
Iy_pol(ix) = 0.d0
|
||||
enddo
|
||||
n_Iy = 0
|
||||
do iy = 0, iorder_p(2)
|
||||
if (abs(P_new(iy,2)) > thresh) then
|
||||
b = P_new(iy,2)
|
||||
do jy = 0, iorder_q(2)
|
||||
e = b*Q_new(jy,2)
|
||||
if (abs(e) < thresh) cycle
|
||||
!DEC$ FORCEINLINE
|
||||
call give_polynom_mult_center_x(P_center(2),Q_center(2),iy,jy,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dy,ny)
|
||||
!DEC$ FORCEINLINE
|
||||
call add_poly_multiply(dy,ny,e,Iy_pol,n_Iy)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
if (n_Iy == -1) then
|
||||
return
|
||||
endif
|
||||
|
||||
iorder = iorder_p(3)+iorder_q(3)+iorder_p(3)+iorder_q(3)
|
||||
do ix=0,iorder
|
||||
Iz_pol(ix) = 0.d0
|
||||
enddo
|
||||
n_Iz = 0
|
||||
do iz = 0, iorder_p(3)
|
||||
if (abs(P_new(iz,3)) > thresh) then
|
||||
c = P_new(iz,3)
|
||||
do jz = 0, iorder_q(3)
|
||||
f = c*Q_new(jz,3)
|
||||
if (abs(f) < thresh) cycle
|
||||
!DEC$ FORCEINLINE
|
||||
call give_polynom_mult_center_x(P_center(3),Q_center(3),iz,jz,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dz,nz)
|
||||
!DEC$ FORCEINLINE
|
||||
call add_poly_multiply(dz,nz,f,Iz_pol,n_Iz)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
if (n_Iz == -1) then
|
||||
return
|
||||
endif
|
||||
|
||||
rho = p*q *pq_inv_2 ! le rho qui va bien
|
||||
dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + &
|
||||
(P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + &
|
||||
(P_center(3) - Q_center(3))*(P_center(3) - Q_center(3))
|
||||
const = dist*rho
|
||||
|
||||
n_pt_tmp = n_Ix+n_Iy
|
||||
do i=0,n_pt_tmp
|
||||
d_poly(i)=0.d0
|
||||
enddo
|
||||
|
||||
!DEC$ FORCEINLINE
|
||||
call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
|
||||
if (n_pt_tmp == -1) then
|
||||
return
|
||||
endif
|
||||
n_pt_out = n_pt_tmp+n_Iz
|
||||
do i=0,n_pt_out
|
||||
d1(i)=0.d0
|
||||
enddo
|
||||
|
||||
!DEC$ FORCEINLINE
|
||||
call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
|
||||
double precision :: rint_sum
|
||||
accu = accu + rint_sum(n_pt_out,const,d1)
|
||||
|
||||
! change p+q in dsqrt
|
||||
general_primitive_integral_erf = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p_plus_q)
|
||||
end
|
||||
|
||||
|
||||
double precision function ERI_erf(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! ATOMIC PRIMTIVE bielectronic integral between the 4 primitives ::
|
||||
! primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2)
|
||||
! primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2)
|
||||
! primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2)
|
||||
! primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2)
|
||||
END_DOC
|
||||
double precision, intent(in) :: delta,gama,alpha,beta
|
||||
integer, intent(in) :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
|
||||
integer :: a_x_2,b_x_2,c_x_2,d_x_2,a_y_2,b_y_2,c_y_2,d_y_2,a_z_2,b_z_2,c_z_2,d_z_2
|
||||
integer :: i,j,k,l,n_pt
|
||||
integer :: n_pt_sup
|
||||
double precision :: p,q,denom,coeff
|
||||
double precision :: I_f
|
||||
integer :: nx,ny,nz
|
||||
include 'Utils/constants.include.F'
|
||||
nx = a_x+b_x+c_x+d_x
|
||||
if(iand(nx,1) == 1) then
|
||||
ERI_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
ny = a_y+b_y+c_y+d_y
|
||||
if(iand(ny,1) == 1) then
|
||||
ERI_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
nz = a_z+b_z+c_z+d_z
|
||||
if(iand(nz,1) == 1) then
|
||||
ERI_erf = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
ASSERT (alpha >= 0.d0)
|
||||
ASSERT (beta >= 0.d0)
|
||||
ASSERT (delta >= 0.d0)
|
||||
ASSERT (gama >= 0.d0)
|
||||
p = alpha + beta
|
||||
q = delta + gama
|
||||
double precision :: p_plus_q
|
||||
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
|
||||
ASSERT (p+q >= 0.d0)
|
||||
n_pt = ishft( nx+ny+nz,1 )
|
||||
|
||||
coeff = pi_5_2 / (p * q * dsqrt(p_plus_q))
|
||||
if (n_pt == 0) then
|
||||
ERI_erf = coeff
|
||||
return
|
||||
endif
|
||||
|
||||
call integrale_new_erf(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt)
|
||||
|
||||
ERI_erf = I_f * coeff
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine integrale_new_erf(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt)
|
||||
BEGIN_DOC
|
||||
! calculate the integral of the polynom ::
|
||||
! I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q)
|
||||
! between ( 0 ; 1)
|
||||
END_DOC
|
||||
|
||||
|
||||
implicit none
|
||||
include 'Utils/constants.include.F'
|
||||
double precision :: p,q
|
||||
integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
|
||||
integer :: i, n_iter, n_pt, j
|
||||
double precision :: I_f, pq_inv, p10_1, p10_2, p01_1, p01_2,rho,pq_inv_2
|
||||
integer :: ix,iy,iz, jx,jy,jz, sx,sy,sz
|
||||
|
||||
j = ishft(n_pt,-1)
|
||||
ASSERT (n_pt > 1)
|
||||
double precision :: p_plus_q
|
||||
p_plus_q = (p+q) * ((p*q)/(p+q) + mu_erf*mu_erf)/(mu_erf*mu_erf)
|
||||
|
||||
pq_inv = 0.5d0/(p_plus_q)
|
||||
pq_inv_2 = pq_inv + pq_inv
|
||||
p10_1 = 0.5d0/p
|
||||
p01_1 = 0.5d0/q
|
||||
p10_2 = 0.5d0 * q /(p * p_plus_q)
|
||||
p01_2 = 0.5d0 * p /(q * p_plus_q)
|
||||
double precision :: B00(n_pt_max_integrals)
|
||||
double precision :: B10(n_pt_max_integrals), B01(n_pt_max_integrals)
|
||||
double precision :: t1(n_pt_max_integrals), t2(n_pt_max_integrals)
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00
|
||||
ix = a_x+b_x
|
||||
jx = c_x+d_x
|
||||
iy = a_y+b_y
|
||||
jy = c_y+d_y
|
||||
iz = a_z+b_z
|
||||
jz = c_z+d_z
|
||||
sx = ix+jx
|
||||
sy = iy+jy
|
||||
sz = iz+jz
|
||||
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i = 1,n_pt
|
||||
B10(i) = p10_1 - gauleg_t2(i,j)* p10_2
|
||||
B01(i) = p01_1 - gauleg_t2(i,j)* p01_2
|
||||
B00(i) = gauleg_t2(i,j)*pq_inv
|
||||
enddo
|
||||
if (sx > 0) then
|
||||
call I_x1_new(ix,jx,B10,B01,B00,t1,n_pt)
|
||||
else
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i = 1,n_pt
|
||||
t1(i) = 1.d0
|
||||
enddo
|
||||
endif
|
||||
if (sy > 0) then
|
||||
call I_x1_new(iy,jy,B10,B01,B00,t2,n_pt)
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i = 1,n_pt
|
||||
t1(i) = t1(i)*t2(i)
|
||||
enddo
|
||||
endif
|
||||
if (sz > 0) then
|
||||
call I_x1_new(iz,jz,B10,B01,B00,t2,n_pt)
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i = 1,n_pt
|
||||
t1(i) = t1(i)*t2(i)
|
||||
enddo
|
||||
endif
|
||||
I_f= 0.d0
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i = 1,n_pt
|
||||
I_f += gauleg_w(i,j)*t1(i)
|
||||
enddo
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value)
|
||||
implicit none
|
||||
use map_module
|
||||
BEGIN_DOC
|
||||
! Parallel client for AO integrals
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: j,l
|
||||
integer,intent(out) :: n_integrals
|
||||
integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num)
|
||||
real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num)
|
||||
|
||||
integer :: i,k
|
||||
double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2
|
||||
double precision :: integral, wall_0
|
||||
double precision :: thr
|
||||
integer :: kk, m, j1, i1
|
||||
|
||||
thr = ao_integrals_threshold
|
||||
|
||||
n_integrals = 0
|
||||
|
||||
j1 = j+ishft(l*l-l,-1)
|
||||
do k = 1, ao_num ! r1
|
||||
i1 = ishft(k*k-k,-1)
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
do i = 1, k
|
||||
i1 += 1
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thr) then
|
||||
cycle
|
||||
endif
|
||||
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thr ) then
|
||||
cycle
|
||||
endif
|
||||
!DIR$ FORCEINLINE
|
||||
integral = ao_bielec_integral_erf(i,k,j,l) ! i,k : r1 j,l : r2
|
||||
if (abs(integral) < thr) then
|
||||
cycle
|
||||
endif
|
||||
n_integrals += 1
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
|
||||
buffer_value(n_integrals) = integral
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
@ -0,0 +1,194 @@
|
||||
subroutine ao_bielec_integrals_erf_in_map_slave_tcp(i)
|
||||
implicit none
|
||||
integer, intent(in) :: i
|
||||
BEGIN_DOC
|
||||
! Computes a buffer of integrals. i is the ID of the current thread.
|
||||
END_DOC
|
||||
call ao_bielec_integrals_erf_in_map_slave(0,i)
|
||||
end
|
||||
|
||||
|
||||
subroutine ao_bielec_integrals_erf_in_map_slave_inproc(i)
|
||||
implicit none
|
||||
integer, intent(in) :: i
|
||||
BEGIN_DOC
|
||||
! Computes a buffer of integrals. i is the ID of the current thread.
|
||||
END_DOC
|
||||
call ao_bielec_integrals_erf_in_map_slave(1,i)
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine ao_bielec_integrals_erf_in_map_slave(thread,iproc)
|
||||
use map_module
|
||||
use f77_zmq
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes a buffer of integrals
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: thread, iproc
|
||||
|
||||
integer :: j,l,n_integrals
|
||||
integer :: rc
|
||||
real(integral_kind), allocatable :: buffer_value(:)
|
||||
integer(key_kind), allocatable :: buffer_i(:)
|
||||
|
||||
integer :: worker_id, task_id
|
||||
character*(512) :: task
|
||||
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_push_socket
|
||||
integer(ZMQ_PTR) :: zmq_socket_push
|
||||
|
||||
character*(64) :: state
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
|
||||
integer, external :: connect_to_taskserver
|
||||
if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
return
|
||||
endif
|
||||
|
||||
zmq_socket_push = new_zmq_push_socket(thread)
|
||||
|
||||
allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
|
||||
|
||||
|
||||
do
|
||||
integer, external :: get_task_from_taskserver
|
||||
if (get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task) == -1) then
|
||||
exit
|
||||
endif
|
||||
if (task_id == 0) exit
|
||||
read(task,*) j, l
|
||||
integer, external :: task_done_to_taskserver
|
||||
call compute_ao_integrals_erf_jl(j,l,n_integrals,buffer_i,buffer_value)
|
||||
if (task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id) == -1) then
|
||||
stop 'Unable to send task_done'
|
||||
endif
|
||||
call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
|
||||
enddo
|
||||
|
||||
integer, external :: disconnect_from_taskserver
|
||||
if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then
|
||||
continue
|
||||
endif
|
||||
deallocate( buffer_i, buffer_value )
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call end_zmq_push_socket(zmq_socket_push,thread)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine ao_bielec_integrals_erf_in_map_collector(zmq_socket_pull)
|
||||
use map_module
|
||||
use f77_zmq
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Collects results from the AO integral calculation
|
||||
END_DOC
|
||||
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
|
||||
integer :: j,l,n_integrals
|
||||
integer :: rc
|
||||
|
||||
real(integral_kind), allocatable :: buffer_value(:)
|
||||
integer(key_kind), allocatable :: buffer_i(:)
|
||||
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_pull_socket
|
||||
|
||||
integer*8 :: control, accu, sze
|
||||
integer :: task_id, more
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
|
||||
sze = ao_num*ao_num
|
||||
allocate ( buffer_i(sze), buffer_value(sze) )
|
||||
|
||||
accu = 0_8
|
||||
more = 1
|
||||
do while (more == 1)
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)
|
||||
if (rc == -1) then
|
||||
n_integrals = 0
|
||||
return
|
||||
endif
|
||||
if (rc /= 4) then
|
||||
print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, n_integrals, 4, 0)'
|
||||
stop 'error'
|
||||
endif
|
||||
|
||||
if (n_integrals >= 0) then
|
||||
|
||||
if (n_integrals > sze) then
|
||||
deallocate (buffer_value, buffer_i)
|
||||
sze = n_integrals
|
||||
allocate (buffer_value(sze), buffer_i(sze))
|
||||
endif
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)
|
||||
if (rc /= key_kind*n_integrals) then
|
||||
print *, rc, key_kind, n_integrals
|
||||
print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_i, key_kind*n_integrals, 0)'
|
||||
stop 'error'
|
||||
endif
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)
|
||||
if (rc /= integral_kind*n_integrals) then
|
||||
print *, irp_here, ': f77_zmq_recv( zmq_socket_pull, buffer_value, integral_kind*n_integrals, 0)'
|
||||
stop 'error'
|
||||
endif
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
|
||||
|
||||
IRP_IF ZMQ_PUSH
|
||||
IRP_ELSE
|
||||
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
|
||||
if (rc /= 4) then
|
||||
print *, irp_here, ' : f77_zmq_send (zmq_socket_pull,...'
|
||||
stop 'error'
|
||||
endif
|
||||
IRP_ENDIF
|
||||
|
||||
|
||||
call insert_into_ao_integrals_erf_map(n_integrals,buffer_i,buffer_value)
|
||||
accu += n_integrals
|
||||
if (task_id /= 0) then
|
||||
integer, external :: zmq_delete_task
|
||||
if (zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) == -1) then
|
||||
stop 'Unable to delete task'
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
deallocate( buffer_i, buffer_value )
|
||||
|
||||
integer (map_size_kind) :: get_ao_erf_map_size
|
||||
control = get_ao_erf_map_size(ao_integrals_erf_map)
|
||||
|
||||
if (control /= accu) then
|
||||
print *, ''
|
||||
print *, irp_here
|
||||
print *, 'Control : ', control
|
||||
print *, 'Accu : ', accu
|
||||
print *, 'Some integrals were lost during the parallel computation.'
|
||||
print *, 'Try to reduce the number of threads.'
|
||||
stop
|
||||
endif
|
||||
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
|
||||
end
|
||||
|
47
src/integrals_bielec_erf/integrals_3_index_erf.irp.f
Normal file
47
src/integrals_bielec_erf/integrals_3_index_erf.irp.f
Normal file
@ -0,0 +1,47 @@
|
||||
BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_erf, (mo_tot_num,mo_tot_num, mo_tot_num)]
|
||||
&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_erf,(mo_tot_num,mo_tot_num, mo_tot_num)]
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
double precision :: get_mo_bielec_integral_erf
|
||||
double precision :: integral
|
||||
|
||||
do k = 1, mo_tot_num
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
l = j
|
||||
integral = get_mo_bielec_integral_erf(i,j,k,l,mo_integrals_erf_map)
|
||||
big_array_coulomb_integrals_erf(j,i,k) = integral
|
||||
l = j
|
||||
integral = get_mo_bielec_integral_erf(i,j,l,k,mo_integrals_erf_map)
|
||||
big_array_exchange_integrals_erf(j,i,k) = integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
!BEGIN_PROVIDER [double precision, big_array_coulomb_integrals_sr, (mo_tot_num,mo_tot_num, mo_tot_num)]
|
||||
!&BEGIN_PROVIDER [double precision, big_array_exchange_integrals_sr,(mo_tot_num,mo_tot_num, mo_tot_num)]
|
||||
!implicit none
|
||||
!integer :: i,j,k,l
|
||||
!double precision :: get_mo_bielec_integral_sr
|
||||
!double precision :: integral
|
||||
|
||||
!do k = 1, mo_tot_num
|
||||
! do i = 1, mo_tot_num
|
||||
! do j = 1, mo_tot_num
|
||||
! l = j
|
||||
! integral = get_mo_bielec_integral_sr(i,j,k,l,mo_integrals_sr_map)
|
||||
! big_array_coulomb_integrals_sr(j,i,k) = integral
|
||||
! l = j
|
||||
! integral = get_mo_bielec_integral_sr(i,j,l,k,mo_integrals_sr_map)
|
||||
! big_array_exchange_integrals_sr(j,i,k) = integral
|
||||
! enddo
|
||||
! enddo
|
||||
!enddo
|
||||
|
||||
|
||||
!ND_PROVIDER
|
||||
|
682
src/integrals_bielec_erf/map_integrals_erf.irp.f
Normal file
682
src/integrals_bielec_erf/map_integrals_erf.irp.f
Normal file
@ -0,0 +1,682 @@
|
||||
use map_module
|
||||
|
||||
!! AO Map
|
||||
!! ======
|
||||
|
||||
BEGIN_PROVIDER [ type(map_type), ao_integrals_erf_map ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! AO integrals
|
||||
END_DOC
|
||||
integer(key_kind) :: key_max
|
||||
integer(map_size_kind) :: sze
|
||||
call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,key_max)
|
||||
sze = key_max
|
||||
call map_init(ao_integrals_erf_map,sze)
|
||||
print*, 'AO map initialized : ', sze
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_min ]
|
||||
&BEGIN_PROVIDER [ integer, ao_integrals_erf_cache_max ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Min and max values of the AOs for which the integrals are in the cache
|
||||
END_DOC
|
||||
ao_integrals_erf_cache_min = max(1,ao_num - 63)
|
||||
ao_integrals_erf_cache_max = ao_num
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_integrals_erf_cache, (0:64*64*64*64) ]
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Cache of AO integrals for fast access
|
||||
END_DOC
|
||||
PROVIDE ao_bielec_integrals_erf_in_map
|
||||
integer :: i,j,k,l,ii
|
||||
integer(key_kind) :: idx
|
||||
real(integral_kind) :: integral
|
||||
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
|
||||
do l=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
|
||||
do k=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
|
||||
do j=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
|
||||
do i=ao_integrals_erf_cache_min,ao_integrals_erf_cache_max
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,idx)
|
||||
!DIR$ FORCEINLINE
|
||||
call map_get(ao_integrals_erf_map,idx,integral)
|
||||
ii = l-ao_integrals_erf_cache_min
|
||||
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
|
||||
ao_integrals_erf_cache(ii) = integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
double precision function get_ao_bielec_integral_erf(i,j,k,l,map) result(result)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Gets one AO bi-electronic integral from the AO map
|
||||
END_DOC
|
||||
integer, intent(in) :: i,j,k,l
|
||||
integer(key_kind) :: idx
|
||||
type(map_type), intent(inout) :: map
|
||||
integer :: ii
|
||||
real(integral_kind) :: tmp
|
||||
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_cache ao_integrals_erf_cache_min
|
||||
!DIR$ FORCEINLINE
|
||||
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
|
||||
tmp = 0.d0
|
||||
else if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < ao_integrals_threshold) then
|
||||
tmp = 0.d0
|
||||
else
|
||||
ii = l-ao_integrals_erf_cache_min
|
||||
ii = ior(ii, k-ao_integrals_erf_cache_min)
|
||||
ii = ior(ii, j-ao_integrals_erf_cache_min)
|
||||
ii = ior(ii, i-ao_integrals_erf_cache_min)
|
||||
if (iand(ii, -64) /= 0) then
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,idx)
|
||||
!DIR$ FORCEINLINE
|
||||
call map_get(map,idx,tmp)
|
||||
tmp = tmp
|
||||
else
|
||||
ii = l-ao_integrals_erf_cache_min
|
||||
ii = ior( ishft(ii,6), k-ao_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), j-ao_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), i-ao_integrals_erf_cache_min)
|
||||
tmp = ao_integrals_erf_cache(ii)
|
||||
endif
|
||||
endif
|
||||
result = tmp
|
||||
end
|
||||
|
||||
|
||||
subroutine get_ao_bielec_integrals_erf(j,k,l,sze,out_val)
|
||||
use map_module
|
||||
BEGIN_DOC
|
||||
! Gets multiple AO bi-electronic integral from the AO map .
|
||||
! All i are retrieved for j,k,l fixed.
|
||||
END_DOC
|
||||
implicit none
|
||||
integer, intent(in) :: j,k,l, sze
|
||||
real(integral_kind), intent(out) :: out_val(sze)
|
||||
|
||||
integer :: i
|
||||
integer(key_kind) :: hash
|
||||
double precision :: thresh
|
||||
PROVIDE ao_bielec_integrals_erf_in_map ao_integrals_erf_map
|
||||
thresh = ao_integrals_threshold
|
||||
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
out_val = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
double precision :: get_ao_bielec_integral_erf
|
||||
do i=1,sze
|
||||
out_val(i) = get_ao_bielec_integral_erf(i,j,k,l,ao_integrals_erf_map)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine get_ao_bielec_integrals_erf_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Gets multiple AO bi-electronic integral from the AO map .
|
||||
! All non-zero i are retrieved for j,k,l fixed.
|
||||
END_DOC
|
||||
integer, intent(in) :: j,k,l, sze
|
||||
real(integral_kind), intent(out) :: out_val(sze)
|
||||
integer, intent(out) :: out_val_index(sze),non_zero_int
|
||||
|
||||
integer :: i
|
||||
integer(key_kind) :: hash
|
||||
double precision :: thresh,tmp
|
||||
PROVIDE ao_bielec_integrals_erf_in_map
|
||||
thresh = ao_integrals_threshold
|
||||
|
||||
non_zero_int = 0
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
out_val = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
non_zero_int = 0
|
||||
do i=1,sze
|
||||
integer, external :: ao_l4
|
||||
double precision, external :: ao_bielec_integral_erf
|
||||
!DIR$ FORCEINLINE
|
||||
if (ao_bielec_integral_erf_schwartz(i,k)*ao_bielec_integral_erf_schwartz(j,l) < thresh) then
|
||||
cycle
|
||||
endif
|
||||
call bielec_integrals_index(i,j,k,l,hash)
|
||||
call map_get(ao_integrals_erf_map, hash,tmp)
|
||||
if (dabs(tmp) < thresh ) cycle
|
||||
non_zero_int = non_zero_int+1
|
||||
out_val_index(non_zero_int) = i
|
||||
out_val(non_zero_int) = tmp
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
function get_ao_erf_map_size()
|
||||
implicit none
|
||||
integer (map_size_kind) :: get_ao_erf_map_size
|
||||
BEGIN_DOC
|
||||
! Returns the number of elements in the AO map
|
||||
END_DOC
|
||||
get_ao_erf_map_size = ao_integrals_erf_map % n_elements
|
||||
end
|
||||
|
||||
subroutine clear_ao_erf_map
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Frees the memory of the AO map
|
||||
END_DOC
|
||||
call map_deinit(ao_integrals_erf_map)
|
||||
FREE ao_integrals_erf_map
|
||||
end
|
||||
|
||||
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
subroutine dump_$ao_integrals(filename)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Save to disk the $ao integrals
|
||||
END_DOC
|
||||
character*(*), intent(in) :: filename
|
||||
integer(cache_key_kind), pointer :: key(:)
|
||||
real(integral_kind), pointer :: val(:)
|
||||
integer*8 :: i,j, n
|
||||
call ezfio_set_work_empty(.False.)
|
||||
open(unit=66,file=filename,FORM='unformatted')
|
||||
write(66) integral_kind, key_kind
|
||||
write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size, &
|
||||
$ao_integrals_map%n_elements
|
||||
do i=0_8,$ao_integrals_map%map_size
|
||||
write(66) $ao_integrals_map%map(i)%sorted, $ao_integrals_map%map(i)%map_size,&
|
||||
$ao_integrals_map%map(i)%n_elements
|
||||
enddo
|
||||
do i=0_8,$ao_integrals_map%map_size
|
||||
key => $ao_integrals_map%map(i)%key
|
||||
val => $ao_integrals_map%map(i)%value
|
||||
n = $ao_integrals_map%map(i)%n_elements
|
||||
write(66) (key(j), j=1,n), (val(j), j=1,n)
|
||||
enddo
|
||||
close(66)
|
||||
|
||||
end
|
||||
|
||||
IRP_IF COARRAY
|
||||
subroutine communicate_$ao_integrals()
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Communicate the $ao integrals with co-array
|
||||
END_DOC
|
||||
integer(cache_key_kind), pointer :: key(:)
|
||||
real(integral_kind), pointer :: val(:)
|
||||
integer*8 :: i,j, k, nmax
|
||||
integer*8, save :: n[*]
|
||||
integer :: copy_n
|
||||
|
||||
real(integral_kind), allocatable :: buffer_val(:)[:]
|
||||
integer(cache_key_kind), allocatable :: buffer_key(:)[:]
|
||||
real(integral_kind), allocatable :: copy_val(:)
|
||||
integer(key_kind), allocatable :: copy_key(:)
|
||||
|
||||
n = 0_8
|
||||
do i=0_8,$ao_integrals_map%map_size
|
||||
n = max(n,$ao_integrals_map%map(i)%n_elements)
|
||||
enddo
|
||||
sync all
|
||||
nmax = 0_8
|
||||
do j=1,num_images()
|
||||
nmax = max(nmax,n[j])
|
||||
enddo
|
||||
allocate( buffer_key(nmax)[*], buffer_val(nmax)[*])
|
||||
allocate( copy_key(nmax), copy_val(nmax))
|
||||
do i=0_8,$ao_integrals_map%map_size
|
||||
key => $ao_integrals_map%map(i)%key
|
||||
val => $ao_integrals_map%map(i)%value
|
||||
n = $ao_integrals_map%map(i)%n_elements
|
||||
do j=1,n
|
||||
buffer_key(j) = key(j)
|
||||
buffer_val(j) = val(j)
|
||||
enddo
|
||||
sync all
|
||||
do j=1,num_images()
|
||||
if (j /= this_image()) then
|
||||
copy_n = n[j]
|
||||
do k=1,copy_n
|
||||
copy_val(k) = buffer_val(k)[j]
|
||||
copy_key(k) = buffer_key(k)[j]
|
||||
copy_key(k) = copy_key(k)+ishft(i,-map_shift)
|
||||
enddo
|
||||
call map_append($ao_integrals_map, copy_key, copy_val, copy_n )
|
||||
endif
|
||||
enddo
|
||||
sync all
|
||||
enddo
|
||||
deallocate( buffer_key, buffer_val, copy_val, copy_key)
|
||||
|
||||
end
|
||||
IRP_ENDIF
|
||||
|
||||
|
||||
integer function load_$ao_integrals(filename)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Read from disk the $ao integrals
|
||||
END_DOC
|
||||
character*(*), intent(in) :: filename
|
||||
integer*8 :: i
|
||||
integer(cache_key_kind), pointer :: key(:)
|
||||
real(integral_kind), pointer :: val(:)
|
||||
integer :: iknd, kknd
|
||||
integer*8 :: n, j
|
||||
load_$ao_integrals = 1
|
||||
open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
|
||||
read(66,err=98,end=98) iknd, kknd
|
||||
if (iknd /= integral_kind) then
|
||||
print *, 'Wrong integrals kind in file :', iknd
|
||||
stop 1
|
||||
endif
|
||||
if (kknd /= key_kind) then
|
||||
print *, 'Wrong key kind in file :', kknd
|
||||
stop 1
|
||||
endif
|
||||
read(66,err=98,end=98) $ao_integrals_map%sorted, $ao_integrals_map%map_size,&
|
||||
$ao_integrals_map%n_elements
|
||||
do i=0_8, $ao_integrals_map%map_size
|
||||
read(66,err=99,end=99) $ao_integrals_map%map(i)%sorted, &
|
||||
$ao_integrals_map%map(i)%map_size, $ao_integrals_map%map(i)%n_elements
|
||||
call cache_map_reallocate($ao_integrals_map%map(i),$ao_integrals_map%map(i)%map_size)
|
||||
enddo
|
||||
do i=0_8, $ao_integrals_map%map_size
|
||||
key => $ao_integrals_map%map(i)%key
|
||||
val => $ao_integrals_map%map(i)%value
|
||||
n = $ao_integrals_map%map(i)%n_elements
|
||||
read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
|
||||
enddo
|
||||
call map_sort($ao_integrals_map)
|
||||
load_$ao_integrals = 0
|
||||
return
|
||||
99 continue
|
||||
call map_deinit($ao_integrals_map)
|
||||
98 continue
|
||||
stop 'Problem reading $ao_integrals_map file in work/'
|
||||
|
||||
end
|
||||
|
||||
SUBST [ ao_integrals_map, ao_integrals, ao_num ]
|
||||
ao_integrals_erf_map ; ao_integrals_erf ; ao_num ;;
|
||||
mo_integrals_erf_map ; mo_integrals_erf ; mo_tot_num;;
|
||||
END_TEMPLATE
|
||||
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ type(map_type), mo_integrals_erf_map ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! MO integrals
|
||||
END_DOC
|
||||
integer(key_kind) :: key_max
|
||||
integer(map_size_kind) :: sze
|
||||
call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max)
|
||||
sze = key_max
|
||||
call map_init(mo_integrals_erf_map,sze)
|
||||
print*, 'MO ERF map initialized'
|
||||
END_PROVIDER
|
||||
|
||||
subroutine insert_into_ao_integrals_erf_map(n_integrals,buffer_i, buffer_values)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Create new entry into AO map
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: n_integrals
|
||||
integer(key_kind), intent(inout) :: buffer_i(n_integrals)
|
||||
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
|
||||
|
||||
call map_append(ao_integrals_erf_map, buffer_i, buffer_values, n_integrals)
|
||||
end
|
||||
|
||||
subroutine insert_into_mo_integrals_erf_map(n_integrals, &
|
||||
buffer_i, buffer_values, thr)
|
||||
use map_module
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Create new entry into MO map, or accumulate in an existing entry
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: n_integrals
|
||||
integer(key_kind), intent(inout) :: buffer_i(n_integrals)
|
||||
real(integral_kind), intent(inout) :: buffer_values(n_integrals)
|
||||
real(integral_kind), intent(in) :: thr
|
||||
call map_update(mo_integrals_erf_map, buffer_i, buffer_values, n_integrals, thr)
|
||||
end
|
||||
|
||||
BEGIN_PROVIDER [ integer, mo_integrals_erf_cache_min ]
|
||||
&BEGIN_PROVIDER [ integer, mo_integrals_erf_cache_max ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Min and max values of the MOs for which the integrals are in the cache
|
||||
END_DOC
|
||||
mo_integrals_erf_cache_min = max(1,elec_alpha_num - 31)
|
||||
mo_integrals_erf_cache_max = min(mo_tot_num,mo_integrals_erf_cache_min+63)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, mo_integrals_erf_cache, (0:64*64*64*64) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Cache of MO integrals for fast access
|
||||
END_DOC
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
integer :: i,j,k,l
|
||||
integer :: ii
|
||||
integer(key_kind) :: idx
|
||||
real(integral_kind) :: integral
|
||||
FREE ao_integrals_erf_cache
|
||||
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
|
||||
do l=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
|
||||
do k=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
|
||||
do j=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
|
||||
do i=mo_integrals_erf_cache_min,mo_integrals_erf_cache_max
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,idx)
|
||||
!DIR$ FORCEINLINE
|
||||
call map_get(mo_integrals_erf_map,idx,integral)
|
||||
ii = l-mo_integrals_erf_cache_min
|
||||
ii = ior( ishft(ii,6), k-mo_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), j-mo_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), i-mo_integrals_erf_cache_min)
|
||||
mo_integrals_erf_cache(ii) = integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
double precision function get_mo_bielec_integral_erf(i,j,k,l,map)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns one integral <ij|kl> in the MO basis
|
||||
END_DOC
|
||||
integer, intent(in) :: i,j,k,l
|
||||
integer(key_kind) :: idx
|
||||
integer :: ii
|
||||
type(map_type), intent(inout) :: map
|
||||
real(integral_kind) :: tmp
|
||||
PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_cache
|
||||
ii = l-mo_integrals_erf_cache_min
|
||||
ii = ior(ii, k-mo_integrals_erf_cache_min)
|
||||
ii = ior(ii, j-mo_integrals_erf_cache_min)
|
||||
ii = ior(ii, i-mo_integrals_erf_cache_min)
|
||||
if (iand(ii, -64) /= 0) then
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,idx)
|
||||
!DIR$ FORCEINLINE
|
||||
call map_get(map,idx,tmp)
|
||||
get_mo_bielec_integral_erf = dble(tmp)
|
||||
else
|
||||
ii = l-mo_integrals_erf_cache_min
|
||||
ii = ior( ishft(ii,6), k-mo_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), j-mo_integrals_erf_cache_min)
|
||||
ii = ior( ishft(ii,6), i-mo_integrals_erf_cache_min)
|
||||
get_mo_bielec_integral_erf = mo_integrals_erf_cache(ii)
|
||||
endif
|
||||
end
|
||||
|
||||
|
||||
double precision function mo_bielec_integral_erf(i,j,k,l)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns one integral <ij|kl> in the MO basis
|
||||
END_DOC
|
||||
integer, intent(in) :: i,j,k,l
|
||||
double precision :: get_mo_bielec_integral_erf
|
||||
PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_cache
|
||||
!DIR$ FORCEINLINE
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
mo_bielec_integral_erf = get_mo_bielec_integral_erf(i,j,k,l,mo_integrals_erf_map)
|
||||
return
|
||||
end
|
||||
|
||||
subroutine get_mo_bielec_integrals_erf(j,k,l,sze,out_val,map)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns multiple integrals <ij|kl> in the MO basis, all
|
||||
! i for j,k,l fixed.
|
||||
END_DOC
|
||||
integer, intent(in) :: j,k,l, sze
|
||||
double precision, intent(out) :: out_val(sze)
|
||||
type(map_type), intent(inout) :: map
|
||||
integer :: i
|
||||
integer(key_kind) :: hash(sze)
|
||||
real(integral_kind) :: tmp_val(sze)
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
|
||||
do i=1,sze
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,hash(i))
|
||||
enddo
|
||||
|
||||
if (key_kind == 8) then
|
||||
call map_get_many(map, hash, out_val, sze)
|
||||
else
|
||||
call map_get_many(map, hash, tmp_val, sze)
|
||||
! Conversion to double precision
|
||||
do i=1,sze
|
||||
out_val(i) = dble(tmp_val(i))
|
||||
enddo
|
||||
endif
|
||||
end
|
||||
|
||||
subroutine get_mo_bielec_integrals_erf_ij(k,l,sze,out_array,map)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns multiple integrals <ij|kl> in the MO basis, all
|
||||
! i(1)j(2) 1/r12 k(1)l(2)
|
||||
! i, j for k,l fixed.
|
||||
END_DOC
|
||||
integer, intent(in) :: k,l, sze
|
||||
double precision, intent(out) :: out_array(sze,sze)
|
||||
type(map_type), intent(inout) :: map
|
||||
integer :: i,j,kk,ll,m
|
||||
integer(key_kind),allocatable :: hash(:)
|
||||
integer ,allocatable :: pairs(:,:), iorder(:)
|
||||
real(integral_kind), allocatable :: tmp_val(:)
|
||||
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), &
|
||||
tmp_val(sze*sze))
|
||||
|
||||
kk=0
|
||||
out_array = 0.d0
|
||||
do j=1,sze
|
||||
do i=1,sze
|
||||
kk += 1
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,hash(kk))
|
||||
pairs(1,kk) = i
|
||||
pairs(2,kk) = j
|
||||
iorder(kk) = kk
|
||||
enddo
|
||||
enddo
|
||||
|
||||
logical :: integral_is_in_map
|
||||
if (key_kind == 8) then
|
||||
call i8radix_sort(hash,iorder,kk,-1)
|
||||
else if (key_kind == 4) then
|
||||
call iradix_sort(hash,iorder,kk,-1)
|
||||
else if (key_kind == 2) then
|
||||
call i2radix_sort(hash,iorder,kk,-1)
|
||||
endif
|
||||
|
||||
call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk)
|
||||
|
||||
do ll=1,kk
|
||||
m = iorder(ll)
|
||||
i=pairs(1,m)
|
||||
j=pairs(2,m)
|
||||
out_array(i,j) = tmp_val(ll)
|
||||
enddo
|
||||
|
||||
deallocate(pairs,hash,iorder,tmp_val)
|
||||
end
|
||||
|
||||
|
||||
subroutine get_mo_bielec_integrals_erf_i1j1(k,l,sze,out_array,map)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns multiple integrals <ik|jl> in the MO basis, all
|
||||
! i(1)j(1) erf(mu_erf * r12) /r12 k(2)l(2)
|
||||
! i, j for k,l fixed.
|
||||
END_DOC
|
||||
integer, intent(in) :: k,l, sze
|
||||
double precision, intent(out) :: out_array(sze,sze)
|
||||
type(map_type), intent(inout) :: map
|
||||
integer :: i,j,kk,ll,m
|
||||
integer(key_kind),allocatable :: hash(:)
|
||||
integer ,allocatable :: pairs(:,:), iorder(:)
|
||||
real(integral_kind), allocatable :: tmp_val(:)
|
||||
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), &
|
||||
tmp_val(sze*sze))
|
||||
|
||||
kk=0
|
||||
out_array = 0.d0
|
||||
do j=1,sze
|
||||
do i=1,sze
|
||||
kk += 1
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,k,j,l,hash(kk))
|
||||
pairs(1,kk) = i
|
||||
pairs(2,kk) = j
|
||||
iorder(kk) = kk
|
||||
enddo
|
||||
enddo
|
||||
|
||||
logical :: integral_is_in_map
|
||||
if (key_kind == 8) then
|
||||
call i8radix_sort(hash,iorder,kk,-1)
|
||||
else if (key_kind == 4) then
|
||||
call iradix_sort(hash,iorder,kk,-1)
|
||||
else if (key_kind == 2) then
|
||||
call i2radix_sort(hash,iorder,kk,-1)
|
||||
endif
|
||||
|
||||
call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk)
|
||||
|
||||
do ll=1,kk
|
||||
m = iorder(ll)
|
||||
i=pairs(1,m)
|
||||
j=pairs(2,m)
|
||||
out_array(i,j) = tmp_val(ll)
|
||||
enddo
|
||||
|
||||
deallocate(pairs,hash,iorder,tmp_val)
|
||||
end
|
||||
|
||||
|
||||
subroutine get_mo_bielec_integrals_erf_coulomb_ii(k,l,sze,out_val,map)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns multiple integrals <ki|li>
|
||||
! k(1)i(2) 1/r12 l(1)i(2) :: out_val(i1)
|
||||
! for k,l fixed.
|
||||
END_DOC
|
||||
integer, intent(in) :: k,l, sze
|
||||
double precision, intent(out) :: out_val(sze)
|
||||
type(map_type), intent(inout) :: map
|
||||
integer :: i
|
||||
integer(key_kind) :: hash(sze)
|
||||
real(integral_kind) :: tmp_val(sze)
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
|
||||
integer :: kk
|
||||
do i=1,sze
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(k,i,l,i,hash(i))
|
||||
enddo
|
||||
|
||||
if (key_kind == 8) then
|
||||
call map_get_many(map, hash, out_val, sze)
|
||||
else
|
||||
call map_get_many(map, hash, tmp_val, sze)
|
||||
! Conversion to double precision
|
||||
do i=1,sze
|
||||
out_val(i) = dble(tmp_val(i))
|
||||
enddo
|
||||
endif
|
||||
end
|
||||
|
||||
subroutine get_mo_bielec_integrals_erf_exch_ii(k,l,sze,out_val,map)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns multiple integrals <ki|il>
|
||||
! k(1)i(2) 1/r12 i(1)l(2) :: out_val(i1)
|
||||
! for k,l fixed.
|
||||
END_DOC
|
||||
integer, intent(in) :: k,l, sze
|
||||
double precision, intent(out) :: out_val(sze)
|
||||
type(map_type), intent(inout) :: map
|
||||
integer :: i
|
||||
integer(key_kind) :: hash(sze)
|
||||
real(integral_kind) :: tmp_val(sze)
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
|
||||
integer :: kk
|
||||
do i=1,sze
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(k,i,i,l,hash(i))
|
||||
enddo
|
||||
|
||||
if (key_kind == 8) then
|
||||
call map_get_many(map, hash, out_val, sze)
|
||||
else
|
||||
call map_get_many(map, hash, tmp_val, sze)
|
||||
! Conversion to double precision
|
||||
do i=1,sze
|
||||
out_val(i) = dble(tmp_val(i))
|
||||
enddo
|
||||
endif
|
||||
end
|
||||
|
||||
|
||||
integer*8 function get_mo_erf_map_size()
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Return the number of elements in the MO map
|
||||
END_DOC
|
||||
get_mo_erf_map_size = mo_integrals_erf_map % n_elements
|
||||
end
|
549
src/integrals_bielec_erf/mo_bi_integrals_erf.irp.f
Normal file
549
src/integrals_bielec_erf/mo_bi_integrals_erf.irp.f
Normal file
@ -0,0 +1,549 @@
|
||||
subroutine mo_bielec_integrals_erf_index(i,j,k,l,i1)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes an unique index for i,j,k,l integrals
|
||||
END_DOC
|
||||
integer, intent(in) :: i,j,k,l
|
||||
integer(key_kind), intent(out) :: i1
|
||||
integer(key_kind) :: p,q,r,s,i2
|
||||
p = min(i,k)
|
||||
r = max(i,k)
|
||||
p = p+ishft(r*r-r,-1)
|
||||
q = min(j,l)
|
||||
s = max(j,l)
|
||||
q = q+ishft(s*s-s,-1)
|
||||
i1 = min(p,q)
|
||||
i2 = max(p,q)
|
||||
i1 = i1+ishft(i2*i2-i2,-1)
|
||||
end
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ logical, mo_bielec_integrals_erf_in_map ]
|
||||
use map_module
|
||||
implicit none
|
||||
integer(bit_kind) :: mask_ijkl(N_int,4)
|
||||
integer(bit_kind) :: mask_ijk(N_int,3)
|
||||
|
||||
BEGIN_DOC
|
||||
! If True, the map of MO bielectronic integrals is provided
|
||||
END_DOC
|
||||
|
||||
real :: map_mb
|
||||
|
||||
mo_bielec_integrals_erf_in_map = .True.
|
||||
if (read_mo_integrals_erf) then
|
||||
print*,'Reading the MO integrals_erf'
|
||||
call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
|
||||
print*, 'MO integrals_erf provided'
|
||||
return
|
||||
else
|
||||
PROVIDE ao_bielec_integrals_erf_in_map
|
||||
endif
|
||||
|
||||
! call four_index_transform_block(ao_integrals_erf_map,mo_integrals_erf_map, &
|
||||
! mo_coef, size(mo_coef,1), &
|
||||
! 1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num, &
|
||||
! 1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num)
|
||||
call add_integrals_to_map_erf(full_ijkl_bitmask_4)
|
||||
integer*8 :: get_mo_erf_map_size, mo_erf_map_size
|
||||
mo_erf_map_size = get_mo_erf_map_size()
|
||||
|
||||
! print*,'Molecular integrals ERF provided:'
|
||||
! print*,' Size of MO ERF map ', map_mb(mo_integrals_erf_map) ,'MB'
|
||||
! print*,' Number of MO ERF integrals: ', mo_erf_map_size
|
||||
if (write_mo_integrals_erf) then
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
|
||||
call ezfio_set_integrals_erf_disk_access_mo_integrals_erf("Read")
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_from_ao, (mo_tot_num,mo_tot_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_exchange_from_ao, (mo_tot_num,mo_tot_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_anti_from_ao, (mo_tot_num,mo_tot_num) ]
|
||||
BEGIN_DOC
|
||||
! mo_bielec_integral_jj_from_ao(i,j) = J_ij
|
||||
! mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij
|
||||
! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j,p,q,r,s
|
||||
double precision :: c
|
||||
real(integral_kind) :: integral
|
||||
integer :: n, pp
|
||||
real(integral_kind), allocatable :: int_value(:)
|
||||
integer, allocatable :: int_idx(:)
|
||||
|
||||
double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:)
|
||||
|
||||
if (.not.do_direct_integrals) then
|
||||
PROVIDE ao_bielec_integrals_erf_in_map mo_coef
|
||||
endif
|
||||
|
||||
mo_bielec_integral_erf_jj_from_ao = 0.d0
|
||||
mo_bielec_integral_erf_jj_exchange_from_ao = 0.d0
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr
|
||||
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, &
|
||||
!$OMP iqrs, iqsr,iqri,iqis) &
|
||||
!$OMP SHARED(mo_tot_num,mo_coef_transp,ao_num,&
|
||||
!$OMP ao_integrals_threshold,do_direct_integrals) &
|
||||
!$OMP REDUCTION(+:mo_bielec_integral_erf_jj_from_ao,mo_bielec_integral_erf_jj_exchange_from_ao)
|
||||
|
||||
allocate( int_value(ao_num), int_idx(ao_num), &
|
||||
iqrs(mo_tot_num,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),&
|
||||
iqsr(mo_tot_num,ao_num) )
|
||||
|
||||
!$OMP DO SCHEDULE (guided)
|
||||
do s=1,ao_num
|
||||
do q=1,ao_num
|
||||
|
||||
do j=1,ao_num
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,mo_tot_num
|
||||
iqrs(i,j) = 0.d0
|
||||
iqsr(i,j) = 0.d0
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if (do_direct_integrals) then
|
||||
double precision :: ao_bielec_integral_erf
|
||||
do r=1,ao_num
|
||||
call compute_ao_bielec_integrals_erf(q,r,s,ao_num,int_value)
|
||||
do p=1,ao_num
|
||||
integral = int_value(p)
|
||||
if (abs(integral) > ao_integrals_threshold) then
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,mo_tot_num
|
||||
iqrs(i,r) += mo_coef_transp(i,p) * integral
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
call compute_ao_bielec_integrals_erf(q,s,r,ao_num,int_value)
|
||||
do p=1,ao_num
|
||||
integral = int_value(p)
|
||||
if (abs(integral) > ao_integrals_threshold) then
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,mo_tot_num
|
||||
iqsr(i,r) += mo_coef_transp(i,p) * integral
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
do r=1,ao_num
|
||||
call get_ao_bielec_integrals_erf_non_zero(q,r,s,ao_num,int_value,int_idx,n)
|
||||
do pp=1,n
|
||||
p = int_idx(pp)
|
||||
integral = int_value(pp)
|
||||
if (abs(integral) > ao_integrals_threshold) then
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,mo_tot_num
|
||||
iqrs(i,r) += mo_coef_transp(i,p) * integral
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
call get_ao_bielec_integrals_erf_non_zero(q,s,r,ao_num,int_value,int_idx,n)
|
||||
do pp=1,n
|
||||
p = int_idx(pp)
|
||||
integral = int_value(pp)
|
||||
if (abs(integral) > ao_integrals_threshold) then
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,mo_tot_num
|
||||
iqsr(i,r) += mo_coef_transp(i,p) * integral
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
iqis = 0.d0
|
||||
iqri = 0.d0
|
||||
do r=1,ao_num
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,mo_tot_num
|
||||
iqis(i) += mo_coef_transp(i,r) * iqrs(i,r)
|
||||
iqri(i) += mo_coef_transp(i,r) * iqsr(i,r)
|
||||
enddo
|
||||
enddo
|
||||
do i=1,mo_tot_num
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do j=1,mo_tot_num
|
||||
c = mo_coef_transp(j,q)*mo_coef_transp(j,s)
|
||||
mo_bielec_integral_erf_jj_from_ao(j,i) += c * iqis(i)
|
||||
mo_bielec_integral_erf_jj_exchange_from_ao(j,i) += c * iqri(i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
deallocate(iqrs,iqsr,int_value,int_idx)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
mo_bielec_integral_erf_jj_anti_from_ao = mo_bielec_integral_erf_jj_from_ao - mo_bielec_integral_erf_jj_exchange_from_ao
|
||||
|
||||
|
||||
! end
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj, (mo_tot_num,mo_tot_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_exchange, (mo_tot_num,mo_tot_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_erf_jj_anti, (mo_tot_num,mo_tot_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! mo_bielec_integral_jj(i,j) = J_ij
|
||||
! mo_bielec_integral_jj_exchange(i,j) = K_ij
|
||||
! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij
|
||||
END_DOC
|
||||
|
||||
integer :: i,j
|
||||
double precision :: get_mo_bielec_integral_erf
|
||||
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
mo_bielec_integral_erf_jj = 0.d0
|
||||
mo_bielec_integral_erf_jj_exchange = 0.d0
|
||||
|
||||
do j=1,mo_tot_num
|
||||
do i=1,mo_tot_num
|
||||
mo_bielec_integral_erf_jj(i,j) = get_mo_bielec_integral_erf(i,j,i,j,mo_integrals_erf_map)
|
||||
mo_bielec_integral_erf_jj_exchange(i,j) = get_mo_bielec_integral_erf(i,j,j,i,mo_integrals_erf_map)
|
||||
mo_bielec_integral_erf_jj_anti(i,j) = mo_bielec_integral_erf_jj(i,j) - mo_bielec_integral_erf_jj_exchange(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine clear_mo_erf_map
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Frees the memory of the MO map
|
||||
END_DOC
|
||||
call map_deinit(mo_integrals_erf_map)
|
||||
FREE mo_integrals_erf_map mo_bielec_integral_erf_jj mo_bielec_integral_erf_jj_anti
|
||||
FREE mo_bielec_integral_Erf_jj_exchange mo_bielec_integrals_erf_in_map
|
||||
|
||||
|
||||
end
|
||||
|
||||
subroutine provide_all_mo_integrals_erf
|
||||
implicit none
|
||||
provide mo_integrals_erf_map mo_bielec_integral_erf_jj mo_bielec_integral_erf_jj_anti
|
||||
provide mo_bielec_integral_erf_jj_exchange mo_bielec_integrals_erf_in_map
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine add_integrals_to_map_erf(mask_ijkl)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Adds integrals to tha MO map according to some bitmask
|
||||
END_DOC
|
||||
|
||||
integer(bit_kind), intent(in) :: mask_ijkl(N_int,4)
|
||||
|
||||
integer :: i,j,k,l
|
||||
integer :: i0,j0,k0,l0
|
||||
double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0
|
||||
|
||||
integer, allocatable :: list_ijkl(:,:)
|
||||
integer :: n_i, n_j, n_k, n_l
|
||||
integer, allocatable :: bielec_tmp_0_idx(:)
|
||||
real(integral_kind), allocatable :: bielec_tmp_0(:,:)
|
||||
double precision, allocatable :: bielec_tmp_1(:)
|
||||
double precision, allocatable :: bielec_tmp_2(:,:)
|
||||
double precision, allocatable :: bielec_tmp_3(:,:,:)
|
||||
!DIR$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3
|
||||
|
||||
integer :: n_integrals
|
||||
integer :: size_buffer
|
||||
integer(key_kind),allocatable :: buffer_i(:)
|
||||
real(integral_kind),allocatable :: buffer_value(:)
|
||||
double precision :: map_mb
|
||||
|
||||
integer :: i1,j1,k1,l1, ii1, kmax, thread_num
|
||||
integer :: i2,i3,i4
|
||||
double precision,parameter :: thr_coef = 1.d-10
|
||||
|
||||
PROVIDE ao_bielec_integrals_in_map mo_coef
|
||||
|
||||
!Get list of MOs for i,j,k and l
|
||||
!-------------------------------
|
||||
|
||||
allocate(list_ijkl(mo_tot_num,4))
|
||||
call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int )
|
||||
call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int )
|
||||
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
|
||||
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
|
||||
character*(2048) :: output(1)
|
||||
print*, 'i'
|
||||
call bitstring_to_str( output(1), mask_ijkl(1,1), N_int )
|
||||
print *, trim(output(1))
|
||||
j = 0
|
||||
do i = 1, N_int
|
||||
j += popcnt(mask_ijkl(i,1))
|
||||
enddo
|
||||
if(j==0)then
|
||||
return
|
||||
endif
|
||||
|
||||
print*, 'j'
|
||||
call bitstring_to_str( output(1), mask_ijkl(1,2), N_int )
|
||||
print *, trim(output(1))
|
||||
j = 0
|
||||
do i = 1, N_int
|
||||
j += popcnt(mask_ijkl(i,2))
|
||||
enddo
|
||||
if(j==0)then
|
||||
return
|
||||
endif
|
||||
|
||||
print*, 'k'
|
||||
call bitstring_to_str( output(1), mask_ijkl(1,3), N_int )
|
||||
print *, trim(output(1))
|
||||
j = 0
|
||||
do i = 1, N_int
|
||||
j += popcnt(mask_ijkl(i,3))
|
||||
enddo
|
||||
if(j==0)then
|
||||
return
|
||||
endif
|
||||
|
||||
print*, 'l'
|
||||
call bitstring_to_str( output(1), mask_ijkl(1,4), N_int )
|
||||
print *, trim(output(1))
|
||||
j = 0
|
||||
do i = 1, N_int
|
||||
j += popcnt(mask_ijkl(i,4))
|
||||
enddo
|
||||
if(j==0)then
|
||||
return
|
||||
endif
|
||||
|
||||
size_buffer = min(ao_num*ao_num*ao_num,16000000)
|
||||
print*, 'Providing the ERF molecular integrals '
|
||||
print*, 'Buffers : ', 8.*(mo_tot_num*(n_j)*(n_k+1) + mo_tot_num+&
|
||||
ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core'
|
||||
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
double precision :: accu_bis
|
||||
accu_bis = 0.d0
|
||||
|
||||
!$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, &
|
||||
!$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,&
|
||||
!$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, &
|
||||
!$OMP wall_0,thread_num,accu_bis) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l, &
|
||||
!$OMP mo_coef_transp, &
|
||||
!$OMP mo_coef_transp_is_built, list_ijkl, &
|
||||
!$OMP mo_coef_is_built, wall_1, &
|
||||
!$OMP mo_coef,mo_integrals_threshold,mo_integrals_erf_map)
|
||||
n_integrals = 0
|
||||
wall_0 = wall_1
|
||||
allocate(bielec_tmp_3(mo_tot_num, n_j, n_k), &
|
||||
bielec_tmp_1(mo_tot_num), &
|
||||
bielec_tmp_0(ao_num,ao_num), &
|
||||
bielec_tmp_0_idx(ao_num), &
|
||||
bielec_tmp_2(mo_tot_num, n_j), &
|
||||
buffer_i(size_buffer), &
|
||||
buffer_value(size_buffer) )
|
||||
|
||||
thread_num = 0
|
||||
!$ thread_num = omp_get_thread_num()
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
do l1 = 1,ao_num
|
||||
bielec_tmp_3 = 0.d0
|
||||
do k1 = 1,ao_num
|
||||
bielec_tmp_2 = 0.d0
|
||||
do j1 = 1,ao_num
|
||||
call get_ao_bielec_integrals_erf(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) ! all integrals for a given l1, k1
|
||||
! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1))
|
||||
enddo
|
||||
do j1 = 1,ao_num
|
||||
kmax = 0
|
||||
do i1 = 1,ao_num
|
||||
c = bielec_tmp_0(i1,j1)
|
||||
if (c == 0.d0) then
|
||||
cycle
|
||||
endif
|
||||
kmax += 1
|
||||
bielec_tmp_0(kmax,j1) = c
|
||||
bielec_tmp_0_idx(kmax) = i1
|
||||
enddo
|
||||
|
||||
if (kmax==0) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
bielec_tmp_1 = 0.d0
|
||||
ii1=1
|
||||
! sum_m c_m^i (m)
|
||||
do ii1 = 1,kmax-4,4
|
||||
i1 = bielec_tmp_0_idx(ii1)
|
||||
i2 = bielec_tmp_0_idx(ii1+1)
|
||||
i3 = bielec_tmp_0_idx(ii1+2)
|
||||
i4 = bielec_tmp_0_idx(ii1+3)
|
||||
do i = list_ijkl(1,1), list_ijkl(n_i,1)
|
||||
bielec_tmp_1(i) = bielec_tmp_1(i) + &
|
||||
mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + &
|
||||
mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + &
|
||||
mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + &
|
||||
mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1)
|
||||
enddo ! i
|
||||
enddo ! ii1
|
||||
|
||||
i2 = ii1
|
||||
do ii1 = i2,kmax
|
||||
i1 = bielec_tmp_0_idx(ii1)
|
||||
do i = list_ijkl(1,1), list_ijkl(n_i,1)
|
||||
bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1)
|
||||
enddo ! i
|
||||
enddo ! ii1
|
||||
c = 0.d0
|
||||
|
||||
do i = list_ijkl(1,1), list_ijkl(n_i,1)
|
||||
c = max(c,abs(bielec_tmp_1(i)))
|
||||
if (c>mo_integrals_threshold) exit
|
||||
enddo
|
||||
if ( c < mo_integrals_threshold ) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
do j0 = 1, n_j
|
||||
j = list_ijkl(j0,2)
|
||||
c = mo_coef_transp(j,j1)
|
||||
if (abs(c) < thr_coef) then
|
||||
cycle
|
||||
endif
|
||||
do i = list_ijkl(1,1), list_ijkl(n_i,1)
|
||||
bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i)
|
||||
enddo ! i
|
||||
enddo ! j
|
||||
enddo !j1
|
||||
if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
|
||||
do k0 = 1, n_k
|
||||
k = list_ijkl(k0,3)
|
||||
c = mo_coef_transp(k,k1)
|
||||
if (abs(c) < thr_coef) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
do j0 = 1, n_j
|
||||
j = list_ijkl(j0,2)
|
||||
do i = list_ijkl(1,1), k
|
||||
bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0)
|
||||
enddo!i
|
||||
enddo !j
|
||||
|
||||
enddo !k
|
||||
enddo !k1
|
||||
|
||||
|
||||
|
||||
do l0 = 1,n_l
|
||||
l = list_ijkl(l0,4)
|
||||
c = mo_coef_transp(l,l1)
|
||||
if (abs(c) < thr_coef) then
|
||||
cycle
|
||||
endif
|
||||
j1 = ishft((l*l-l),-1)
|
||||
do j0 = 1, n_j
|
||||
j = list_ijkl(j0,2)
|
||||
if (j > l) then
|
||||
exit
|
||||
endif
|
||||
j1 += 1
|
||||
do k0 = 1, n_k
|
||||
k = list_ijkl(k0,3)
|
||||
i1 = ishft((k*k-k),-1)
|
||||
if (i1<=j1) then
|
||||
continue
|
||||
else
|
||||
exit
|
||||
endif
|
||||
bielec_tmp_1 = 0.d0
|
||||
do i0 = 1, n_i
|
||||
i = list_ijkl(i0,1)
|
||||
if (i>k) then
|
||||
exit
|
||||
endif
|
||||
bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0)
|
||||
! i1+=1
|
||||
enddo
|
||||
|
||||
do i0 = 1, n_i
|
||||
i = list_ijkl(i0,1)
|
||||
if(i> min(k,j1-i1+list_ijkl(1,1)-1))then
|
||||
exit
|
||||
endif
|
||||
if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then
|
||||
cycle
|
||||
endif
|
||||
n_integrals += 1
|
||||
buffer_value(n_integrals) = bielec_tmp_1(i)
|
||||
!DIR$ FORCEINLINE
|
||||
call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
|
||||
if (n_integrals == size_buffer) then
|
||||
call insert_into_mo_integrals_erf_map(n_integrals,buffer_i,buffer_value,&
|
||||
real(mo_integrals_threshold,integral_kind))
|
||||
n_integrals = 0
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall_2)
|
||||
if (thread_num == 0) then
|
||||
if (wall_2 - wall_0 > 1.d0) then
|
||||
wall_0 = wall_2
|
||||
print*, 100.*float(l1)/float(ao_num), '% in ', &
|
||||
wall_2-wall_1, 's', map_mb(mo_integrals_erf_map) ,'MB'
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3)
|
||||
|
||||
integer :: index_needed
|
||||
|
||||
call insert_into_mo_integrals_erf_map(n_integrals,buffer_i,buffer_value,&
|
||||
real(mo_integrals_threshold,integral_kind))
|
||||
deallocate(buffer_i, buffer_value)
|
||||
!$OMP END PARALLEL
|
||||
call map_merge(mo_integrals_erf_map)
|
||||
|
||||
call wall_time(wall_2)
|
||||
call cpu_time(cpu_2)
|
||||
integer*8 :: get_mo_erf_map_size, mo_map_size
|
||||
mo_map_size = get_mo_erf_map_size()
|
||||
|
||||
deallocate(list_ijkl)
|
||||
|
||||
|
||||
print*,'Molecular ERF integrals provided:'
|
||||
print*,' Size of MO ERF map ', map_mb(mo_integrals_erf_map) ,'MB'
|
||||
print*,' Number of MO ERF integrals: ', mo_map_size
|
||||
print*,' cpu time :',cpu_2 - cpu_1, 's'
|
||||
print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
|
||||
|
||||
end
|
||||
|
||||
|
125
src/integrals_bielec_erf/providers_ao_erf.irp.f
Normal file
125
src/integrals_bielec_erf/providers_ao_erf.irp.f
Normal file
@ -0,0 +1,125 @@
|
||||
|
||||
BEGIN_PROVIDER [ logical, ao_bielec_integrals_erf_in_map ]
|
||||
implicit none
|
||||
use f77_zmq
|
||||
use map_module
|
||||
BEGIN_DOC
|
||||
! Map of Atomic integrals
|
||||
! i(r1) j(r2) 1/r12 k(r1) l(r2)
|
||||
END_DOC
|
||||
|
||||
integer :: i,j,k,l
|
||||
double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2
|
||||
double precision :: integral, wall_0
|
||||
include 'Utils/constants.include.F'
|
||||
|
||||
! For integrals file
|
||||
integer(key_kind),allocatable :: buffer_i(:)
|
||||
integer,parameter :: size_buffer = 1024*64
|
||||
real(integral_kind),allocatable :: buffer_value(:)
|
||||
|
||||
integer :: n_integrals, rc
|
||||
integer :: kk, m, j1, i1, lmax
|
||||
character*(64) :: fmt
|
||||
|
||||
integral = ao_bielec_integral_erf(1,1,1,1)
|
||||
|
||||
double precision :: map_mb
|
||||
PROVIDE read_ao_integrals_erf disk_access_ao_integrals_erf
|
||||
if (read_ao_integrals_erf) then
|
||||
print*,'Reading the AO ERF integrals'
|
||||
call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
|
||||
print*, 'AO ERF integrals provided'
|
||||
ao_bielec_integrals_erf_in_map = .True.
|
||||
return
|
||||
endif
|
||||
|
||||
print*, 'Providing the AO ERF integrals'
|
||||
call wall_time(wall_0)
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
|
||||
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals_erf')
|
||||
|
||||
character(len=:), allocatable :: task
|
||||
allocate(character(len=ao_num*12) :: task)
|
||||
write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))'
|
||||
do l=1,ao_num
|
||||
write(task,fmt) (i,l, i=1,l)
|
||||
integer, external :: add_task_to_taskserver
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then
|
||||
stop 'Unable to add task to server'
|
||||
endif
|
||||
enddo
|
||||
deallocate(task)
|
||||
|
||||
integer, external :: zmq_set_running
|
||||
if (zmq_set_running(zmq_to_qp_run_socket) == -1) then
|
||||
print *, irp_here, ': Failed in zmq_set_running'
|
||||
endif
|
||||
|
||||
PROVIDE nproc
|
||||
!$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1)
|
||||
i = omp_get_thread_num()
|
||||
if (i==0) then
|
||||
call ao_bielec_integrals_erf_in_map_collector(zmq_socket_pull)
|
||||
else
|
||||
call ao_bielec_integrals_erf_in_map_slave_inproc(i)
|
||||
endif
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals_erf')
|
||||
|
||||
|
||||
print*, 'Sorting the map'
|
||||
call map_sort(ao_integrals_erf_map)
|
||||
call cpu_time(cpu_2)
|
||||
call wall_time(wall_2)
|
||||
integer(map_size_kind) :: get_ao_erf_map_size, ao_erf_map_size
|
||||
ao_erf_map_size = get_ao_erf_map_size()
|
||||
|
||||
print*, 'AO ERF integrals provided:'
|
||||
print*, ' Size of AO ERF map : ', map_mb(ao_integrals_erf_map) ,'MB'
|
||||
print*, ' Number of AO ERF integrals :', ao_erf_map_size
|
||||
print*, ' cpu time :',cpu_2 - cpu_1, 's'
|
||||
print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
|
||||
|
||||
ao_bielec_integrals_erf_in_map = .True.
|
||||
|
||||
if (write_ao_integrals_erf) then
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
|
||||
call ezfio_set_integrals_erf_disk_access_ao_integrals_erf("Read")
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_bielec_integral_erf_schwartz,(ao_num,ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Needed to compute Schwartz inequalities
|
||||
END_DOC
|
||||
|
||||
integer :: i,k
|
||||
double precision :: ao_bielec_integral_erf,cpu_1,cpu_2, wall_1, wall_2
|
||||
|
||||
ao_bielec_integral_erf_schwartz(1,1) = ao_bielec_integral_erf(1,1,1,1)
|
||||
!$OMP PARALLEL DO PRIVATE(i,k) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP SHARED (ao_num,ao_bielec_integral_erf_schwartz) &
|
||||
!$OMP SCHEDULE(dynamic)
|
||||
do i=1,ao_num
|
||||
do k=1,i
|
||||
ao_bielec_integral_erf_schwartz(i,k) = dsqrt(ao_bielec_integral_erf(i,k,i,k))
|
||||
ao_bielec_integral_erf_schwartz(k,i) = ao_bielec_integral_erf_schwartz(i,k)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
32
src/integrals_bielec_erf/qp_ao_erf_ints.irp.f
Normal file
32
src/integrals_bielec_erf/qp_ao_erf_ints.irp.f
Normal file
@ -0,0 +1,32 @@
|
||||
program qp_ao_ints
|
||||
use omp_lib
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Increments a running calculation to compute AO integral_erfs
|
||||
END_DOC
|
||||
integer :: i
|
||||
|
||||
call switch_qp_run_to_master
|
||||
|
||||
zmq_context = f77_zmq_ctx_new ()
|
||||
|
||||
! Set the state of the ZMQ
|
||||
zmq_state = 'ao_integral_erfs'
|
||||
|
||||
! Provide everything needed
|
||||
double precision :: integral_erf, ao_bielec_integral_erf
|
||||
integral_erf = ao_bielec_integral_erf(1,1,1,1)
|
||||
|
||||
character*(64) :: state
|
||||
call wait_for_state(zmq_state,state)
|
||||
do while (state /= 'Stopped')
|
||||
!$OMP PARALLEL DEFAULT(PRIVATE) PRIVATE(i)
|
||||
i = omp_get_thread_num()
|
||||
call ao_bielec_integrals_erf_in_map_slave_tcp(i)
|
||||
!$OMP END PARALLEL
|
||||
call wait_for_state(zmq_state,state)
|
||||
enddo
|
||||
|
||||
print *, 'Done'
|
||||
end
|
||||
|
47
src/integrals_bielec_erf/read_write.irp.f
Normal file
47
src/integrals_bielec_erf/read_write.irp.f
Normal file
@ -0,0 +1,47 @@
|
||||
BEGIN_PROVIDER [ logical, read_ao_integrals_erf ]
|
||||
&BEGIN_PROVIDER [ logical, read_mo_integrals_erf ]
|
||||
&BEGIN_PROVIDER [ logical, write_ao_integrals_erf ]
|
||||
&BEGIN_PROVIDER [ logical, write_mo_integrals_erf ]
|
||||
|
||||
BEGIN_DOC
|
||||
! One level of abstraction for disk_access_ao_integrals_erf and disk_access_mo_integrals_erf
|
||||
END_DOC
|
||||
implicit none
|
||||
|
||||
if (disk_access_ao_integrals_erf.EQ.'Read') then
|
||||
read_ao_integrals_erf = .True.
|
||||
write_ao_integrals_erf = .False.
|
||||
|
||||
else if (disk_access_ao_integrals_erf.EQ.'Write') then
|
||||
read_ao_integrals_erf = .False.
|
||||
write_ao_integrals_erf = .True.
|
||||
|
||||
else if (disk_access_ao_integrals_erf.EQ.'None') then
|
||||
read_ao_integrals_erf = .False.
|
||||
write_ao_integrals_erf = .False.
|
||||
|
||||
else
|
||||
print *, 'bielec_integrals_erf/disk_access_ao_integrals_erf has a wrong type'
|
||||
stop 1
|
||||
|
||||
endif
|
||||
|
||||
if (disk_access_mo_integrals_erf.EQ.'Read') then
|
||||
read_mo_integrals_erf = .True.
|
||||
write_mo_integrals_erf = .False.
|
||||
|
||||
else if (disk_access_mo_integrals_erf.EQ.'Write') then
|
||||
read_mo_integrals_erf = .False.
|
||||
write_mo_integrals_erf = .True.
|
||||
|
||||
else if (disk_access_mo_integrals_erf.EQ.'None') then
|
||||
read_mo_integrals_erf = .False.
|
||||
write_mo_integrals_erf = .False.
|
||||
|
||||
else
|
||||
print *, 'bielec_integrals_erf/disk_access_mo_integrals_erf has a wrong type'
|
||||
stop 1
|
||||
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
36
src/integrals_bielec_erf/routines_save_integrals_erf.irp.f
Normal file
36
src/integrals_bielec_erf/routines_save_integrals_erf.irp.f
Normal file
@ -0,0 +1,36 @@
|
||||
subroutine save_erf_bi_elec_integrals_mo
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
|
||||
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
|
||||
end
|
||||
|
||||
subroutine save_erf_bi_elec_integrals_ao
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
PROVIDE ao_bielec_integrals_erf_in_map
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
|
||||
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
|
||||
end
|
||||
|
||||
subroutine save_erf_bi_elec_integrals_mo_into_regular_integrals_mo
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
PROVIDE mo_bielec_integrals_erf_in_map
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_erf_map)
|
||||
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
|
||||
end
|
||||
|
||||
subroutine save_erf_bi_elec_integrals_ao_into_regular_integrals_ao
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
PROVIDE ao_bielec_integrals_erf_in_map
|
||||
call ezfio_set_work_empty(.False.)
|
||||
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_erf_map)
|
||||
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
|
||||
end
|
||||
|
20
src/integrals_bielec_erf/write_integrals_erf.irp.f
Normal file
20
src/integrals_bielec_erf/write_integrals_erf.irp.f
Normal file
@ -0,0 +1,20 @@
|
||||
program write_integrals
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! This program saves the bielec erf integrals into the EZFIO folder
|
||||
END_DOC
|
||||
disk_access_mo_integrals = "None"
|
||||
touch disk_access_mo_integrals
|
||||
disk_access_ao_integrals = "None"
|
||||
touch disk_access_ao_integrals
|
||||
call routine
|
||||
|
||||
end
|
||||
|
||||
subroutine routine
|
||||
implicit none
|
||||
call save_erf_bi_elec_integrals_ao
|
||||
call save_erf_bi_elec_integrals_mo
|
||||
|
||||
end
|
||||
|
@ -0,0 +1,21 @@
|
||||
program write_integrals_erf_to_regular_integrals
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! This program saves the bielec erf integrals into the EZFIO folder but at the regular bielec integrals place.
|
||||
! Threfore, if the user runs a future calculation with a reading of the integrals, the calculation will be performed with the erf bielec integrals instead of the regular bielec integrals
|
||||
END_DOC
|
||||
disk_access_mo_integrals = "None"
|
||||
touch disk_access_mo_integrals
|
||||
disk_access_ao_integrals = "None"
|
||||
touch disk_access_ao_integrals
|
||||
call routine
|
||||
|
||||
end
|
||||
|
||||
subroutine routine
|
||||
implicit none
|
||||
call save_erf_bi_elec_integrals_mo_into_regular_integrals_mo
|
||||
call save_erf_bi_elec_integrals_ao_into_regular_integrals_ao
|
||||
|
||||
end
|
||||
|
55
src/mo_basis/mos_in_r.irp.f
Normal file
55
src/mo_basis/mos_in_r.irp.f
Normal file
@ -0,0 +1,55 @@
|
||||
|
||||
subroutine give_all_mos_at_r(r,mos_array)
|
||||
implicit none
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: mos_array(mo_tot_num)
|
||||
double precision :: aos_array(ao_num)
|
||||
call give_all_aos_at_r(r,aos_array)
|
||||
call dgemv('N',mo_tot_num,ao_num,1.d0,mo_coef_transp,mo_tot_num,aos_array,1,0.d0,mos_array,1)
|
||||
end
|
||||
|
||||
subroutine give_all_mos_and_grad_at_r(r,mos_array,mos_grad_array)
|
||||
implicit none
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: mos_array(mo_tot_num)
|
||||
double precision, intent(out) :: mos_grad_array(mo_tot_num,3)
|
||||
integer :: i,j,k
|
||||
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3)
|
||||
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
|
||||
mos_array=0d0
|
||||
mos_grad_array=0d0
|
||||
do j = 1, mo_tot_num
|
||||
do k=1, ao_num
|
||||
mos_array(j) += mo_coef(k,j)*aos_array(k)
|
||||
mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1)
|
||||
mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
|
||||
mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_lapl_array)
|
||||
implicit none
|
||||
double precision, intent(in) :: r(3)
|
||||
double precision, intent(out) :: mos_array(mo_tot_num)
|
||||
double precision, intent(out) :: mos_grad_array(mo_tot_num,3),mos_lapl_array(mo_tot_num,3)
|
||||
integer :: i,j,k
|
||||
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3)
|
||||
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
|
||||
mos_array=0d0
|
||||
mos_grad_array=0d0
|
||||
mos_lapl_array=0d0
|
||||
do j = 1, mo_tot_num
|
||||
do k=1, ao_num
|
||||
mos_array(j) += mo_coef(k,j)*aos_array(k)
|
||||
mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1)
|
||||
mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
|
||||
mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
|
||||
mos_lapl_array(j,1) += mo_coef(k,j)*aos_lapl_array(k,1)
|
||||
mos_lapl_array(j,2) += mo_coef(k,j)*aos_lapl_array(k,2)
|
||||
mos_lapl_array(j,3) += mo_coef(k,j)*aos_lapl_array(k,3)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
Loading…
Reference in New Issue
Block a user