PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-26 15:12:19 +02:00
Code Releases Activity

added tc_keywords and three_body_ints

Browse Source
This commit is contained in:
eginer 2022-10-05 17:29:26 +02:00
parent 57e592870c
commit ad0203c959
19 changed files with 2501 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
src/some_mu_of_r/EZFIO.cfg Normal file
View File

@ -0,0 +1,59 @@
[constant_mu]
type: logical
doc: If |true|, the mu(r) is constant and set to mu_erf
interface: ezfio,provider,ocaml
default: True
[mu_of_r_tc_ints]
type: character*(32)
doc: type of mu(r) for the TC Hamiltonian : can be [ basis| rsc | lda ]
interface: ezfio, provider, ocaml
default: lda
[damped_mu_of_r]
type: logical
doc: If |true|, the mu(r) is damped by an error function to take a minimal value of mu_erf
interface: ezfio,provider,ocaml
default: False
[mu_of_r_min]
type: double precision
doc: minimal value of mu(r)
interface: ezfio,provider,ocaml
default: 0.001
ezfio_name: mu_of_r_min
[ampl_cos]
type: double precision
doc: amplitude of the cos for mu_test(r)
interface: ezfio,provider,ocaml
default: 0.1
ezfio_name: ampl_cos
[omega_cos]
type: double precision
doc: pulsation of the cos for mu_test(r)
interface: ezfio,provider,ocaml
default: 0.1
ezfio_name: omega_cos
[dexp_gauss]
type: double precision
doc: pulsation of the cos for mu_test(r)
interface: ezfio,provider,ocaml
default: 10.0
ezfio_name: dexp_gauss
[mu_test_choice]
type: character*(32)
doc: type of mu(r) for the TC Hamiltonian : can be [ cos | gauss ]
interface: ezfio, provider, ocaml
default: cos
[rescaled_on_top_mu]
type: logical
doc: If |true|, the mu(r) is rescaled by the ratio of ontop and density at HF level
interface: ezfio,provider,ocaml
default: False

2
src/some_mu_of_r/NEED Normal file
View File

@ -0,0 +1,2 @@
dft_utils_func
ao_two_e_erf_ints

144
src/some_mu_of_r/mu_grad_n.irp.f Normal file
View File

@ -0,0 +1,144 @@
BEGIN_PROVIDER [double precision, mu_of_r_extra_grid_grad_n , (n_points_extra_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
double precision :: weight, rho_a_hf, rho_b_hf, g0,rho_hf
double precision :: rs,grad_n
double precision :: g0_UEG_mu_inf
double precision :: cst_rs,alpha_rs
double precision :: elec_a,elec_b
double precision :: r(3),dx,mu_grad_n,mu_tmp,mu,damped_mu, mu_min
dx = 1.d-5
elec_a = 0.d0
elec_b = 0.d0
mu_min = mu_erf
do ipoint = 1, n_points_extra_final_grid
weight = final_weight_at_r_vector_extra(ipoint)
r(:) = final_grid_points_extra(:,ipoint)
mu_tmp = mu_grad_n(r)
if(damped_mu_of_r)then
mu = damped_mu(mu_tmp,mu_min)
else
mu = max(mu_tmp,mu_of_r_min)
endif
mu_of_r_extra_grid_grad_n(ipoint,1) = mu
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, average_mu_grad_n ]
&BEGIN_PROVIDER [double precision, mu_of_r_grad_n , (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_mu_of_r_grad_n , (3,n_points_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
double precision :: weight, rho_a_hf, rho_b_hf, g0,rho_hf
double precision :: rs,grad_n
double precision :: g0_UEG_mu_inf
double precision :: cst_rs,alpha_rs
double precision :: elec_a,elec_b,grad_mu(3),mu_grad_n
double precision :: mu, r(3),dx,damped_mu,mu_min,mu_tmp
mu_min = mu_erf
dx = 1.d-5
average_mu_grad_n = 0.d0
elec_a = 0.d0
elec_b = 0.d0
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
r(:) = final_grid_points(:,ipoint)
mu_tmp = mu_grad_n(r)
if(damped_mu_of_r)then
mu = damped_mu(mu_tmp,mu_min)
mu = max(mu_tmp,mu_of_r_min)
else
mu = mu_tmp
endif
mu_of_r_grad_n(ipoint,1) = mu
average_mu_grad_n += mu_of_r_grad_n(ipoint,1) * weight * rho_hf
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
r(:) = final_grid_points(:,ipoint)
if(damped_mu_of_r)then
call get_grad_damped_mu_grad_n(r,dx,mu_min,grad_mu)
else
call get_grad_mu_grad_n(r,dx,grad_mu)
if(mu_tmp.lt.mu)then
grad_mu = 0.d0
endif
endif
grad_mu_of_r_grad_n(:,ipoint,1) = grad_mu(:)
enddo
average_mu_grad_n = average_mu_grad_n / dble(elec_a+ elec_b)
END_PROVIDER
double precision function mu_grad_n(r)
implicit none
double precision, intent(in) :: r(3)
double precision, allocatable :: aos_array(:),grad_aos_array(:,:)
double precision, allocatable :: dm_a(:),dm_b(:), dm_a_grad(:,:), dm_b_grad(:,:)
double precision :: grad_n, rho
integer :: m
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))
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)
grad_n = 0.D0
do m = 1, 3
grad_n += dm_a_grad(m,1)**2 + dm_b_grad(m,1)**2 + 2.d0 * dm_a_grad(m,1)*dm_b_grad(m,1)
enddo
rho = dm_a(1) + dm_b(1)
grad_n = dsqrt(grad_n)
if(dabs(rho).gt.1.d-20)then
mu_grad_n = grad_n/(4.d0 * rho)
else
mu_grad_n = 0.d0
endif
end
double precision function mu_grad_n_damped(r,mu_min)
implicit none
double precision, intent(in) :: r(3), mu_min
double precision :: mu_grad_n, damped_mu, mu_tmp
mu_tmp = mu_grad_n(r)
mu_grad_n_damped = damped_mu(mu_tmp,mu_min)
end
subroutine get_grad_mu_grad_n(r,dx,grad_mu)
implicit none
double precision, intent(in) :: r(3),dx
double precision, intent(out):: grad_mu(3)
double precision :: r1(3),mu_plus,mu_minus,mu_grad_n
integer :: m
do m = 1, 3 ! compute grad mu
r1 = r
r1(m) += dx
mu_plus = mu_grad_n(r1)
r1 = r
r1(m) -= dx
mu_minus = mu_grad_n(r1)
grad_mu(m) = (mu_plus - mu_minus)/(2.d0 * dx)
enddo
end
subroutine get_grad_damped_mu_grad_n(r,dx,mu_min,grad_mu)
implicit none
double precision, intent(in) :: r(3),dx,mu_min
double precision, intent(out):: grad_mu(3)
double precision :: r1(3),mu_plus,mu_minus,mu_grad_n_damped
integer :: m
do m = 1, 3 ! compute grad mu
r1 = r
r1(m) += dx
mu_plus = mu_grad_n_damped(r1,mu_min)
r1 = r
r1(m) -= dx
mu_minus = mu_grad_n_damped(r1,mu_min)
grad_mu(m) = (mu_plus - mu_minus)/(2.d0 * dx)
enddo
end

208
src/some_mu_of_r/mu_lda.irp.f Normal file
View File

@ -0,0 +1,208 @@
BEGIN_PROVIDER [double precision, mu_of_r_extra_grid_lda , (n_points_extra_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
double precision :: weight, rho_a_hf, rho_b_hf, g0,rho_hf
double precision :: rs,grad_n
double precision :: g0_UEG_mu_inf
double precision :: cst_rs,alpha_rs
double precision :: elec_a,elec_b
double precision :: r(3),dx,mu_lda,mu_tmp,mu,damped_mu, mu_min
dx = 1.d-5
elec_a = 0.d0
elec_b = 0.d0
mu_min = mu_erf
do ipoint = 1, n_points_extra_final_grid
weight = final_weight_at_r_vector_extra(ipoint)
r(:) = final_grid_points_extra(:,ipoint)
call dm_dft_alpha_beta_at_r(r,rho_a_hf,rho_b_hf)
mu_tmp = mu_lda(rho_a_hf,rho_b_hf)
if(damped_mu_of_r)then
mu = damped_mu(mu_tmp,mu_min)
else
mu = max(mu_tmp,mu_of_r_min)
endif
if(rescaled_on_top_mu)then
mu = mu * dsqrt(0.25d0 * (rho_a_hf+rho_b_hf)**2/(rho_a_hf*rho_b_hf+1.d-12))
else
mu=mu
endif
mu_of_r_extra_grid_lda(ipoint,1) = mu
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, average_mu_lda ]
&BEGIN_PROVIDER [double precision, mu_of_r_lda , (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_mu_of_r_lda , (3,n_points_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
double precision :: weight, rho_a_hf, rho_b_hf, g0,rho_hf
double precision :: rs,grad_n
double precision :: g0_UEG_mu_inf
double precision :: cst_rs,alpha_rs
double precision :: elec_a,elec_b,grad_mu(3),mu_lda
double precision :: mu, r(3),dx,damped_mu,mu_min,mu_tmp
mu_min = mu_erf
dx = 1.d-5
average_mu_lda = 0.d0
elec_a = 0.d0
elec_b = 0.d0
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
rho_a_hf = one_e_dm_and_grad_alpha_in_r(4,ipoint,1)
rho_b_hf = one_e_dm_and_grad_beta_in_r(4,ipoint,1)
rho_hf = rho_a_hf + rho_b_hf
mu_tmp = mu_lda(rho_a_hf,rho_b_hf)
if(damped_mu_of_r)then
mu = damped_mu(mu_tmp,mu_min)
mu = max(mu_tmp,mu_of_r_min)
else
mu = mu_tmp
endif
if(rescaled_on_top_mu)then
mu = mu * dsqrt(0.25d0 * (rho_a_hf+rho_b_hf)**2/(rho_a_hf*rho_b_hf+1.d-12))
else
mu=mu
endif
mu_of_r_lda(ipoint,1) = mu
average_mu_lda += mu_of_r_lda(ipoint,1) * weight * rho_hf
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
r(:) = final_grid_points(:,ipoint)
if(damped_mu_of_r)then
call get_grad_damped_mu_lda(r,dx,mu_min,grad_mu)
else
call get_grad_mu_lda(r,dx,grad_mu)
if(mu_tmp.lt.mu)then
grad_mu = 0.d0
endif
endif
grad_mu_of_r_lda(:,ipoint,1) = grad_mu(:)
enddo
average_mu_lda = average_mu_lda / dble(elec_a+ elec_b)
END_PROVIDER
double precision function mu_lda(rho_a,rho_b)
implicit none
double precision, intent(in) :: rho_a,rho_b
include 'constants.include.F'
double precision :: g0,g0_UEG_mu_inf
g0 = g0_UEG_mu_inf(rho_a,rho_b)
mu_lda = - 1.d0 / (dlog(2.d0 * g0) * sqpi)
end
double precision function damped_mu(mu,mu_min)
implicit none
double precision, intent(in) :: mu, mu_min
damped_mu = mu_min * (1.d0 - derf(mu)) + derf(mu)*mu
end
double precision function mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
implicit none
double precision, intent(in) :: rho_b_hf,rho_a_hf,mu_min
double precision :: mu_lda,mu,damped_mu
mu = mu_lda(rho_a_hf,rho_b_hf)
mu_lda_damped = damped_mu(mu,mu_min)
end
subroutine get_grad_mu_lda(r,dx,grad_mu)
implicit none
double precision, intent(in) :: r(3),dx
double precision, intent(out):: grad_mu(3)
double precision :: r1(3),rho_a_hf,rho_b_hf,mu_plus,mu_minus,mu_lda
integer :: m
do m = 1, 3 ! compute grad mu
r1 = r
r1(m) += dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_plus = mu_lda(rho_a_hf,rho_b_hf)
r1 = r
r1(m) -= dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_minus = mu_lda(rho_a_hf,rho_b_hf)
grad_mu(m) = (mu_plus - mu_minus)/(2.d0 * dx)
enddo
end
subroutine get_grad_damped_mu_lda(r,dx,mu_min,grad_mu)
implicit none
double precision, intent(in) :: r(3),dx,mu_min
double precision, intent(out):: grad_mu(3)
double precision :: r1(3),rho_a_hf,rho_b_hf,mu_plus,mu_minus,mu_lda_damped
integer :: m
do m = 1, 3 ! compute grad mu
r1 = r
r1(m) += dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_plus = mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
r1 = r
r1(m) -= dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_minus = mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
grad_mu(m) = (mu_plus - mu_minus)/(2.d0 * dx)
enddo
end
double precision function grad_mu_lda_comp(r,dx,mu_min,m)
implicit none
double precision, intent(in) :: r(3),dx,mu_min
integer, intent(in) :: m
double precision :: r1(3),rho_a_hf,rho_b_hf,mu_plus,mu_minus,mu_lda_damped
r1 = r
r1(m) += dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_plus = mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
r1 = r
r1(m) -= dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_minus = mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
grad_mu_lda_comp = (mu_plus - mu_minus)/(2.d0 * dx)
end
subroutine get_lapl_mu_lda(r,dx,mu_min,lapl_mu)
implicit none
double precision, intent(in) :: r(3),dx,mu_min
double precision, intent(out):: lapl_mu(3)
double precision :: r1(3),rho_a_hf,rho_b_hf,mu_plus,mu_minus,mu_lda_damped,mu
integer :: m
do m = 1, 3
r1 = r
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu = mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
r1 = r
r1(m) += 2.d0 * dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_plus = mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
r1 = r
r1(m) -= 2.d0 * dx
call dm_dft_alpha_beta_at_r(r1,rho_a_hf,rho_b_hf)
mu_minus = mu_lda_damped(rho_a_hf,rho_b_hf,mu_min)
lapl_mu(m) = (mu_plus + mu_minus - 2.d0 * mu)/(4.d0 * dx * dx)
enddo
lapl_mu = 0.d0
end
BEGIN_PROVIDER [ double precision, average_mu_rs_c_lda]
implicit none
average_mu_rs_c_lda = 0.5d0 * (average_mu_rs_c + average_mu_lda)
END_PROVIDER

150
src/some_mu_of_r/mu_of_r_ints.irp.f Normal file
View File

@ -0,0 +1,150 @@
BEGIN_PROVIDER [double precision, average_mu_of_r_for_ints, (N_states) ]
&BEGIN_PROVIDER [double precision, average_grad_mu_of_r, (N_states) ]
&BEGIN_PROVIDER [double precision, av_grad_inv_mu_mu_of_r, (N_states) ]
&BEGIN_PROVIDER [double precision, inv_2_mu_of_r_for_ints, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, mu_of_r_for_ints, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, inv_4_mu_of_r_for_ints, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_mu_of_r_for_ints, (3,n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_mu_of_r_transp_for_ints, (n_points_final_grid,N_states,3) ]
&BEGIN_PROVIDER [double precision, grad_sq_mu_of_r_for_ints, (n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
!
! mu(r) and its gradient for evaluation f integrals for the TC hamiltonian
END_DOC
integer :: ipoint,istate,mm
double precision :: wall0,wall1,mu_max
mu_max = 1.d-4
print*,'providing mu_of_r ...'
call wall_time(wall0)
! if(.not.constant_mu)then
! do istate = 1, N_states
! do ipoint = 1, n_points_final_grid
! if(mu_of_r_tc_ints.EQ."basis")then
! mu_of_r_for_ints(ipoint,istate) = mu_of_r_basis_hf(ipoint)
! grad_mu_of_r_for_ints(:,ipoint,istate) = grad_mu_of_r_basis_hf(:,ipoint)
! if(mu_of_r_tc_ints.EQ."rsc")then
! mu_of_r_for_ints(ipoint,istate) = mu_of_r_rs_c(ipoint,istate)
! grad_mu_of_r_for_ints(:,ipoint,istate) = grad_mu_of_r_rs_c(:,ipoint,istate)
! else if(mu_of_r_tc_ints.EQ."lda")then
! mu_of_r_for_ints(ipoint,istate) = mu_of_r_lda(ipoint,istate)
! grad_mu_of_r_for_ints(:,ipoint,istate) = grad_mu_of_r_lda(:,ipoint,istate)
! else if(mu_of_r_tc_ints.EQ."rsc_lda")then
! mu_of_r_for_ints(ipoint,istate) = 0.5d0 * (mu_of_r_lda(ipoint,istate) + mu_of_r_rs_c(ipoint,istate))
! grad_mu_of_r_for_ints(:,ipoint,istate) = 0.5d0 * (grad_mu_of_r_lda(:,ipoint,istate) + grad_mu_of_r_rs_c(:,ipoint,istate))
! else if(mu_of_r_tc_ints.EQ."grad_n")then
! mu_of_r_for_ints(ipoint,istate) = mu_of_r_grad_n(ipoint,istate)
! grad_mu_of_r_for_ints(:,ipoint,istate) = grad_mu_of_r_grad_n(:,ipoint,istate)
! else if(mu_of_r_tc_ints.EQ."mu_test")then
! mu_of_r_for_ints(ipoint,istate) = mu_of_r_test_func(ipoint,istate)
! grad_mu_of_r_for_ints(:,ipoint,istate) = grad_mu_of_r_test_func(:,ipoint,istate)
! else
! print*,'you requested the following mu_of_r_tc_ints'
! print*,mu_of_r_tc_ints
! print*,'which does not correspond to any of the options for such keyword'
! stop
! endif
! enddo
! enddo
! else
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
mu_of_r_for_ints(ipoint,istate) = mu_erf
grad_mu_of_r_for_ints(:,ipoint,istate) = 0.d0
grad_sq_mu_of_r_for_ints(ipoint,istate) = 0.d0
enddo
enddo
! endif
! do istate = 1, N_states
! do ipoint = 1, n_points_final_grid
! mu_of_r_for_ints(ipoint,istate) = max(mu_of_r_for_ints(ipoint,istate),mu_max)
! inv_2_mu_of_r_for_ints(ipoint,istate) = 1.d0/(mu_of_r_for_ints(ipoint,istate))**2
! inv_4_mu_of_r_for_ints(ipoint,istate) = 1.d0/(mu_of_r_for_ints(ipoint,istate))**4
! do mm = 1, 3
! grad_mu_of_r_transp_for_ints(ipoint,istate,mm) = grad_mu_of_r_for_ints(mm,ipoint,istate)
! enddo
!
! grad_sq_mu_of_r_for_ints(ipoint,istate) = 0.d0
! do mm = 1, 3
! grad_sq_mu_of_r_for_ints(ipoint,istate) += grad_mu_of_r_for_ints(mm,ipoint,istate)**2.d0
! enddo
! enddo
! enddo
double precision :: elec_tot,dm,weight,grad_mu_sq
average_grad_mu_of_r = 0.d0
average_mu_of_r_for_ints = 0.d0
av_grad_inv_mu_mu_of_r = 0.d0
! do istate = 1, N_states
! elec_tot = 0.d0
! do ipoint = 1, n_points_final_grid
! dm = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! weight = final_weight_at_r_vector_extra(ipoint)
! average_grad_mu_of_r(istate) += dsqrt(grad_sq_mu_of_r_for_ints(ipoint,istate)) * dm * weight
! average_mu_of_r_for_ints(istate) += mu_of_r_for_ints(ipoint,istate) * dm * weight
! av_grad_inv_mu_mu_of_r(istate) += mu_of_r_for_ints(ipoint,istate)**(-2) * dsqrt(grad_sq_mu_of_r_for_ints(ipoint,istate)) * dm * weight
! elec_tot += dm * weight
! enddo
! average_mu_of_r_for_ints(istate) = average_mu_of_r_for_ints(istate) / elec_tot
! average_grad_mu_of_r(istate) = average_grad_mu_of_r(istate) / elec_tot
! av_grad_inv_mu_mu_of_r(istate) = av_grad_inv_mu_mu_of_r(istate)/elec_tot
! enddo
call wall_time(wall1)
print*,'Time to provide mu_of_r_for_ints = ',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_of_r_extra_grid_for_ints, (n_points_extra_final_grid,N_states) ]
implicit none
BEGIN_DOC
!
! mu(r) and its gradient for evaluation f integrals for the TC hamiltonian
END_DOC
integer :: ipoint,istate,mm,mu_tmp
double precision :: wall0,wall1,mu_max
mu_max = 1.d-4
print*,'providing mu_of_r_extra_grid ...'
call wall_time(wall0)
! if(.not.constant_mu)then
! do istate = 1, N_states
! do ipoint = 1, n_points_extra_final_grid
! if(mu_of_r_tc_ints.EQ."basis")then
! mu_of_r_extra_grid_for_ints(ipoint,istate) = mu_of_r_extra_grid_basis_hf(ipoint)
! if(mu_of_r_tc_ints.EQ."rsc")then
! mu_of_r_extra_grid_for_ints(ipoint,istate) = mu_of_r_extra_grid_rs_c(ipoint,istate)
! else if(mu_of_r_tc_ints.EQ."lda")then
! mu_of_r_extra_grid_for_ints(ipoint,istate) = mu_of_r_extra_grid_lda(ipoint,istate)
! else if(mu_of_r_tc_ints.EQ."rsc_lda")then
! mu_of_r_extra_grid_for_ints(ipoint,istate) = 0.5d0 * (mu_of_r_extra_grid_lda(ipoint,istate) + mu_of_r_extra_grid_rs_c(ipoint,istate))
! else if(mu_of_r_tc_ints.EQ."grad_n")then
! mu_of_r_extra_grid_for_ints(ipoint,istate) = mu_of_r_extra_grid_grad_n(ipoint,istate)
! else if(mu_of_r_tc_ints.EQ."mu_test")then
! mu_of_r_extra_grid_for_ints(ipoint,istate) = mu_of_r_extra_grid_test_func(ipoint,istate)
! else
! print*,'you requested the following mu_of_r_extra_grid_for_ints'
! print*,mu_of_r_tc_ints
! print*,'which does not correspond to any of the options for such keyword'
! stop
! endif
! mu_of_r_extra_grid_for_ints(ipoint,istate) = max(mu_of_r_extra_grid_for_ints(ipoint,istate),mu_max)
! enddo
! enddo
! else
do istate = 1, N_states
do ipoint = 1, n_points_extra_final_grid
mu_of_r_extra_grid_for_ints(ipoint,istate) = mu_erf
enddo
enddo
! endif
call wall_time(wall1)
print*,'Time to provide mu_of_r_extra_grid_for_ints = ',wall1-wall0
END_PROVIDER

162
src/some_mu_of_r/mu_rsc.irp.f Normal file
View File

@ -0,0 +1,162 @@
BEGIN_PROVIDER [double precision, average_mu_rs_c ]
&BEGIN_PROVIDER [double precision, mu_of_r_rs_c, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_mu_of_r_rs_c, (3,n_points_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
include 'constants.include.F'
double precision :: weight, rho_a_hf, rho_b_hf, rho_hf, mu_rs_c,r(3)
average_mu_rs_c = 0.d0
double precision :: elec_a,elec_b,grad_mu(3)
double precision :: dx,mu_lda, mu_tmp,mu_min,mu,damped_mu
dx = 1.d-5
elec_a = 0.d0
elec_b = 0.d0
mu_min = mu_erf
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
r(:) = final_grid_points(:,ipoint)
rho_a_hf = one_e_dm_and_grad_alpha_in_r(4,ipoint,1)
rho_b_hf = one_e_dm_and_grad_beta_in_r(4,ipoint,1)
rho_hf = rho_a_hf + rho_b_hf
rho_hf = max(rho_hf,1.D-10)
mu_tmp = mu_rs_c(rho_hf)
if(damped_mu_of_r)then
mu = damped_mu(mu_tmp,mu_min)
mu = max(mu_tmp,mu_of_r_min)
else
mu = mu_tmp
endif
if(rescaled_on_top_mu)then
mu = mu *dsqrt(0.25d0 * (rho_hf+rho_hf)**2/(rho_a_hf*rho_b_hf+1.d-12))
else
mu=mu
endif
mu_of_r_rs_c(ipoint,1) = mu
average_mu_rs_c += mu_of_r_rs_c(ipoint,1) * rho_hf * weight
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
if(damped_mu_of_r)then
call get_grad_damped_mu_rsc(r,mu_min,dx,grad_mu)
else
call get_grad_mu_rsc(r,dx,grad_mu)
endif
grad_mu_of_r_rs_c(:,ipoint,1) = grad_mu(:)
enddo
average_mu_rs_c = average_mu_rs_c / dble(elec_a+ elec_b)
END_PROVIDER
!double precision function mu_rs_c(rho)
! implicit none
! double precision, intent(in) :: rho
! include 'constants.include.F'
! double precision :: cst_rs,alpha_rs,rs
! cst_rs = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0)
! alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi
!
! rs = cst_rs * rho**(-1.d0/3.d0)
! mu_rs_c = alpha_rs/dsqrt(rs)
!
!end
double precision function damped_mu_rs_c(rho,mu_min)
implicit none
double precision, intent(in) :: rho,mu_min
double precision :: mu_rs_c,mu_tmp,damped_mu
mu_tmp = mu_rs_c(rho)
damped_mu_rs_c = damped_mu(mu_tmp, mu_min)
end
subroutine get_grad_mu_rsc(r,dx,grad_mu)
implicit none
double precision, intent(in) :: r(3),dx
double precision, intent(out):: grad_mu(3)
integer :: m
double precision :: r_tmp(3),mu_plus, mu_minus,mu_rs_c,rho_a_hf,rho_b_hf,rho_hf
do m = 1, 3
r_tmp = r
r_tmp(m) += dx
call dm_dft_alpha_beta_at_r(r_tmp,rho_a_hf,rho_b_hf)
rho_hf = rho_a_hf + rho_b_hf
rho_hf = max(rho_hf,1.D-10)
mu_plus = mu_rs_c(rho_hf)
r_tmp = r
r_tmp(m) -= dx
call dm_dft_alpha_beta_at_r(r_tmp,rho_a_hf,rho_b_hf)
rho_hf = rho_a_hf + rho_b_hf
rho_hf = max(rho_hf,1.D-10)
mu_minus = mu_rs_c(rho_hf)
grad_mu(m) = (mu_plus - mu_minus)/(2.d0 * dx)
enddo
end
subroutine get_grad_damped_mu_rsc(r,mu_min,dx,grad_mu)
implicit none
double precision, intent(in) :: r(3),mu_min,dx
double precision, intent(out):: grad_mu(3)
integer :: m
double precision :: r_tmp(3),mu_plus, mu_minus,mu_rs_c,mu_tmp
double precision :: rho_a_hf,rho_b_hf,rho_hf,damped_mu,rho,damped_mu_rs_c
do m = 1, 3
r_tmp = r
r_tmp(m) += dx
call dm_dft_alpha_beta_at_r(r_tmp,rho_a_hf,rho_b_hf)
rho_hf = rho_a_hf + rho_b_hf
rho_hf = max(rho_hf,1.D-10)
mu_plus = damped_mu_rs_c(rho_hf,mu_min)
r_tmp = r
r_tmp(m) -= dx
call dm_dft_alpha_beta_at_r(r_tmp,rho_a_hf,rho_b_hf)
rho_hf = rho_a_hf + rho_b_hf
rho_hf = max(rho_hf,1.D-10)
mu_minus = damped_mu_rs_c(rho_hf,mu_min)
grad_mu(m) = (mu_plus - mu_minus)/(2.d0 * dx)
enddo
end
BEGIN_PROVIDER [double precision, mu_of_r_extra_grid_rs_c, (n_points_extra_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
include 'constants.include.F'
double precision :: weight, rho_a_hf, rho_b_hf, rho_hf, mu_rs_c
average_mu_rs_c = 0.d0
double precision :: elec_a,elec_b,r(3)
double precision :: mu_tmp ,damped_mu, mu ,dx,mu_min
dx = 1.d-5
elec_a = 0.d0
elec_b = 0.d0
mu_min = mu_erf
do ipoint = 1, n_points_extra_final_grid
weight = final_weight_at_r_vector_extra(ipoint)
r(:) = final_grid_points_extra(:,ipoint)
call dm_dft_alpha_beta_at_r(r,rho_a_hf,rho_b_hf)
rho_hf = rho_a_hf + rho_b_hf
rho_hf = max(rho_hf,1.D-10)
mu_tmp = mu_rs_c(rho_hf)
if(damped_mu_of_r)then
mu = damped_mu(mu_tmp,mu_min)
else
mu = max(mu_tmp,mu_of_r_min)
endif
if(rescaled_on_top_mu)then
mu = mu *dsqrt(0.25d0 * (rho_a_hf+rho_b_hf)**2/(rho_a_hf*rho_b_hf+1.d-12))
else
mu=mu
endif
mu_of_r_extra_grid_rs_c(ipoint,1) = mu
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
enddo
END_PROVIDER

119
src/some_mu_of_r/mu_test.irp.f Normal file
View File

@ -0,0 +1,119 @@
BEGIN_PROVIDER [double precision, mu_of_r_extra_grid_test_func , (n_points_extra_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
double precision :: weight, rho_a_hf, rho_b_hf, g0,rho_hf
double precision :: rs,grad_n
double precision :: g0_UEG_mu_inf
double precision :: cst_rs,alpha_rs
double precision :: elec_a,elec_b
double precision :: r(3),dx,mu_test_func,mu_tmp,mu,damped_mu, mu_min
dx = 1.d-5
elec_a = 0.d0
elec_b = 0.d0
mu_min = mu_erf
do ipoint = 1, n_points_extra_final_grid
weight = final_weight_at_r_vector_extra(ipoint)
r(:) = final_grid_points_extra(:,ipoint)
mu = mu_test_func(r)
mu_of_r_extra_grid_test_func(ipoint,1) = mu
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, average_mu_test_func ]
&BEGIN_PROVIDER [double precision, mu_of_r_test_func , (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_mu_of_r_test_func , (3,n_points_final_grid,N_states) ]
implicit none
integer :: ipoint,i,m
double precision :: weight, rho_a_hf, rho_b_hf, g0,rho_hf
double precision :: rs,grad_n
double precision :: g0_UEG_mu_inf
double precision :: cst_rs,alpha_rs
double precision :: elec_a,elec_b,grad_mu(3),mu_test_func
double precision :: mu, r(3),dx,damped_mu,mu_min,mu_tmp
mu_min = mu_erf
dx = 1.d-5
average_mu_test_func = 0.d0
elec_a = 0.d0
elec_b = 0.d0
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
r(:) = final_grid_points(:,ipoint)
mu_tmp = mu_test_func(r)
mu_of_r_test_func(ipoint,1) = mu_tmp
average_mu_test_func += mu_of_r_test_func(ipoint,1) * weight * rho_hf
elec_a += rho_a_hf * weight
elec_b += rho_b_hf * weight
call get_grad_mu_test_func(r,dx,grad_mu)
if(mu_tmp.lt.mu)then
grad_mu = 0.d0
endif
grad_mu_of_r_test_func(:,ipoint,1) = grad_mu(:)
enddo
average_mu_test_func = average_mu_test_func / dble(elec_a+ elec_b)
END_PROVIDER
double precision function mu_test_func(r)
implicit none
double precision, intent(in) :: r(3)
double precision :: x,y,z
x = r(1)
y = r(2)
z = r(3)
if(mu_test_choice == "cos")then
mu_test_func = mu_erf + ampl_cos * dcos(omega_cos * x) * dcos(omega_cos * y) * dcos(omega_cos * z) !&
! * dexp(-dexp_gauss * (x**2 + y**2 + z**2))
else if(mu_test_choice== "gauss" )then
mu_test_func = mu_erf + ampl_cos * dexp(-dexp_gauss * (x**2 + y**2 + z**2))
endif
mu_test_func = max(mu_of_r_min,mu_test_func)
end
subroutine get_grad_mu_test_func(r,dx,grad_mu)
implicit none
double precision, intent(in) :: r(3),dx
double precision, intent(out):: grad_mu(3)
double precision :: x,y,z
if(mu_test_choice == "cos")then
x = r(1)
y = r(2)
z = r(3)
grad_mu(1) = -ampl_cos * omega_cos * dsin(omega_cos * x) * dcos(omega_cos * y) * dcos(omega_cos * z)
grad_mu(2) = -ampl_cos * omega_cos * dcos(omega_cos * x) * dsin(omega_cos * y) * dcos(omega_cos * z)
grad_mu(3) = -ampl_cos * omega_cos * dcos(omega_cos * x) * dcos(omega_cos * y) * dsin(omega_cos * z)
else
call get_grad_mu_test_num(r,dx,grad_mu)
endif
end
subroutine get_grad_mu_test_num(r,dx,grad_mu)
implicit none
double precision, intent(in) :: r(3), dx
double precision, intent(out):: grad_mu(3)
double precision :: r1(3),mu_plus,mu_minus,mu_test_func
integer :: m
do m = 1, 3 ! compute grad mu
r1 = r
r1(m) += dx
mu_plus = mu_test_func(r1)
r1 = r
r1(m) -= dx
mu_minus = mu_test_func(r1)
grad_mu(m) = (mu_plus - mu_minus)/(2.d0 * dx)
enddo
end

7
src/some_mu_of_r/some_mu_of_r.irp.f Normal file
View File

@ -0,0 +1,7 @@
program some_mu_of_r
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
end

107
src/tc_keywords/EZFIO.cfg Normal file
View File

@ -0,0 +1,107 @@
[read_rl_eigv]
type: logical
doc: If |true|, read the right/left eigenvectors from ezfio
interface: ezfio,provider,ocaml
default: False
[comp_left_eigv]
type: logical
doc: If |true|, computes also the left-eigenvector
interface: ezfio,provider,ocaml
default: False
[three_body_h_tc]
type: logical
doc: If |true|, three-body terms are included
interface: ezfio,provider,ocaml
default: True
[pure_three_body_h_tc]
type: logical
doc: If |true|, pure triple excitation three-body terms are included
interface: ezfio,provider,ocaml
default: False
[double_normal_ord]
type: logical
doc: If |true|, contracted double excitation three-body terms are included
interface: ezfio,provider,ocaml
default: False
[core_tc_op]
type: logical
doc: If |true|, takes the usual Hamiltonian for core orbitals (assumed to be doubly occupied)
interface: ezfio,provider,ocaml
default: False
[full_tc_h_solver]
type: logical
doc: If |true|, you diagonalize the full TC H matrix
interface: ezfio,provider,ocaml
default: False
[thresh_it_dav]
type: Threshold
doc: Thresholds on the energy for iterative Davidson used in TC
interface: ezfio,provider,ocaml
default: 1.e-5
[max_it_dav]
type: integer
doc: nb max of iteration in Davidson used in TC
interface: ezfio,provider,ocaml
default: 1000
[thresh_psi_r]
type: Threshold
doc: Thresholds on the coefficients of the right-eigenvector. Used for PT2 computation.
interface: ezfio,provider,ocaml
default: 0.000005
[thresh_psi_r_norm]
type: logical
doc: If |true|, you prune the WF to compute the PT1 coef based on the norm. If False, the pruning is done through the amplitude on the right-coefficient.
interface: ezfio,provider,ocaml
default: False
[state_following_tc]
type: logical
doc: If |true|, the states are re-ordered to match the input states
default: False
interface: ezfio,provider,ocaml
[bi_ortho]
type: logical
doc: If |true|, the MO basis is assumed to be bi-orthonormal
interface: ezfio,provider,ocaml
default: True
[symetric_fock_tc]
type: logical
doc: If |true|, using F+F^\dagger as Fock TC
interface: ezfio,provider,ocaml
default: False
[thresh_tcscf]
type: Threshold
doc: Threshold on the convergence of the Hartree Fock energy.
interface: ezfio,provider,ocaml
default: 1.e-10
[n_it_tcscf_max]
type: Strictly_positive_int
doc: Maximum number of SCF iterations
interface: ezfio,provider,ocaml
default: 500
[max_ov_tc_scf]
type: logical
doc: If |true|, the TC-SCF is done with a kind of maximum overlap with RHF MOs
interface: ezfio,provider,ocaml
default: True
[selection_tc]
type: integer
doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative
interface: ezfio,provider,ocaml
default: 0

1
src/tc_keywords/NEED Normal file
View File

@ -0,0 +1 @@
ezfio_files

7
src/tc_keywords/tc_keywords.irp.f Normal file
View File

@ -0,0 +1,7 @@
program tc_keywords
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
end

20
src/three_body_ints/EZFIO.cfg Normal file
View File

@ -0,0 +1,20 @@
[io_three_body_ints]
type: Disk_access
doc: Read/Write the 6 index tensor three-body terms from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[symm_3_body_tensor]
type: logical
doc: If |true|, you have a symmetrized two body tensor
interface: ezfio,provider,ocaml
default: False
[read_3_body_tc_ints]
type: logical
doc: If |true|, you read the 3 body integrals from an FCIDUMP like file
interface: ezfio,provider,ocaml
default: False

1
src/three_body_ints/NEED Normal file
View File

@ -0,0 +1 @@
bi_ort_ints

63
src/three_body_ints/io_6_index_tensor.irp.f Normal file
View File

@ -0,0 +1,63 @@
subroutine write_array_6_index_tensor(n_orb,array_tmp,name_file)
implicit none
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,n_orb,n_orb)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'W')
write(i_unit_output)array_tmp
close(unit=i_unit_output)
end
subroutine read_array_6_index_tensor(n_orb,array_tmp,name_file)
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,n_orb,n_orb)
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'R')
read(i_unit_output)array_tmp
close(unit=i_unit_output)
end
subroutine read_fcidump_3_tc(array)
implicit none
double precision, intent(out) :: array(mo_num, mo_num, mo_num, mo_num, mo_num, mo_num)
integer :: i,j,k,l,m,n,i_mo, Reason
double precision :: integral
print*,'Reading the THREE-body integrals from a TC FCIDUMP'
open (unit=15, file="TCDUMP-nosym", status='old', &
access='sequential', action='read' )
read(15,*)i_mo
if(i_mo.ne.mo_num)then
print*,'Something went wrong in the read_fcidump_3_tc !'
print*,'i_mo.ne.mo_num !'
print*,i_mo,mo_num
stop
endif
do
read(15,*,IOSTAT=Reason)integral,i, j, m, k, l, n
if(Reason > 0)then
print*,'Something went wrong in the I/O of read_fcidump_3_tc'
stop
else if(Reason < 0)then
exit
else
! 1 2 3 1 2 3
! <i j m|k l n>
! (ik|jl|mn)
! integral = integral * 1.d0/3.d0 !!!! For NECI convention
array(i,j,m,k,l,n) = integral * 3.d0
endif
enddo
end

207
src/three_body_ints/semi_num_ints_mo.irp.f Normal file
View File

@ -0,0 +1,207 @@
BEGIN_PROVIDER [ double precision, mo_v_ij_erf_rk_cst_mu_naive, ( mo_num, mo_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1 )/(2|r - R|) on the MO basis
END_DOC
integer :: i,j,k,l,ipoint
do ipoint = 1, n_points_final_grid
mo_v_ij_erf_rk_cst_mu_naive(:,:,ipoint) = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, ao_num
do l = 1, ao_num
mo_v_ij_erf_rk_cst_mu_naive(j,i,ipoint) += mo_coef(l,j) * 0.5d0 * v_ij_erf_rk_cst_mu(l,k,ipoint) * mo_coef(k,i)
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_v_ij_erf_rk_cst_mu, ( mo_num, mo_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the MO basis
END_DOC
integer :: ipoint
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint) &
!$OMP SHARED (n_points_final_grid,v_ij_erf_rk_cst_mu,mo_v_ij_erf_rk_cst_mu)
!$OMP DO SCHEDULE (dynamic)
do ipoint = 1, n_points_final_grid
call ao_to_mo(v_ij_erf_rk_cst_mu(1,1,ipoint),size(v_ij_erf_rk_cst_mu,1),mo_v_ij_erf_rk_cst_mu(1,1,ipoint),size(mo_v_ij_erf_rk_cst_mu,1))
enddo
!$OMP END DO
!$OMP END PARALLEL
mo_v_ij_erf_rk_cst_mu = mo_v_ij_erf_rk_cst_mu * 0.5d0
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_v_ij_erf_rk_cst_mu_transp, ( n_points_final_grid,mo_num, mo_num)]
implicit none
BEGIN_DOC
! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the MO basis
END_DOC
integer :: ipoint,i,j
do i = 1, mo_num
do j = 1, mo_num
do ipoint = 1, n_points_final_grid
mo_v_ij_erf_rk_cst_mu_transp(ipoint,j,i) = mo_v_ij_erf_rk_cst_mu(j,i,ipoint)
enddo
enddo
enddo
FREE mo_v_ij_erf_rk_cst_mu
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_x_v_ij_erf_rk_cst_mu_naive, ( mo_num, mo_num,3,n_points_final_grid)]
implicit none
BEGIN_DOC
! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1 )/|r - R| on the MO basis
END_DOC
integer :: i,j,k,l,ipoint,m
do ipoint = 1, n_points_final_grid
mo_x_v_ij_erf_rk_cst_mu_naive(:,:,:,ipoint) = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do m = 1, 3
do k = 1, ao_num
do l = 1, ao_num
mo_x_v_ij_erf_rk_cst_mu_naive(j,i,m,ipoint) += mo_coef(l,j) * 0.5d0 * x_v_ij_erf_rk_cst_mu_transp(l,k,m,ipoint) * mo_coef(k,i)
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_x_v_ij_erf_rk_cst_mu, ( mo_num, mo_num,3,n_points_final_grid)]
implicit none
BEGIN_DOC
! int dr x * phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/2|r - R| on the MO basis
END_DOC
integer :: ipoint,m
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,m) &
!$OMP SHARED (n_points_final_grid,x_v_ij_erf_rk_cst_mu_transp,mo_x_v_ij_erf_rk_cst_mu)
!$OMP DO SCHEDULE (dynamic)
do ipoint = 1, n_points_final_grid
do m = 1, 3
call ao_to_mo(x_v_ij_erf_rk_cst_mu_transp(1,1,m,ipoint),size(x_v_ij_erf_rk_cst_mu_transp,1),mo_x_v_ij_erf_rk_cst_mu(1,1,m,ipoint),size(mo_x_v_ij_erf_rk_cst_mu,1))
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
mo_x_v_ij_erf_rk_cst_mu = 0.5d0 * mo_x_v_ij_erf_rk_cst_mu
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_x_v_ij_erf_rk_cst_mu_transp, (n_points_final_grid,3, mo_num, mo_num)]
implicit none
integer :: i,j,m,ipoint
do i = 1, mo_num
do j = 1, mo_num
do m = 1, 3
do ipoint = 1, n_points_final_grid
mo_x_v_ij_erf_rk_cst_mu_transp(ipoint,m,j,i) = mo_x_v_ij_erf_rk_cst_mu(j,i,m,ipoint)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, x_W_ij_erf_rk, ( n_points_final_grid,3,mo_num, mo_num)]
implicit none
BEGIN_DOC
! W_mn^X(R) = \int dr phi_m(r) phi_n(r) (1 - erf(mu |r-R|)) (x-X)
END_DOC
include 'constants.include.F'
integer :: ipoint,m,i,j
double precision :: xyz,cst
double precision :: wall0, wall1
cst = 0.5d0 * inv_sq_pi
print*,'providing x_W_ij_erf_rk ...'
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,m,i,j,xyz) &
!$OMP SHARED (x_W_ij_erf_rk,n_points_final_grid,mo_x_v_ij_erf_rk_cst_mu_transp,mo_v_ij_erf_rk_cst_mu_transp,mo_num,final_grid_points)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = 1, 3
do ipoint = 1, n_points_final_grid
xyz = final_grid_points(m,ipoint)
x_W_ij_erf_rk(ipoint,m,j,i) = mo_x_v_ij_erf_rk_cst_mu_transp(ipoint,m,j,i) - xyz * mo_v_ij_erf_rk_cst_mu_transp(ipoint,j,i)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
FREE mo_v_ij_erf_rk_cst_mu_transp
FREE mo_x_v_ij_erf_rk_cst_mu_transp
call wall_time(wall1)
print*,'time to provide x_W_ij_erf_rk = ',wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, sqrt_weight_at_r, (n_points_final_grid)]
implicit none
integer :: ipoint
do ipoint = 1, n_points_final_grid
sqrt_weight_at_r(ipoint) = dsqrt(final_weight_at_r_vector(ipoint))
enddo
END_PROVIDER
!BEGIN_PROVIDER [ double precision, mos_in_r_array_transp_sq_weight, (n_points_final_grid,mo_num)]
!BEGIN_PROVIDER [ double precision, gauss_ij_rk_transp, (ao_num, ao_num, n_points_final_grid) ]
! implicit none
! integer :: i,j,ipoint
! do ipoint = 1, n_points_final_grid
! do j = 1, ao_num
! do i = 1, ao_num
! gauss_ij_rk_transp(i,j,ipoint) = gauss_ij_rk(ipoint,i,j)
! enddo
! enddo
! enddo
!END_PROVIDER
!
!
!BEGIN_PROVIDER [ double precision, mo_gauss_ij_rk, ( mo_num, mo_num,n_points_final_grid)]
! implicit none
! integer :: ipoint
! !$OMP PARALLEL &
! !$OMP DEFAULT (NONE) &
! !$OMP PRIVATE (ipoint) &
! !$OMP SHARED (n_points_final_grid,gauss_ij_rk_transp,mo_gauss_ij_rk)
! !$OMP DO SCHEDULE (dynamic)
! do ipoint = 1, n_points_final_grid
! call ao_to_mo(gauss_ij_rk_transp(1,1,ipoint),size(gauss_ij_rk_transp,1),mo_gauss_ij_rk(1,1,ipoint),size(mo_gauss_ij_rk,1))
! enddo
! !$OMP END DO
! !$OMP END PARALLEL
!
!END_PROVIDER
!
!BEGIN_PROVIDER [ double precision, mo_gauss_ij_rk_transp, (n_points_final_grid, mo_num, mo_num)]
! implicit none
! integer :: i,j,ipoint
! do ipoint = 1, n_points_final_grid
! do j = 1, mo_num
! do i = 1, mo_num
! mo_gauss_ij_rk_transp(ipoint,i,j) = mo_gauss_ij_rk(i,j,ipoint)
! enddo
! enddo
! enddo
!
!END_PROVIDER
!

106
src/three_body_ints/three_body_tensor.irp.f Normal file
View File

@ -0,0 +1,106 @@
BEGIN_PROVIDER [ double precision, three_body_ints, (mo_num, mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! matrix element of the -L three-body operator
!
! notice the -1 sign: in this way three_body_ints can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_ints = 0.d0
print*,'Providing the three_body_ints ...'
call wall_time(wall0)
name_file = 'six_index_tensor'
if(read_three_body_ints)then
call read_fcidump_3_tc(three_body_ints)
else
if(read_three_body_ints)then
print*,'Reading three_body_ints from disk ...'
call read_array_6_index_tensor(mo_num,three_body_ints,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_ints)
!$OMP DO SCHEDULE (dynamic)
do n = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
do m = n, mo_num
do j = l, mo_num
do i = k, mo_num
!! if(i>=j)then
integral = 0.d0
call give_integrals_3_body(i,j,m,k,l,n,integral)
three_body_ints(i,j,m,k,l,n) = -1.d0 * integral
! permutation with k,i
three_body_ints(k,j,m,i,l,n) = -1.d0 * integral ! i,k
! two permutations with k,i
three_body_ints(k,l,m,i,j,n) = -1.d0 * integral
three_body_ints(k,j,n,i,l,m) = -1.d0 * integral
! three permutations with k,i
three_body_ints(k,l,n,i,j,m) = -1.d0 * integral
! permutation with l,j
three_body_ints(i,l,m,k,j,n) = -1.d0 * integral ! j,l
! two permutations with l,j
three_body_ints(k,l,m,i,j,n) = -1.d0 * integral
three_body_ints(i,l,n,k,j,m) = -1.d0 * integral
! two permutations with l,j
!!!! three_body_ints(k,l,n,i,j,m) = -1.d0 * integral
! permutation with m,n
three_body_ints(i,j,n,k,l,m) = -1.d0 * integral ! m,n
! two permutations with m,n
three_body_ints(k,j,n,i,l,m) = -1.d0 * integral ! m,n
three_body_ints(i,l,n,k,j,m) = -1.d0 * integral ! m,n
! three permutations with k,i
!!!! three_body_ints(k,l,n,i,j,m) = -1.d0 * integral ! m,n
!! endif
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
endif
call wall_time(wall1)
print*,'wall time for three_body_ints',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_ints on disk ...'
call write_array_6_index_tensor(mo_num,three_body_ints,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
subroutine give_integrals_3_body(i,j,m,k,l,n,integral)
implicit none
double precision, intent(out) :: integral
integer, intent(in) :: i,j,m,k,l,n
double precision :: weight
BEGIN_DOC
! <ijm|L|kln>
END_DOC
integer :: ipoint,mm
integral = 0.d0
do mm = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
integral += weight * mos_in_r_array_transp(ipoint,i) * mos_in_r_array_transp(ipoint,k) * x_W_ij_erf_rk(ipoint,mm,m,n) * x_W_ij_erf_rk(ipoint,mm,j,l)
integral += weight * mos_in_r_array_transp(ipoint,j) * mos_in_r_array_transp(ipoint,l) * x_W_ij_erf_rk(ipoint,mm,m,n) * x_W_ij_erf_rk(ipoint,mm,i,k)
integral += weight * mos_in_r_array_transp(ipoint,m) * mos_in_r_array_transp(ipoint,n) * x_W_ij_erf_rk(ipoint,mm,j,l) * x_W_ij_erf_rk(ipoint,mm,i,k)
enddo
enddo
end

338
src/three_body_ints/three_e_3_idx.irp.f Normal file
View File

@ -0,0 +1,338 @@
BEGIN_PROVIDER [ double precision, three_body_3_index, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 3 index matrix element of the -L three-body operator
!
! three_body_3_index(k,l,n) = < phi_k phi_l phi_n | phi_k phi_l phi_n >
!
! notice the -1 sign: in this way three_body_3_index can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
print*,'Providing the three_body_3_index ...'
name_file = 'three_body_3_index'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_ints from disk ...'
call read_array_3_index_tensor(mo_num,three_body_3_index,name_file)
else
provide x_W_ij_erf_rk
three_body_3_index = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_body_3_index)
!$OMP DO SCHEDULE (guided) COLLAPSE(3)
do m = 1, mo_num ! 3
do j = 1, mo_num ! 2
do i = 1, mo_num ! 1
integral = 0.d0
! 1 2 3 1 2 3
call give_integrals_3_body(i,j,m,i,j,m,integral)
three_body_3_index(i,j,m) = -1.d0 * integral
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_3_index',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_3_index on disk ...'
call write_array_3_index_tensor(mo_num,three_body_3_index,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_3_index_exch_12, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 3 index matrix EXCHANGE element of the -L three-body operator
!
! three_body_3_index_exch_12(k,l,n) = < phi_k phi_l phi_n | phi_l phi_k phi_n >
!
! notice the -1 sign: in this way three_body_3_index_exch_12 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
name_file = 'three_body_3_index_exch_12'
print*,'Providing the three_body_3_index_exch_12 ...'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_ints from disk ...'
call read_array_3_index_tensor(mo_num,three_body_3_index_exch_12,name_file)
else
provide x_W_ij_erf_rk
three_body_3_index_exch_12 = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_body_3_index_exch_12)
!$OMP DO SCHEDULE (guided) COLLAPSE(3)
do m = 1, mo_num ! 3
do j = 1, mo_num ! 2
do i = 1, mo_num ! 1
integral = 0.d0
! 1 2 3 1 2 3
call give_integrals_3_body(i,j,m,j,i,m,integral)
three_body_3_index_exch_12(i,j,m) = -1.d0 * integral
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_3_index_exch_12',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_3_index_exch_12 on disk ...'
call write_array_3_index_tensor(mo_num,three_body_3_index_exch_12,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_3_index_exch_23, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 3 index matrix EXCHANGE element of the -L three-body operator
!
! three_body_3_index_exch_12(k,l,n) = < phi_k phi_l phi_n | phi_k phi_n phi_l >
!
! notice the -1 sign: in this way three_body_3_index_exch_12 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
print*,'Providing the three_body_3_index_exch_23 ...'
call wall_time(wall0)
name_file = 'three_body_3_index_exch_23'
if(read_three_body_ints)then
print*,'Reading three_body_ints from disk ...'
call read_array_3_index_tensor(mo_num,three_body_3_index_exch_23,name_file)
else
provide x_W_ij_erf_rk
three_body_3_index_exch_23 = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_body_3_index_exch_23)
!$OMP DO SCHEDULE (guided) COLLAPSE(3)
do m = 1, mo_num ! 3
do j = 1, mo_num ! 2
do i = 1, mo_num ! 1
integral = 0.d0
! 1 2 3 1 2 3
call give_integrals_3_body(i,j,m,i,m,j,integral)
three_body_3_index_exch_23(i,j,m) = -1.d0 * integral
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
endif
print*,'wall time for three_body_3_index_exch_23',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_3_index_exch_23 on disk ...'
call write_array_3_index_tensor(mo_num,three_body_3_index_exch_23,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_3_index_exch_13, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 3 index matrix EXCHANGE element of the -L three-body operator
!
! three_body_3_index_exch_12(k,l,n) = < phi_k phi_l phi_n | phi_k phi_n phi_l >
!
! notice the -1 sign: in this way three_body_3_index_exch_12 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
print*,'Providing the three_body_3_index_exch_13 ...'
call wall_time(wall0)
name_file = 'three_body_3_index_exch_13'
if(read_three_body_ints)then
print*,'Reading three_body_ints from disk ...'
call read_array_3_index_tensor(mo_num,three_body_3_index_exch_13,name_file)
else
provide x_W_ij_erf_rk
three_body_3_index_exch_13 = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_body_3_index_exch_13)
!$OMP DO SCHEDULE (guided)
do m = 1, mo_num ! 3
do j = 1, mo_num ! 2
do i = 1, mo_num ! 1
integral = 0.d0
! 1 2 3 1 2 3
call give_integrals_3_body(i,j,m,m,j,i,integral)
three_body_3_index_exch_13(i,j,m) = -1.d0 * integral
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_3_index_exch_13',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_3_index_exch_13 on disk ...'
call write_array_3_index_tensor(mo_num,three_body_3_index_exch_13,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_3_index_exch_231, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 3 index matrix element of the -L three-body operator
!
! three_body_3_index_exch_231(k,l,n) = < phi_k phi_l phi_n | phi_l phi_n phi_k >
!
! notice the -1 sign: in this way three_body_3_index can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
print*,'Providing the three_body_3_index_231 ...'
call wall_time(wall0)
name_file = 'three_body_3_index_exch_231'
if(read_three_body_ints)then
print*,'Reading three_body_ints from disk ...'
call read_array_3_index_tensor(mo_num,three_body_3_index_exch_231,name_file)
else
provide x_W_ij_erf_rk
three_body_3_index_exch_231 = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_body_3_index_exch_231)
!$OMP DO SCHEDULE (guided) COLLAPSE(3)
do m = 1, mo_num ! 3
do j = 1, mo_num ! 2
do i = 1, mo_num ! 1
integral = 0.d0
! 1 2 3 1 2 3
call give_integrals_3_body(i,j,m,j,m,i,integral)
three_body_3_index_exch_231(i,j,m) = -1.d0 * integral
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_3_index_exch_231 ',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_3_index_exch_231 on disk ...'
call write_array_3_index_tensor(mo_num,three_body_3_index_exch_231,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_3_index_exch_312, (mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 3 index matrix element of the -L three-body operator
!
! three_body_3_index(k,l,n) = < phi_k phi_l phi_n | phi_l phi_n phi_k >
!
! notice the -1 sign: in this way three_body_3_index can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,m
double precision :: integral, wall1, wall0
character*(128) :: name_file
print*,'Providing the three_body_3_index_312 ...'
call wall_time(wall0)
name_file = 'three_body_3_index_exch_312'
if(read_three_body_ints)then
print*,'Reading three_body_ints from disk ...'
call read_array_3_index_tensor(mo_num,three_body_3_index_exch_312,name_file)
else
provide x_W_ij_erf_rk
three_body_3_index_exch_312 = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_body_3_index_exch_312)
!$OMP DO SCHEDULE (guided) COLLAPSE(3)
do m = 1, mo_num ! 3
do j = 1, mo_num ! 2
do i = 1, mo_num ! 1
integral = 0.d0
! 1 2 3 1 2 3
call give_integrals_3_body(i,j,m,m,i,j,integral)
three_body_3_index_exch_312(i,j,m) = -1.d0 * integral
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_3_index_312',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_3_index_exch_312 on disk ...'
call write_array_3_index_tensor(mo_num,three_body_3_index_exch_312,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
subroutine write_array_3_index_tensor(n_orb,array_tmp,name_file)
implicit none
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'W')
write(i_unit_output)array_tmp
close(unit=i_unit_output)
end
subroutine read_array_3_index_tensor(n_orb,array_tmp,name_file)
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb)
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'R')
read(i_unit_output)array_tmp
close(unit=i_unit_output)
end

347
src/three_body_ints/three_e_4_idx.irp.f Normal file
View File

@ -0,0 +1,347 @@
BEGIN_PROVIDER [ double precision, three_body_4_index, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 4 index matrix direct element of the -L three-body operator
!
! three_body_4_index(j,m,k,i) = < phi_j phi_m phi_k | phi_j phi_m phi_i >
!
! notice the -1 sign: in this way three_body_4_index can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_4_index = 0.d0
print*,'Providing the three_body_4_index ...'
call wall_time(wall0)
name_file = 'three_body_4_index'
if(read_three_body_ints)then
print*,'Reading three_body_4_index from disk ...'
call read_array_4_index_tensor(mo_num,three_body_4_index,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,k,integral) &
!$OMP SHARED (mo_num,three_body_4_index)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
call give_integrals_3_body(i,j,m,k,j,m,integral)
three_body_4_index(j,m,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_4_index',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_4_index on disk ...'
call write_array_4_index_tensor(mo_num,three_body_4_index,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_4_index_exch_12, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 4 index matrix EXCHANGE element of the -L three-body operator
!
! three_body_4_index_exch_12(j,m,k,i) = < phi_m phi_j phi_i | phi_j phi_m phi_k >
!
! notice the -1 sign: in this way three_body_3_index_exch_12 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_4_index_exch_12 = 0.d0
print*,'Providing the three_body_4_index_exch_12 ...'
call wall_time(wall0)
name_file = 'three_body_4_index_exch_12'
if(read_three_body_ints)then
print*,'Reading three_body_4_index_exch_12 from disk ...'
call read_array_4_index_tensor(mo_num,three_body_4_index_exch_12,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,k,integral) &
!$OMP SHARED (mo_num,three_body_4_index_exch_12)
!$OMP DO SCHEDULE (guided) COLLAPSE(4)
do i = 1, mo_num
do k = 1, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
call give_integrals_3_body(i,m,j,k,j,m,integral)
three_body_4_index_exch_12(j,m,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_4_index_exch_12',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_4_index_exch_12 on disk ...'
call write_array_4_index_tensor(mo_num,three_body_4_index_exch_12,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_4_index_exch_12_part, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 4 index matrix EXCHANGE element of the -L three-body operator
!
! three_body_4_index_exch_12_part(j,m,k,i) = < phi_m phi_j phi_i | phi_m phi_k phi_j >
!
! notice the -1 sign: in this way three_body_3_index_exch_12 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_4_index_exch_12_part = 0.d0
print*,'Providing the three_body_4_index_exch_12_part ...'
call wall_time(wall0)
name_file = 'three_body_4_index_exch_12_part'
if(read_three_body_ints)then
print*,'Reading three_body_4_index_exch_12_part from disk ...'
call read_array_4_index_tensor(mo_num,three_body_4_index_exch_12_part,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,k,integral) &
!$OMP SHARED (mo_num,three_body_4_index_exch_12_part)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
!
call give_integrals_3_body(i,j,m,j,k,m,integral)
three_body_4_index_exch_12_part(j,m,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
endif
print*,'wall time for three_body_4_index_exch_12_part',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_4_index_exch_12_part on disk ...'
call write_array_4_index_tensor(mo_num,three_body_4_index_exch_12_part,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_4_index_exch_12_part_bis, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 4 index matrix EXCHANGE element of the -L three-body operator
!
! three_body_4_index_exch_12_part_bis(j,m,k,i) = < phi_m phi_j phi_i | phi_m phi_k phi_j >
!
! notice the -1 sign: in this way three_body_3_index_exch_12 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_4_index_exch_12_part_bis = 0.d0
print*,'Providing the three_body_4_index_exch_12_part_bis ...'
call wall_time(wall0)
name_file = 'three_body_4_index_exch_12_part_bis'
if(read_three_body_ints)then
print*,'Reading three_body_4_index_exch_12_part_bisfrom disk ...'
call read_array_4_index_tensor(mo_num,three_body_4_index_exch_12_part_bis,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,k,integral) &
!$OMP SHARED (mo_num,three_body_4_index_exch_12_part_bis)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
!
call give_integrals_3_body(i,j,m,m,j,k,integral)
three_body_4_index_exch_12_part_bis(j,m,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_4_index_exch_12_part_bis',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_4_index_exch_12_part_bis on disk ...'
call write_array_4_index_tensor(mo_num,three_body_4_index_exch_12_part_bis,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_4_index_exch_231, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 4 index matrix direct element of the -L three-body operator
!
! three_body_4_index_exch_231(j,m,k,i) = < phi_j phi_m phi_k | phi_j phi_m phi_i >
!
! notice the -1 sign: in this way three_body_4_index_exch_231 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_4_index_exch_231 = 0.d0
print*,'Providing the three_body_4_index_exch_231 ...'
call wall_time(wall0)
name_file = 'three_body_4_index_exch_231'
if(read_three_body_ints)then
print*,'Reading three_body_4_index_exch_231 from disk ...'
call read_array_4_index_tensor(mo_num,three_body_4_index_exch_231,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,k,integral) &
!$OMP SHARED (mo_num,three_body_4_index_exch_231)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
call give_integrals_3_body(i,j,m,j,m,k,integral)
three_body_4_index_exch_231(j,m,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_4_index_exch_231',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_4_index_exch_231 on disk ...'
call write_array_4_index_tensor(mo_num,three_body_4_index_exch_231,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_4_index_exch_312, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 4 index matrix direct element of the -L three-body operator
!
! three_body_4_index_exch_312(j,m,k,i) = < phi_j phi_m phi_k | phi_j phi_m phi_i >
!
! notice the -1 sign: in this way three_body_4_index_exch_312 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_4_index_exch_312 = 0.d0
print*,'Providing the three_body_4_index_exch_312 ...'
call wall_time(wall0)
name_file = 'three_body_4_index_exch_312'
if(read_three_body_ints)then
print*,'Reading three_body_4_index_exch_312 from disk ...'
call read_array_4_index_tensor(mo_num,three_body_4_index_exch_312,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,k,integral) &
!$OMP SHARED (mo_num,three_body_4_index_exch_312)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
call give_integrals_3_body(i,j,m,m,k,j,integral)
three_body_4_index_exch_312(j,m,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_4_index_exch_312',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_4_index_exch_312 on disk ...'
call write_array_4_index_tensor(mo_num,three_body_4_index_exch_312,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
subroutine write_array_4_index_tensor(n_orb,array_tmp,name_file)
implicit none
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'W')
write(i_unit_output)array_tmp
close(unit=i_unit_output)
end
subroutine read_array_4_index_tensor(n_orb,array_tmp,name_file)
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb)
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'R')
read(i_unit_output)array_tmp
close(unit=i_unit_output)
end

453
src/three_body_ints/three_e_5_idx.irp.f Normal file
View File

@ -0,0 +1,453 @@
BEGIN_PROVIDER [ double precision, three_body_5_index, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 5 index matrix element of the -L three-body operator
!
! three_body_5_index(i,j,m,l,n) = < phi_i phi_j phi_m | phi_i phi_l phi_n >
!
! notice the -1 sign: in this way three_body_5_index can be directly used to compute Slater rules :)
END_DOC
integer :: j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_5_index(1:mo_num, 1:mo_num, 1:mo_num, 1:mo_num, 1:mo_num) = 0.d0
print*,'Providing the three_body_5_index ...'
name_file = 'three_body_5_index'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_5_index from disk ...'
call read_array_5_index_tensor(mo_num,three_body_5_index,name_file)
else
provide x_W_ij_erf_rk
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_5_index)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do n = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
call give_integrals_3_body(j,m,k,l,n,k,integral)
three_body_5_index(k,j,m,l,n) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_5_index',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_5_index on disk ...'
call write_array_5_index_tensor(mo_num,three_body_5_index,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
! do n = 1, mo_num
! do l = 1, mo_num
! do k = 1, mo_num
! do m = 1, n-1
! do j = 1, l-1
! three_body_5_index(k,j,m,l,n) = three_body_5_index(k,l,n,j,m)
! three_body_5_index(k,j,m,l,n)
! enddo
! enddo
! enddo
! enddo
! enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_5_index_exch_13, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 5 index matrix element of the -L three-body operator
!
! three_body_5_index_exch_13(k,j,m,l,n) = < phi_j phi_m phi_k | phi_k phi_n phi_l >
!
! notice the -1 sign: in this way three_body_5_index_exch_13 can be directly used to compute Slater rules :)
END_DOC
integer :: j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_5_index_exch_13 = 0.d0
name_file = 'three_body_5_index_exch_13'
print*,'Providing the three_body_5_index_exch_13 ...'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_5_index_exch_13 from disk ...'
call read_array_5_index_tensor(mo_num,three_body_5_index_exch_13,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_5_index_exch_13)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do n = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
!! j,m,k,l,n,k : direct (case 2)
call give_integrals_3_body(j,m,k,k,n,l,integral)
!! j,m,k,k,n,l : exchange 1 3
three_body_5_index_exch_13(k,j,m,l,n) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_5_index_exch_13',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_5_index_exch_13 on disk ...'
call write_array_5_index_tensor(mo_num,three_body_5_index_exch_13,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
! do n = 1, mo_num
! do l = 1, mo_num
! do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
! three_body_5_index_exch_13(k,l,n,j,m) = three_body_5_index_exch_13(k,j,m,l,n)
! enddo
! enddo
! enddo
! enddo
! enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_5_index_exch_32, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 5 index matrix element of the -L three-body operator
!
! three_body_5_index_exch_32(i,j,m,l,n) = < phi_i phi_j phi_m | phi_i phi_l phi_n >
!
! notice the -1 sign: in this way three_body_5_index_exch_32 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(328) :: name_file
three_body_5_index_exch_32 = 0.d0
name_file = 'three_body_5_index_exch_32'
print*,'Providing the three_body_5_index_exch_32 ...'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_5_index_exch_32 from disk ...'
call read_array_5_index_tensor(mo_num,three_body_5_index_exch_32,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_5_index_exch_32)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do n = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
!! j,m,k,l,n,k : direct (case 3)
call give_integrals_3_body(j,m,k,l,k,n,integral)
!! j,m,k,l,k,n : exchange 2 3
three_body_5_index_exch_32(k,j,m,l,n) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_5_index_exch_32',wall1 - wall0
if(write_three_body_ints)then
print*,'Writing three_body_5_index_exch_32 on disk ...'
call write_array_5_index_tensor(mo_num,three_body_5_index_exch_32,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
! do n = 1, mo_num
! do l = 1, mo_num
! do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
! three_body_5_index_exch_32(k,l,n,j,m) = three_body_5_index_exch_32(k,j,m,l,n)
! enddo
! enddo
! enddo
! enddo
! enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_5_index_exch_12, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 5 index matrix element of the -L three-body operator
!
! three_body_5_index_exch_12(i,j,m,l,n) = < phi_i phi_j phi_m | phi_i phi_l phi_n >
!
! notice the -1 sign: in this way three_body_5_index_exch_12 can be directly used to compute Slater rules :)
END_DOC
integer :: i,j,k,l,m,n
double precision :: integral, wall1, wall0
character*(328) :: name_file
three_body_5_index_exch_12 = 0.d0
name_file = 'three_body_5_index_exch_12'
print*,'Providing the three_body_5_index_exch_12 ...'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_5_index_exch_12 from disk ...'
call read_array_5_index_tensor(mo_num,three_body_5_index_exch_12,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_5_index_exch_12)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do n = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
!! j,m,k,l,n,k : direct (case 1)
call give_integrals_3_body(j,m,k,n,l,k,integral)
!! j,m,k,l,k,n : exchange 2 3
three_body_5_index_exch_12(k,j,m,l,n) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_5_index_exch_12',wall1 - wall0
! do n = 1, mo_num
! do l = 1, mo_num
! do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
! three_body_5_index_exch_12(k,l,n,j,m) = three_body_5_index_exch_12(k,j,m,l,n)
! enddo
! enddo
! enddo
! enddo
! enddo
if(write_three_body_ints)then
print*,'Writing three_body_5_index_exch_12 on disk ...'
call write_array_5_index_tensor(mo_num,three_body_5_index_exch_12,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_5_index_312, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 5 index matrix element of the -L three-body operator
!
! three_body_5_index_312(i,j,m,l,n) = < phi_i phi_j phi_m | phi_i phi_l phi_n >
!
! notice the -1 sign: in this way three_body_5_index_312 can be directly used to compute Slater rules :)
END_DOC
integer :: j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_5_index_312 = 0.d0
name_file = 'three_body_5_index_312'
print*,'Providing the three_body_5_index_312 ...'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_5_index_312 from disk ...'
call read_array_5_index_tensor(mo_num,three_body_5_index_312,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_5_index_312)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do n = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
! <j m k|l n k> - > <j m k|n k l>
call give_integrals_3_body(j,m,k,n,k,l,integral)
three_body_5_index_312(k,j,m,l,n) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_5_index_312',wall1 - wall0
! do n = 1, mo_num
! do l = 1, mo_num
! do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
! three_body_5_index_312(k,l,n,j,m) = three_body_5_index_312(k,j,m,l,n)
! enddo
! enddo
! enddo
! enddo
! enddo
if(write_three_body_ints)then
print*,'Writing three_body_5_index_312 on disk ...'
call write_array_5_index_tensor(mo_num,three_body_5_index_312,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_body_5_index_132, (mo_num, mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! 5 index matrix element of the -L three-body operator
!
! three_body_5_index_132(i,j,m,l,n) = < phi_i phi_j phi_m | phi_i phi_l phi_n >
!
! notice the -1 sign: in this way three_body_5_index_132 can be directly used to compute Slater rules :)
END_DOC
integer :: j,k,l,m,n
double precision :: integral, wall1, wall0
character*(128) :: name_file
three_body_5_index_132 = 0.d0
name_file = 'three_body_5_index_132'
print*,'Providing the three_body_5_index_132 ...'
call wall_time(wall0)
if(read_three_body_ints)then
print*,'Reading three_body_5_index_132 from disk ...'
call read_array_5_index_tensor(mo_num,three_body_5_index_132,name_file)
else
provide x_W_ij_erf_rk
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j,k,l,m,n,integral) &
!$OMP SHARED (mo_num,three_body_5_index_132)
!$OMP DO SCHEDULE (guided) COLLAPSE(2)
do n = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
do m = 1, mo_num
do j = 1, mo_num
integral = 0.d0
! <j m k|l n k> - > <j m k|k l n>
call give_integrals_3_body(j,m,k,k,l,n,integral)
three_body_5_index_132(k,j,m,l,n) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call wall_time(wall1)
print*,'wall time for three_body_5_index_132',wall1 - wall0
! do n = 1, mo_num
! do l = 1, mo_num
! do k = 1, mo_num
! do m = n, mo_num
! do j = l, mo_num
! three_body_5_index_132(k,l,n,j,m) = three_body_5_index_132(k,j,m,l,n)
! enddo
! enddo
! enddo
! enddo
! enddo
if(write_three_body_ints)then
print*,'Writing three_body_5_index_132 on disk ...'
call write_array_5_index_tensor(mo_num,three_body_5_index_132,name_file)
call ezfio_set_three_body_ints_io_three_body_ints("Read")
endif
END_PROVIDER
subroutine write_array_5_index_tensor(n_orb,array_tmp,name_file)
implicit none
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,n_orb)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'W')
write(i_unit_output)array_tmp
close(unit=i_unit_output)
end
subroutine read_array_5_index_tensor(n_orb,array_tmp,name_file)
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer, intent(in) :: n_orb
character*(128), intent(in) :: name_file
double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,n_orb)
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'R')
read(i_unit_output)array_tmp
close(unit=i_unit_output)
end