mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-18 11:23:38 +01:00
added three_body_ints
This commit is contained in:
eginer
parent
4472a6d9be
commit
3a68b36515
20
src/three_body_ints/EZFIO.cfg
Normal file
20
src/three_body_ints/EZFIO.cfg
Normal 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
1
src/three_body_ints/NEED
Normal file
@ -0,0 +1 @@
|
|||||||
|
bi_ort_ints
|
||||||
63
src/three_body_ints/io_6_index_tensor.irp.f
Normal file
63
src/three_body_ints/io_6_index_tensor.irp.f
Normal 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
207
src/three_body_ints/semi_num_ints_mo.irp.f
Normal 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
106
src/three_body_ints/three_body_tensor.irp.f
Normal 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
338
src/three_body_ints/three_e_3_idx.irp.f
Normal 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
347
src/three_body_ints/three_e_4_idx.irp.f
Normal 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
453
src/three_body_ints/three_e_5_idx.irp.f
Normal 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
|
||||||
Loading…
Reference in New Issue
Block a user