LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-10-07 00:25:57 +02:00
Code Releases Activity

fixed bug in AO Fock UHF

Browse Source
This commit is contained in:
AbdAmmar 2022-12-21 14:10:39 +01:00
parent 76d502bd35
commit 8eecb342c9
2 changed files with 139 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

132
src/tc_scf/fock_3e_bi_ortho_uhf.irp.f
View File

@ -244,6 +244,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)]
double precision :: dm_dk_a, dm_dk_b, dm_dk double precision :: dm_dk_a, dm_dk_b, dm_dk
double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu
double precision :: ti, tf double precision :: ti, tf
double precision, allocatable :: f_tmp(:,:)
print *, ' PROVIDING fock_3e_uhf_ao_a ...' print *, ' PROVIDING fock_3e_uhf_ao_a ...'
call wall_time(ti) call wall_time(ti)
@ -251,37 +252,37 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)]
fock_3e_uhf_ao_a = 0.d0 fock_3e_uhf_ao_a = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, & !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp, &
!$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
!$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a) !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a)
allocate(f_tmp(ao_num,ao_num))
f_tmp = 0.d0
!$OMP DO !$OMP DO
do g = 1, ao_num do g = 1, ao_num
do e = 1, ao_num do e = 1, ao_num
dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
dm_ge = dm_ge_a + dm_ge_b dm_ge = dm_ge_a + dm_ge_b
do d = 1, ao_num do d = 1, ao_num
do k = 1, ao_num do k = 1, ao_num
dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
dm_dk = dm_dk_a + dm_dk_b dm_dk = dm_dk_a + dm_dk_b
do mu = 1, ao_num do mu = 1, ao_num
do nu = 1, ao_num do nu = 1, ao_num
call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
f_tmp(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek &
fock_3e_uhf_ao_a(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek &
+ dm_ge_a * dm_dk_a * i_mugd_eknu & + dm_ge_a * dm_dk_a * i_mugd_eknu &
+ dm_ge_a * dm_dk_a * i_mugd_knue & + dm_ge_a * dm_dk_a * i_mugd_knue &
- dm_ge * dm_dk_a * i_mugd_kenu &
- dm_ge_a * dm_dk * i_mugd_enuk & - dm_ge_a * dm_dk * i_mugd_enuk &
- dm_ge * dm_dk_a * i_mugd_kenu &
- dm_ge_a * dm_dk_a * i_mugd_nuke & - dm_ge_a * dm_dk_a * i_mugd_nuke &
- dm_ge_b * dm_dk_b * i_mugd_nuke ) - dm_ge_b * dm_dk_b * i_mugd_nuke )
enddo enddo
@ -290,9 +291,59 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)]
enddo enddo
enddo enddo
enddo enddo
!$OMP END DO !$OMP END DO NOWAIT
!$OMP CRITICAL
do mu = 1, ao_num
do nu = 1, ao_num
fock_3e_uhf_ao_a(mu,nu) += f_tmp(mu,nu)
enddo
enddo
!$OMP END CRITICAL
deallocate(f_tmp)
!$OMP END PARALLEL !$OMP END PARALLEL
! TODO
! !$OMP PARALLEL DEFAULT (NONE) &
! !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, &
! !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
! !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a)
! !$OMP DO
! do g = 1, ao_num
! do e = 1, ao_num
! dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
! dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
! dm_ge = dm_ge_a + dm_ge_b
! do d = 1, ao_num
! do k = 1, ao_num
! dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
! dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
! dm_dk = dm_dk_a + dm_dk_b
! do mu = 1, ao_num
! do nu = 1, ao_num
! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
! fock_3e_uhf_ao_a(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek &
! + dm_ge_a * dm_dk_a * i_mugd_eknu &
! + dm_ge_a * dm_dk_a * i_mugd_knue &
! - dm_ge_a * dm_dk * i_mugd_enuk &
! - dm_ge * dm_dk_a * i_mugd_kenu &
! - dm_ge_a * dm_dk_a * i_mugd_nuke &
! - dm_ge_b * dm_dk_b * i_mugd_nuke )
! enddo
! enddo
! enddo
! enddo
! enddo
! enddo
! !$OMP END DO
! !$OMP END PARALLEL
call wall_time(tf) call wall_time(tf)
print *, ' total Wall time for fock_3e_uhf_ao_a =', tf - ti print *, ' total Wall time for fock_3e_uhf_ao_a =', tf - ti
@ -319,6 +370,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)]
double precision :: dm_dk_a, dm_dk_b, dm_dk double precision :: dm_dk_a, dm_dk_b, dm_dk
double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu
double precision :: ti, tf double precision :: ti, tf
double precision, allocatable :: f_tmp(:,:)
print *, ' PROVIDING fock_3e_uhf_ao_b ...' print *, ' PROVIDING fock_3e_uhf_ao_b ...'
call wall_time(ti) call wall_time(ti)
@ -326,37 +378,37 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)]
fock_3e_uhf_ao_b = 0.d0 fock_3e_uhf_ao_b = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, & !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp, &
!$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
!$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b) !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b)
allocate(f_tmp(ao_num,ao_num))
f_tmp = 0.d0
!$OMP DO !$OMP DO
do g = 1, ao_num do g = 1, ao_num
do e = 1, ao_num do e = 1, ao_num
dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
dm_ge = dm_ge_a + dm_ge_b dm_ge = dm_ge_a + dm_ge_b
do d = 1, ao_num do d = 1, ao_num
do k = 1, ao_num do k = 1, ao_num
dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
dm_dk = dm_dk_a + dm_dk_b dm_dk = dm_dk_a + dm_dk_b
do mu = 1, ao_num do mu = 1, ao_num
do nu = 1, ao_num do nu = 1, ao_num
call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
f_tmp(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek &
fock_3e_uhf_ao_b(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek &
+ dm_ge_b * dm_dk_b * i_mugd_eknu & + dm_ge_b * dm_dk_b * i_mugd_eknu &
+ dm_ge_b * dm_dk_b * i_mugd_knue & + dm_ge_b * dm_dk_b * i_mugd_knue &
- dm_ge * dm_dk_b * i_mugd_kenu &
- dm_ge_b * dm_dk * i_mugd_enuk & - dm_ge_b * dm_dk * i_mugd_enuk &
- dm_ge * dm_dk_b * i_mugd_kenu &
- dm_ge_b * dm_dk_b * i_mugd_nuke & - dm_ge_b * dm_dk_b * i_mugd_nuke &
- dm_ge_a * dm_dk_a * i_mugd_nuke ) - dm_ge_a * dm_dk_a * i_mugd_nuke )
enddo enddo
@ -365,9 +417,57 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)]
enddo enddo
enddo enddo
enddo enddo
!$OMP END DO !$OMP END DO NOWAIT
!$OMP CRITICAL
do mu = 1, ao_num
do nu = 1, ao_num
fock_3e_uhf_ao_b(mu,nu) += f_tmp(mu,nu)
enddo
enddo
!$OMP END CRITICAL
deallocate(f_tmp)
!$OMP END PARALLEL !$OMP END PARALLEL
! TODO
! !$OMP PARALLEL DO DEFAULT (NONE) &
! !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, &
! !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
! !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b)
! do g = 1, ao_num
! do e = 1, ao_num
! dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
! dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
! dm_ge = dm_ge_a + dm_ge_b
! do d = 1, ao_num
! do k = 1, ao_num
! dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
! dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
! dm_dk = dm_dk_a + dm_dk_b
! do mu = 1, ao_num
! do nu = 1, ao_num
! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
! call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
! fock_3e_uhf_ao_b(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek &
! + dm_ge_b * dm_dk_b * i_mugd_eknu &
! + dm_ge_b * dm_dk_b * i_mugd_knue &
! - dm_ge_b * dm_dk * i_mugd_enuk &
! - dm_ge * dm_dk_b * i_mugd_kenu &
! - dm_ge_b * dm_dk_b * i_mugd_nuke &
! - dm_ge_a * dm_dk_a * i_mugd_nuke )
! enddo
! enddo
! enddo
! enddo
! enddo
! enddo
! !$OMP END PARALLEL DO
call wall_time(tf) call wall_time(tf)
print *, ' total Wall time for fock_3e_uhf_ao_b =', tf - ti print *, ' total Wall time for fock_3e_uhf_ao_b =', tf - ti

8
src/tc_scf/test_int.irp.f
View File

@ -25,7 +25,7 @@ program test_ints
!call routine_int2_grad1u2_grad2u2_j1b2 !call routine_int2_grad1u2_grad2u2_j1b2
!call test_fock_3e_uhf_ao() call test_fock_3e_uhf_ao()
call test_fock_3e_uhf_mo() call test_fock_3e_uhf_mo()
end end
@ -364,11 +364,10 @@ subroutine test_fock_3e_uhf_ao()
thr_ih = 1d-7 thr_ih = 1d-7
PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth
PROVIDE fock_3e_uhf_ao_a fock_3e_uhf_ao_b
! --- ! ---
PROVIDE fock_3e_uhf_ao_a
allocate(fock_3e_uhf_ao_a_mo(mo_num,mo_num)) allocate(fock_3e_uhf_ao_a_mo(mo_num,mo_num))
call ao_to_mo_bi_ortho( fock_3e_uhf_ao_a , size(fock_3e_uhf_ao_a , 1) & call ao_to_mo_bi_ortho( fock_3e_uhf_ao_a , size(fock_3e_uhf_ao_a , 1) &
, fock_3e_uhf_ao_a_mo, size(fock_3e_uhf_ao_a_mo, 1) ) , fock_3e_uhf_ao_a_mo, size(fock_3e_uhf_ao_a_mo, 1) )
@ -397,8 +396,6 @@ subroutine test_fock_3e_uhf_ao()
! --- ! ---
PROVIDE fock_3e_uhf_ao_b
allocate(fock_3e_uhf_ao_b_mo(mo_num,mo_num)) allocate(fock_3e_uhf_ao_b_mo(mo_num,mo_num))
call ao_to_mo_bi_ortho( fock_3e_uhf_ao_b , size(fock_3e_uhf_ao_b , 1) & call ao_to_mo_bi_ortho( fock_3e_uhf_ao_b , size(fock_3e_uhf_ao_b , 1) &
, fock_3e_uhf_ao_b_mo, size(fock_3e_uhf_ao_b_mo, 1) ) , fock_3e_uhf_ao_b_mo, size(fock_3e_uhf_ao_b_mo, 1) )
@ -421,6 +418,7 @@ subroutine test_fock_3e_uhf_ao()
enddo enddo
enddo enddo
print *, ' diff on F_b = ', diff_tot/norm print *, ' diff on F_b = ', diff_tot/norm
print *, ' '
deallocate(fock_3e_uhf_ao_b_mo) deallocate(fock_3e_uhf_ao_b_mo)