LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-09-08 08:54:44 +02:00
Code Releases Activity

Merge pull request #305 from AbdAmmar/dev-stable-tc-scf

Browse Source 
Dev stable tc scf
This commit is contained in:
AbdAmmar 2023-09-19 14:15:38 +02:00 committed by GitHub
commit a3e8b4732a
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
55 changed files with 8346 additions and 3344 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

154
src/ao_many_one_e_ints/ao_erf_gauss.irp.f
View File

@ -1245,3 +1245,157 @@ end subroutine NAI_pol_x2_mult_erf_ao
! ---
subroutine NAI_pol_012_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints)
BEGIN_DOC
!
! Computes the following integral :
!
! ints(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
! ints(2) = $\int_{-\infty}^{infty} dr x^1 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! ints(3) = $\int_{-\infty}^{infty} dr y^1 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! ints(4) = $\int_{-\infty}^{infty} dr z^1 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
! ints(5) = $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! ints(6) = $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! ints(7) = $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3)
double precision, intent(out) :: ints(7)
integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, m
integer :: power_A1(3), power_A2(3)
double precision :: Ai_center(3), Aj_center(3), alphai, alphaj, coef, coefi
double precision :: integral0, integral1, integral2
double precision, external :: NAI_pol_mult_erf_with1s
ASSERT(beta .ge. 0.d0)
if(beta .lt. 1d-10) then
call NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
return
endif
ints = 0.d0
power_Ai(1:3) = ao_power(i_ao,1:3)
power_Aj(1:3) = ao_power(j_ao,1:3)
Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alphai = ao_expo_ordered_transp (i,i_ao)
coefi = ao_coef_normalized_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
alphaj = ao_expo_ordered_transp (j,j_ao)
coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao)
integral0 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_Ai, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in)
ints(1) += coef * integral0
do m = 1, 3
power_A1 = power_Ai
power_A1(m) += 1
integral1 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A1, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in)
ints(1+m) += coef * (integral1 + Ai_center(m)*integral0)
power_A2 = power_Ai
power_A2(m) += 2
integral2 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A2, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in)
ints(4+m) += coef * (integral2 + Ai_center(m) * (2.d0*integral1 + Ai_center(m)*integral0))
enddo
enddo
enddo
end subroutine NAI_pol_012_mult_erf_ao_with1s
! ---
subroutine NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
BEGIN_DOC
!
! Computes the following integral :
!
! int(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
! int(2) = $\int_{-\infty}^{infty} dr x^1 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! int(3) = $\int_{-\infty}^{infty} dr y^1 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! int(4) = $\int_{-\infty}^{infty} dr z^1 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
! int(5) = $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! int(6) = $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! int(7) = $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: mu_in, C_center(3)
double precision, intent(out) :: ints(7)
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, m
integer :: power_A1(3), power_A2(3)
double precision :: A_center(3), B_center(3), alpha, beta, coef
double precision :: integral0, integral1, integral2
double precision :: NAI_pol_mult_erf
ints = 0.d0
num_A = ao_nucl(i_ao)
power_A(1:3) = ao_power(i_ao,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j_ao)
power_B(1:3) = ao_power(j_ao,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
n_pt_in = n_pt_max_integrals
do i = 1, ao_prim_num(i_ao)
alpha = ao_expo_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
beta = ao_expo_ordered_transp(j,j_ao)
coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
integral0 = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
ints(1) += coef * integral0
do m = 1, 3
power_A1 = power_A
power_A1(m) += 1
integral1 = NAI_pol_mult_erf(A_center, B_center, power_A1, power_B, alpha, beta, C_center, n_pt_in, mu_in)
ints(1+m) += coef * (integral1 + A_center(m)*integral0)
power_A2 = power_A
power_A2(m) += 2
integral2 = NAI_pol_mult_erf(A_center, B_center, power_A2, power_B, alpha, beta, C_center, n_pt_in, mu_in)
ints(4+m) += coef * (integral2 + A_center(m) * (2.d0*integral1 + A_center(m)*integral0))
enddo
enddo
enddo
end subroutine NAI_pol_012_mult_erf_ao
! ---

131
src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f
View File

@ -299,15 +299,12 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
!
! TODO
! one subroutine for all integrals
!
END_DOC
include 'constants.include.F'
@ -325,7 +322,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin
double precision, external :: overlap_gauss_r12_ao_with1s
double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_an ...'
print*, ' providing v_ij_u_cst_mu_j1b_an_old ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
@ -333,7 +330,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin
ct = inv_sq_pi_2 / mu_erf
v_ij_u_cst_mu_j1b_an = 0.d0
v_ij_u_cst_mu_j1b_an_old = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
@ -342,7 +339,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, &
!$OMP final_grid_points, mu_erf, ct, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an)
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an_old)
!$OMP DO
do ipoint = 1, n_points_final_grid
@ -413,6 +410,125 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin
! ---
v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = tmp
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = v_ij_u_cst_mu_j1b_an_old(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_an_old', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
!
END_DOC
include 'constants.include.F'
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3), r1_2
double precision :: int_o
double precision :: int_c(7), int_e(7)
double precision :: coef, beta, B_center(3)
double precision :: tmp, ct
double precision :: wall0, wall1
double precision, external :: overlap_gauss_r12_ao_with1s
double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_an ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
ct = inv_sq_pi_2 / mu_erf
v_ij_u_cst_mu_j1b_an = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP r1_2, tmp, int_c, int_e, int_o) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, &
!$OMP final_grid_points, mu_erf, ct, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
r1_2 = 0.5d0 * (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
do i = 1, ao_num
do j = i, ao_num
! ---
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e)
int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j)
tmp = coef &
* ( r1_2 * (int_c(1) - int_e(1)) &
- r(1) * (int_c(2) - int_e(2)) - r(2) * (int_c(3) - int_e(3)) - r(3) * (int_c(4) - int_e(4)) &
+ 0.5d0 * (int_c(5) + int_c(6) + int_c(7) - int_e(5) - int_e(6) - int_e(7)) &
- ct * int_o &
)
! ---
do i_1s = 2, List_all_comb_b2_size
coef = List_all_comb_b2_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e)
int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j)
tmp = tmp + coef &
* ( r1_2 * (int_c(1) - int_e(1)) &
- r(1) * (int_c(2) - int_e(2)) - r(2) * (int_c(3) - int_e(3)) - r(3) * (int_c(4) - int_e(4)) &
+ 0.5d0 * (int_c(5) + int_c(6) + int_c(7) - int_e(5) - int_e(6) - int_e(7)) &
- ct * int_o &
)
enddo
! ---
v_ij_u_cst_mu_j1b_an(j,i,ipoint) = tmp
enddo
enddo
@ -434,4 +550,3 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin
END_PROVIDER
! ---

4
src/becke_numerical_grid/EZFIO.cfg
View File

@ -33,6 +33,10 @@ doc: Number of angular grid points given from input. Warning, this number cannot
interface: ezfio,provider,ocaml
default: 1202
[n_points_extra_final_grid]
type: integer
doc: Total number of extra_grid points
interface: ezfio
[extra_grid_type_sgn]
type: integer

4
src/becke_numerical_grid/extra_grid.irp.f
View File

@ -14,7 +14,7 @@
implicit none
if(.not.my_extra_grid_becke)then
if(.not. my_extra_grid_becke) then
select case (extra_grid_type_sgn)
case(0)
n_points_extra_radial_grid = 23
@ -33,7 +33,7 @@
stop
end select
else
n_points_extra_radial_grid = my_n_pt_r_extra_grid
n_points_extra_radial_grid = my_n_pt_r_extra_grid
n_points_extra_integration_angular = my_n_pt_a_extra_grid
endif

34
src/becke_numerical_grid/extra_grid_vector.irp.f
View File

@ -23,29 +23,33 @@ BEGIN_PROVIDER [integer, n_points_extra_final_grid]
enddo
enddo
print*,'n_points_extra_final_grid = ',n_points_extra_final_grid
print*,'n max point = ',n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1)
! call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid)
print*, ' n_points_extra_final_grid = ', n_points_extra_final_grid
print*, ' n max point = ', n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1)
call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid)
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, final_grid_points_extra, (3,n_points_extra_final_grid)]
&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_extra, (n_points_extra_final_grid) ]
&BEGIN_PROVIDER [integer, index_final_points_extra, (3,n_points_extra_final_grid) ]
&BEGIN_PROVIDER [integer, index_final_points_extra_reverse, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ]
implicit none
&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_extra, (n_points_extra_final_grid)]
&BEGIN_PROVIDER [integer, index_final_points_extra, (3,n_points_extra_final_grid)]
&BEGIN_PROVIDER [integer, index_final_points_extra_reverse, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num)]
BEGIN_DOC
! final_grid_points_extra(1:3,j) = (/ x, y, z /) of the jth grid point
!
! final_weight_at_r_vector_extra(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
!
! index_final_points_extra(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
!
! index_final_points_extra_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
! final_grid_points_extra(1:3,j) = (/ x, y, z /) of the jth grid point
!
! final_weight_at_r_vector_extra(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
!
! index_final_points_extra(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
!
! index_final_points_extra_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
END_DOC
implicit none
integer :: i,j,k,l,i_count
double precision :: r(3)
i_count = 0
do j = 1, nucl_num
do i = 1, n_points_extra_radial_grid -1
@ -67,3 +71,5 @@ END_PROVIDER
enddo
END_PROVIDER

5
src/becke_numerical_grid/grid_becke.irp.f
View File

@ -14,7 +14,7 @@
implicit none
if(.not.my_grid_becke)then
if(.not. my_grid_becke) then
select case (grid_type_sgn)
case(0)
n_points_radial_grid = 23
@ -37,6 +37,9 @@
n_points_integration_angular = my_n_pt_a_grid
endif
print*, " n_points_radial_grid = ", n_points_radial_grid
print*, " n_points_integration_angular = ", n_points_integration_angular
END_PROVIDER
! ---

35
src/bi_ort_ints/bi_ort_ints.irp.f
View File

@ -18,10 +18,11 @@ program bi_ort_ints
! call test_5idx
! call test_5idx2
call test_4idx()
call test_4idx_n4()
!call test_4idx_n4()
!call test_4idx2()
!call test_5idx2
!call test_5idx
end
subroutine test_5idx2
@ -340,7 +341,7 @@ subroutine test_4idx()
implicit none
integer :: i, j, k, l
double precision :: accu, contrib, new, ref, thr
double precision :: accu, contrib, new, ref, thr, norm
thr = 1d-10
@ -348,6 +349,7 @@ subroutine test_4idx()
PROVIDE three_e_4_idx_direct_bi_ort
accu = 0.d0
norm = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
@ -356,7 +358,6 @@ subroutine test_4idx()
new = three_e_4_idx_direct_bi_ort (l,k,j,i)
ref = three_e_4_idx_direct_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_direct_bi_ort'
print*, l, k, j, i
@ -364,11 +365,14 @@ subroutine test_4idx()
stop
endif
accu += contrib
norm += dabs(ref)
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_direct_bi_ort = ', accu / dble(mo_num)**4
print*, ' accu on three_e_4_idx_direct_bi_ort (%) = ', 100.d0 * accu / norm
! ---
@ -376,6 +380,7 @@ subroutine test_4idx()
PROVIDE three_e_4_idx_exch13_bi_ort
accu = 0.d0
norm = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
@ -384,7 +389,6 @@ subroutine test_4idx()
new = three_e_4_idx_exch13_bi_ort (l,k,j,i)
ref = three_e_4_idx_exch13_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_exch13_bi_ort'
print*, l, k, j, i
@ -392,11 +396,14 @@ subroutine test_4idx()
stop
endif
accu += contrib
norm += dabs(ref)
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_exch13_bi_ort = ', accu / dble(mo_num)**4
print*, ' accu on three_e_4_idx_exch13_bi_ort (%) = ', 100.d0 * accu / norm
! ---
@ -404,6 +411,7 @@ subroutine test_4idx()
PROVIDE three_e_4_idx_cycle_1_bi_ort
accu = 0.d0
norm = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
@ -412,7 +420,6 @@ subroutine test_4idx()
new = three_e_4_idx_cycle_1_bi_ort (l,k,j,i)
ref = three_e_4_idx_cycle_1_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_cycle_1_bi_ort'
print*, l, k, j, i
@ -420,11 +427,14 @@ subroutine test_4idx()
stop
endif
accu += contrib
norm += dabs(ref)
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_cycle_1_bi_ort = ', accu / dble(mo_num)**4
print*, ' accu on three_e_4_idx_cycle_1_bi_ort (%) = ', 100.d0 * accu / norm
! ---
@ -432,6 +442,7 @@ subroutine test_4idx()
PROVIDE three_e_4_idx_exch23_bi_ort
accu = 0.d0
norm = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
@ -440,7 +451,6 @@ subroutine test_4idx()
new = three_e_4_idx_exch23_bi_ort (l,k,j,i)
ref = three_e_4_idx_exch23_bi_ort_old(l,k,j,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. thr) then
print*, ' problem in three_e_4_idx_exch23_bi_ort'
print*, l, k, j, i
@ -448,13 +458,18 @@ subroutine test_4idx()
stop
endif
accu += contrib
norm += dabs(ref)
enddo
enddo
enddo
enddo
print*, ' accu on three_e_4_idx_exch23_bi_ort = ', accu / dble(mo_num)**4
print*, ' accu on three_e_4_idx_exch23_bi_ort (%) = ', 100.d0 * accu / norm
! ---
return
end

1610
src/bi_ort_ints/no_dressing.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

66
src/bi_ort_ints/no_dressing_energy.irp.f Normal file
View File

@ -0,0 +1,66 @@
! ---
BEGIN_PROVIDER [double precision, energy_1e_noL_HF]
implicit none
integer :: i
PROVIDE mo_bi_ortho_tc_one_e
energy_1e_noL_HF = 0.d0
do i = 1, elec_beta_num
energy_1e_noL_HF += mo_bi_ortho_tc_one_e(i,i)
enddo
do i = 1, elec_alpha_num
energy_1e_noL_HF += mo_bi_ortho_tc_one_e(i,i)
enddo
print*, "energy_1e_noL_HF = ", energy_1e_noL_HF
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, energy_2e_noL_HF]
implicit none
integer :: i, j
PROVIDE mo_bi_ortho_tc_two_e
energy_2e_noL_HF = 0.d0
! down-down & down-down
do i = 1, elec_beta_num
do j = 1, elec_beta_num
energy_2e_noL_HF += (mo_bi_ortho_tc_two_e(i,j,i,j) - mo_bi_ortho_tc_two_e(j,i,i,j))
enddo
enddo
! down-down & up-up
do i = 1, elec_beta_num
do j = 1, elec_alpha_num
energy_2e_noL_HF += mo_bi_ortho_tc_two_e(i,j,i,j)
enddo
enddo
! up-up & down-down
do i = 1, elec_alpha_num
do j = 1, elec_beta_num
energy_2e_noL_HF += mo_bi_ortho_tc_two_e(i,j,i,j)
enddo
enddo
! up-up & up-up
do i = 1, elec_alpha_num
do j = 1, elec_alpha_num
energy_2e_noL_HF += (mo_bi_ortho_tc_two_e(i,j,i,j) - mo_bi_ortho_tc_two_e(j,i,i,j))
enddo
enddo
! 0.5 x is in the Slater-Condon rules and not in the integrals
energy_2e_noL_HF = 0.5d0 * energy_2e_noL_HF
print*, "energy_2e_noL_HF = ", energy_2e_noL_HF
END_PROVIDER
! ---

512
src/bi_ort_ints/no_dressing_naive.irp.f Normal file
View File

@ -0,0 +1,512 @@
! ---
BEGIN_PROVIDER [double precision, noL_0e_naive]
implicit none
integer :: ii, jj, kk
integer :: i, j, k
double precision :: sigma_i, sigma_j, sigma_k
double precision :: I_ijk_ijk, I_ijk_kij, I_ijk_jki, I_ijk_jik, I_ijk_kji, I_ijk_ikj
double precision :: t0, t1
double precision, allocatable :: tmp(:)
print*, " Providing noL_0e_naive ..."
call wall_time(t0)
allocate(tmp(elec_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ii, i, sigma_i, jj, j, sigma_j, kk, k, sigma_k, &
!$OMP I_ijk_ijk, I_ijk_kij, I_ijk_jki, I_ijk_jik, &
!$OMP I_ijk_kji, I_ijk_ikj) &
!$OMP SHARED (elec_beta_num, elec_num, tmp)
!$OMP DO
do ii = 1, elec_num
if(ii .le. elec_beta_num) then
i = ii
sigma_i = -1.d0
else
i = ii - elec_beta_num
sigma_i = +1.d0
endif
tmp(ii) = 0.d0
do jj = 1, elec_num
if(jj .le. elec_beta_num) then
j = jj
sigma_j = -1.d0
else
j = jj - elec_beta_num
sigma_j = +1.d0
endif
do kk = 1, elec_num
if(kk .le. elec_beta_num) then
k = kk
sigma_k = -1.d0
else
k = kk - elec_beta_num
sigma_k = +1.d0
endif
call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k &
, i, sigma_i, j, sigma_j, k, sigma_k &
, I_ijk_ijk)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k &
, k, sigma_k, i, sigma_i, j, sigma_j &
, I_ijk_kij)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k &
, j, sigma_j, k, sigma_k, i, sigma_i &
, I_ijk_jki)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k &
, j, sigma_j, i, sigma_i, k, sigma_k &
, I_ijk_jik)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k &
, k, sigma_k, j, sigma_j, i, sigma_i &
, I_ijk_kji)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k &
, i, sigma_i, k, sigma_k, j, sigma_j &
, I_ijk_ikj)
tmp(ii) = tmp(ii) + I_ijk_ijk + I_ijk_kij + I_ijk_jki - I_ijk_jik - I_ijk_kji - I_ijk_ikj
! = tmp(ii) + I_ijk_ijk + 2.d0 * I_ijk_kij - 3.d0 * I_ijk_jik
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
noL_0e_naive = -1.d0 * (sum(tmp)) / 6.d0
deallocate(tmp)
call wall_time(t1)
print*, " Wall time for noL_0e_naive (min) = ", (t1 - t0)/60.d0
print*, " noL_0e_naive = ", noL_0e_naive
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, noL_1e_naive, (mo_num, mo_num)]
BEGIN_DOC
!
! < p | H(1) | s > is dressed with noL_1e_naive(p,s)
!
END_DOC
implicit none
integer :: ii, jj
integer :: i, j, p, s
double precision :: sigma_i, sigma_j, sigma_p, sigma_s
double precision :: I_pij_sji, I_pij_sij, I_pij_jis, I_pij_ijs, I_pij_isj, I_pij_jsi
double precision :: t0, t1
print*, " Providing noL_1e_naive ..."
call wall_time(t0)
! ----
! up-up part
sigma_p = +1.d0
sigma_s = +1.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ii, i, sigma_i, jj, j, sigma_j, &
!$OMP I_pij_sji, I_pij_sij, I_pij_jis, &
!$OMP I_pij_ijs, I_pij_isj, I_pij_jsi ) &
!$OMP SHARED (mo_num, elec_beta_num, elec_num, &
!$OMP sigma_p, sigma_s, noL_1e_naive)
!$OMP DO COLLAPSE (2)
do s = 1, mo_num
do p = 1, mo_num
noL_1e_naive(p,s) = 0.d0
do ii = 1, elec_num
if(ii .le. elec_beta_num) then
i = ii
sigma_i = -1.d0
else
i = ii - elec_beta_num
sigma_i = +1.d0
endif
do jj = 1, elec_num
if(jj .le. elec_beta_num) then
j = jj
sigma_j = -1.d0
else
j = jj - elec_beta_num
sigma_j = +1d0
endif
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, s, sigma_s, j, sigma_j, i, sigma_i &
, I_pij_sji)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, s, sigma_s, i, sigma_i, j, sigma_j &
, I_pij_sij)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, j, sigma_j, i, sigma_i, s, sigma_s &
, I_pij_jis)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, i, sigma_i, j, sigma_j, s, sigma_s &
, I_pij_ijs)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, i, sigma_i, s, sigma_s, j, sigma_j &
, I_pij_isj)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, j, sigma_j, s, sigma_s, i, sigma_i &
, I_pij_jsi)
! x 0.5 because we consider 0.5 (up + down)
noL_1e_naive(p,s) = noL_1e_naive(p,s) - 0.25d0 * (I_pij_sji - I_pij_sij + I_pij_jis - I_pij_ijs + I_pij_isj - I_pij_jsi)
enddo ! j
enddo ! i
enddo ! s
enddo ! p
!$OMP END DO
!$OMP END PARALLEL
! ----
! down-down part
sigma_p = -1.d0
sigma_s = -1.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ii, i, sigma_i, jj, j, sigma_j, &
!$OMP I_pij_sji, I_pij_sij, I_pij_jis, &
!$OMP I_pij_ijs, I_pij_isj, I_pij_jsi ) &
!$OMP SHARED (mo_num, elec_beta_num, elec_num, &
!$OMP sigma_p, sigma_s, noL_1e_naive)
!$OMP DO COLLAPSE (2)
do s = 1, mo_num
do p = 1, mo_num
do ii = 1, elec_num
if(ii .le. elec_beta_num) then
i = ii
sigma_i = -1.d0
else
i = ii - elec_beta_num
sigma_i = +1.d0
endif
do jj = 1, elec_num
if(jj .le. elec_beta_num) then
j = jj
sigma_j = -1.d0
else
j = jj - elec_beta_num
sigma_j = +1d0
endif
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, s, sigma_s, j, sigma_j, i, sigma_i &
, I_pij_sji)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, s, sigma_s, i, sigma_i, j, sigma_j &
, I_pij_sij)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, j, sigma_j, i, sigma_i, s, sigma_s &
, I_pij_jis)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, i, sigma_i, j, sigma_j, s, sigma_s &
, I_pij_ijs)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, i, sigma_i, s, sigma_s, j, sigma_j &
, I_pij_isj)
call give_integrals_3_body_bi_ort_spin( p, sigma_p, i, sigma_i, j, sigma_j &
, j, sigma_j, s, sigma_s, i, sigma_i &
, I_pij_jsi)
! x 0.5 because we consider 0.5 (up + down)
noL_1e_naive(p,s) = noL_1e_naive(p,s) - 0.25d0 * (I_pij_sji - I_pij_sij + I_pij_jis - I_pij_ijs + I_pij_isj - I_pij_jsi)
enddo ! j
enddo ! i
enddo ! s
enddo ! p
!$OMP END DO
!$OMP END PARALLEL
! ---
call wall_time(t1)
print*, " Wall time for noL_1e_naive (min) = ", (t1 - t0)/60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, noL_2e_naive, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! < p q | H(2) | s t > is dressed with noL_2e_naive(p,q,s,t)
!
END_DOC
implicit none
integer :: ii
integer :: i, p, q, s, t
double precision :: sigma_i, sigma_p, sigma_q, sigma_s, sigma_t
double precision :: I_ipq_ist, I_ipq_sit, I_ipq_tsi
double precision :: t0, t1
print*, " Providing noL_2e_naive ..."
call wall_time(t0)
! ----
! up-up & up-up part
sigma_p = +1.d0
sigma_s = +1.d0
sigma_q = +1.d0
sigma_t = +1.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ii, i, sigma_i, p, q, s, t, &
!$OMP I_ipq_ist, I_ipq_sit, I_ipq_tsi) &
!$OMP SHARED (mo_num, elec_beta_num, elec_num, &
!$OMP sigma_p, sigma_q, sigma_s, sigma_t, &
!$OMP noL_2e_naive)
!$OMP DO COLLAPSE (4)
do t = 1, mo_num
do s = 1, mo_num
do q = 1, mo_num
do p = 1, mo_num
noL_2e_naive(p,q,s,t) = 0.d0
do ii = 1, elec_num
if(ii .le. elec_beta_num) then
i = ii
sigma_i = -1.d0
else
i = ii - elec_beta_num
sigma_i = +1.d0
endif
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, i, sigma_i, s, sigma_s, t, sigma_t &
, I_ipq_ist)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, s, sigma_s, i, sigma_i, t, sigma_t &
, I_ipq_sit)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, t, sigma_t, s, sigma_s, i, sigma_i &
, I_ipq_tsi)
! x 0.25 because we consider 0.25 (up-up + up-down + down-up + down-down)
noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
enddo ! i
enddo ! p
enddo ! q
enddo ! s
enddo ! t
!$OMP END DO
!$OMP END PARALLEL
! ----
! up-up & down-down part
sigma_p = +1.d0
sigma_s = +1.d0
sigma_q = -1.d0
sigma_t = -1.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ii, i, sigma_i, p, q, s, t, &
!$OMP I_ipq_ist, I_ipq_sit, I_ipq_tsi) &
!$OMP SHARED (mo_num, elec_beta_num, elec_num, &
!$OMP sigma_p, sigma_q, sigma_s, sigma_t, &
!$OMP noL_2e_naive)
!$OMP DO COLLAPSE (4)
do t = 1, mo_num
do s = 1, mo_num
do q = 1, mo_num
do p = 1, mo_num
do ii = 1, elec_num
if(ii .le. elec_beta_num) then
i = ii
sigma_i = -1.d0
else
i = ii - elec_beta_num
sigma_i = +1.d0
endif
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, i, sigma_i, s, sigma_s, t, sigma_t &
, I_ipq_ist)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, s, sigma_s, i, sigma_i, t, sigma_t &
, I_ipq_sit)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, t, sigma_t, s, sigma_s, i, sigma_i &
, I_ipq_tsi)
! x 0.25 because we consider 0.25 (up-up + up-down + down-up + down-down)
noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
enddo ! i
enddo ! p
enddo ! q
enddo ! s
enddo ! t
!$OMP END DO
!$OMP END PARALLEL
! ----
! down-down & up-up part
sigma_p = -1.d0
sigma_s = -1.d0
sigma_q = +1.d0
sigma_t = +1.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ii, i, sigma_i, p, q, s, t, &
!$OMP I_ipq_ist, I_ipq_sit, I_ipq_tsi) &
!$OMP SHARED (mo_num, elec_beta_num, elec_num, &
!$OMP sigma_p, sigma_q, sigma_s, sigma_t, &
!$OMP noL_2e_naive)
!$OMP DO COLLAPSE (4)
do t = 1, mo_num
do s = 1, mo_num
do q = 1, mo_num
do p = 1, mo_num
do ii = 1, elec_num
if(ii .le. elec_beta_num) then
i = ii
sigma_i = -1.d0
else
i = ii - elec_beta_num
sigma_i = +1.d0
endif
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, i, sigma_i, s, sigma_s, t, sigma_t &
, I_ipq_ist)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, s, sigma_s, i, sigma_i, t, sigma_t &
, I_ipq_sit)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, t, sigma_t, s, sigma_s, i, sigma_i &
, I_ipq_tsi)
! x 0.25 because we consider 0.25 (up-up + up-down + down-up + down-down)
noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
enddo ! i
enddo ! p
enddo ! q
enddo ! s
enddo ! t
!$OMP END DO
!$OMP END PARALLEL
! ----
! down-down & down-down part
sigma_p = -1.d0
sigma_s = -1.d0
sigma_q = -1.d0
sigma_t = -1.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ii, i, sigma_i, p, q, s, t, &
!$OMP I_ipq_ist, I_ipq_sit, I_ipq_tsi) &
!$OMP SHARED (mo_num, elec_beta_num, elec_num, &
!$OMP sigma_p, sigma_q, sigma_s, sigma_t, &
!$OMP noL_2e_naive)
!$OMP DO COLLAPSE (4)
do t = 1, mo_num
do s = 1, mo_num
do q = 1, mo_num
do p = 1, mo_num
do ii = 1, elec_num
if(ii .le. elec_beta_num) then
i = ii
sigma_i = -1.d0
else
i = ii - elec_beta_num
sigma_i = +1.d0
endif
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, i, sigma_i, s, sigma_s, t, sigma_t &
, I_ipq_ist)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, s, sigma_s, i, sigma_i, t, sigma_t &
, I_ipq_sit)
call give_integrals_3_body_bi_ort_spin( i, sigma_i, p, sigma_p, q, sigma_q &
, t, sigma_t, s, sigma_s, i, sigma_i &
, I_ipq_tsi)
! x 0.25 because we consider 0.25 (up-up + up-down + down-up + down-down)
noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
enddo ! i
enddo ! p
enddo ! q
enddo ! s
enddo ! t
!$OMP END DO
!$OMP END PARALLEL
call wall_time(t1)
print*, " Wall time for noL_2e_naive (min) = ", (t1 - t0)/60.d0
END_PROVIDER
! ---

21
src/bi_ort_ints/one_e_bi_ort.irp.f
View File

@ -29,7 +29,7 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_one_e, (mo_num, mo_num)]
BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_one_e, (mo_num, mo_num)]
BEGIN_DOC
!
@ -41,6 +41,11 @@ BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_one_e, (mo_num, mo_num)]
call ao_to_mo_bi_ortho(ao_one_e_integrals_tc_tot, ao_num, mo_bi_ortho_tc_one_e, mo_num)
if(noL_standard) then
PROVIDE noL_1e
mo_bi_ortho_tc_one_e = mo_bi_ortho_tc_one_e + noL_1e
endif
END_PROVIDER
! ---
@ -48,12 +53,14 @@ END_PROVIDER
BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_x , (mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_y , (mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_z , (mo_num,mo_num)]
BEGIN_DOC
! array of the integrals of Left MO_i * x Right MO_j
! array of the integrals of Left MO_i * y Right MO_j
! array of the integrals of Left MO_i * z Right MO_j
END_DOC
implicit none
BEGIN_DOC
! array of the integrals of Left MO_i * x Right MO_j
! array of the integrals of Left MO_i * y Right MO_j
! array of the integrals of Left MO_i * z Right MO_j
END_DOC
implicit none
call ao_to_mo_bi_ortho( &
ao_dipole_x, &

183
src/bi_ort_ints/three_body_ijmk.irp.f
View File