mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-07 14:03:37 +01:00
TC-SCF CLEANED
This commit is contained in:
Abdallah Ammar
parent
c50018e8bd
commit
da8eac81e0
@ -322,6 +322,12 @@ END_PROVIDER
|
|||||||
|
|
||||||
BEGIN_PROVIDER [double precision, noL_0e]
|
BEGIN_PROVIDER [double precision, noL_0e]
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! < Phi_left | L | Phi_right >
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i, j, k, ipoint
|
integer :: i, j, k, ipoint
|
||||||
double precision :: t0, t1
|
double precision :: t0, t1
|
||||||
@ -330,7 +336,6 @@ BEGIN_PROVIDER [double precision, noL_0e]
|
|||||||
double precision, allocatable :: tmp_M(:,:), tmp_S(:), tmp_O(:), tmp_J(:,:)
|
double precision, allocatable :: tmp_M(:,:), tmp_S(:), tmp_O(:), tmp_J(:,:)
|
||||||
double precision, allocatable :: tmp_M_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:)
|
double precision, allocatable :: tmp_M_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:)
|
||||||
|
|
||||||
|
|
||||||
call wall_time(t0)
|
call wall_time(t0)
|
||||||
print*, " Providing noL_0e ..."
|
print*, " Providing noL_0e ..."
|
||||||
|
|
||||||
|
|||||||
@ -9,3 +9,39 @@ doc: If |true|, tc-scf has converged
|
|||||||
interface: ezfio,provider,ocaml
|
interface: ezfio,provider,ocaml
|
||||||
default: False
|
default: False
|
||||||
|
|
||||||
|
[max_dim_diis_tcscf]
|
||||||
|
type: integer
|
||||||
|
doc: Maximum size of the DIIS extrapolation procedure
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 15
|
||||||
|
|
||||||
|
[level_shift_tcscf]
|
||||||
|
type: Positive_float
|
||||||
|
doc: Energy shift on the virtual MOs to improve TCSCF convergence
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 0.
|
||||||
|
|
||||||
|
[im_thresh_tcscf]
|
||||||
|
type: Threshold
|
||||||
|
doc: Thresholds on the Imag part of energy
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 1.e-7
|
||||||
|
|
||||||
|
[thresh_tcscf]
|
||||||
|
type: Threshold
|
||||||
|
doc: Threshold on the convergence of the Hartree Fock energy.
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 1.e-8
|
||||||
|
|
||||||
|
[n_it_tcscf_max]
|
||||||
|
type: Strictly_positive_int
|
||||||
|
doc: Maximum number of SCF iterations
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 50
|
||||||
|
|
||||||
|
[tc_Brillouin_Right]
|
||||||
|
type: logical
|
||||||
|
doc: If |true|, impose only right-Brillouin condition
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: False
|
||||||
|
|
||||||
|
|||||||
@ -1,107 +0,0 @@
|
|||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, good_hermit_tc_fock_mat, (mo_num, mo_num)]
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! good_hermit_tc_fock_mat = Hermitian Upper triangular Fock matrix
|
|
||||||
!
|
|
||||||
! The converged eigenvectors of such matrix yield to orthonormal vectors satisfying the left Brillouin theorem
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: i, j
|
|
||||||
|
|
||||||
good_hermit_tc_fock_mat = Fock_matrix_tc_mo_tot
|
|
||||||
do j = 1, mo_num
|
|
||||||
do i = 1, j-1
|
|
||||||
good_hermit_tc_fock_mat(i,j) = Fock_matrix_tc_mo_tot(j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, hermit_average_tc_fock_mat, (mo_num, mo_num)]
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! hermit_average_tc_fock_mat = (F + F^\dagger)/2
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: i, j
|
|
||||||
|
|
||||||
hermit_average_tc_fock_mat = Fock_matrix_tc_mo_tot
|
|
||||||
do j = 1, mo_num
|
|
||||||
do i = 1, mo_num
|
|
||||||
hermit_average_tc_fock_mat(i,j) = 0.5d0 * (Fock_matrix_tc_mo_tot(j,i) + Fock_matrix_tc_mo_tot(i,j))
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
! ---
|
|
||||||
BEGIN_PROVIDER [ double precision, grad_hermit]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! square of gradient of the energy
|
|
||||||
END_DOC
|
|
||||||
if(symetric_fock_tc)then
|
|
||||||
grad_hermit = grad_hermit_average_tc_fock_mat
|
|
||||||
else
|
|
||||||
grad_hermit = grad_good_hermit_tc_fock_mat
|
|
||||||
endif
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, grad_good_hermit_tc_fock_mat]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! grad_good_hermit_tc_fock_mat = norm of gradients of the upper triangular TC fock
|
|
||||||
END_DOC
|
|
||||||
integer :: i, j
|
|
||||||
grad_good_hermit_tc_fock_mat = 0.d0
|
|
||||||
do i = 1, elec_alpha_num
|
|
||||||
do j = elec_alpha_num+1, mo_num
|
|
||||||
grad_good_hermit_tc_fock_mat += dabs(good_hermit_tc_fock_mat(i,j))
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, grad_hermit_average_tc_fock_mat]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! grad_hermit_average_tc_fock_mat = norm of gradients of the upper triangular TC fock
|
|
||||||
END_DOC
|
|
||||||
integer :: i, j
|
|
||||||
grad_hermit_average_tc_fock_mat = 0.d0
|
|
||||||
do i = 1, elec_alpha_num
|
|
||||||
do j = elec_alpha_num+1, mo_num
|
|
||||||
grad_hermit_average_tc_fock_mat += dabs(hermit_average_tc_fock_mat(i,j))
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine save_good_hermit_tc_eigvectors()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: sign
|
|
||||||
character*(64) :: label
|
|
||||||
logical :: output
|
|
||||||
|
|
||||||
sign = 1
|
|
||||||
label = "Canonical"
|
|
||||||
output = .False.
|
|
||||||
|
|
||||||
if(symetric_fock_tc)then
|
|
||||||
call mo_as_eigvectors_of_mo_matrix(hermit_average_tc_fock_mat, mo_num, mo_num, label, sign, output)
|
|
||||||
else
|
|
||||||
call mo_as_eigvectors_of_mo_matrix(good_hermit_tc_fock_mat, mo_num, mo_num, label, sign, output)
|
|
||||||
endif
|
|
||||||
end subroutine save_good_hermit_tc_eigvectors
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -110,23 +110,14 @@ BEGIN_PROVIDER [double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num)]
|
|||||||
double precision :: t0, t1, tt0, tt1
|
double precision :: t0, t1, tt0, tt1
|
||||||
double precision, allocatable :: tmp(:,:)
|
double precision, allocatable :: tmp(:,:)
|
||||||
|
|
||||||
if(bi_ortho) then
|
PROVIDE mo_l_coef mo_r_coef
|
||||||
|
|
||||||
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) )
|
||||||
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
|
|
||||||
PROVIDE fock_3e_mo_a
|
|
||||||
Fock_matrix_tc_mo_alpha += fock_3e_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) )
|
|
||||||
|
|
||||||
|
if(three_body_h_tc) then
|
||||||
|
PROVIDE fock_3e_mo_a
|
||||||
|
Fock_matrix_tc_mo_alpha += fock_3e_mo_a
|
||||||
endif
|
endif
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
@ -142,21 +133,12 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
|
|||||||
implicit none
|
implicit none
|
||||||
double precision, allocatable :: tmp(:,:)
|
double precision, allocatable :: tmp(:,:)
|
||||||
|
|
||||||
if(bi_ortho) then
|
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
|
||||||
|
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
|
||||||
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
|
|
||||||
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
|
|
||||||
|
|
||||||
if(three_body_h_tc) then
|
|
||||||
PROVIDE fock_3e_mo_b
|
|
||||||
Fock_matrix_tc_mo_beta += fock_3e_mo_b
|
|
||||||
endif
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
|
|
||||||
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
|
|
||||||
|
|
||||||
|
if(three_body_h_tc) then
|
||||||
|
PROVIDE fock_3e_mo_b
|
||||||
|
Fock_matrix_tc_mo_beta += fock_3e_mo_b
|
||||||
endif
|
endif
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|||||||
@ -132,7 +132,7 @@
|
|||||||
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
|
||||||
@ -153,8 +153,21 @@
|
|||||||
enddo
|
enddo
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if(.not.bi_ortho .and. three_body_h_tc)then
|
if(tc_Brillouin_Right) then
|
||||||
Fock_matrix_tc_mo_tot += fock_3_mat
|
|
||||||
|
double precision, allocatable :: tmp(:,:)
|
||||||
|
allocate(tmp(mo_num,mo_num))
|
||||||
|
|
||||||
|
tmp = Fock_matrix_tc_mo_tot
|
||||||
|
do j = 1, mo_num
|
||||||
|
do i = 1, j-1
|
||||||
|
tmp(i,j) = Fock_matrix_tc_mo_tot(j,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
Fock_matrix_tc_mo_tot = tmp
|
||||||
|
deallocate(tmp)
|
||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|||||||
@ -1,771 +0,0 @@
|
|||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i,j
|
|
||||||
double precision :: contrib
|
|
||||||
|
|
||||||
fock_3_mat = 0.d0
|
|
||||||
if(.not.bi_ortho .and. three_body_h_tc) then
|
|
||||||
|
|
||||||
call give_fock_ia_three_e_total(1, 1, contrib)
|
|
||||||
!! !$OMP PARALLEL &
|
|
||||||
!! !$OMP DEFAULT (NONE) &
|
|
||||||
!! !$OMP PRIVATE (i,j,m,integral) &
|
|
||||||
!! !$OMP SHARED (mo_num,three_body_3_index)
|
|
||||||
!! !$OMP DO SCHEDULE (guided) COLLAPSE(3)
|
|
||||||
do i = 1, mo_num
|
|
||||||
do j = 1, mo_num
|
|
||||||
call give_fock_ia_three_e_total(j,i,contrib)
|
|
||||||
fock_3_mat(j,i) = -contrib
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!else if(bi_ortho.and.three_body_h_tc) then
|
|
||||||
!! !$OMP END DO
|
|
||||||
!! !$OMP END PARALLEL
|
|
||||||
!! do i = 1, mo_num
|
|
||||||
!! do j = 1, i-1
|
|
||||||
!! mat_three(j,i) = mat_three(i,j)
|
|
||||||
!! enddo
|
|
||||||
!! enddo
|
|
||||||
endif
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
subroutine give_fock_ia_three_e_total(i,a,contrib)
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! contrib is the TOTAL (same spins / opposite spins) contribution from the three body term to the Fock operator
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
integer, intent(in) :: i,a
|
|
||||||
double precision, intent(out) :: contrib
|
|
||||||
double precision :: int_1, int_2, int_3
|
|
||||||
double precision :: mos_i, mos_a, w_ia
|
|
||||||
double precision :: mos_ia, weight
|
|
||||||
|
|
||||||
integer :: mm, ipoint,k,l
|
|
||||||
|
|
||||||
int_1 = 0.d0
|
|
||||||
int_2 = 0.d0
|
|
||||||
int_3 = 0.d0
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
mos_i = mos_in_r_array_transp(ipoint,i)
|
|
||||||
mos_a = mos_in_r_array_transp(ipoint,a)
|
|
||||||
mos_ia = mos_a * mos_i
|
|
||||||
w_ia = x_W_ij_erf_rk(ipoint,mm,i,a)
|
|
||||||
|
|
||||||
int_1 += weight * fock_3_w_kk_sum(ipoint,mm) * (4.d0 * fock_3_rho_beta(ipoint) * w_ia &
|
|
||||||
+ 2.0d0 * mos_ia * fock_3_w_kk_sum(ipoint,mm) &
|
|
||||||
- 2.0d0 * fock_3_w_ki_mos_k(ipoint,mm,i) * mos_a &
|
|
||||||
- 2.0d0 * fock_3_w_ki_mos_k(ipoint,mm,a) * mos_i )
|
|
||||||
int_2 += weight * (-1.d0) * ( 2.0d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * w_ia &
|
|
||||||
+ 2.0d0 * fock_3_rho_beta(ipoint) * fock_3_w_ki_wk_a(ipoint,mm,i,a) &
|
|
||||||
+ 1.0d0 * mos_ia * fock_3_trace_w_tilde(ipoint,mm) )
|
|
||||||
|
|
||||||
int_3 += weight * 1.d0 * (fock_3_w_kl_wla_phi_k(ipoint,mm,i) * mos_a + fock_3_w_kl_wla_phi_k(ipoint,mm,a) * mos_i &
|
|
||||||
+fock_3_w_ki_mos_k(ipoint,mm,i) * fock_3_w_ki_mos_k(ipoint,mm,a) )
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
contrib = int_1 + int_2 + int_3
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, diag_three_elem_hf]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, k, ipoint, mm
|
|
||||||
double precision :: contrib, weight, four_third, one_third, two_third, exchange_int_231
|
|
||||||
double precision :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb
|
|
||||||
double precision, allocatable :: tmp(:)
|
|
||||||
double precision, allocatable :: tmp_L(:,:), tmp_R(:,:)
|
|
||||||
double precision, allocatable :: tmp_M(:,:), tmp_S(:), tmp_O(:), tmp_J(:,:)
|
|
||||||
double precision, allocatable :: tmp_M_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:)
|
|
||||||
|
|
||||||
PROVIDE mo_l_coef mo_r_coef
|
|
||||||
|
|
||||||
!print *, ' providing diag_three_elem_hf'
|
|
||||||
|
|
||||||
if(.not. three_body_h_tc) then
|
|
||||||
|
|
||||||
if(noL_standard) then
|
|
||||||
PROVIDE noL_0e
|
|
||||||
diag_three_elem_hf = noL_0e
|
|
||||||
else
|
|
||||||
diag_three_elem_hf = 0.d0
|
|
||||||
endif
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
if(.not. bi_ortho) then
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
one_third = 1.d0/3.d0
|
|
||||||
two_third = 2.d0/3.d0
|
|
||||||
four_third = 4.d0/3.d0
|
|
||||||
diag_three_elem_hf = 0.d0
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
call give_integrals_3_body(k, j, i, j, i, k, exchange_int_231)
|
|
||||||
diag_three_elem_hf += two_third * exchange_int_231
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
contrib = 3.d0 * fock_3_w_kk_sum(ipoint,mm) * fock_3_rho_beta(ipoint) * fock_3_w_kk_sum(ipoint,mm) &
|
|
||||||
- 2.d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * fock_3_w_kk_sum(ipoint,mm) &
|
|
||||||
- 1.d0 * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm)
|
|
||||||
contrib *= four_third
|
|
||||||
contrib += -two_third * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) &
|
|
||||||
-four_third * fock_3_w_kk_sum(ipoint,mm) * fock_3_w_kl_mo_k_mo_l(ipoint,mm)
|
|
||||||
diag_three_elem_hf += weight * contrib
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
diag_three_elem_hf = - diag_three_elem_hf
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
! ------------
|
|
||||||
! SLOW VERSION
|
|
||||||
! ------------
|
|
||||||
|
|
||||||
!call give_aaa_contrib(integral_aaa)
|
|
||||||
!call give_aab_contrib(integral_aab)
|
|
||||||
!call give_abb_contrib(integral_abb)
|
|
||||||
!call give_bbb_contrib(integral_bbb)
|
|
||||||
!diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb
|
|
||||||
|
|
||||||
! ------------
|
|
||||||
! ------------
|
|
||||||
|
|
||||||
PROVIDE int2_grad1_u12_bimo_t
|
|
||||||
PROVIDE mos_l_in_r_array_transp
|
|
||||||
PROVIDE mos_r_in_r_array_transp
|
|
||||||
|
|
||||||
if(elec_alpha_num .eq. elec_beta_num) then
|
|
||||||
|
|
||||||
allocate(tmp(elec_beta_num))
|
|
||||||
allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
|
||||||
!$OMP DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R) &
|
|
||||||
!$OMP SHARED(elec_beta_num, n_points_final_grid, &
|
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
|
||||||
!$OMP int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
|
|
||||||
tmp_L = 0.d0
|
|
||||||
tmp_R = 0.d0
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
|
|
||||||
tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
tmp(j) = 0.d0
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
|
||||||
enddo
|
|
||||||
enddo ! j
|
|
||||||
!$OMP END DO
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
diag_three_elem_hf = -2.d0 * sum(tmp)
|
|
||||||
|
|
||||||
deallocate(tmp)
|
|
||||||
deallocate(tmp_L, tmp_R)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
|
|
||||||
tmp_O = 0.d0
|
|
||||||
tmp_J = 0.d0
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
|
||||||
!$OMP DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv) &
|
|
||||||
!$OMP SHARED(elec_beta_num, n_points_final_grid, &
|
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
|
||||||
!$OMP int2_grad1_u12_bimo_t, tmp_O, tmp_J)
|
|
||||||
|
|
||||||
allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
|
|
||||||
tmp_O_priv = 0.d0
|
|
||||||
tmp_J_priv = 0.d0
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp_O_priv(ipoint) = tmp_O_priv(ipoint) + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
|
|
||||||
tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
|
|
||||||
tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP CRITICAL
|
|
||||||
tmp_O = tmp_O + tmp_O_priv
|
|
||||||
tmp_J = tmp_J + tmp_J_priv
|
|
||||||
!$OMP END CRITICAL
|
|
||||||
|
|
||||||
deallocate(tmp_O_priv, tmp_J_priv)
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
|
|
||||||
tmp_M = 0.d0
|
|
||||||
tmp_S = 0.d0
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
|
||||||
!$OMP DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv) &
|
|
||||||
!$OMP SHARED(elec_beta_num, n_points_final_grid, &
|
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
|
||||||
!$OMP int2_grad1_u12_bimo_t, tmp_M, tmp_S)
|
|
||||||
|
|
||||||
allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
|
|
||||||
tmp_M_priv = 0.d0
|
|
||||||
tmp_S_priv = 0.d0
|
|
||||||
|
|
||||||
!$OMP DO COLLAPSE(2)
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
|
|
||||||
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
|
||||||
+ int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
|
||||||
+ int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP CRITICAL
|
|
||||||
tmp_M = tmp_M + tmp_M_priv
|
|
||||||
tmp_S = tmp_S + tmp_S_priv
|
|
||||||
!$OMP END CRITICAL
|
|
||||||
|
|
||||||
deallocate(tmp_M_priv, tmp_S_priv)
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
allocate(tmp(n_points_final_grid))
|
|
||||||
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
|
|
||||||
|
|
||||||
tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint) &
|
|
||||||
- 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
|
|
||||||
+ tmp_J(ipoint,2) * tmp_M(ipoint,2) &
|
|
||||||
+ tmp_J(ipoint,3) * tmp_M(ipoint,3)))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
diag_three_elem_hf = diag_three_elem_hf -2.d0 * (sum(tmp))
|
|
||||||
|
|
||||||
deallocate(tmp)
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
allocate(tmp(elec_alpha_num))
|
|
||||||
allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
|
||||||
!$OMP DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R) &
|
|
||||||
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
|
||||||
!$OMP int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
|
|
||||||
tmp_L = 0.d0
|
|
||||||
tmp_R = 0.d0
|
|
||||||
do i = elec_beta_num+1, elec_alpha_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_L(ipoint,1) = tmp_L(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,2) = tmp_L(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,3) = tmp_L(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
|
|
||||||
tmp_R(ipoint,1) = tmp_R(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,2) = tmp_R(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,3) = tmp_R(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
tmp(j) = 0.d0
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
|
|
||||||
tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
|
||||||
enddo
|
|
||||||
enddo ! j
|
|
||||||
!$OMP END DO
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
|
||||||
!$OMP DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R) &
|
|
||||||
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
|
||||||
!$OMP int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do j = elec_beta_num+1, elec_alpha_num
|
|
||||||
|
|
||||||
tmp_L = 0.d0
|
|
||||||
tmp_R = 0.d0
|
|
||||||
do i = 1, elec_alpha_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
|
||||||
|
|
||||||
tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
tmp(j) = 0.d0
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp(j) = tmp(j) + 0.5d0 * final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
|
||||||
enddo
|
|
||||||
enddo ! j
|
|
||||||
!$OMP END DO
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
diag_three_elem_hf = -2.d0 * sum(tmp)
|
|
||||||
|
|
||||||
deallocate(tmp)
|
|
||||||
deallocate(tmp_L, tmp_R)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
|
|
||||||
tmp_O = 0.d0
|
|
||||||
tmp_J = 0.d0
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
|
||||||
!$OMP DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv) &
|
|
||||||
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
|
||||||
!$OMP int2_grad1_u12_bimo_t, tmp_O, tmp_J)
|
|
||||||
|
|
||||||
allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
|
|
||||||
tmp_O_priv = 0.d0
|
|
||||||
tmp_J_priv = 0.d0
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp_O_priv(ipoint) = tmp_O_priv(ipoint) + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
|
|
||||||
tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
|
|
||||||
tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do i = elec_beta_num+1, elec_alpha_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp_O_priv(ipoint) = tmp_O_priv(ipoint) + 0.5d0 * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,i)
|
|
||||||
tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,i)
|
|
||||||
tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP CRITICAL
|
|
||||||
tmp_O = tmp_O + tmp_O_priv
|
|
||||||
tmp_J = tmp_J + tmp_J_priv
|
|
||||||
!$OMP END CRITICAL
|
|
||||||
|
|
||||||
deallocate(tmp_O_priv, tmp_J_priv)
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
|
|
||||||
tmp_M = 0.d0
|
|
||||||
tmp_S = 0.d0
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
|
||||||
!$OMP DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv) &
|
|
||||||
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
|
||||||
!$OMP int2_grad1_u12_bimo_t, tmp_M, tmp_S)
|
|
||||||
|
|
||||||
allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
|
|
||||||
tmp_M_priv = 0.d0
|
|
||||||
tmp_S_priv = 0.d0
|
|
||||||
|
|
||||||
!$OMP DO COLLAPSE(2)
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
|
|
||||||
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
|
||||||
+ int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
|
||||||
+ int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP DO COLLAPSE(2)
|
|
||||||
do i = elec_beta_num+1, elec_alpha_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
|
|
||||||
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
|
|
||||||
|
|
||||||
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
|
||||||
+ int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
|
||||||
+ int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP DO COLLAPSE(2)
|
|
||||||
do i = elec_beta_num+1, elec_alpha_num
|
|
||||||
do j = elec_beta_num+1, elec_alpha_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
|
|
||||||
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
|
||||||
+ 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
|
||||||
+ 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP CRITICAL
|
|
||||||
tmp_M = tmp_M + tmp_M_priv
|
|
||||||
tmp_S = tmp_S + tmp_S_priv
|
|
||||||
!$OMP END CRITICAL
|
|
||||||
|
|
||||||
deallocate(tmp_M_priv, tmp_S_priv)
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
allocate(tmp(n_points_final_grid))
|
|
||||||
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
|
|
||||||
tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
|
|
||||||
|
|
||||||
tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint) &
|
|
||||||
- 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
|
|
||||||
+ tmp_J(ipoint,2) * tmp_M(ipoint,2) &
|
|
||||||
+ tmp_J(ipoint,3) * tmp_M(ipoint,3)))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
diag_three_elem_hf = diag_three_elem_hf - 2.d0 * (sum(tmp))
|
|
||||||
|
|
||||||
deallocate(tmp)
|
|
||||||
|
|
||||||
endif
|
|
||||||
|
|
||||||
|
|
||||||
endif
|
|
||||||
|
|
||||||
endif
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_mat_a_op_sh, (mo_num, mo_num)]
|
|
||||||
implicit none
|
|
||||||
integer :: h,p,i,j
|
|
||||||
double precision :: direct_int, exch_int, exchange_int_231, exchange_int_312
|
|
||||||
double precision :: exchange_int_23, exchange_int_12, exchange_int_13
|
|
||||||
|
|
||||||
fock_3_mat_a_op_sh = 0.d0
|
|
||||||
do h = 1, mo_num
|
|
||||||
do p = 1, mo_num
|
|
||||||
!F_a^{ab}(h,p)
|
|
||||||
do i = 1, elec_beta_num ! beta
|
|
||||||
do j = elec_beta_num+1, elec_alpha_num ! alpha
|
|
||||||
call give_integrals_3_body(h,j,i,p,j,i,direct_int) ! <hji|pji>
|
|
||||||
call give_integrals_3_body(h,j,i,j,p,i,exch_int)
|
|
||||||
fock_3_mat_a_op_sh(h,p) -= direct_int - exch_int
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!F_a^{aa}(h,p)
|
|
||||||
do i = 1, elec_beta_num ! alpha
|
|
||||||
do j = elec_beta_num+1, elec_alpha_num ! alpha
|
|
||||||
call give_integrals_3_body(h,j,i,p,j,i,direct_int)
|
|
||||||
call give_integrals_3_body(h,j,i,i,p,j,exchange_int_231)
|
|
||||||
call give_integrals_3_body(h,j,i,j,i,p,exchange_int_312)
|
|
||||||
call give_integrals_3_body(h,j,i,p,i,j,exchange_int_23)
|
|
||||||
call give_integrals_3_body(h,j,i,i,j,p,exchange_int_12)
|
|
||||||
call give_integrals_3_body(h,j,i,j,p,i,exchange_int_13)
|
|
||||||
fock_3_mat_a_op_sh(h,p) -= ( direct_int + exchange_int_231 + exchange_int_312 &
|
|
||||||
- exchange_int_23 & ! i <-> j
|
|
||||||
- exchange_int_12 & ! p <-> j
|
|
||||||
- exchange_int_13 )! p <-> i
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
! symmetrized
|
|
||||||
! do p = 1, elec_beta_num
|
|
||||||
! do h = elec_alpha_num +1, mo_num
|
|
||||||
! fock_3_mat_a_op_sh(h,p) = fock_3_mat_a_op_sh(p,h)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
|
|
||||||
! do h = elec_beta_num+1, elec_alpha_num
|
|
||||||
! do p = elec_alpha_num +1, mo_num
|
|
||||||
! !F_a^{bb}(h,p)
|
|
||||||
! do i = 1, elec_beta_num
|
|
||||||
! do j = i+1, elec_beta_num
|
|
||||||
! call give_integrals_3_body(h,j,i,p,j,i,direct_int)
|
|
||||||
! call give_integrals_3_body(h,j,i,p,i,j,exch_int)
|
|
||||||
! fock_3_mat_a_op_sh(h,p) -= direct_int - exch_int
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_mat_b_op_sh, (mo_num, mo_num)]
|
|
||||||
implicit none
|
|
||||||
integer :: h,p,i,j
|
|
||||||
double precision :: direct_int, exch_int
|
|
||||||
fock_3_mat_b_op_sh = 0.d0
|
|
||||||
do h = 1, elec_beta_num
|
|
||||||
do p = elec_alpha_num +1, mo_num
|
|
||||||
!F_b^{aa}(h,p)
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = elec_beta_num+1, elec_alpha_num
|
|
||||||
call give_integrals_3_body(h,j,i,p,j,i,direct_int)
|
|
||||||
call give_integrals_3_body(h,j,i,p,i,j,exch_int)
|
|
||||||
fock_3_mat_b_op_sh(h,p) += direct_int - exch_int
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!F_b^{ab}(h,p)
|
|
||||||
do i = elec_beta_num+1, elec_beta_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
call give_integrals_3_body(h,j,i,p,j,i,direct_int)
|
|
||||||
call give_integrals_3_body(h,j,i,j,p,i,exch_int)
|
|
||||||
fock_3_mat_b_op_sh(h,p) += direct_int - exch_int
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_w_kk_sum, (n_points_final_grid,3)]
|
|
||||||
implicit none
|
|
||||||
integer :: mm, ipoint,k
|
|
||||||
double precision :: w_kk
|
|
||||||
fock_3_w_kk_sum = 0.d0
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
w_kk = x_W_ij_erf_rk(ipoint,mm,k,k)
|
|
||||||
fock_3_w_kk_sum(ipoint,mm) += w_kk
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_w_ki_mos_k, (n_points_final_grid,3,mo_num)]
|
|
||||||
implicit none
|
|
||||||
integer :: mm, ipoint,k,i
|
|
||||||
double precision :: w_ki, mo_k
|
|
||||||
fock_3_w_ki_mos_k = 0.d0
|
|
||||||
do i = 1, mo_num
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
w_ki = x_W_ij_erf_rk(ipoint,mm,k,i)
|
|
||||||
mo_k = mos_in_r_array(k,ipoint)
|
|
||||||
fock_3_w_ki_mos_k(ipoint,mm,i) += w_ki * mo_k
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_w_kl_w_kl, (n_points_final_grid,3)]
|
|
||||||
implicit none
|
|
||||||
integer :: k,j,ipoint,mm
|
|
||||||
double precision :: w_kj
|
|
||||||
fock_3_w_kl_w_kl = 0.d0
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
w_kj = x_W_ij_erf_rk(ipoint,mm,k,j)
|
|
||||||
fock_3_w_kl_w_kl(ipoint,mm) += w_kj * w_kj
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_rho_beta, (n_points_final_grid)]
|
|
||||||
implicit none
|
|
||||||
integer :: ipoint,k
|
|
||||||
fock_3_rho_beta = 0.d0
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
fock_3_rho_beta(ipoint) += mos_in_r_array(k,ipoint) * mos_in_r_array(k,ipoint)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_w_kl_mo_k_mo_l, (n_points_final_grid,3)]
|
|
||||||
implicit none
|
|
||||||
integer :: ipoint,k,l,mm
|
|
||||||
double precision :: mos_k, mos_l, w_kl
|
|
||||||
fock_3_w_kl_mo_k_mo_l = 0.d0
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
do l = 1, elec_beta_num
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
mos_k = mos_in_r_array_transp(ipoint,k)
|
|
||||||
mos_l = mos_in_r_array_transp(ipoint,l)
|
|
||||||
w_kl = x_W_ij_erf_rk(ipoint,mm,l,k)
|
|
||||||
fock_3_w_kl_mo_k_mo_l(ipoint,mm) += w_kl * mos_k * mos_l
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_w_ki_wk_a, (n_points_final_grid,3,mo_num, mo_num)]
|
|
||||||
implicit none
|
|
||||||
integer :: ipoint,i,a,k,mm
|
|
||||||
double precision :: w_ki,w_ka
|
|
||||||
fock_3_w_ki_wk_a = 0.d0
|
|
||||||
do i = 1, mo_num
|
|
||||||
do a = 1, mo_num
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
w_ki = x_W_ij_erf_rk(ipoint,mm,k,i)
|
|
||||||
w_ka = x_W_ij_erf_rk(ipoint,mm,k,a)
|
|
||||||
fock_3_w_ki_wk_a(ipoint,mm,a,i) += w_ki * w_ka
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_trace_w_tilde, (n_points_final_grid,3)]
|
|
||||||
implicit none
|
|
||||||
integer :: ipoint,k,mm
|
|
||||||
fock_3_trace_w_tilde = 0.d0
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
fock_3_trace_w_tilde(ipoint,mm) += fock_3_w_ki_wk_a(ipoint,mm,k,k)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, fock_3_w_kl_wla_phi_k, (n_points_final_grid,3,mo_num)]
|
|
||||||
implicit none
|
|
||||||
integer :: ipoint,a,k,mm,l
|
|
||||||
double precision :: w_kl,w_la, mo_k
|
|
||||||
fock_3_w_kl_wla_phi_k = 0.d0
|
|
||||||
do a = 1, mo_num
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
do l = 1, elec_beta_num
|
|
||||||
do mm = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
w_kl = x_W_ij_erf_rk(ipoint,mm,l,k)
|
|
||||||
w_la = x_W_ij_erf_rk(ipoint,mm,l,a)
|
|
||||||
mo_k = mos_in_r_array_transp(ipoint,k)
|
|
||||||
fock_3_w_kl_wla_phi_k(ipoint,mm,a) += w_kl * w_la * mo_k
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -1,391 +0,0 @@
|
|||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, tc_scf_dm_in_r, (n_points_final_grid) ]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j
|
|
||||||
|
|
||||||
tc_scf_dm_in_r = 0.d0
|
|
||||||
do i = 1, n_points_final_grid
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
tc_scf_dm_in_r(i) += mos_r_in_r_array(j,i) * mos_l_in_r_array(j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, w_sum_in_r, (n_points_final_grid, 3)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: ipoint, j, xi
|
|
||||||
|
|
||||||
w_sum_in_r = 0.d0
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
!w_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,j)
|
|
||||||
w_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, ww_sum_in_r, (n_points_final_grid, 3)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: ipoint, j, xi
|
|
||||||
double precision :: tmp
|
|
||||||
|
|
||||||
ww_sum_in_r = 0.d0
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
tmp = x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
|
|
||||||
ww_sum_in_r(ipoint,xi) += tmp * tmp
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, W1_r_in_r, (n_points_final_grid, 3, mo_num)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, xi, ipoint
|
|
||||||
|
|
||||||
! TODO: call lapack
|
|
||||||
|
|
||||||
W1_r_in_r = 0.d0
|
|
||||||
do i = 1, mo_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
W1_r_in_r(ipoint,xi,i) += mos_r_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, W1_l_in_r, (n_points_final_grid, 3, mo_num)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, xi, ipoint
|
|
||||||
|
|
||||||
! TODO: call lapack
|
|
||||||
|
|
||||||
W1_l_in_r = 0.d0
|
|
||||||
do i = 1, mo_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
W1_l_in_r(ipoint,xi,i) += mos_l_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, W1_in_r, (n_points_final_grid, 3)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: j, xi, ipoint
|
|
||||||
|
|
||||||
! TODO: call lapack
|
|
||||||
|
|
||||||
W1_in_r = 0.d0
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
W1_in_r(ipoint,xi) += W1_l_in_r(ipoint,xi,j) * mos_r_in_r_array_transp(ipoint,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, W1_diag_in_r, (n_points_final_grid, 3)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: j, xi, ipoint
|
|
||||||
|
|
||||||
! TODO: call lapack
|
|
||||||
|
|
||||||
W1_diag_in_r = 0.d0
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
W1_diag_in_r(ipoint,xi) += mos_r_in_r_array_transp(ipoint,j) * mos_l_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, v_sum_in_r, (n_points_final_grid, 3)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, xi, ipoint
|
|
||||||
|
|
||||||
! TODO: call lapack
|
|
||||||
v_sum_in_r = 0.d0
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
v_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, W1_W1_r_in_r, (n_points_final_grid, 3, mo_num)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, m, xi, ipoint
|
|
||||||
|
|
||||||
! TODO: call lapack
|
|
||||||
|
|
||||||
W1_W1_r_in_r = 0.d0
|
|
||||||
do i = 1, mo_num
|
|
||||||
do m = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
W1_W1_r_in_r(ipoint,xi,i) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,m,i) * W1_r_in_r(ipoint,xi,m)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, W1_W1_l_in_r, (n_points_final_grid, 3, mo_num)]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, xi, ipoint
|
|
||||||
|
|
||||||
! TODO: call lapack
|
|
||||||
|
|
||||||
W1_W1_l_in_r = 0.d0
|
|
||||||
do i = 1, mo_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
W1_W1_l_in_r(ipoint,xi,i) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j) * W1_l_in_r(ipoint,xi,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine direct_term_imj_bi_ortho(a, i, integral)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! computes sum_(j,m = 1, elec_beta_num) < a m j | i m j > with bi ortho mos
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i, a
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
|
|
||||||
integer :: ipoint, xi
|
|
||||||
double precision :: weight, tmp
|
|
||||||
|
|
||||||
integral = 0.d0
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
!integral += ( mos_l_in_r_array(a,ipoint) * mos_r_in_r_array(i,ipoint) * w_sum_in_r(ipoint,xi) * w_sum_in_r(ipoint,xi) &
|
|
||||||
! + 2.d0 * tc_scf_dm_in_r(ipoint) * w_sum_in_r(ipoint,xi) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) ) * weight
|
|
||||||
|
|
||||||
tmp = w_sum_in_r(ipoint,xi)
|
|
||||||
|
|
||||||
integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * tmp * tmp &
|
|
||||||
+ 2.d0 * tc_scf_dm_in_r(ipoint) * tmp * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) &
|
|
||||||
) * weight
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine exch_term_jmi_bi_ortho(a, i, integral)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! computes sum_(j,m = 1, elec_beta_num) < a m j | j m i > with bi ortho mos
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i, a
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
|
|
||||||
integer :: ipoint, xi, j
|
|
||||||
double precision :: weight, tmp
|
|
||||||
|
|
||||||
integral = 0.d0
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
|
|
||||||
tmp = 0.d0
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
tmp = tmp + x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_r_in_r(ipoint,xi,i) * w_sum_in_r(ipoint,xi) &
|
|
||||||
+ tc_scf_dm_in_r(ipoint) * tmp &
|
|
||||||
+ mos_r_in_r_array_transp(ipoint,i) * W1_l_in_r(ipoint,xi,a) * w_sum_in_r(ipoint,xi) &
|
|
||||||
) * weight
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine exch_term_ijm_bi_ortho(a, i, integral)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! computes sum_(j,m = 1, elec_beta_num) < a m j | i j m > with bi ortho mos
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i, a
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
|
|
||||||
integer :: ipoint, xi
|
|
||||||
double precision :: weight
|
|
||||||
|
|
||||||
integral = 0.d0
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
|
|
||||||
integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * v_sum_in_r(ipoint,xi) &
|
|
||||||
+ 2.d0 * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) * W1_in_r(ipoint,xi) &
|
|
||||||
) * weight
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine direct_term_ijj_bi_ortho(a, i, integral)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! computes sum_(j = 1, elec_beta_num) < a j j | i j j > with bi ortho mos
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i, a
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
|
|
||||||
integer :: ipoint, xi
|
|
||||||
double precision :: weight
|
|
||||||
|
|
||||||
integral = 0.d0
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
|
|
||||||
integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * ww_sum_in_r(ipoint,xi) &
|
|
||||||
+ 2.d0 * W1_diag_in_r(ipoint, xi) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) &
|
|
||||||
) * weight
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine cyclic_term_jim_bi_ortho(a, i, integral)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! computes sum_(j,m = 1, elec_beta_num) < a m j | j i m > with bi ortho mos
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i, a
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
|
|
||||||
integer :: ipoint, xi
|
|
||||||
double precision :: weight
|
|
||||||
|
|
||||||
integral = 0.d0
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
|
|
||||||
integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_W1_r_in_r(ipoint,xi,i) &
|
|
||||||
+ W1_W1_l_in_r(ipoint,xi,a) * mos_r_in_r_array_transp(ipoint,i) &
|
|
||||||
+ W1_l_in_r(ipoint,xi,a) * W1_r_in_r(ipoint,xi,i) &
|
|
||||||
) * weight
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine cyclic_term_mji_bi_ortho(a, i, integral)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! computes sum_(j,m = 1, elec_beta_num) < a m j | m j i > with bi ortho mos
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i, a
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
|
|
||||||
integer :: ipoint, xi
|
|
||||||
double precision :: weight
|
|
||||||
|
|
||||||
integral = 0.d0
|
|
||||||
do xi = 1, 3
|
|
||||||
do ipoint = 1, n_points_final_grid
|
|
||||||
weight = final_weight_at_r_vector(ipoint)
|
|
||||||
|
|
||||||
integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_W1_r_in_r(ipoint,xi,i) &
|
|
||||||
+ W1_l_in_r(ipoint,xi,a) * W1_r_in_r(ipoint,xi,i) &
|
|
||||||
+ W1_W1_l_in_r(ipoint,xi,a) * mos_r_in_r_array_transp(ipoint,i) &
|
|
||||||
) * weight
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -1,318 +0,0 @@
|
|||||||
BEGIN_PROVIDER [integer , m_max_sm_7]
|
|
||||||
&BEGIN_PROVIDER [integer , n_max_sm_7]
|
|
||||||
&BEGIN_PROVIDER [integer , o_max_sm_7]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! maximum value of the "m", "n" and "o" integer in the Jastrow function as in Eq. (4)
|
|
||||||
! of Schmidt,Moskowitz, JCP, 93, 4172 (1990) for the SM_7 version of Table IV
|
|
||||||
END_DOC
|
|
||||||
m_max_sm_7 = 4
|
|
||||||
n_max_sm_7 = 0
|
|
||||||
o_max_sm_7 = 4
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [integer , m_max_sm_9]
|
|
||||||
&BEGIN_PROVIDER [integer , n_max_sm_9]
|
|
||||||
&BEGIN_PROVIDER [integer , o_max_sm_9]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! maximum value of the "m", "n" and "o" integer in the Jastrow function as in Eq. (4)
|
|
||||||
! of Schmidt,Moskowitz, JCP, 93, 4172 (1990) for the SM_9 version of Table IV
|
|
||||||
END_DOC
|
|
||||||
m_max_sm_9 = 4
|
|
||||||
n_max_sm_9 = 2
|
|
||||||
o_max_sm_9 = 4
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [integer , m_max_sm_17]
|
|
||||||
&BEGIN_PROVIDER [integer , n_max_sm_17]
|
|
||||||
&BEGIN_PROVIDER [integer , o_max_sm_17]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! maximum value of the "m", "n" and "o" integer in the Jastrow function as in Eq. (4)
|
|
||||||
! of Schmidt,Moskowitz, JCP, 93, 4172 (1990) for the SM_17 version of Table IV
|
|
||||||
END_DOC
|
|
||||||
m_max_sm_17 = 6
|
|
||||||
n_max_sm_17 = 2
|
|
||||||
o_max_sm_17 = 6
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, c_mn_o_sm_7, (0:m_max_sm_7,0:n_max_sm_7,0:o_max_sm_7,2:10)]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
!
|
|
||||||
!c_mn_o_7(0:4,0:4,2:10) = coefficient for the SM_7 correlation factor as given is Table IV of
|
|
||||||
! Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
! the first index (0:4) is the "m" integer for the 1e part
|
|
||||||
! the second index(0:0) is the "n" integer for the 1e part WHICH IS ALWAYS SET TO 0 FOR SM_7
|
|
||||||
! the third index (0:4) is the "o" integer for the 2e part
|
|
||||||
! the fourth index (2:10) is the nuclear charge of the atom
|
|
||||||
END_DOC
|
|
||||||
c_mn_o_sm_7 = 0.d0
|
|
||||||
integer :: i
|
|
||||||
do i = 2, 10 ! loop over nuclear charge
|
|
||||||
c_mn_o_sm_7(0,0,1,i) = 0.5d0 ! all the linear terms are set to 1/2 to satisfy the anti-parallel spin condition
|
|
||||||
enddo
|
|
||||||
! He atom
|
|
||||||
! two electron terms
|
|
||||||
c_mn_o_sm_7(0,0,2,2) = 0.50516d0
|
|
||||||
c_mn_o_sm_7(0,0,3,2) = -0.19313d0
|
|
||||||
c_mn_o_sm_7(0,0,4,2) = 0.30276d0
|
|
||||||
! one-electron terms
|
|
||||||
c_mn_o_sm_7(2,0,0,2) = -0.16995d0
|
|
||||||
c_mn_o_sm_7(3,0,0,2) = -0.34505d0
|
|
||||||
c_mn_o_sm_7(4,0,0,2) = -0.54777d0
|
|
||||||
! Ne atom
|
|
||||||
! two electron terms
|
|
||||||
c_mn_o_sm_7(0,0,2,10) = -0.792d0
|
|
||||||
c_mn_o_sm_7(0,0,3,10) = 1.05232d0
|
|
||||||
c_mn_o_sm_7(0,0,4,10) = -0.65615d0
|
|
||||||
! one-electron terms
|
|
||||||
c_mn_o_sm_7(2,0,0,10) = -0.13312d0
|
|
||||||
c_mn_o_sm_7(3,0,0,10) = -0.00131d0
|
|
||||||
c_mn_o_sm_7(4,0,0,10) = 0.09083d0
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, c_mn_o_sm_9, (0:m_max_sm_9,0:n_max_sm_9,0:o_max_sm_9,2:10)]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
!
|
|
||||||
!c_mn_o_9(0:4,0:4,2:10) = coefficient for the SM_9 correlation factor as given is Table IV of
|
|
||||||
! Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
! the first index (0:4) is the "m" integer for the 1e part
|
|
||||||
! the second index(0:0) is the "n" integer for the 1e part WHICH IS ALWAYS SET TO 0 FOR SM_9
|
|
||||||
! the third index (0:4) is the "o" integer for the 2e part
|
|
||||||
! the fourth index (2:10) is the nuclear charge of the atom
|
|
||||||
END_DOC
|
|
||||||
c_mn_o_sm_9 = 0.d0
|
|
||||||
integer :: i
|
|
||||||
do i = 2, 10 ! loop over nuclear charge
|
|
||||||
c_mn_o_sm_9(0,0,1,i) = 0.5d0 ! all the linear terms are set to 1/2 to satisfy the anti-parallel spin condition
|
|
||||||
enddo
|
|
||||||
! He atom
|
|
||||||
! two electron terms
|
|
||||||
c_mn_o_sm_9(0,0,2,2) = 0.50516d0
|
|
||||||
c_mn_o_sm_9(0,0,3,2) = -0.19313d0
|
|
||||||
c_mn_o_sm_9(0,0,4,2) = 0.30276d0
|
|
||||||
! one-electron terms
|
|
||||||
c_mn_o_sm_9(2,0,0,2) = -0.16995d0
|
|
||||||
c_mn_o_sm_9(3,0,0,2) = -0.34505d0
|
|
||||||
c_mn_o_sm_9(4,0,0,2) = -0.54777d0
|
|
||||||
! Ne atom
|
|
||||||
! two electron terms
|
|
||||||
c_mn_o_sm_9(0,0,2,10) = -0.792d0
|
|
||||||
c_mn_o_sm_9(0,0,3,10) = 1.05232d0
|
|
||||||
c_mn_o_sm_9(0,0,4,10) = -0.65615d0
|
|
||||||
! one-electron terms
|
|
||||||
c_mn_o_sm_9(2,0,0,10) = -0.13312d0
|
|
||||||
c_mn_o_sm_9(3,0,0,10) = -0.00131d0
|
|
||||||
c_mn_o_sm_9(4,0,0,10) = 0.09083d0
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, c_mn_o_sm_17, (0:m_max_sm_17,0:n_max_sm_17,0:o_max_sm_17,2:10)]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
!
|
|
||||||
!c_mn_o_17(0:4,0:4,2:10) = coefficient for the SM_17 correlation factor as given is Table IV of
|
|
||||||
! Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
! the first index (0:4) is the "m" integer for the 1e part
|
|
||||||
! the second index(0:0) is the "n" integer for the 1e part WHICH IS ALWAYS SET TO 0 FOR SM_17
|
|
||||||
! the third index (0:4) is the "o" integer for the 2e part
|
|
||||||
! the fourth index (2:10) is the nuclear charge of the atom
|
|
||||||
END_DOC
|
|
||||||
c_mn_o_sm_17 = 0.d0
|
|
||||||
integer :: i
|
|
||||||
do i = 2, 10 ! loop over nuclear charge
|
|
||||||
c_mn_o_sm_17(0,0,1,i) = 0.5d0 ! all the linear terms are set to 1/2 to satisfy the anti-parallel spin condition
|
|
||||||
enddo
|
|
||||||
! He atom
|
|
||||||
! two electron terms
|
|
||||||
c_mn_o_sm_17(0,0,2,2) = 0.09239d0
|
|
||||||
c_mn_o_sm_17(0,0,3,2) = -0.38664d0
|
|
||||||
c_mn_o_sm_17(0,0,4,2) = 0.95764d0
|
|
||||||
! one-electron terms
|
|
||||||
c_mn_o_sm_17(2,0,0,2) = 0.23208d0
|
|
||||||
c_mn_o_sm_17(3,0,0,2) = -0.45032d0
|
|
||||||
c_mn_o_sm_17(4,0,0,2) = 0.82777d0
|
|
||||||
c_mn_o_sm_17(2,2,0,2) = -4.15388d0
|
|
||||||
! ee-n terms
|
|
||||||
c_mn_o_sm_17(2,0,2,2) = 0.80622d0
|
|
||||||
c_mn_o_sm_17(2,2,2,2) = 10.19704d0
|
|
||||||
c_mn_o_sm_17(4,0,2,2) = -4.96259d0
|
|
||||||
c_mn_o_sm_17(2,0,4,2) = -1.35647d0
|
|
||||||
c_mn_o_sm_17(4,2,2,2) = -5.90907d0
|
|
||||||
c_mn_o_sm_17(6,0,2,2) = 0.90343d0
|
|
||||||
c_mn_o_sm_17(4,0,4,2) = 5.50739d0
|
|
||||||
c_mn_o_sm_17(2,2,4,2) = -0.03154d0
|
|
||||||
c_mn_o_sm_17(2,0,6,2) = -1.1051860
|
|
||||||
|
|
||||||
|
|
||||||
! Ne atom
|
|
||||||
! two electron terms
|
|
||||||
c_mn_o_sm_17(0,0,2,10) = -0.80909d0
|
|
||||||
c_mn_o_sm_17(0,0,3,10) = -0.00219d0
|
|
||||||
c_mn_o_sm_17(0,0,4,10) = 0.59188d0
|
|
||||||
! one-electron terms
|
|
||||||
c_mn_o_sm_17(2,0,0,10) = -0.00567d0
|
|
||||||
c_mn_o_sm_17(3,0,0,10) = 0.14011d0
|
|
||||||
c_mn_o_sm_17(4,0,0,10) = -0.05671d0
|
|
||||||
c_mn_o_sm_17(2,2,0,10) = -3.33767d0
|
|
||||||
! ee-n terms
|
|
||||||
c_mn_o_sm_17(2,0,2,10) = 1.95067d0
|
|
||||||
c_mn_o_sm_17(2,2,2,10) = 6.83340d0
|
|
||||||
c_mn_o_sm_17(4,0,2,10) = -3.29231d0
|
|
||||||
c_mn_o_sm_17(2,0,4,10) = -2.44998d0
|
|
||||||
c_mn_o_sm_17(4,2,2,10) = -2.13029d0
|
|
||||||
c_mn_o_sm_17(6,0,2,10) = 2.25768d0
|
|
||||||
c_mn_o_sm_17(4,0,4,10) = 1.97951d0
|
|
||||||
c_mn_o_sm_17(2,2,4,10) = -2.0924160
|
|
||||||
c_mn_o_sm_17(2,0,6,10) = 0.35493d0
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, b_I_sm_90,(2:10)]
|
|
||||||
&BEGIN_PROVIDER [ double precision, d_I_sm_90,(2:10)]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! "b_I" and "d_I" parameters of Eqs. (4) and (5) of Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
END_DOC
|
|
||||||
b_I_sm_90 = 1.d0
|
|
||||||
d_I_sm_90 = 1.d0
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
subroutine get_full_sm_90_jastrow(r1,r2,rI,sm_j,i_charge, j_1e,j_2e,j_een,j_tot)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(in) :: r1(3),r2(3),rI(3)
|
|
||||||
integer, intent(in) :: sm_j, i_charge
|
|
||||||
double precision, intent(out):: j_1e,j_2e,j_een,j_tot
|
|
||||||
BEGIN_DOC
|
|
||||||
! Jastrow function as in Eq. (4) of Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
! the i_charge variable is the integer specifying the charge of the atom for the Jastrow
|
|
||||||
! the sm_j integer variable represents the "quality" of the jastrow : sm_j = 7, 9, 17
|
|
||||||
END_DOC
|
|
||||||
double precision :: r_inucl,r_jnucl,r_ij,b_I, d_I
|
|
||||||
b_I = b_I_sm_90(i_charge)
|
|
||||||
d_I = d_I_sm_90(i_charge)
|
|
||||||
call get_rescaled_variables_j_sm_90(r1,r2,rI,b_I,d_I,r_inucl,r_jnucl,r_ij)
|
|
||||||
call jastrow_func_sm_90(r_inucl,r_jnucl,r_ij,sm_j,i_charge, j_1e,j_2e,j_een,j_tot)
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine get_rescaled_variables_j_sm_90(r1,r2,rI,b_I,d_I,r_inucl,r_jnucl,r_ij)
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! rescaled variables of Eq. (5) and (6) of Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
! the "b_I" and "d_I" parameters are the same as in Eqs. (5) and (6)
|
|
||||||
END_DOC
|
|
||||||
double precision, intent(in) :: r1(3),r2(3),rI(3)
|
|
||||||
double precision, intent(in) :: b_I, d_I
|
|
||||||
double precision, intent(out):: r_inucl,r_jnucl,r_ij
|
|
||||||
double precision :: rin, rjn, rij
|
|
||||||
integer :: i
|
|
||||||
rin = 0.d0
|
|
||||||
rjn = 0.d0
|
|
||||||
rij = 0.d0
|
|
||||||
do i = 1,3
|
|
||||||
rin += (r1(i) - rI(i)) * (r1(i) - rI(i))
|
|
||||||
rjn += (r2(i) - rI(i)) * (r2(i) - rI(i))
|
|
||||||
rij += (r2(i) - r1(i)) * (r2(i) - r1(i))
|
|
||||||
enddo
|
|
||||||
rin = dsqrt(rin)
|
|
||||||
rjn = dsqrt(rjn)
|
|
||||||
rij = dsqrt(rij)
|
|
||||||
r_inucl = b_I * rin/(1.d0 + b_I * rin)
|
|
||||||
r_jnucl = b_I * rjn/(1.d0 + b_I * rjn)
|
|
||||||
r_ij = d_I * rij/(1.d0 + b_I * rij)
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine jastrow_func_sm_90(r_inucl,r_jnucl,r_ij,sm_j,i_charge, j_1e,j_2e,j_een,j_tot)
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! Jastrow function as in Eq. (4) of Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
! Here the r_inucl, r_jnucl are the rescaled variables as defined in Eq. (5) with "b_I"
|
|
||||||
! r_ij is the rescaled variable as defined in Eq. (6) with "d_I"
|
|
||||||
! the i_charge variable is the integer specifying the charge of the atom for the Jastrow
|
|
||||||
! the sm_j integer variable represents the "quality" of the jastrow : sm_j = 7, 9, 17
|
|
||||||
!
|
|
||||||
! it returns the j_1e : sum of terms with "o" = "n" = 0, "m" /= 0,
|
|
||||||
! j_2e : sum of terms with "m" = "n" = 0, "o" /= 0,
|
|
||||||
! j_een : sum of terms with "m" /=0, "n" /= 0, "o" /= 0,
|
|
||||||
! j_tot : the total sum
|
|
||||||
END_DOC
|
|
||||||
double precision, intent(in) :: r_inucl,r_jnucl,r_ij
|
|
||||||
integer, intent(in) :: sm_j,i_charge
|
|
||||||
double precision, intent(out):: j_1e,j_2e,j_een,j_tot
|
|
||||||
j_1e = 0.D0
|
|
||||||
j_2e = 0.D0
|
|
||||||
j_een = 0.D0
|
|
||||||
double precision :: delta_mn,jastrow_sm_90_atomic
|
|
||||||
integer :: m,n,o
|
|
||||||
BEGIN_TEMPLATE
|
|
||||||
! pure 2e part
|
|
||||||
n = 0
|
|
||||||
m = 0
|
|
||||||
if(sm_j == $X )then
|
|
||||||
do o = 1, o_max_sm_$X
|
|
||||||
if(dabs(c_mn_o_sm_$X(m,n,o,i_charge)).lt.1.d-10)cycle
|
|
||||||
j_2e += c_mn_o_sm_$X(m,n,o,i_charge) * jastrow_sm_90_atomic(m,n,o,i_charge,r_inucl,r_jnucl,r_ij)
|
|
||||||
enddo
|
|
||||||
! else
|
|
||||||
! print*,'sm_j = ',sm_j
|
|
||||||
! print*,'not implemented, stop'
|
|
||||||
! stop
|
|
||||||
endif
|
|
||||||
! pure one-e part
|
|
||||||
o = 0
|
|
||||||
if(sm_j == $X)then
|
|
||||||
do n = 2, n_max_sm_$X
|
|
||||||
do m = 2, m_max_sm_$X
|
|
||||||
j_1e += c_mn_o_sm_$X(m,n,o,i_charge) * jastrow_sm_90_atomic(m,n,o,i_charge,r_inucl,r_jnucl,r_ij)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
! else
|
|
||||||
! print*,'sm_j = ',sm_j
|
|
||||||
! print*,'not implemented, stop'
|
|
||||||
! stop
|
|
||||||
endif
|
|
||||||
! e-e-n part
|
|
||||||
if(sm_j == $X)then
|
|
||||||
do o = 1, o_max_sm_$X
|
|
||||||
do m = 2, m_max_sm_$X
|
|
||||||
do n = 2, n_max_sm_$X
|
|
||||||
j_een += c_mn_o_sm_$X(m,n,o,i_charge) * jastrow_sm_90_atomic(m,n,o,i_charge,r_inucl,r_jnucl,r_ij)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
else
|
|
||||||
! print*,'sm_j = ',sm_j
|
|
||||||
! print*,'not implemented, stop'
|
|
||||||
! stop
|
|
||||||
endif
|
|
||||||
j_tot = j_1e + j_2e + j_een
|
|
||||||
SUBST [ X]
|
|
||||||
7 ;;
|
|
||||||
9 ;;
|
|
||||||
17 ;;
|
|
||||||
END_TEMPLATE
|
|
||||||
end
|
|
||||||
|
|
||||||
double precision function jastrow_sm_90_atomic(m,n,o,i_charge,r_inucl,r_jnucl,r_ij)
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! contribution to the function of Eq. (4) of Schmidt,Moskowitz, JCP, 93, 4172 (1990)
|
|
||||||
! for a given m,n,o and atom
|
|
||||||
END_DOC
|
|
||||||
double precision, intent(in) :: r_inucl,r_jnucl,r_ij
|
|
||||||
integer , intent(in) :: m,n,o,i_charge
|
|
||||||
double precision :: delta_mn
|
|
||||||
if(m==n)then
|
|
||||||
delta_mn = 0.5d0
|
|
||||||
else
|
|
||||||
delta_mn = 1.D0
|
|
||||||
endif
|
|
||||||
jastrow_sm_90_atomic = delta_mn * (r_inucl**m * r_jnucl**n + r_jnucl**m * r_inucl**n)*r_ij**o
|
|
||||||
end
|
|
||||||
@ -1,69 +0,0 @@
|
|||||||
program plot_j
|
|
||||||
implicit none
|
|
||||||
double precision :: r1(3),rI(3),r2(3)
|
|
||||||
double precision :: r12,dx,xmax, j_1e,j_2e,j_een,j_tot
|
|
||||||
double precision :: j_mu_F_x_j
|
|
||||||
integer :: i,nx,m,i_charge,sm_j
|
|
||||||
|
|
||||||
character*(128) :: output
|
|
||||||
integer :: i_unit_output_He_sm_7,i_unit_output_Ne_sm_7
|
|
||||||
integer :: i_unit_output_He_sm_17,i_unit_output_Ne_sm_17
|
|
||||||
integer :: getUnitAndOpen
|
|
||||||
output='J_SM_7_He'
|
|
||||||
i_unit_output_He_sm_7 = getUnitAndOpen(output,'w')
|
|
||||||
output='J_SM_7_Ne'
|
|
||||||
i_unit_output_Ne_sm_7 = getUnitAndOpen(output,'w')
|
|
||||||
|
|
||||||
output='J_SM_17_He'
|
|
||||||
i_unit_output_He_sm_17 = getUnitAndOpen(output,'w')
|
|
||||||
output='J_SM_17_Ne'
|
|
||||||
i_unit_output_Ne_sm_17 = getUnitAndOpen(output,'w')
|
|
||||||
|
|
||||||
rI = 0.d0
|
|
||||||
r1 = 0.d0
|
|
||||||
r2 = 0.d0
|
|
||||||
r1(1) = 1.5d0
|
|
||||||
xmax = 20.d0
|
|
||||||
r2(1) = -xmax*0.5d0
|
|
||||||
nx = 1000
|
|
||||||
dx = xmax/dble(nx)
|
|
||||||
do i = 1, nx
|
|
||||||
r12 = 0.d0
|
|
||||||
do m = 1, 3
|
|
||||||
r12 += (r1(m) - r2(m))*(r1(m) - r2(m))
|
|
||||||
enddo
|
|
||||||
r12 = dsqrt(r12)
|
|
||||||
double precision :: jmu,env_nucl,jmu_env,jmu_scaled, jmu_scaled_env
|
|
||||||
double precision :: b_I,d_I,r_inucl,r_jnucl,r_ij
|
|
||||||
b_I = 1.D0
|
|
||||||
d_I = 1.D0
|
|
||||||
call get_rescaled_variables_j_sm_90(r1,r2,rI,b_I,d_I,r_inucl,r_jnucl,r_ij)
|
|
||||||
jmu=j_mu_F_x_j(r12)
|
|
||||||
jmu_scaled=j_mu_F_x_j(r_ij)
|
|
||||||
jmu_env = jmu * env_nucl(r1) * env_nucl(r2)
|
|
||||||
! jmu_scaled_env= jmu_scaled * (1.d0 - env_coef(1) * dexp(-env_expo(1)*r_inucl**2)) * (1.d0 - env_coef(1) * dexp(-env_expo(1)*r_jnucl**2))
|
|
||||||
jmu_scaled_env= jmu_scaled * env_nucl(r1) * env_nucl(r2)
|
|
||||||
! He
|
|
||||||
i_charge = 2
|
|
||||||
! SM 7 Jastrow
|
|
||||||
sm_j = 7
|
|
||||||
call get_full_sm_90_jastrow(r1,r2,rI,sm_j,i_charge, j_1e,j_2e,j_een,j_tot)
|
|
||||||
write(i_unit_output_He_sm_7,'(100(F16.10,X))')r2(1),r12,j_mu_F_x_j(r12), j_1e,j_2e,j_een,j_tot,jmu_env,jmu_scaled,jmu_scaled_env
|
|
||||||
! SM 17 Jastrow
|
|
||||||
sm_j = 17
|
|
||||||
call get_full_sm_90_jastrow(r1,r2,rI,sm_j,i_charge, j_1e,j_2e,j_een,j_tot)
|
|
||||||
write(i_unit_output_He_sm_17,'(100(F16.10,X))')r2(1),r12,j_mu_F_x_j(r12), j_1e,j_2e,j_een,j_tot,jmu_env,jmu_scaled,jmu_scaled_env
|
|
||||||
! Ne
|
|
||||||
i_charge = 10
|
|
||||||
! SM 7 Jastrow
|
|
||||||
sm_j = 7
|
|
||||||
call get_full_sm_90_jastrow(r1,r2,rI,sm_j,i_charge, j_1e,j_2e,j_een,j_tot)
|
|
||||||
write(i_unit_output_Ne_sm_7,'(100(F16.10,X))')r2(1),r12,j_mu_F_x_j(r12), j_1e,j_2e,j_een,j_tot,jmu_env,jmu_scaled,jmu_scaled_env
|
|
||||||
! SM 17 Jastrow
|
|
||||||
sm_j = 17
|
|
||||||
call get_full_sm_90_jastrow(r1,r2,rI,sm_j,i_charge, j_1e,j_2e,j_een,j_tot)
|
|
||||||
write(i_unit_output_Ne_sm_17,'(100(F16.10,X))')r2(1),r12,j_mu_F_x_j(r12), j_1e,j_2e,j_een,j_tot,jmu_env,jmu_scaled,jmu_scaled_env
|
|
||||||
r2(1) += dx
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
@ -1,59 +0,0 @@
|
|||||||
program print_fit_param
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! TODO : Put the documentation of the program here
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
my_grid_becke = .True.
|
|
||||||
PROVIDE tc_grid1_a tc_grid1_r
|
|
||||||
my_n_pt_r_grid = tc_grid1_r
|
|
||||||
my_n_pt_a_grid = tc_grid1_a
|
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
|
||||||
|
|
||||||
!call create_guess
|
|
||||||
!call orthonormalize_mos
|
|
||||||
|
|
||||||
call main()
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine main()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i
|
|
||||||
|
|
||||||
mu_erf = 1.d0
|
|
||||||
touch mu_erf
|
|
||||||
|
|
||||||
print *, ' fit for (1 - erf(x))^2'
|
|
||||||
do i = 1, n_max_fit_slat
|
|
||||||
print*, expo_gauss_1_erf_x_2(i), coef_gauss_1_erf_x_2(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]'
|
|
||||||
do i = 1, n_max_fit_slat
|
|
||||||
print *, expo_gauss_j_mu_x(i), 2.d0 * coef_gauss_j_mu_x(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]^2'
|
|
||||||
do i = 1, n_max_fit_slat
|
|
||||||
print *, expo_gauss_j_mu_x_2(i), 4.d0 * coef_gauss_j_mu_x_2(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)] x [1 - erf(mu * r12)]'
|
|
||||||
do i = 1, n_max_fit_slat
|
|
||||||
print *, expo_gauss_j_mu_1_erf(i), 4.d0 * coef_gauss_j_mu_1_erf(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
return
|
|
||||||
end subroutine main
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -1,55 +0,0 @@
|
|||||||
program print_tcscf_energy
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! TODO : Put the documentation of the program here
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
print *, 'Hello world'
|
|
||||||
my_grid_becke = .True.
|
|
||||||
PROVIDE tc_grid1_a tc_grid1_r
|
|
||||||
my_n_pt_r_grid = tc_grid1_r
|
|
||||||
my_n_pt_a_grid = tc_grid1_a
|
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
|
||||||
|
|
||||||
call main()
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine main()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
double precision :: etc_tot, etc_1e, etc_2e, etc_3e
|
|
||||||
|
|
||||||
PROVIDE j2e_type mu_erf
|
|
||||||
PROVIDE j1e_type j1e_coef j1e_expo
|
|
||||||
PROVIDE env_type env_coef env_expo
|
|
||||||
|
|
||||||
print*, ' j2e_type = ', j2e_type
|
|
||||||
print*, ' j1e_type = ', j1e_type
|
|
||||||
print*, ' env_type = ', env_type
|
|
||||||
|
|
||||||
print*, ' mu_erf = ', mu_erf
|
|
||||||
|
|
||||||
etc_tot = TC_HF_energy
|
|
||||||
etc_1e = TC_HF_one_e_energy
|
|
||||||
etc_2e = TC_HF_two_e_energy
|
|
||||||
etc_3e = 0.d0
|
|
||||||
if(three_body_h_tc) then
|
|
||||||
!etc_3e = diag_three_elem_hf
|
|
||||||
etc_3e = tcscf_energy_3e_naive
|
|
||||||
endif
|
|
||||||
|
|
||||||
print *, " E_TC = ", etc_tot
|
|
||||||
print *, " E_1e = ", etc_1e
|
|
||||||
print *, " E_2e = ", etc_2e
|
|
||||||
print *, " E_3e = ", etc_3e
|
|
||||||
|
|
||||||
return
|
|
||||||
end subroutine main
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -1,129 +0,0 @@
|
|||||||
! ---
|
|
||||||
|
|
||||||
subroutine rh_tcscf_simple()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, it, dim_DIIS
|
|
||||||
double precision :: t0, t1
|
|
||||||
double precision :: e_save, e_delta, rho_delta
|
|
||||||
double precision :: etc_tot, etc_1e, etc_2e, etc_3e, tc_grad
|
|
||||||
double precision :: er_DIIS
|
|
||||||
double precision, allocatable :: rho_old(:,:), rho_new(:,:)
|
|
||||||
|
|
||||||
allocate(rho_old(ao_num,ao_num), rho_new(ao_num,ao_num))
|
|
||||||
|
|
||||||
it = 0
|
|
||||||
e_save = 0.d0
|
|
||||||
dim_DIIS = 0
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
if(.not. bi_ortho) then
|
|
||||||
print *, ' grad_hermit = ', grad_hermit
|
|
||||||
call save_good_hermit_tc_eigvectors
|
|
||||||
TOUCH mo_coef
|
|
||||||
call save_mos
|
|
||||||
endif
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
if(bi_ortho) then
|
|
||||||
|
|
||||||
PROVIDE level_shift_tcscf
|
|
||||||
PROVIDE mo_l_coef mo_r_coef
|
|
||||||
|
|
||||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
|
||||||
'====', '================', '================', '================', '================', '================' &
|
|
||||||
, '================', '================', '================', '====', '========'
|
|
||||||
|
|
||||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
|
||||||
' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' &
|
|
||||||
, ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)'
|
|
||||||
|
|
||||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
|
||||||
'====', '================', '================', '================', '================', '================' &
|
|
||||||
, '================', '================', '================', '====', '========'
|
|
||||||
|
|
||||||
|
|
||||||
! first iteration (HF orbitals)
|
|
||||||
call wall_time(t0)
|
|
||||||
|
|
||||||
etc_tot = TC_HF_energy
|
|
||||||
etc_1e = TC_HF_one_e_energy
|
|
||||||
etc_2e = TC_HF_two_e_energy
|
|
||||||
etc_3e = 0.d0
|
|
||||||
if(three_body_h_tc) then
|
|
||||||
etc_3e = diag_three_elem_hf
|
|
||||||
endif
|
|
||||||
tc_grad = grad_non_hermit
|
|
||||||
er_DIIS = maxval(abs(FQS_SQF_mo))
|
|
||||||
e_delta = dabs(etc_tot - e_save)
|
|
||||||
e_save = etc_tot
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
|
||||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
|
||||||
|
|
||||||
do while(tc_grad .gt. dsqrt(thresh_tcscf))
|
|
||||||
call wall_time(t0)
|
|
||||||
|
|
||||||
it += 1
|
|
||||||
if(it > n_it_tcscf_max) then
|
|
||||||
print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
mo_l_coef = fock_tc_leigvec_ao
|
|
||||||
mo_r_coef = fock_tc_reigvec_ao
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
|
||||||
TOUCH mo_l_coef mo_r_coef
|
|
||||||
|
|
||||||
etc_tot = TC_HF_energy
|
|
||||||
etc_1e = TC_HF_one_e_energy
|
|
||||||
etc_2e = TC_HF_two_e_energy
|
|
||||||
etc_3e = 0.d0
|
|
||||||
if(three_body_h_tc) then
|
|
||||||
etc_3e = diag_three_elem_hf
|
|
||||||
endif
|
|
||||||
tc_grad = grad_non_hermit
|
|
||||||
er_DIIS = maxval(abs(FQS_SQF_mo))
|
|
||||||
e_delta = dabs(etc_tot - e_save)
|
|
||||||
e_save = etc_tot
|
|
||||||
|
|
||||||
call ezfio_set_tc_scf_tcscf_energy(etc_tot)
|
|
||||||
|
|
||||||
call wall_time(t1)
|
|
||||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
|
||||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
|
||||||
enddo
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
do while( (grad_hermit.gt.dsqrt(thresh_tcscf)) .and. (it.lt.n_it_tcscf_max) )
|
|
||||||
print*,'grad_hermit = ',grad_hermit
|
|
||||||
it += 1
|
|
||||||
print *, 'iteration = ', it
|
|
||||||
print *, '***'
|
|
||||||
print *, 'TC HF total energy = ', TC_HF_energy
|
|
||||||
print *, 'TC HF 1 e energy = ', TC_HF_one_e_energy
|
|
||||||
print *, 'TC HF 2 e energy = ', TC_HF_two_e_energy
|
|
||||||
print *, 'TC HF 3 body = ', diag_three_elem_hf
|
|
||||||
print *, '***'
|
|
||||||
print *, ''
|
|
||||||
call save_good_hermit_tc_eigvectors
|
|
||||||
TOUCH mo_coef
|
|
||||||
call save_mos
|
|
||||||
enddo
|
|
||||||
|
|
||||||
endif
|
|
||||||
|
|
||||||
print *, ' TCSCF Simple converged !'
|
|
||||||
!call print_energy_and_mos(good_angles)
|
|
||||||
|
|
||||||
deallocate(rho_old, rho_new)
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -1,369 +0,0 @@
|
|||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
program rotate_tcscf_orbitals
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! TODO : Rotate the bi-orthonormal orbitals in order to minimize left-right angles when degenerate
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
my_grid_becke = .True.
|
|
||||||
PROVIDE tc_grid1_a tc_grid1_r
|
|
||||||
my_n_pt_r_grid = tc_grid1_r
|
|
||||||
my_n_pt_a_grid = tc_grid1_a
|
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
|
||||||
|
|
||||||
bi_ortho = .True.
|
|
||||||
touch bi_ortho
|
|
||||||
|
|
||||||
call minimize_tc_orb_angles()
|
|
||||||
!call maximize_overlap()
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine maximize_overlap()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, m, n
|
|
||||||
double precision :: accu_d, accu_nd
|
|
||||||
double precision, allocatable :: C(:,:), R(:,:), L(:,:), W(:,:), e(:)
|
|
||||||
double precision, allocatable :: S(:,:)
|
|
||||||
|
|
||||||
n = ao_num
|
|
||||||
m = mo_num
|
|
||||||
|
|
||||||
allocate(L(n,m), R(n,m), C(n,m), W(n,n), e(m))
|
|
||||||
L = mo_l_coef
|
|
||||||
R = mo_r_coef
|
|
||||||
C = mo_coef
|
|
||||||
W = ao_overlap
|
|
||||||
|
|
||||||
print*, ' fock matrix diag elements'
|
|
||||||
do i = 1, m
|
|
||||||
e(i) = Fock_matrix_tc_mo_tot(i,i)
|
|
||||||
print*, e(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
print *, ' overlap before :'
|
|
||||||
print *, ' '
|
|
||||||
|
|
||||||
allocate(S(m,m))
|
|
||||||
|
|
||||||
call LTxSxR(n, m, L, W, R, S)
|
|
||||||
!print*, " L.T x R"
|
|
||||||
!do i = 1, m
|
|
||||||
! write(*, '(100(F16.10,X))') S(i,i)
|
|
||||||
!enddo
|
|
||||||
call LTxSxR(n, m, L, W, C, S)
|
|
||||||
print*, " L.T x C"
|
|
||||||
do i = 1, m
|
|
||||||
write(*, '(100(F16.10,X))') S(i,:)
|
|
||||||
enddo
|
|
||||||
call LTxSxR(n, m, C, W, R, S)
|
|
||||||
print*, " C.T x R"
|
|
||||||
do i = 1, m
|
|
||||||
write(*, '(100(F16.10,X))') S(i,:)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
deallocate(S)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
call rotate_degen_eigvec_to_maximize_overlap(n, m, e, C, W, L, R)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
print *, ' overlap after :'
|
|
||||||
print *, ' '
|
|
||||||
|
|
||||||
allocate(S(m,m))
|
|
||||||
|
|
||||||
call LTxSxR(n, m, L, W, R, S)
|
|
||||||
!print*, " L.T x R"
|
|
||||||
!do i = 1, m
|
|
||||||
! write(*, '(100(F16.10,X))') S(i,i)
|
|
||||||
!enddo
|
|
||||||
call LTxSxR(n, m, L, W, C, S)
|
|
||||||
print*, " L.T x C"
|
|
||||||
do i = 1, m
|
|
||||||
write(*, '(100(F16.10,X))') S(i,:)
|
|
||||||
enddo
|
|
||||||
call LTxSxR(n, m, C, W, R, S)
|
|
||||||
print*, " C.T x R"
|
|
||||||
do i = 1, m
|
|
||||||
write(*, '(100(F16.10,X))') S(i,:)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
deallocate(S)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
mo_l_coef = L
|
|
||||||
mo_r_coef = R
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
deallocate(L, R, C, W, e)
|
|
||||||
|
|
||||||
end subroutine maximize_overlap
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine rotate_degen_eigvec_to_maximize_overlap(n, m, e0, C0, W0, L0, R0)
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer, intent(in) :: n, m
|
|
||||||
double precision, intent(in) :: e0(m), W0(n,n), C0(n,m)
|
|
||||||
double precision, intent(inout) :: L0(n,m), R0(n,m)
|
|
||||||
|
|
||||||
|
|
||||||
integer :: i, j, k, kk, mm, id1, tot_deg
|
|
||||||
double precision :: ei, ej, de, de_thr
|
|
||||||
integer, allocatable :: deg_num(:)
|
|
||||||
double precision, allocatable :: L(:,:), R(:,:), C(:,:), Lnew(:,:), Rnew(:,:), tmp(:,:)
|
|
||||||
!double precision, allocatable :: S(:,:), Snew(:,:), T(:,:), Ttmp(:,:), Stmp(:,:)
|
|
||||||
double precision, allocatable :: S(:,:), Snew(:,:), T(:,:), Ttmp(:,:), Stmp(:,:)
|
|
||||||
!real*8 :: S(m,m), Snew(m,m), T(m,m)
|
|
||||||
|
|
||||||
id1 = 700
|
|
||||||
allocate(S(id1,id1), Snew(id1,id1), T(id1,id1))
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
allocate( deg_num(m) )
|
|
||||||
do i = 1, m
|
|
||||||
deg_num(i) = 1
|
|
||||||
enddo
|
|
||||||
|
|
||||||
de_thr = thr_degen_tc
|
|
||||||
|
|
||||||
do i = 1, m-1
|
|
||||||
ei = e0(i)
|
|
||||||
|
|
||||||
! already considered in degen vectors
|
|
||||||
if(deg_num(i).eq.0) cycle
|
|
||||||
|
|
||||||
do j = i+1, m
|
|
||||||
ej = e0(j)
|
|
||||||
de = dabs(ei - ej)
|
|
||||||
|
|
||||||
if(de .lt. de_thr) then
|
|
||||||
deg_num(i) = deg_num(i) + 1
|
|
||||||
deg_num(j) = 0
|
|
||||||
endif
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
tot_deg = 0
|
|
||||||
do i = 1, m
|
|
||||||
if(deg_num(i).gt.1) then
|
|
||||||
print *, ' degen on', i, deg_num(i)
|
|
||||||
tot_deg = tot_deg + 1
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
if(tot_deg .eq. 0) then
|
|
||||||
print *, ' no degen'
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
do i = 1, m
|
|
||||||
mm = deg_num(i)
|
|
||||||
|
|
||||||
if(mm .gt. 1) then
|
|
||||||
|
|
||||||
allocate(L(n,mm), R(n,mm), C(n,mm))
|
|
||||||
do j = 1, mm
|
|
||||||
L(1:n,j) = L0(1:n,i+j-1)
|
|
||||||
R(1:n,j) = R0(1:n,i+j-1)
|
|
||||||
C(1:n,j) = C0(1:n,i+j-1)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
! C.T x W0 x R
|
|
||||||
allocate(tmp(mm,n), Stmp(mm,mm))
|
|
||||||
call dgemm( 'T', 'N', mm, n, n, 1.d0 &
|
|
||||||
, C, size(C, 1), W0, size(W0, 1) &
|
|
||||||
, 0.d0, tmp, size(tmp, 1) )
|
|
||||||
call dgemm( 'N', 'N', mm, mm, n, 1.d0 &
|
|
||||||
, tmp, size(tmp, 1), R, size(R, 1) &
|
|
||||||
, 0.d0, Stmp, size(Stmp, 1) )
|
|
||||||
deallocate(C, tmp)
|
|
||||||
|
|
||||||
S = 0.d0
|
|
||||||
do k = 1, mm
|
|
||||||
do kk = 1, mm
|
|
||||||
S(kk,k) = Stmp(kk,k)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
deallocate(Stmp)
|
|
||||||
|
|
||||||
!print*, " overlap bef"
|
|
||||||
!do k = 1, mm
|
|
||||||
! write(*, '(100(F16.10,X))') (S(k,kk), kk=1, mm)
|
|
||||||
!enddo
|
|
||||||
|
|
||||||
T = 0.d0
|
|
||||||
Snew = 0.d0
|
|
||||||
call maxovl(mm, mm, S, T, Snew)
|
|
||||||
|
|
||||||
!print*, " overlap aft"
|
|
||||||
!do k = 1, mm
|
|
||||||
! write(*, '(100(F16.10,X))') (Snew(k,kk), kk=1, mm)
|
|
||||||
!enddo
|
|
||||||
|
|
||||||
allocate(Ttmp(mm,mm))
|
|
||||||
Ttmp(1:mm,1:mm) = T(1:mm,1:mm)
|
|
||||||
|
|
||||||
allocate(Lnew(n,mm), Rnew(n,mm))
|
|
||||||
call dgemm( 'N', 'N', n, mm, mm, 1.d0 &
|
|
||||||
, R, size(R, 1), Ttmp(1,1), size(Ttmp, 1) &
|
|
||||||
, 0.d0, Rnew, size(Rnew, 1) )
|
|
||||||
call dgemm( 'N', 'N', n, mm, mm, 1.d0 &
|
|
||||||
, L, size(L, 1), Ttmp(1,1), size(Ttmp, 1) &
|
|
||||||
, 0.d0, Lnew, size(Lnew, 1) )
|
|
||||||
|
|
||||||
deallocate(L, R)
|
|
||||||
deallocate(Ttmp)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
do j = 1, mm
|
|
||||||
L0(1:n,i+j-1) = Lnew(1:n,j)
|
|
||||||
R0(1:n,i+j-1) = Rnew(1:n,j)
|
|
||||||
enddo
|
|
||||||
deallocate(Lnew, Rnew)
|
|
||||||
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
deallocate(S, Snew, T)
|
|
||||||
|
|
||||||
end subroutine rotate_degen_eigvec_to_maximize_overlap
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine fix_right_to_one()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, m, n, mm, tot_deg
|
|
||||||
double precision :: accu_d, accu_nd
|
|
||||||
double precision :: de_thr, ei, ej, de
|
|
||||||
integer, allocatable :: deg_num(:)
|
|
||||||
double precision, allocatable :: R0(:,:), L0(:,:), W(:,:), e0(:)
|
|
||||||
double precision, allocatable :: R(:,:), L(:,:), S(:,:), Stmp(:,:), tmp(:,:)
|
|
||||||
|
|
||||||
n = ao_num
|
|
||||||
m = mo_num
|
|
||||||
|
|
||||||
allocate(L0(n,m), R0(n,m), W(n,n), e0(m))
|
|
||||||
L0 = mo_l_coef
|
|
||||||
R0 = mo_r_coef
|
|
||||||
W = ao_overlap
|
|
||||||
|
|
||||||
print*, ' fock matrix diag elements'
|
|
||||||
do i = 1, m
|
|
||||||
e0(i) = Fock_matrix_tc_mo_tot(i,i)
|
|
||||||
print*, e0(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
allocate( deg_num(m) )
|
|
||||||
do i = 1, m
|
|
||||||
deg_num(i) = 1
|
|
||||||
enddo
|
|
||||||
|
|
||||||
de_thr = 1d-6
|
|
||||||
|
|
||||||
do i = 1, m-1
|
|
||||||
ei = e0(i)
|
|
||||||
|
|
||||||
! already considered in degen vectors
|
|
||||||
if(deg_num(i).eq.0) cycle
|
|
||||||
|
|
||||||
do j = i+1, m
|
|
||||||
ej = e0(j)
|
|
||||||
de = dabs(ei - ej)
|
|
||||||
|
|
||||||
if(de .lt. de_thr) then
|
|
||||||
deg_num(i) = deg_num(i) + 1
|
|
||||||
deg_num(j) = 0
|
|
||||||
endif
|
|
||||||
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
deallocate(e0)
|
|
||||||
|
|
||||||
tot_deg = 0
|
|
||||||
do i = 1, m
|
|
||||||
if(deg_num(i).gt.1) then
|
|
||||||
print *, ' degen on', i, deg_num(i)
|
|
||||||
tot_deg = tot_deg + 1
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
if(tot_deg .eq. 0) then
|
|
||||||
print *, ' no degen'
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
do i = 1, m
|
|
||||||
mm = deg_num(i)
|
|
||||||
|
|
||||||
if(mm .gt. 1) then
|
|
||||||
|
|
||||||
allocate(L(n,mm), R(n,mm))
|
|
||||||
do j = 1, mm
|
|
||||||
L(1:n,j) = L0(1:n,i+j-1)
|
|
||||||
R(1:n,j) = R0(1:n,i+j-1)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
call impose_weighted_orthog_svd(n, mm, W, R)
|
|
||||||
call impose_weighted_biorthog_qr(n, mm, thresh_biorthog_diag, thresh_biorthog_nondiag, R, W, L)
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
do j = 1, mm
|
|
||||||
L0(1:n,i+j-1) = L(1:n,j)
|
|
||||||
R0(1:n,i+j-1) = R(1:n,j)
|
|
||||||
enddo
|
|
||||||
deallocate(L, R)
|
|
||||||
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
call check_weighted_biorthog_binormalize(n, m, L0, W, R0, thresh_biorthog_diag, thresh_biorthog_nondiag, .true.)
|
|
||||||
|
|
||||||
deallocate(W, deg_num)
|
|
||||||
|
|
||||||
mo_l_coef = L0
|
|
||||||
mo_r_coef = R0
|
|
||||||
deallocate(L0, R0)
|
|
||||||
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
|
||||||
print *, ' orbitals are rotated '
|
|
||||||
|
|
||||||
return
|
|
||||||
end subroutine fix_right_to_one
|
|
||||||
|
|
||||||
! ---
|
|
||||||
@ -1,91 +0,0 @@
|
|||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
program tc_petermann_factor
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! TODO : Put the documentation of the program here
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
my_grid_becke = .True.
|
|
||||||
PROVIDE tc_grid1_a tc_grid1_r
|
|
||||||
my_n_pt_r_grid = tc_grid1_r
|
|
||||||
my_n_pt_a_grid = tc_grid1_a
|
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
|
||||||
|
|
||||||
call main()
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine main()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j
|
|
||||||
double precision :: Pf_diag_av
|
|
||||||
double precision, allocatable :: Sl(:,:), Sr(:,:), Pf(:,:)
|
|
||||||
|
|
||||||
allocate(Sl(mo_num,mo_num), Sr(mo_num,mo_num), Pf(mo_num,mo_num))
|
|
||||||
|
|
||||||
|
|
||||||
call LTxSxR(ao_num, mo_num, mo_l_coef, ao_overlap, mo_r_coef, Sl)
|
|
||||||
!call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
|
|
||||||
! , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
|
|
||||||
! , 0.d0, Sl, size(Sl, 1) )
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' left-right orthog matrix:'
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(100(F8.4,X))') Sl(:,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
call LTxSxR(ao_num, mo_num, mo_l_coef, ao_overlap, mo_l_coef, Sl)
|
|
||||||
!call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
|
|
||||||
! , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
|
|
||||||
! , 0.d0, Sl, size(Sl, 1) )
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' left-orthog matrix:'
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(100(F8.4,X))') Sl(:,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
call LTxSxR(ao_num, mo_num, mo_r_coef, ao_overlap, mo_r_coef, Sr)
|
|
||||||
! call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
|
|
||||||
! , mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
|
|
||||||
! , 0.d0, Sr, size(Sr, 1) )
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' right-orthog matrix:'
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(100(F8.4,X))') Sr(:,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' Petermann matrix:'
|
|
||||||
do i = 1, mo_num
|
|
||||||
do j = 1, mo_num
|
|
||||||
Pf(j,i) = Sl(j,i) * Sr(j,i)
|
|
||||||
enddo
|
|
||||||
write(*,'(100(F8.4,X))') Pf(:,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
Pf_diag_av = 0.d0
|
|
||||||
do i = 1, mo_num
|
|
||||||
Pf_diag_av = Pf_diag_av + Pf(i,i)
|
|
||||||
enddo
|
|
||||||
Pf_diag_av = Pf_diag_av / dble(mo_num)
|
|
||||||
|
|
||||||
print *, ''
|
|
||||||
print *, ' mean of the diagonal Petermann factor = ', Pf_diag_av
|
|
||||||
|
|
||||||
deallocate(Sl, Sr, Pf)
|
|
||||||
|
|
||||||
return
|
|
||||||
end subroutine
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -10,13 +10,10 @@ program tc_scf
|
|||||||
integer :: i
|
integer :: i
|
||||||
logical :: good_angles
|
logical :: good_angles
|
||||||
|
|
||||||
PROVIDE j1e_type
|
|
||||||
PROVIDE j2e_type
|
|
||||||
PROVIDE tcscf_algorithm
|
|
||||||
|
|
||||||
print *, ' TC-SCF with:'
|
print *, ' TC-SCF with:'
|
||||||
print *, ' j1e_type = ', j1e_type
|
|
||||||
print *, ' j2e_type = ', j2e_type
|
print *, ' j2e_type = ', j2e_type
|
||||||
|
print *, ' j1e_type = ', j1e_type
|
||||||
|
print *, ' env_type = ', env_type
|
||||||
|
|
||||||
write(json_unit,json_array_open_fmt) 'tc-scf'
|
write(json_unit,json_array_open_fmt) 'tc-scf'
|
||||||
|
|
||||||
@ -29,7 +26,6 @@ program tc_scf
|
|||||||
call write_int(6, my_n_pt_r_grid, 'radial external grid over')
|
call write_int(6, my_n_pt_r_grid, 'radial external grid over')
|
||||||
call write_int(6, my_n_pt_a_grid, 'angular external grid over')
|
call write_int(6, my_n_pt_a_grid, 'angular external grid over')
|
||||||
|
|
||||||
|
|
||||||
if(tc_integ_type .eq. "numeric") then
|
if(tc_integ_type .eq. "numeric") then
|
||||||
my_extra_grid_becke = .True.
|
my_extra_grid_becke = .True.
|
||||||
PROVIDE tc_grid2_a tc_grid2_r
|
PROVIDE tc_grid2_a tc_grid2_r
|
||||||
@ -41,17 +37,7 @@ program tc_scf
|
|||||||
call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
|
call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
|
||||||
endif
|
endif
|
||||||
|
|
||||||
!call create_guess()
|
call rh_tcscf_diis()
|
||||||
!call orthonormalize_mos()
|
|
||||||
|
|
||||||
if(tcscf_algorithm == 'DIIS') then
|
|
||||||
call rh_tcscf_diis()
|
|
||||||
elseif(tcscf_algorithm == 'Simple') then
|
|
||||||
call rh_tcscf_simple()
|
|
||||||
else
|
|
||||||
print *, ' not implemented yet', tcscf_algorithm
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
PROVIDE Fock_matrix_tc_diag_mo_tot
|
PROVIDE Fock_matrix_tc_diag_mo_tot
|
||||||
print*, ' Eigenvalues:'
|
print*, ' Eigenvalues:'
|
||||||
@ -59,14 +45,11 @@ program tc_scf
|
|||||||
print*, i, Fock_matrix_tc_diag_mo_tot(i)
|
print*, i, Fock_matrix_tc_diag_mo_tot(i)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
! TODO
|
if(minimize_lr_angles) then
|
||||||
! rotate angles in separate code only if necessary
|
|
||||||
if(minimize_lr_angles)then
|
|
||||||
call minimize_tc_orb_angles()
|
call minimize_tc_orb_angles()
|
||||||
endif
|
endif
|
||||||
call print_energy_and_mos(good_angles)
|
call print_energy_and_mos(good_angles)
|
||||||
|
|
||||||
|
|
||||||
write(json_unit,json_array_close_fmtx)
|
write(json_unit,json_array_close_fmtx)
|
||||||
call json_close
|
call json_close
|
||||||
|
|
||||||
|
|||||||
@ -10,16 +10,8 @@ BEGIN_PROVIDER [double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num)
|
|||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
if(bi_ortho) then
|
PROVIDE mo_l_coef mo_r_coef
|
||||||
|
TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta
|
||||||
PROVIDE mo_l_coef mo_r_coef
|
|
||||||
TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
TCSCF_density_matrix_ao_beta = SCF_density_matrix_ao_beta
|
|
||||||
|
|
||||||
endif
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -35,16 +27,8 @@ BEGIN_PROVIDER [double precision, TCSCF_density_matrix_ao_alpha, (ao_num, ao_num
|
|||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
if(bi_ortho) then
|
PROVIDE mo_l_coef mo_r_coef
|
||||||
|
TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha
|
||||||
PROVIDE mo_l_coef mo_r_coef
|
|
||||||
TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
TCSCF_density_matrix_ao_alpha = SCF_density_matrix_ao_alpha
|
|
||||||
|
|
||||||
endif
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|||||||
@ -34,3 +34,426 @@ END_PROVIDER
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, diag_three_elem_hf]
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! < Phi_left | L | Phi_right >
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! if three_body_h_tc == false and noL_standard == true ==> do a normal ordering
|
||||||
|
!
|
||||||
|
! todo
|
||||||
|
! this should be equivalent to
|
||||||
|
! three_body_h_tc == true and noL_standard == false
|
||||||
|
!
|
||||||
|
! if three_body_h_tc == false and noL_standard == false ==> this is equal to 0
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer :: i, j, k, ipoint, mm
|
||||||
|
double precision :: contrib, weight, four_third, one_third, two_third, exchange_int_231
|
||||||
|
double precision :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb
|
||||||
|
double precision, allocatable :: tmp(:)
|
||||||
|
double precision, allocatable :: tmp_L(:,:), tmp_R(:,:)
|
||||||
|
double precision, allocatable :: tmp_M(:,:), tmp_S(:), tmp_O(:), tmp_J(:,:)
|
||||||
|
double precision, allocatable :: tmp_M_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:)
|
||||||
|
|
||||||
|
PROVIDE mo_l_coef mo_r_coef
|
||||||
|
|
||||||
|
if(.not. three_body_h_tc) then
|
||||||
|
|
||||||
|
if(noL_standard) then
|
||||||
|
PROVIDE noL_0e
|
||||||
|
diag_three_elem_hf = noL_0e
|
||||||
|
else
|
||||||
|
diag_three_elem_hf = 0.d0
|
||||||
|
endif
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
PROVIDE int2_grad1_u12_bimo_t
|
||||||
|
PROVIDE mos_l_in_r_array_transp
|
||||||
|
PROVIDE mos_r_in_r_array_transp
|
||||||
|
|
||||||
|
if(elec_alpha_num .eq. elec_beta_num) then
|
||||||
|
|
||||||
|
allocate(tmp(elec_beta_num))
|
||||||
|
allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R) &
|
||||||
|
!$OMP SHARED(elec_beta_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do j = 1, elec_beta_num
|
||||||
|
|
||||||
|
tmp_L = 0.d0
|
||||||
|
tmp_R = 0.d0
|
||||||
|
do i = 1, elec_beta_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
|
||||||
|
tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
tmp(j) = 0.d0
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
||||||
|
enddo
|
||||||
|
enddo ! j
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
diag_three_elem_hf = -2.d0 * sum(tmp)
|
||||||
|
|
||||||
|
deallocate(tmp)
|
||||||
|
deallocate(tmp_L, tmp_R)
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
|
||||||
|
tmp_O = 0.d0
|
||||||
|
tmp_J = 0.d0
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv) &
|
||||||
|
!$OMP SHARED(elec_beta_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, tmp_O, tmp_J)
|
||||||
|
|
||||||
|
allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
|
||||||
|
tmp_O_priv = 0.d0
|
||||||
|
tmp_J_priv = 0.d0
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do i = 1, elec_beta_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp_O_priv(ipoint) = tmp_O_priv(ipoint) + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
|
||||||
|
tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
|
||||||
|
tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
|
!$OMP CRITICAL
|
||||||
|
tmp_O = tmp_O + tmp_O_priv
|
||||||
|
tmp_J = tmp_J + tmp_J_priv
|
||||||
|
!$OMP END CRITICAL
|
||||||
|
|
||||||
|
deallocate(tmp_O_priv, tmp_J_priv)
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
|
||||||
|
tmp_M = 0.d0
|
||||||
|
tmp_S = 0.d0
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv) &
|
||||||
|
!$OMP SHARED(elec_beta_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, tmp_M, tmp_S)
|
||||||
|
|
||||||
|
allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
|
||||||
|
tmp_M_priv = 0.d0
|
||||||
|
tmp_S_priv = 0.d0
|
||||||
|
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = 1, elec_beta_num
|
||||||
|
do j = 1, elec_beta_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
|
||||||
|
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
|
!$OMP CRITICAL
|
||||||
|
tmp_M = tmp_M + tmp_M_priv
|
||||||
|
tmp_S = tmp_S + tmp_S_priv
|
||||||
|
!$OMP END CRITICAL
|
||||||
|
|
||||||
|
deallocate(tmp_M_priv, tmp_S_priv)
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
allocate(tmp(n_points_final_grid))
|
||||||
|
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
|
||||||
|
|
||||||
|
tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint) &
|
||||||
|
- 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
|
||||||
|
+ tmp_J(ipoint,2) * tmp_M(ipoint,2) &
|
||||||
|
+ tmp_J(ipoint,3) * tmp_M(ipoint,3)))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
diag_three_elem_hf = diag_three_elem_hf -2.d0 * (sum(tmp))
|
||||||
|
|
||||||
|
deallocate(tmp)
|
||||||
|
|
||||||
|
else ! elec_alpha_num .neq. elec_beta_num
|
||||||
|
|
||||||
|
allocate(tmp(elec_alpha_num))
|
||||||
|
allocate(tmp_L(n_points_final_grid,3), tmp_R(n_points_final_grid,3))
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R) &
|
||||||
|
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do j = 1, elec_beta_num
|
||||||
|
|
||||||
|
tmp_L = 0.d0
|
||||||
|
tmp_R = 0.d0
|
||||||
|
do i = elec_beta_num+1, elec_alpha_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_L(ipoint,1) = tmp_L(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,2) = tmp_L(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,3) = tmp_L(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
|
||||||
|
tmp_R(ipoint,1) = tmp_R(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,2) = tmp_R(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,3) = tmp_R(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
tmp(j) = 0.d0
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do i = 1, elec_beta_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
|
||||||
|
tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp(j) = tmp(j) + final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
||||||
|
enddo
|
||||||
|
enddo ! j
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(j, i, ipoint, tmp_L, tmp_R) &
|
||||||
|
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, tmp, final_weight_at_r_vector)
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do j = elec_beta_num+1, elec_alpha_num
|
||||||
|
|
||||||
|
tmp_L = 0.d0
|
||||||
|
tmp_R = 0.d0
|
||||||
|
do i = 1, elec_alpha_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp_L(ipoint,1) = tmp_L(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,2) = tmp_L(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_L(ipoint,3) = tmp_L(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i)
|
||||||
|
|
||||||
|
tmp_R(ipoint,1) = tmp_R(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,2) = tmp_R(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_R(ipoint,3) = tmp_R(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
tmp(j) = 0.d0
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp(j) = tmp(j) + 0.5d0 * final_weight_at_r_vector(ipoint) * (tmp_L(ipoint,1)*tmp_R(ipoint,1) + tmp_L(ipoint,2)*tmp_R(ipoint,2) + tmp_L(ipoint,3)*tmp_R(ipoint,3))
|
||||||
|
enddo
|
||||||
|
enddo ! j
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
diag_three_elem_hf = -2.d0 * sum(tmp)
|
||||||
|
|
||||||
|
deallocate(tmp)
|
||||||
|
deallocate(tmp_L, tmp_R)
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
allocate(tmp_O(n_points_final_grid), tmp_J(n_points_final_grid,3))
|
||||||
|
tmp_O = 0.d0
|
||||||
|
tmp_J = 0.d0
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(i, ipoint, tmp_O_priv, tmp_J_priv) &
|
||||||
|
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, tmp_O, tmp_J)
|
||||||
|
|
||||||
|
allocate(tmp_O_priv(n_points_final_grid), tmp_J_priv(n_points_final_grid,3))
|
||||||
|
tmp_O_priv = 0.d0
|
||||||
|
tmp_J_priv = 0.d0
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do i = 1, elec_beta_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp_O_priv(ipoint) = tmp_O_priv(ipoint) + mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,i,i)
|
||||||
|
tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,i,i)
|
||||||
|
tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,i,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
|
!$OMP DO
|
||||||
|
do i = elec_beta_num+1, elec_alpha_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp_O_priv(ipoint) = tmp_O_priv(ipoint) + 0.5d0 * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_J_priv(ipoint,1) = tmp_J_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,i)
|
||||||
|
tmp_J_priv(ipoint,2) = tmp_J_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,i)
|
||||||
|
tmp_J_priv(ipoint,3) = tmp_J_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
|
!$OMP CRITICAL
|
||||||
|
tmp_O = tmp_O + tmp_O_priv
|
||||||
|
tmp_J = tmp_J + tmp_J_priv
|
||||||
|
!$OMP END CRITICAL
|
||||||
|
|
||||||
|
deallocate(tmp_O_priv, tmp_J_priv)
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
allocate(tmp_M(n_points_final_grid,3), tmp_S(n_points_final_grid))
|
||||||
|
tmp_M = 0.d0
|
||||||
|
tmp_S = 0.d0
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE(i, j, ipoint, tmp_M_priv, tmp_S_priv) &
|
||||||
|
!$OMP SHARED(elec_beta_num, elec_alpha_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, tmp_M, tmp_S)
|
||||||
|
|
||||||
|
allocate(tmp_M_priv(n_points_final_grid,3), tmp_S_priv(n_points_final_grid))
|
||||||
|
tmp_M_priv = 0.d0
|
||||||
|
tmp_S_priv = 0.d0
|
||||||
|
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = 1, elec_beta_num
|
||||||
|
do j = 1, elec_beta_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
|
||||||
|
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = elec_beta_num+1, elec_alpha_num
|
||||||
|
do j = 1, elec_beta_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
|
||||||
|
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i)
|
||||||
|
|
||||||
|
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do i = elec_beta_num+1, elec_alpha_num
|
||||||
|
do j = elec_beta_num+1, elec_alpha_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_M_priv(ipoint,1) = tmp_M_priv(ipoint,1) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,2) = tmp_M_priv(ipoint,2) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
tmp_M_priv(ipoint,3) = tmp_M_priv(ipoint,3) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,j,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j)
|
||||||
|
|
||||||
|
tmp_S_priv(ipoint) = tmp_S_priv(ipoint) + 0.5d0 * int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
|
||||||
|
+ 0.5d0 * int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
|
||||||
|
+ 0.5d0 * int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO NOWAIT
|
||||||
|
|
||||||
|
!$OMP CRITICAL
|
||||||
|
tmp_M = tmp_M + tmp_M_priv
|
||||||
|
tmp_S = tmp_S + tmp_S_priv
|
||||||
|
!$OMP END CRITICAL
|
||||||
|
|
||||||
|
deallocate(tmp_M_priv, tmp_S_priv)
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
allocate(tmp(n_points_final_grid))
|
||||||
|
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
|
tmp_S(ipoint) = 2.d0 * (tmp_J(ipoint,1)*tmp_J(ipoint,1) + tmp_J(ipoint,2)*tmp_J(ipoint,2) + tmp_J(ipoint,3)*tmp_J(ipoint,3)) - tmp_S(ipoint)
|
||||||
|
|
||||||
|
tmp(ipoint) = final_weight_at_r_vector(ipoint) * ( tmp_O(ipoint) * tmp_S(ipoint) &
|
||||||
|
- 2.d0 * ( tmp_J(ipoint,1) * tmp_M(ipoint,1) &
|
||||||
|
+ tmp_J(ipoint,2) * tmp_M(ipoint,2) &
|
||||||
|
+ tmp_J(ipoint,3) * tmp_M(ipoint,3)))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
diag_three_elem_hf = diag_three_elem_hf - 2.d0 * (sum(tmp))
|
||||||
|
|
||||||
|
deallocate(tmp)
|
||||||
|
|
||||||
|
endif ! alpha/beta condition
|
||||||
|
|
||||||
|
endif ! three_body_h_tc
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|||||||
@ -1,80 +0,0 @@
|
|||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, tcscf_energy_3e_naive]
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, k
|
|
||||||
integer :: neu, ned, D(elec_num)
|
|
||||||
integer :: ii, jj, kk
|
|
||||||
integer :: si, sj, sk
|
|
||||||
double precision :: I_ijk, I_jki, I_kij, I_jik, I_ikj, I_kji
|
|
||||||
double precision :: I_tot
|
|
||||||
|
|
||||||
PROVIDE mo_l_coef mo_r_coef
|
|
||||||
|
|
||||||
neu = elec_alpha_num
|
|
||||||
ned = elec_beta_num
|
|
||||||
if (neu > 0) D(1:neu) = [(2*i-1, i = 1, neu)]
|
|
||||||
if (ned > 0) D(neu+1:neu+ned) = [(2*i, i = 1, ned)]
|
|
||||||
|
|
||||||
!print*, "D = "
|
|
||||||
!do i = 1, elec_num
|
|
||||||
! ii = (D(i) - 1) / 2 + 1
|
|
||||||
! si = mod(D(i), 2)
|
|
||||||
! print*, i, D(i), ii, si
|
|
||||||
!enddo
|
|
||||||
|
|
||||||
tcscf_energy_3e_naive = 0.d0
|
|
||||||
|
|
||||||
do i = 1, elec_num - 2
|
|
||||||
ii = (D(i) - 1) / 2 + 1
|
|
||||||
si = mod(D(i), 2)
|
|
||||||
|
|
||||||
do j = i + 1, elec_num - 1
|
|
||||||
jj = (D(j) - 1) / 2 + 1
|
|
||||||
sj = mod(D(j), 2)
|
|
||||||
|
|
||||||
do k = j + 1, elec_num
|
|
||||||
kk = (D(k) - 1) / 2 + 1
|
|
||||||
sk = mod(D(k), 2)
|
|
||||||
|
|
||||||
call give_integrals_3_body_bi_ort(ii, jj, kk, ii, jj, kk, I_ijk)
|
|
||||||
I_tot = I_ijk
|
|
||||||
|
|
||||||
if(sj==si .and. sk==sj) then
|
|
||||||
call give_integrals_3_body_bi_ort(ii, jj, kk, jj, kk, ii, I_jki)
|
|
||||||
I_tot += I_jki
|
|
||||||
endif
|
|
||||||
|
|
||||||
if(sk==si .and. si==sj) then
|
|
||||||
call give_integrals_3_body_bi_ort(ii, jj, kk, kk, ii, jj, I_kij)
|
|
||||||
I_tot += I_kij
|
|
||||||
endif
|
|
||||||
|
|
||||||
if(sj==si) then
|
|
||||||
call give_integrals_3_body_bi_ort(ii, jj, kk, jj, ii, kk, I_jik)
|
|
||||||
I_tot -= I_jik
|
|
||||||
endif
|
|
||||||
|
|
||||||
if(sk==sj) then
|
|
||||||
call give_integrals_3_body_bi_ort(ii, jj, kk, ii, kk, jj, I_ikj)
|
|
||||||
I_tot -= I_ikj
|
|
||||||
endif
|
|
||||||
|
|
||||||
if(sk==si) then
|
|
||||||
call give_integrals_3_body_bi_ort(ii, jj, kk, kk, jj, ii, I_kji)
|
|
||||||
I_tot -= I_kji
|
|
||||||
endif
|
|
||||||
|
|
||||||
tcscf_energy_3e_naive += I_tot
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
tcscf_energy_3e_naive = -tcscf_energy_3e_naive
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -1,189 +0,0 @@
|
|||||||
|
|
||||||
subroutine contrib_3e_diag_sss(i, j, k, integral)
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
! returns the pure same spin contribution to diagonal matrix element of 3e term
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i, j, k
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
|
|
||||||
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j >
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int) ! < i k j | j i k >
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i >
|
|
||||||
integral = direct_int + c_3_int + c_minus_3_int
|
|
||||||
|
|
||||||
! negative terms :: exchange contrib
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, k, i, j, exch_12_int)!!! < i k j | k i j > : E_12
|
|
||||||
|
|
||||||
integral += - exch_13_int - exch_23_int - exch_12_int
|
|
||||||
integral = -integral
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine contrib_3e_diag_soo(i,j,k,integral)
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: i,j,k
|
|
||||||
BEGIN_DOC
|
|
||||||
! returns the pure same spin contribution to diagonal matrix element of 3e term
|
|
||||||
END_DOC
|
|
||||||
double precision, intent(out) :: integral
|
|
||||||
double precision :: direct_int, exch_23_int
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int) ! < i k j | i k j >
|
|
||||||
call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)! < i k j | i j k > : E_23
|
|
||||||
integral = direct_int - exch_23_int
|
|
||||||
integral = -integral
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
subroutine give_aaa_contrib_bis(integral_aaa)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(out) :: integral_aaa
|
|
||||||
double precision :: integral
|
|
||||||
integer :: i,j,k
|
|
||||||
integral_aaa = 0.d0
|
|
||||||
do i = 1, elec_alpha_num
|
|
||||||
do j = i+1, elec_alpha_num
|
|
||||||
do k = j+1, elec_alpha_num
|
|
||||||
call contrib_3e_diag_sss(i,j,k,integral)
|
|
||||||
integral_aaa += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine give_aaa_contrib(integral_aaa)
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, k
|
|
||||||
double precision :: integral
|
|
||||||
double precision, intent(out) :: integral_aaa
|
|
||||||
|
|
||||||
integral_aaa = 0.d0
|
|
||||||
do i = 1, elec_alpha_num
|
|
||||||
do j = 1, elec_alpha_num
|
|
||||||
do k = 1, elec_alpha_num
|
|
||||||
call contrib_3e_diag_sss(i, j, k, integral)
|
|
||||||
integral_aaa += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
integral_aaa *= 1.d0/6.d0
|
|
||||||
|
|
||||||
return
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine give_aab_contrib(integral_aab)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(out) :: integral_aab
|
|
||||||
double precision :: integral
|
|
||||||
integer :: i,j,k
|
|
||||||
integral_aab = 0.d0
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = 1, elec_alpha_num
|
|
||||||
do k = 1, elec_alpha_num
|
|
||||||
call contrib_3e_diag_soo(i,j,k,integral)
|
|
||||||
integral_aab += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
integral_aab *= 0.5d0
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
subroutine give_aab_contrib_bis(integral_aab)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(out) :: integral_aab
|
|
||||||
double precision :: integral
|
|
||||||
integer :: i,j,k
|
|
||||||
integral_aab = 0.d0
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = 1, elec_alpha_num
|
|
||||||
do k = j+1, elec_alpha_num
|
|
||||||
call contrib_3e_diag_soo(i,j,k,integral)
|
|
||||||
integral_aab += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
subroutine give_abb_contrib(integral_abb)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(out) :: integral_abb
|
|
||||||
double precision :: integral
|
|
||||||
integer :: i,j,k
|
|
||||||
integral_abb = 0.d0
|
|
||||||
do i = 1, elec_alpha_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
call contrib_3e_diag_soo(i,j,k,integral)
|
|
||||||
integral_abb += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
integral_abb *= 0.5d0
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine give_abb_contrib_bis(integral_abb)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(out) :: integral_abb
|
|
||||||
double precision :: integral
|
|
||||||
integer :: i,j,k
|
|
||||||
integral_abb = 0.d0
|
|
||||||
do i = 1, elec_alpha_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do k = j+1, elec_beta_num
|
|
||||||
call contrib_3e_diag_soo(i,j,k,integral)
|
|
||||||
integral_abb += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine give_bbb_contrib_bis(integral_bbb)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(out) :: integral_bbb
|
|
||||||
double precision :: integral
|
|
||||||
integer :: i,j,k
|
|
||||||
integral_bbb = 0.d0
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = i+1, elec_beta_num
|
|
||||||
do k = j+1, elec_beta_num
|
|
||||||
call contrib_3e_diag_sss(i,j,k,integral)
|
|
||||||
integral_bbb += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine give_bbb_contrib(integral_bbb)
|
|
||||||
implicit none
|
|
||||||
double precision, intent(out) :: integral_bbb
|
|
||||||
double precision :: integral
|
|
||||||
integer :: i,j,k
|
|
||||||
integral_bbb = 0.d0
|
|
||||||
do i = 1, elec_beta_num
|
|
||||||
do j = 1, elec_beta_num
|
|
||||||
do k = 1, elec_beta_num
|
|
||||||
call contrib_3e_diag_sss(i,j,k,integral)
|
|
||||||
integral_bbb += integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
integral_bbb *= 1.d0/6.d0
|
|
||||||
end
|
|
||||||
|
|
||||||
|
|
||||||
@ -4,11 +4,9 @@ program write_ao_2e_tc_integ
|
|||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
PROVIDE j1e_type
|
|
||||||
PROVIDE j2e_type
|
|
||||||
|
|
||||||
print *, ' j1e_type = ', j1e_type
|
|
||||||
print *, ' j2e_type = ', j2e_type
|
print *, ' j2e_type = ', j2e_type
|
||||||
|
print *, ' j1e_type = ', j1e_type
|
||||||
|
print *, ' env_type = ', env_type
|
||||||
|
|
||||||
my_grid_becke = .True.
|
my_grid_becke = .True.
|
||||||
PROVIDE tc_grid1_a tc_grid1_r
|
PROVIDE tc_grid1_a tc_grid1_r
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user