diff --git a/src/three_body_ints/EZFIO.cfg b/src/three_body_ints/EZFIO.cfg new file mode 100644 index 00000000..9624c161 --- /dev/null +++ b/src/three_body_ints/EZFIO.cfg @@ -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 + + diff --git a/src/three_body_ints/NEED b/src/three_body_ints/NEED new file mode 100644 index 00000000..ad7b6bf8 --- /dev/null +++ b/src/three_body_ints/NEED @@ -0,0 +1 @@ +bi_ort_ints diff --git a/src/three_body_ints/io_6_index_tensor.irp.f b/src/three_body_ints/io_6_index_tensor.irp.f new file mode 100644 index 00000000..dd654f7e --- /dev/null +++ b/src/three_body_ints/io_6_index_tensor.irp.f @@ -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 + ! + ! (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 diff --git a/src/three_body_ints/semi_num_ints_mo.irp.f b/src/three_body_ints/semi_num_ints_mo.irp.f new file mode 100644 index 00000000..831ceb9b --- /dev/null +++ b/src/three_body_ints/semi_num_ints_mo.irp.f @@ -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 +! diff --git a/src/three_body_ints/three_body_tensor.irp.f b/src/three_body_ints/three_body_tensor.irp.f new file mode 100644 index 00000000..2b65a925 --- /dev/null +++ b/src/three_body_ints/three_body_tensor.irp.f @@ -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 +! + 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 + diff --git a/src/three_body_ints/three_e_3_idx.irp.f b/src/three_body_ints/three_e_3_idx.irp.f new file mode 100644 index 00000000..13210f00 --- /dev/null +++ b/src/three_body_ints/three_e_3_idx.irp.f @@ -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 diff --git a/src/three_body_ints/three_e_4_idx.irp.f b/src/three_body_ints/three_e_4_idx.irp.f new file mode 100644 index 00000000..0c6743f0 --- /dev/null +++ b/src/three_body_ints/three_e_4_idx.irp.f @@ -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 diff --git a/src/three_body_ints/three_e_5_idx.irp.f b/src/three_body_ints/three_e_5_idx.irp.f new file mode 100644 index 00000000..914601ff --- /dev/null +++ b/src/three_body_ints/three_e_5_idx.irp.f @@ -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 + + ! - > + 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 + + ! - > + 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