PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-19 04:22:32 +01:00
Code Releases Activity

fixed bug in TC integrals

Browse Source
This commit is contained in:
AbdAmmar 2023-04-26 08:52:06 +02:00
parent e1dc8b9ec3
commit 207e52d220
15 changed files with 453 additions and 364 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -12,13 +12,16 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
do i = 1, ao_num print *, ' do things properly !'
do j = 1, ao_num stop
ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) &
+ j1b_gauss_hermII (j,i) & !do i = 1, ao_num
+ j1b_gauss_nonherm(j,i) ) ! do j = 1, ao_num
enddo ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) &
enddo ! + j1b_gauss_hermII (j,i) &
! + j1b_gauss_nonherm(j,i) )
! enddo
!enddo
endif endif

86
src/bi_ort_ints/semi_num_ints_mo.irp.f
View File

@ -110,27 +110,36 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
print *, ' providing int2_grad1_u12_ao_transp ...' print *, ' providing int2_grad1_u12_ao_transp ...'
call wall_time(wall0) call wall_time(wall0)
if(test_cycle_tc)then if(test_cycle_tc) then
do ipoint = 1, n_points_final_grid
do i = 1, ao_num PROVIDE int2_grad1_u12_ao_test
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1) do ipoint = 1, n_points_final_grid
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2) do i = 1, ao_num
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3) do j = 1, ao_num
enddo int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1)
enddo int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2)
enddo int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3)
enddo
enddo
enddo
else else
do ipoint = 1, n_points_final_grid
do i = 1, ao_num PROVIDE int2_grad1_u12_ao
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1) do ipoint = 1, n_points_final_grid
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2) do i = 1, ao_num
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3) do j = 1, ao_num
enddo int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1)
enddo int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2)
enddo int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3)
enddo
enddo
enddo
endif endif
call wall_time(wall1) call wall_time(wall1)
print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0 print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0
@ -144,9 +153,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num,
integer :: ipoint integer :: ipoint
double precision :: wall0, wall1 double precision :: wall0, wall1
!print *, ' providing int2_grad1_u12_bimo_transp' PROVIDE mo_l_coef mo_r_coef
PROVIDE int2_grad1_u12_ao_transp
!print *, ' providing int2_grad1_u12_bimo_transp'
!call wall_time(wall0)
call wall_time(wall0)
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint) & !$OMP PRIVATE (ipoint) &
@ -163,25 +175,31 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num,
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
call wall_time(wall1) !call wall_time(wall1)
!print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0 !print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,3, mo_num, mo_num )] BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)]
implicit none
integer :: i, j, ipoint implicit none
do ipoint = 1, n_points_final_grid integer :: i, j, ipoint
do i = 1, mo_num
do j = 1, mo_num PROVIDE mo_l_coef mo_r_coef
int2_grad1_u12_bimo_t(ipoint,1,j,i) = int2_grad1_u12_bimo_transp(j,i,1,ipoint) PROVIDE int2_grad1_u12_bimo_transp
int2_grad1_u12_bimo_t(ipoint,2,j,i) = int2_grad1_u12_bimo_transp(j,i,2,ipoint)
int2_grad1_u12_bimo_t(ipoint,3,j,i) = int2_grad1_u12_bimo_transp(j,i,3,ipoint) do ipoint = 1, n_points_final_grid
enddo do i = 1, mo_num
enddo do j = 1, mo_num
enddo int2_grad1_u12_bimo_t(ipoint,1,j,i) = int2_grad1_u12_bimo_transp(j,i,1,ipoint)
int2_grad1_u12_bimo_t(ipoint,2,j,i) = int2_grad1_u12_bimo_transp(j,i,2,ipoint)
int2_grad1_u12_bimo_t(ipoint,3,j,i) = int2_grad1_u12_bimo_transp(j,i,3,ipoint)
enddo
enddo
enddo
END_PROVIDER END_PROVIDER
! --- ! ---

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

@ -81,21 +81,24 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
integer :: ipoint integer :: ipoint
double precision :: weight double precision :: weight
PROVIDE mo_l_coef mo_r_coef
PROVIDE int2_grad1_u12_bimo_t
integral = 0.d0 integral = 0.d0
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint) weight = final_weight_at_r_vector(ipoint)
integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
* ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j) & * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j) &
+ int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j) & + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j) &
+ int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,l,j) ) + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,l,j) )
integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
* ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,k,i) ) + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
* ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & * ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,3,l,j) * int2_grad1_u12_bimo_t(ipoint,3,k,i) ) + int2_grad1_u12_bimo_t(ipoint,3,l,j) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
enddo enddo

3
src/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -20,6 +20,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_vartc_tot, (ao_num, ao_num, ao_num, a
enddo enddo
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_num) ] BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_num) ]
@ -66,6 +67,8 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
else else
PROVIDE ao_tc_int_chemist
do j = 1, ao_num do j = 1, ao_num
do l = 1, ao_num do l = 1, ao_num
do i = 1, ao_num do i = 1, ao_num

4
src/bi_ortho_mos/mos_rl.irp.f
View File

@ -17,6 +17,8 @@ subroutine ao_to_mo_bi_ortho(A_ao, LDA_ao, A_mo, LDA_mo)
double precision, intent(out) :: A_mo(LDA_mo,mo_num) double precision, intent(out) :: A_mo(LDA_mo,mo_num)
double precision, allocatable :: T(:,:) double precision, allocatable :: T(:,:)
PROVIDE mo_l_coef mo_r_coef
allocate ( T(ao_num,mo_num) ) allocate ( T(ao_num,mo_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
@ -54,6 +56,8 @@ subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao)
double precision, intent(out) :: A_ao(LDA_ao,ao_num) double precision, intent(out) :: A_ao(LDA_ao,ao_num)
double precision, allocatable :: tmp_1(:,:), tmp_2(:,:) double precision, allocatable :: tmp_1(:,:), tmp_2(:,:)
PROVIDE mo_l_coef mo_r_coef
! ao_overlap x mo_r_coef ! ao_overlap x mo_r_coef
allocate( tmp_1(ao_num,mo_num) ) allocate( tmp_1(ao_num,mo_num) )
call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0 & call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0 &

8
src/non_h_ints_mu/grad_squared.irp.f
View File

@ -352,7 +352,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
BEGIN_DOC BEGIN_DOC
! !
! tc_grad_square_ao(k,i,l,j) = 1/2 <kl | |\grad_1 u(r1,r2)|^2 + |\grad_2 u(r1,r2)|^2 | ij> ! tc_grad_square_ao(k,i,l,j) = -1/2 <kl | |\grad_1 u(r1,r2)|^2 + |\grad_2 u(r1,r2)|^2 | ij>
! !
END_DOC END_DOC
@ -373,6 +373,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
else else
PROVIDE int2_grad1_u12_square_ao
allocate(b_mat(n_points_final_grid,ao_num,ao_num)) allocate(b_mat(n_points_final_grid,ao_num,ao_num))
b_mat = 0.d0 b_mat = 0.d0
@ -392,8 +394,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
!$OMP END PARALLEL !$OMP END PARALLEL
tc_grad_square_ao = 0.d0 tc_grad_square_ao = 0.d0
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, int2_grad1_u12_squared_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & , int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
, 0.d0, tc_grad_square_ao, ao_num*ao_num) , 0.d0, tc_grad_square_ao, ao_num*ao_num)
deallocate(b_mat) deallocate(b_mat)

137
src/non_h_ints_mu/new_grad_tc.irp.f
View File

@ -1,68 +1,68 @@
! --- ! ---
BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)] !BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)]
!
BEGIN_DOC ! BEGIN_DOC
! ! !
! int1_grad2_u12_ao(:,i,j,ipoint) = \int dr1 [-1 * \grad_r2 J(r1,r2)] \phi_i(r1) \phi_j(r1) ! ! int1_grad2_u12_ao(:,i,j,ipoint) = \int dr1 [-1 * \grad_r2 J(r1,r2)] \phi_i(r1) \phi_j(r1)
! ! !
! where r1 = r(ipoint) ! ! where r1 = r(ipoint)
! ! !
! if J(r1,r2) = u12: ! ! if J(r1,r2) = u12:
! ! !
! int1_grad2_u12_ao(:,i,j,ipoint) = +0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r1) \phi_j(r1) ! ! int1_grad2_u12_ao(:,i,j,ipoint) = +0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r1) \phi_j(r1)
! = -0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ] ! ! = -0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ]
! = -int2_grad1_u12_ao(i,j,ipoint,:) ! ! = -int2_grad1_u12_ao(i,j,ipoint,:)
! ! !
! if J(r1,r2) = u12 x v1 x v2 ! ! if J(r1,r2) = u12 x v1 x v2
! ! !
! int1_grad2_u12_ao(:,i,j,ipoint) = v2 x [ 0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] v1 \phi_i(r1) \phi_j(r1) ] ! ! int1_grad2_u12_ao(:,i,j,ipoint) = v2 x [ 0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] v1 \phi_i(r1) \phi_j(r1) ]
! - \grad_2 v2 x [ \int dr1 u12 v1 \phi_i(r1) \phi_j(r1) ] ! ! - \grad_2 v2 x [ \int dr1 u12 v1 \phi_i(r1) \phi_j(r1) ]
! = -0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) ! ! = -0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:)
! + 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) ! ! + 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:)
! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) ! ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
! ! !
! ! !
END_DOC ! END_DOC
!
implicit none ! implicit none
integer :: ipoint, i, j ! integer :: ipoint, i, j
double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 ! double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
!
PROVIDE j1b_type ! PROVIDE j1b_type
!
if(j1b_type .eq. 3) then ! if(j1b_type .eq. 3) then
!
do ipoint = 1, n_points_final_grid ! do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint) ! x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint) ! y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint) ! z = final_grid_points(3,ipoint)
!
tmp0 = 0.5d0 * v_1b(ipoint) ! tmp0 = 0.5d0 * v_1b(ipoint)
tmp_x = v_1b_grad(1,ipoint) ! tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint) ! tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint) ! tmp_z = v_1b_grad(3,ipoint)
!
do j = 1, ao_num ! do j = 1, ao_num
do i = 1, ao_num ! do i = 1, ao_num
!
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) ! tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) ! tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
!
int1_grad2_u12_ao(1,i,j,ipoint) = -tmp1 * x + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x ! int1_grad2_u12_ao(1,i,j,ipoint) = -tmp1 * x + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y ! int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z ! int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
enddo ! enddo
enddo ! enddo
enddo ! enddo
!
else ! else
!
int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao ! int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao
!
endif ! endif
!
END_PROVIDER !END_PROVIDER
! --- ! ---
@ -192,7 +192,10 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
! !
! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) . \grad_1 | ij > ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) . \grad_1 | ij >
! !
! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) ! = -1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2)
! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 (-1) \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2)
!
! -1 in \int dr2
! !
! This is obtained by integration by parts. ! This is obtained by integration by parts.
! !
@ -207,8 +210,6 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
print*, ' providing tc_grad_and_lapl_ao ...' print*, ' providing tc_grad_and_lapl_ao ...'
call wall_time(time0) call wall_time(time0)
PROVIDE int2_grad1_u12_ao
if(read_tc_integ) then if(read_tc_integ) then
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao', action="read") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao', action="read")
@ -217,6 +218,8 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
else else
PROVIDE int2_grad1_u12_ao
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3)) allocate(b_mat(n_points_final_grid,ao_num,ao_num,3))
b_mat = 0.d0 b_mat = 0.d0
@ -247,7 +250,7 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
do m = 1, 3 do m = 1, 3
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & , int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
, 0.d0, tc_grad_and_lapl_ao, ao_num*ao_num) , 1.d0, tc_grad_and_lapl_ao, ao_num*ao_num)
enddo enddo
deallocate(b_mat) deallocate(b_mat)

246
src/non_h_ints_mu/tc_integ.irp.f
View File

@ -1,8 +1,7 @@
! --- ! ---
BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)] BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)]
&BEGIN_PROVIDER [double precision, int2_grad1_u12_squared_ao, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
@ -23,10 +22,6 @@
! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) ! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:)
! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
! !
!
!
! int2_grad1_u12_squared_ao = int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2
!
END_DOC END_DOC
implicit none implicit none
@ -34,25 +29,22 @@
double precision :: time0, time1 double precision :: time0, time1
double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
print*, ' providing int2_grad1_u12_ao & int2_grad1_u12_squared_ao ...' print*, ' providing int2_grad1_u12_ao ...'
call wall_time(time0)
PROVIDE j1b_type PROVIDE j1b_type
if(read_tc_integ) then if(read_tc_integ) then
call wall_time(time0)
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read") open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read")
read(11) int2_grad1_u12_ao read(11) int2_grad1_u12_ao
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
endif endif
if(j1b_type .eq. 3) then if(j1b_type .eq. 3) then
! ---
if(.not.read_tc_integ) then if(.not.read_tc_integ) then
call wall_time(time0)
PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
int2_grad1_u12_ao = 0.d0 int2_grad1_u12_ao = 0.d0
! TODO OPENMP ! TODO OPENMP
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -73,42 +65,14 @@
enddo enddo
enddo enddo
enddo enddo
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
endif endif
! ---
call wall_time(time0)
PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
int2_grad1_u12_squared_ao = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, ipoint) &
!$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
!$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
do i = 1, ao_num
int2_grad1_u12_squared_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_squared_ao =', time1-time0
! ---
elseif(j1b_type .ge. 100) then elseif(j1b_type .ge. 100) then
! ---
call wall_time(time0)
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
PROVIDE grad1_u12_squared_num
PROVIDE grad1_u12_num PROVIDE grad1_u12_num
double precision, allocatable :: tmp(:,:,:) double precision, allocatable :: tmp(:,:,:)
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
tmp = 0.d0 tmp = 0.d0
@ -130,34 +94,14 @@
if(.not.read_tc_integ) then if(.not.read_tc_integ) then
int2_grad1_u12_ao = 0.d0 int2_grad1_u12_ao = 0.d0
do m = 1, 3 do m = 1, 3
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, 1.d0 & !call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0 &
! this work also because of the symmetry in K(1,2) and sign compensation in L(1,2,3)
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -1.d0 &
, tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid & , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid &
, 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num) , 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num)
enddo enddo
!! DEBUG !! these dgemm are equivalen to
!PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
!int2_grad1_u12_ao = 0.d0
!do ipoint = 1, n_points_final_grid
! x = final_grid_points(1,ipoint)
! y = final_grid_points(2,ipoint)
! z = final_grid_points(3,ipoint)
! tmp0 = 0.5d0 * v_1b(ipoint)
! tmp_x = v_1b_grad(1,ipoint)
! tmp_y = v_1b_grad(2,ipoint)
! tmp_z = v_1b_grad(3,ipoint)
! do j = 1, ao_num
! do i = 1, ao_num
! tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
! tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
! int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
! int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
! int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
! enddo
! enddo
!enddo
! these dgemm are equivalen to
!!$OMP PARALLEL & !!$OMP PARALLEL &
!!$OMP DEFAULT (NONE) & !!$OMP DEFAULT (NONE) &
!!$OMP PRIVATE (j, i, ipoint, jpoint, w) & !!$OMP PRIVATE (j, i, ipoint, jpoint, w) &
@ -169,7 +113,9 @@
! do j = 1, ao_num ! do j = 1, ao_num
! do i = 1, ao_num ! do i = 1, ao_num
! do jpoint = 1, n_points_extra_final_grid ! do jpoint = 1, n_points_extra_final_grid
! w = tmp(jpoint,i,j) ! w = -tmp(jpoint,i,j)
! !w = tmp(jpoint,i,j) this work also because of the symmetry in K(1,2)
! ! and sign compensation in L(1,2,3)
! int2_grad1_u12_ao(i,j,ipoint,1) += w * grad1_u12_num(jpoint,ipoint,1) ! int2_grad1_u12_ao(i,j,ipoint,1) += w * grad1_u12_num(jpoint,ipoint,1)
! int2_grad1_u12_ao(i,j,ipoint,2) += w * grad1_u12_num(jpoint,ipoint,2) ! int2_grad1_u12_ao(i,j,ipoint,2) += w * grad1_u12_num(jpoint,ipoint,2)
! int2_grad1_u12_ao(i,j,ipoint,3) += w * grad1_u12_num(jpoint,ipoint,3) ! int2_grad1_u12_ao(i,j,ipoint,3) += w * grad1_u12_num(jpoint,ipoint,3)
@ -179,67 +125,15 @@
!enddo !enddo
!!$OMP END DO !!$OMP END DO
!!$OMP END PARALLEL !!$OMP END PARALLEL
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
endif endif
! ---
call wall_time(time0)
int2_grad1_u12_squared_ao = 0.d0
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 &
, tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid &
, 0.d0, int2_grad1_u12_squared_ao(1,1,1), ao_num*ao_num)
!! DEBUG
!PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
!int2_grad1_u12_squared_ao = 0.d0
!!$OMP PARALLEL &
!!$OMP DEFAULT (NONE) &
!!$OMP PRIVATE (i, j, ipoint) &
!!$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
!!$OMP DO SCHEDULE (static)
!do ipoint = 1, n_points_final_grid
! do j = 1, ao_num
! do i = 1, ao_num
! int2_grad1_u12_squared_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
! enddo
! enddo
!enddo
!!$OMP END DO
!!$OMP END PARALLEL
!call wall_time(time1)
!! this dgemm is equivalen to
!!$OMP PARALLEL &
!!$OMP DEFAULT (NONE) &
!!$OMP PRIVATE (i, j, ipoint, jpoint, w) &
!!$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, &
!!$OMP n_points_extra_final_grid, final_weight_at_r_vector_extra, &
!!$OMP aos_in_r_array_extra_transp, grad1_u12_squared_num, tmp)
!!$OMP DO SCHEDULE (static)
!do ipoint = 1, n_points_final_grid
! do j = 1, ao_num
! do i = 1, ao_num
! do jpoint = 1, n_points_extra_final_grid
! w = 0.5d0 * tmp(jpoint,i,j)
! int2_grad1_u12_squared_ao(i,j,ipoint) += w * grad1_u12_squared_num(jpoint,ipoint)
! enddo
! enddo
! enddo
!enddo
!!$OMP END DO
!!$OMP END PARALLEL
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_squared_ao =', time1-time0
! ---
deallocate(tmp) deallocate(tmp)
else else
print *, ' j1b_type = ', j1b_type, 'not implemented yet' print *, ' j1b_type = ', j1b_type, 'not implemented yet'
stop stop
endif endif
if(write_tc_integ.and.mpi_master) then if(write_tc_integ.and.mpi_master) then
@ -250,6 +144,112 @@
call ezfio_set_tc_keywords_io_tc_integ('Read') call ezfio_set_tc_keywords_io_tc_integ('Read')
endif endif
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int2_grad1_u12_square_ao = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2
!
END_DOC
implicit none
integer :: ipoint, i, j, m, jpoint
double precision :: time0, time1
double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
print*, ' providing int2_grad1_u12_square_ao ...'
call wall_time(time0)
PROVIDE j1b_type
if(j1b_type .eq. 3) then
PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
int2_grad1_u12_square_ao = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, ipoint) &
!$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
!$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
do i = 1, ao_num
int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
elseif(j1b_type .ge. 100) then
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
PROVIDE grad1_u12_squared_num
double precision, allocatable :: tmp(:,:,:)
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
tmp = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j, i, jpoint) &
!$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
!$OMP DO SCHEDULE (static)
do j = 1, ao_num
do i = 1, ao_num
do jpoint = 1, n_points_extra_final_grid
tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
int2_grad1_u12_square_ao = 0.d0
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 &
, tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid &
, 0.d0, int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num)
!! this dgemm is equivalen to
!!$OMP PARALLEL &
!!$OMP DEFAULT (NONE) &
!!$OMP PRIVATE (i, j, ipoint, jpoint, w) &
!!$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, &
!!$OMP n_points_extra_final_grid, final_weight_at_r_vector_extra, &
!!$OMP aos_in_r_array_extra_transp, grad1_u12_squared_num, tmp)
!!$OMP DO SCHEDULE (static)
!do ipoint = 1, n_points_final_grid
! do j = 1, ao_num
! do i = 1, ao_num
! do jpoint = 1, n_points_extra_final_grid
! w = -0.5d0 * tmp(jpoint,i,j)
! int2_grad1_u12_square_ao(i,j,ipoint) += w * grad1_u12_squared_num(jpoint,ipoint)
! enddo
! enddo
! enddo
!enddo
!!$OMP END DO
!!$OMP END PARALLEL
deallocate(tmp)
else
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
stop
endif
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_square_ao =', time1-time0
END_PROVIDER END_PROVIDER
! --- ! ---

30
src/non_h_ints_mu/total_tc_int.irp.f
View File

@ -48,9 +48,14 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
print *, ' providing ao_tc_int_chemist ...' print *, ' providing ao_tc_int_chemist ...'
call wall_time(wall0) call wall_time(wall0)
if(test_cycle_tc)then if(test_cycle_tc) then
ao_tc_int_chemist = ao_tc_int_chemist_test ao_tc_int_chemist = ao_tc_int_chemist_test
else else
PROVIDE tc_grad_square_ao tc_grad_and_lapl_ao ao_two_e_coul
do j = 1, ao_num do j = 1, ao_num
do l = 1, ao_num do l = 1, ao_num
do i = 1, ao_num do i = 1, ao_num
@ -68,27 +73,34 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, ao_num, ao_num)]
! --- ! ---
BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, ao_num, ao_num)]
implicit none implicit none
integer :: i, j, k, l integer :: i, j, k, l
double precision :: wall1, wall0 double precision :: wall1, wall0
print *, ' providing ao_tc_int_chemist_no_cycle ...' print *, ' providing ao_tc_int_chemist_no_cycle ...'
call wall_time(wall0) call wall_time(wall0)
do j = 1, ao_num
do l = 1, ao_num do j = 1, ao_num
do i = 1, ao_num do l = 1, ao_num
do k = 1, ao_num do i = 1, ao_num
ao_tc_int_chemist_no_cycle(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) do k = 1, ao_num
! ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j) ao_tc_int_chemist_no_cycle(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
enddo !ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j)
enddo enddo
enddo enddo
enddo enddo
enddo
call wall_time(wall1) call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0 print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)] BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)]
implicit none implicit none

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

@ -92,17 +92,22 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
allocate(F(ao_num,ao_num)) allocate(F(ao_num,ao_num))
if(var_tc) then if(var_tc) then
do i = 1, ao_num do i = 1, ao_num
do j = 1, ao_num do j = 1, ao_num
F(j,i) = Fock_matrix_vartc_ao_tot(j,i) F(j,i) = Fock_matrix_vartc_ao_tot(j,i)
enddo enddo
enddo enddo
else else
PROVIDE Fock_matrix_tc_ao_tot
do i = 1, ao_num do i = 1, ao_num
do j = 1, ao_num do j = 1, ao_num
F(j,i) = Fock_matrix_tc_ao_tot(j,i) F(j,i) = Fock_matrix_tc_ao_tot(j,i)
enddo enddo
enddo enddo
endif endif
allocate(tmp(ao_num,ao_num)) allocate(tmp(ao_num,ao_num))
@ -139,6 +144,9 @@ BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)]
implicit none implicit none
PROVIDE mo_r_coef mo_l_coef
PROVIDE FQS_SQF_ao
call ao_to_mo_bi_ortho( FQS_SQF_ao, size(FQS_SQF_ao, 1) & call ao_to_mo_bi_ortho( FQS_SQF_ao, size(FQS_SQF_ao, 1) &
, FQS_SQF_mo, size(FQS_SQF_mo, 1) ) , FQS_SQF_mo, size(FQS_SQF_mo, 1) )

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

@ -50,19 +50,23 @@ END_PROVIDER
BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
BEGIN_DOC BEGIN_DOC
! ALPHA part of the Fock matrix from three-electron terms !
! ! ALPHA part of the Fock matrix from three-electron terms
! WARNING :: non hermitian if bi-ortho MOS used !
! WARNING :: non hermitian if bi-ortho MOS used
!
END_DOC END_DOC
implicit none implicit none
integer :: a, b, i, j, o integer :: a, b, i, j, o
double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia
double precision :: ti, tf double precision :: ti, tf
PROVIDE mo_l_coef mo_r_coef PROVIDE mo_l_coef mo_r_coef
PROVIDE fock_3e_uhf_mo_cs
!print *, ' PROVIDING fock_3e_uhf_mo_a ...' !print *, ' PROVIDING fock_3e_uhf_mo_a ...'
call wall_time(ti) !call wall_time(ti)
o = elec_beta_num + 1 o = elec_beta_num + 1
@ -142,7 +146,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
enddo enddo
enddo enddo
call wall_time(tf) !call wall_time(tf)
!print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti !print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti
END_PROVIDER END_PROVIDER

14
src/tc_scf/fock_tc.irp.f
View File

@ -18,6 +18,8 @@
double precision :: density, density_a, density_b double precision :: density, density_a, density_b
double precision :: t0, t1 double precision :: t0, t1
PROVIDE ao_two_e_tc_tot
!print*, ' providing two_e_tc_non_hermit_integral_seq ...' !print*, ' providing two_e_tc_non_hermit_integral_seq ...'
!call wall_time(t0) !call wall_time(t0)
@ -80,6 +82,10 @@ END_PROVIDER
double precision :: t0, t1 double precision :: t0, t1
double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
PROVIDE ao_two_e_tc_tot
PROVIDE mo_l_coef mo_r_coef
PROVIDE TCSCF_density_matrix_ao_alpha TCSCF_density_matrix_ao_beta
!print*, ' providing two_e_tc_non_hermit_integral ...' !print*, ' providing two_e_tc_non_hermit_integral ...'
!call wall_time(t0) !call wall_time(t0)
@ -142,7 +148,7 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_alpha, (ao_num, ao_num)]
implicit none implicit none
Fock_matrix_tc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_tc_non_hermit_integral_alpha Fock_matrix_tc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_tc_non_hermit_integral_alpha
END_PROVIDER END_PROVIDER
@ -181,14 +187,17 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
!call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1)) !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1))
!deallocate(tmp) !deallocate(tmp)
PROVIDE mo_l_coef mo_r_coef
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
if(three_body_h_tc) then if(three_body_h_tc) then
!Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
PROVIDE fock_3e_uhf_mo_a
Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a
endif endif
else else
call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
@ -276,6 +285,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
implicit none implicit none
PROVIDE mo_l_coef mo_r_coef
PROVIDE Fock_matrix_tc_mo_tot
call mo_to_ao_bi_ortho( Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) & call mo_to_ao_bi_ortho( Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) &
, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) ) , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) )

218
src/tc_scf/fock_tc_mo_tot.irp.f
View File

@ -1,107 +1,120 @@
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ] BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)] &BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)]
implicit none
BEGIN_DOC
! TC-Fock matrix on the MO basis. WARNING !!! NON HERMITIAN !!!
! For open shells, the ROHF Fock Matrix is ::
!
! | F-K | F + K/2 | F |
! |---------------------------------|
! | F + K/2 | F | F - K/2 |
! |---------------------------------|
! | F | F - K/2 | F + K |
!
!
! F = 1/2 (Fa + Fb)
!
! K = Fb - Fa
!
END_DOC
integer :: i,j,n
if (elec_alpha_num == elec_beta_num) then
Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha
else
do j=1,elec_beta_num BEGIN_DOC
! F-K ! TC-Fock matrix on the MO basis. WARNING !!! NON HERMITIAN !!!
do i=1,elec_beta_num !CC ! For open shells, the ROHF Fock Matrix is ::
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& !
! | F-K | F + K/2 | F |
! |---------------------------------|
! | F + K/2 | F | F - K/2 |
! |---------------------------------|
! | F | F - K/2 | F + K |
!
!
! F = 1/2 (Fa + Fb)
!
! K = Fb - Fa
!
END_DOC
implicit none
integer :: i, j, n
if(elec_alpha_num == elec_beta_num) then
PROVIDE Fock_matrix_tc_mo_alpha
Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha
else
PROVIDE Fock_matrix_tc_mo_beta Fock_matrix_tc_mo_alpha
do j = 1, elec_beta_num
! F-K
do i = 1, elec_beta_num !CC
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
- (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo enddo
! F+K/2 ! F+K/2
do i=elec_beta_num+1,elec_alpha_num !CA do i = elec_beta_num+1, elec_alpha_num !CA
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+ 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo enddo
! F ! F
do i=elec_alpha_num+1, mo_num !CV do i = elec_alpha_num+1, mo_num !CV
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
enddo enddo
enddo enddo
do j=elec_beta_num+1,elec_alpha_num do j = elec_beta_num+1, elec_alpha_num
! F+K/2 ! F+K/2
do i=1,elec_beta_num !AC do i = 1, elec_beta_num !AC
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+ 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo enddo
! F ! F
do i=elec_beta_num+1,elec_alpha_num !AA do i = elec_beta_num+1, elec_alpha_num !AA
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
enddo enddo
! F-K/2 ! F-K/2
do i=elec_alpha_num+1, mo_num !AV do i = elec_alpha_num+1, mo_num !AV
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
- 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo enddo
enddo enddo
do j=elec_alpha_num+1, mo_num do j = elec_alpha_num+1, mo_num
! F ! F
do i=1,elec_beta_num !VC do i = 1, elec_beta_num !VC
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
enddo enddo
! F-K/2 ! F-K/2
do i=elec_beta_num+1,elec_alpha_num !VA do i = elec_beta_num+1, elec_alpha_num !VA
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
- 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo enddo
! F+K ! F+K
do i=elec_alpha_num+1,mo_num !VV do i = elec_alpha_num+1, mo_num !VV
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) & Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) &
+ (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo enddo
enddo enddo
if(three_body_h_tc)then
if(three_body_h_tc) then
PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth
! C-O ! C-O
do j = 1, elec_beta_num do j = 1, elec_beta_num
do i = elec_beta_num+1, elec_alpha_num do i = elec_beta_num+1, elec_alpha_num
Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo enddo
enddo enddo
! C-V ! C-V
do j = 1, elec_beta_num do j = 1, elec_beta_num
do i = elec_alpha_num+1, mo_num do i = elec_alpha_num+1, mo_num
Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo enddo
enddo enddo
! O-V ! O-V
do j = elec_beta_num+1, elec_alpha_num do j = elec_beta_num+1, elec_alpha_num
do i = elec_alpha_num+1, mo_num do i = elec_alpha_num+1, mo_num
Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo enddo
enddo enddo
endif endif
endif endif
do i = 1, mo_num do i = 1, mo_num
Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i) Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i)
enddo enddo
if(frozen_orb_scf)then if(frozen_orb_scf)then
@ -116,28 +129,29 @@
enddo enddo
endif endif
if(no_oa_or_av_opt)then if(no_oa_or_av_opt)then
do i = 1, n_act_orb do i = 1, n_act_orb
iorb = list_act(i) iorb = list_act(i)
do j = 1, n_inact_orb do j = 1, n_inact_orb
jorb = list_inact(j) jorb = list_inact(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo enddo
do j = 1, n_virt_orb do j = 1, n_virt_orb
jorb = list_virt(j) jorb = list_virt(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo enddo
do j = 1, n_core_orb do j = 1, n_core_orb
jorb = list_core(j) jorb = list_core(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo enddo
enddo enddo
endif endif
if(.not.bi_ortho .and. three_body_h_tc)then if(.not.bi_ortho .and. three_body_h_tc)then
Fock_matrix_tc_mo_tot += fock_3_mat Fock_matrix_tc_mo_tot += fock_3_mat
endif endif
END_PROVIDER END_PROVIDER

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

@ -4,14 +4,16 @@
BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)]
BEGIN_DOC BEGIN_DOC
! Alpha part of the Fock matrix from three-electron terms ! Alpha part of the Fock matrix from three-electron terms
! !
! WARNING :: non hermitian if bi-ortho MOS used ! WARNING :: non hermitian if bi-ortho MOS used
END_DOC END_DOC
implicit none implicit none
integer :: i, a integer :: i, a
PROVIDE mo_l_coef mo_r_coef PROVIDE mo_l_coef mo_r_coef
PROVIDE fock_cs_3e_bi_orth fock_a_tmp1_bi_ortho fock_a_tmp2_bi_ortho
fock_a_tot_3e_bi_orth = 0.d0 fock_a_tot_3e_bi_orth = 0.d0

1
src/tc_scf/tc_scf_energy.irp.f
View File

@ -11,6 +11,7 @@
integer :: i, j integer :: i, j
PROVIDE mo_l_coef mo_r_coef PROVIDE mo_l_coef mo_r_coef
PROVIDE two_e_tc_non_hermit_integral_alpha two_e_tc_non_hermit_integral_beta
TC_HF_energy = nuclear_repulsion TC_HF_energy = nuclear_repulsion
TC_HF_one_e_energy = 0.d0 TC_HF_one_e_energy = 0.d0