LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-17 10:05:18 +02:00
Code Releases Activity

Merge pull request #279 from AbdAmmar/dev-stable-tc-scf
Some checks failed
continuous-integration/drone/push Build is failing
Details

Browse Source 
fixed bug in TC integrals
This commit is contained in:
Emmanuel Giner 2023-04-26 10:05:53 +02:00 committed by GitHub
commit aa97c943b6
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
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

@ -11,14 +11,17 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
provide j1b_type
if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
print *, ' do things properly !'
stop
do i = 1, ao_num
do j = 1, ao_num
ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) &
+ j1b_gauss_hermII (j,i) &
+ j1b_gauss_nonherm(j,i) )
enddo
enddo
!do i = 1, ao_num
! do j = 1, ao_num
! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) &
! + j1b_gauss_hermII (j,i) &
! + j1b_gauss_nonherm(j,i) )
! enddo
!enddo
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 ...'
call wall_time(wall0)
if(test_cycle_tc)then
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3)
enddo
enddo
enddo
if(test_cycle_tc) then
PROVIDE int2_grad1_u12_ao_test
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3)
enddo
enddo
enddo
else
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3)
enddo
enddo
enddo
PROVIDE int2_grad1_u12_ao
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, ao_num
int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1)
int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2)
int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3)
enddo
enddo
enddo
endif
call wall_time(wall1)
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
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 DEFAULT (NONE) &
!$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 PARALLEL
call wall_time(wall1)
!call wall_time(wall1)
!print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,3, mo_num, mo_num )]
implicit none
integer :: i, j, ipoint
do ipoint = 1, n_points_final_grid
do i = 1, mo_num
do j = 1, mo_num
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
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)]
implicit none
integer :: i, j, ipoint
PROVIDE mo_l_coef mo_r_coef
PROVIDE int2_grad1_u12_bimo_transp
do ipoint = 1, n_points_final_grid
do i = 1, mo_num
do j = 1, mo_num
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
! ---

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
double precision :: weight
PROVIDE mo_l_coef mo_r_coef
PROVIDE int2_grad1_u12_bimo_t
integral = 0.d0
do ipoint = 1, n_points_final_grid
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) &
* ( 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) &
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,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) )
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,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
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,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) )
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,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
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,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) )
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
END_PROVIDER
! ---
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
PROVIDE ao_tc_int_chemist
do j = 1, ao_num
do l = 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, allocatable :: T(:,:)
PROVIDE mo_l_coef mo_r_coef
allocate ( T(ao_num,mo_num) )
!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, allocatable :: tmp_1(:,:), tmp_2(:,:)
PROVIDE mo_l_coef mo_r_coef
! ao_overlap x mo_r_coef
allocate( tmp_1(ao_num,mo_num) )
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
!
! 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
@ -373,6 +373,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
else
PROVIDE int2_grad1_u12_square_ao
allocate(b_mat(n_points_final_grid,ao_num,ao_num))
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
tc_grad_square_ao = 0.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 &
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, 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)
deallocate(b_mat)

139
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_DOC
!
! 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)
!
! 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)
! = -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,:)
!
! 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) ]
! - \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) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:)
! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
!
!
END_DOC
implicit none
integer :: ipoint, i, j
double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
PROVIDE j1b_type
if(j1b_type .eq. 3) then
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)
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(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
else
int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao
endif
END_PROVIDER
!BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)]
!
! BEGIN_DOC
! !
! ! 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)
! !
! ! 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)
! ! = -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,:)
! !
! ! 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) ]
! ! - \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) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:)
! ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
! !
! !
! END_DOC
!
! implicit none
! integer :: ipoint, i, j
! double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
!
! PROVIDE j1b_type
!
! if(j1b_type .eq. 3) then
!
! 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)
!
! 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(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
!
! else
!
! int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao
!
! endif
!
!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 >
!
! = 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.
!
@ -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 ...'
call wall_time(time0)
PROVIDE int2_grad1_u12_ao
if(read_tc_integ) then
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
PROVIDE int2_grad1_u12_ao
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3))
b_mat = 0.d0
@ -247,10 +250,10 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
do m = 1, 3
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 &
, 0.d0, tc_grad_and_lapl_ao, ao_num*ao_num)
, 1.d0, tc_grad_and_lapl_ao, ao_num*ao_num)
enddo
deallocate(b_mat)
call sum_A_At(tc_grad_and_lapl_ao(1,1,1,1), ao_num*ao_num)
! !$OMP PARALLEL &
! !$OMP DEFAULT (NONE) &

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_squared_ao, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_DOC
!
@ -23,10 +22,6 @@
! - 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)
!
!
!
! int2_grad1_u12_squared_ao = int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2
!
END_DOC
implicit none
@ -34,25 +29,22 @@
double precision :: time0, time1
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
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")
read(11) int2_grad1_u12_ao
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
endif
if(j1b_type .eq. 3) 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
! TODO OPENMP
do ipoint = 1, n_points_final_grid
@ -73,42 +65,14 @@
enddo
enddo
enddo
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
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
! ---
call wall_time(time0)
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
PROVIDE grad1_u12_squared_num
PROVIDE grad1_u12_num
double precision, allocatable :: tmp(:,:,:)
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
tmp = 0.d0
@ -130,34 +94,14 @@
if(.not.read_tc_integ) then
int2_grad1_u12_ao = 0.d0
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 &
, 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num)
enddo
!! DEBUG
!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
!! these dgemm are equivalen to
!!$OMP PARALLEL &
!!$OMP DEFAULT (NONE) &
!!$OMP PRIVATE (j, i, ipoint, jpoint, w) &
@ -169,7 +113,9 @@
! do j = 1, ao_num
! do i = 1, ao_num
! 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,2) += w * grad1_u12_num(jpoint,ipoint,2)
! int2_grad1_u12_ao(i,j,ipoint,3) += w * grad1_u12_num(jpoint,ipoint,3)
@ -179,67 +125,15 @@
!enddo
!!$OMP END DO
!!$OMP END PARALLEL
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao =', time1-time0
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)
else
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
stop
endif
if(write_tc_integ.and.mpi_master) then
@ -250,6 +144,112 @@
call ezfio_set_tc_keywords_io_tc_integ('Read')
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
! ---

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 ...'
call wall_time(wall0)
if(test_cycle_tc)then
if(test_cycle_tc) then
ao_tc_int_chemist = ao_tc_int_chemist_test
else
PROVIDE tc_grad_square_ao tc_grad_and_lapl_ao ao_two_e_coul
do j = 1, ao_num
do l = 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
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
integer :: i, j, k, l
double precision :: wall1, wall0
print *, ' providing ao_tc_int_chemist_no_cycle ...'
call wall_time(wall0)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 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)
! ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j)
enddo
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 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)
!ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j)
enddo
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)]
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))
if(var_tc) then
do i = 1, ao_num
do j = 1, ao_num
F(j,i) = Fock_matrix_vartc_ao_tot(j,i)
enddo
enddo
else
PROVIDE Fock_matrix_tc_ao_tot
do i = 1, ao_num
do j = 1, ao_num
F(j,i) = Fock_matrix_tc_ao_tot(j,i)
enddo
enddo
endif
allocate(tmp(ao_num,ao_num))
@ -139,6 +144,9 @@ BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)]
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) &
, 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_DOC
! ALPHA part of the Fock matrix from three-electron terms
!
! WARNING :: non hermitian if bi-ortho MOS used
!
! ALPHA part of the Fock matrix from three-electron terms
!
! WARNING :: non hermitian if bi-ortho MOS used
!
END_DOC
implicit none
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 :: ti, tf
PROVIDE mo_l_coef mo_r_coef
PROVIDE fock_3e_uhf_mo_cs
!print *, ' PROVIDING fock_3e_uhf_mo_a ...'
call wall_time(ti)
!call wall_time(ti)
o = elec_beta_num + 1
@ -142,7 +146,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
enddo
enddo
call wall_time(tf)
!call wall_time(tf)
!print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti
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 :: t0, t1
PROVIDE ao_two_e_tc_tot
!print*, ' providing two_e_tc_non_hermit_integral_seq ...'
!call wall_time(t0)
@ -80,6 +82,10 @@ END_PROVIDER
double precision :: t0, t1
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 ...'
!call wall_time(t0)
@ -142,7 +148,7 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_alpha, (ao_num, ao_num)]
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
@ -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))
!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) &
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
if(three_body_h_tc) then
!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
endif
else
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) )
@ -276,6 +285,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
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) &
, 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_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
! 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))&
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
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))
enddo
! F+K/2
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))&
enddo
! F+K/2
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))&
+ 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
! F
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))
enddo
enddo
enddo
! F
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))
enddo
enddo
do j=elec_beta_num+1,elec_alpha_num
! F+K/2
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))&
do j = elec_beta_num+1, elec_alpha_num
! F+K/2
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))&
+ 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
! F
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))
enddo
! F-K/2
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))&
enddo
! F
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))
enddo
! F-K/2
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))&
- 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
! F
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))
enddo
! F-K/2
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))&
do j = elec_alpha_num+1, mo_num
! F
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))
enddo
! F-K/2
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))&
- 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
! F+K
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)) &
enddo
! F+K
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_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
enddo
if(three_body_h_tc)then
enddo
enddo
if(three_body_h_tc) then
PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth
! C-O
do j = 1, elec_beta_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(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
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(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
enddo
! C-V
do j = 1, elec_beta_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(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
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(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
enddo
! O-V
do j = elec_beta_num+1, elec_alpha_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(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
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(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
enddo
enddo
endif
endif
endif
endif
do i = 1, mo_num
Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i)
enddo
do i = 1, mo_num
Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i)
enddo
if(frozen_orb_scf)then
@ -116,28 +129,29 @@
enddo
endif
if(no_oa_or_av_opt)then
do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
do j = 1, n_core_orb
jorb = list_core(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
enddo
endif
if(no_oa_or_av_opt)then
do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
do j = 1, n_core_orb
jorb = list_core(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
enddo
endif
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
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_DOC
! Alpha part of the Fock matrix from three-electron terms
!
! WARNING :: non hermitian if bi-ortho MOS used
! Alpha part of the Fock matrix from three-electron terms
!
! WARNING :: non hermitian if bi-ortho MOS used
END_DOC
implicit none
integer :: i, a
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

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

@ -11,6 +11,7 @@
integer :: i, j
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_one_e_energy = 0.d0