PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 20:34:58 +01:00
Code Releases Activity

TC-SCF CLEANED

Browse Source
This commit is contained in:
Abdallah Ammar 2024-05-01 21:52:00 +02:00
parent c50018e8bd
commit d43d960b1a
27 changed files with 94 additions and 2796 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
plugins/local/bi_ort_ints/no_dressing.irp.f
View File

@ -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 ..."

2
plugins/local/non_hermit_dav/biorthog.irp.f
View File

@ -43,7 +43,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
! track & sort the real eigenvalues ! track & sort the real eigenvalues
n_good = 0 n_good = 0
thr = Im_thresh_tcscf thr = Im_thresh_tc
do i = 1, n do i = 1, n
if(dabs(WI(i)) .lt. thr) then if(dabs(WI(i)) .lt. thr) then
n_good += 1 n_good += 1

1
plugins/local/slater_tc/NEED
View File

@ -5,3 +5,4 @@ bi_ortho_mos
tc_keywords tc_keywords
non_hermit_dav non_hermit_dav
dav_general_mat dav_general_mat
tc_scf

0
plugins/local/tc_bi_ortho/symmetrized_3_e_int_prov.irp.f → plugins/local/slater_tc/symmetrized_3_e_int_prov.irp.f
View File

33
plugins/local/tc_bi_ortho/test_tc_fock.irp.f
View File

@ -24,44 +24,12 @@ program test_tc_fock
!call routine_2 !call routine_2
! call routine_3() ! call routine_3()
! call test_3e
call routine_tot call routine_tot
end end
! --- ! ---
subroutine test_3e
implicit none
double precision :: integral_aaa,integral_aab,integral_abb,integral_bbb,accu
double precision :: hmono, htwoe, hthree, htot
call htilde_mu_mat_bi_ortho_slow(ref_bitmask, ref_bitmask, N_int, hmono, htwoe, hthree, htot)
print*,'hmono = ',hmono
print*,'htwoe = ',htwoe
print*,'hthree= ',hthree
print*,'htot = ',htot
print*,''
print*,''
print*,'TC_one= ',tc_hf_one_e_energy
print*,'TC_two= ',TC_HF_two_e_energy
print*,'TC_3e = ',diag_three_elem_hf
print*,'TC_tot= ',TC_HF_energy
print*,''
print*,''
call give_aaa_contrib(integral_aaa)
print*,'integral_aaa = ',integral_aaa
call give_aab_contrib(integral_aab)
print*,'integral_aab = ',integral_aab
call give_abb_contrib(integral_abb)
print*,'integral_abb = ',integral_abb
call give_bbb_contrib(integral_bbb)
print*,'integral_bbb = ',integral_bbb
accu = integral_aaa + integral_aab + integral_abb + integral_bbb
print*,'accu = ',accu
print*,'delta = ',hthree - accu
end
subroutine routine_3() subroutine routine_3()
use bitmasks ! you need to include the bitmasks_module.f90 features use bitmasks ! you need to include the bitmasks_module.f90 features
@ -86,7 +54,6 @@ subroutine routine_3()
do i = 1, elec_num_tab(s1) do i = 1, elec_num_tab(s1)
do a = elec_num_tab(s1)+1, mo_num ! virtual do a = elec_num_tab(s1)+1, mo_num ! virtual
det_i = ref_bitmask det_i = ref_bitmask
call do_single_excitation(det_i, i, a, s1, i_ok) call do_single_excitation(det_i, i, a, s1, i_ok)
if(i_ok == -1) then if(i_ok == -1) then

48
plugins/local/tc_keywords/EZFIO.cfg
View File

@ -100,30 +100,12 @@ doc: If |true|, the states are re-ordered to match the input states
default: False default: False
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
[bi_ortho]
type: logical
doc: If |true|, the MO basis is assumed to be bi-orthonormal
interface: ezfio,provider,ocaml
default: True
[symetric_fock_tc] [symetric_fock_tc]
type: logical type: logical
doc: If |true|, using F+F^t as Fock TC doc: If |true|, using F+F^t as Fock TC
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[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
[selection_tc] [selection_tc]
type: integer type: integer
doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative
@ -160,30 +142,6 @@ doc: If |true|, maximize the overlap between orthogonalized left- and right eige
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.
[tcscf_algorithm]
type: character*(32)
doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS]
interface: ezfio,provider,ocaml
default: DIIS
[im_thresh_tcscf]
type: Threshold
doc: Thresholds on the Imag part of energy
interface: ezfio,provider,ocaml
default: 1.e-7
[test_cycle_tc] [test_cycle_tc]
type: logical type: logical
doc: If |true|, the integrals of the three-body jastrow are computed with cycles doc: If |true|, the integrals of the three-body jastrow are computed with cycles
@ -304,3 +262,9 @@ doc: If |true|, more calc but less mem
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[im_thresh_tc]
type: Threshold
doc: Thresholds on the Imag part of TC energy
interface: ezfio,provider,ocaml
default: 1.e-7

7
plugins/local/tc_keywords/tc_keywords.irp.f
View File

@ -1,7 +0,0 @@
program tc_keywords
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
end

30
plugins/local/tc_scf/EZFIO.cfg
View File

@ -9,3 +9,33 @@ 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.
[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

107
plugins/local/tc_scf/fock_hermit.irp.f
View File

@ -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
! ---

18
plugins/local/tc_scf/fock_tc.irp.f
View File

@ -110,8 +110,6 @@ 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) & call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
@ -122,13 +120,6 @@ BEGIN_PROVIDER [double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num)]
Fock_matrix_tc_mo_alpha += fock_3e_mo_a Fock_matrix_tc_mo_alpha += fock_3e_mo_a
endif 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) )
endif
END_PROVIDER END_PROVIDER
! --- ! ---
@ -142,8 +133,6 @@ 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) & 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) ) , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
@ -152,13 +141,6 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
Fock_matrix_tc_mo_beta += fock_3e_mo_b Fock_matrix_tc_mo_beta += fock_3e_mo_b
endif 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) )
endif
END_PROVIDER END_PROVIDER
! --- ! ---

19
plugins/local/tc_scf/fock_tc_mo_tot.irp.f
View File

@ -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

771
plugins/local/tc_scf/fock_three_hermit.irp.f
View File

@ -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

391
plugins/local/tc_scf/integrals_in_r_stuff.irp.f
View File

@ -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
! ---

318
plugins/local/tc_scf/jast_schmos_90.irp.f
View File

@ -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

69
plugins/local/tc_scf/plot_j_schMos.irp.f
View File

@ -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

59
plugins/local/tc_scf/print_fit_param.irp.f
View File

@ -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
! ---

55
plugins/local/tc_scf/print_tcscf_energy.irp.f
View File

@ -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
! ---

4
plugins/local/tc_scf/rh_tcscf_diis.irp.f
View File

@ -61,7 +61,7 @@ subroutine rh_tcscf_diis()
etc_tot = TC_HF_energy etc_tot = TC_HF_energy
etc_1e = TC_HF_one_e_energy etc_1e = TC_HF_one_e_energy
etc_2e = TC_HF_two_e_energy etc_2e = TC_HF_two_e_energy
etc_3e = diag_three_elem_hf etc_3e = TC_HF_three_e_energy
!tc_grad = grad_non_hermit !tc_grad = grad_non_hermit
er_DIIS = maxval(abs(FQS_SQF_mo)) er_DIIS = maxval(abs(FQS_SQF_mo))
e_delta = dabs(etc_tot - e_save) e_delta = dabs(etc_tot - e_save)
@ -189,7 +189,7 @@ subroutine rh_tcscf_diis()
etc_tot = TC_HF_energy etc_tot = TC_HF_energy
etc_1e = TC_HF_one_e_energy etc_1e = TC_HF_one_e_energy
etc_2e = TC_HF_two_e_energy etc_2e = TC_HF_two_e_energy
etc_3e = diag_three_elem_hf etc_3e = TC_HF_three_e_energy
!tc_grad = grad_non_hermit !tc_grad = grad_non_hermit
er_DIIS = maxval(abs(FQS_SQF_mo)) er_DIIS = maxval(abs(FQS_SQF_mo))
e_delta = dabs(etc_tot - e_save) e_delta = dabs(etc_tot - e_save)

129
plugins/local/tc_scf/rh_tcscf_simple.irp.f
View File

@ -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
! ---

369
plugins/local/tc_scf/rotate_tcscf_orbitals.irp.f
View File

@ -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
! ---

91
plugins/local/tc_scf/tc_petermann_factor.irp.f
View File

@ -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
! ---

23
plugins/local/tc_scf/tc_scf.irp.f
View File

@ -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 orthonormalize_mos()
if(tcscf_algorithm == 'DIIS') then
call rh_tcscf_diis() 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

16
plugins/local/tc_scf/tc_scf_dm.irp.f
View File

@ -10,17 +10,9 @@ 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 PROVIDE mo_l_coef mo_r_coef
TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta 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,17 +27,9 @@ 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 PROVIDE mo_l_coef mo_r_coef
TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha 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

16
plugins/local/tc_scf/tc_scf_energy.irp.f
View File

@ -1,7 +1,8 @@
BEGIN_PROVIDER [ double precision, TC_HF_energy ] BEGIN_PROVIDER [double precision, TC_HF_energy ]
&BEGIN_PROVIDER [ double precision, TC_HF_one_e_energy] &BEGIN_PROVIDER [double precision, TC_HF_one_e_energy ]
&BEGIN_PROVIDER [ double precision, TC_HF_two_e_energy] &BEGIN_PROVIDER [double precision, TC_HF_two_e_energy ]
&BEGIN_PROVIDER [double precision, TC_HF_three_e_energy]
BEGIN_DOC BEGIN_DOC
! TC Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. ! TC Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
@ -27,8 +28,13 @@
enddo enddo
enddo enddo
TC_HF_energy += TC_HF_one_e_energy + TC_HF_two_e_energy if((three_body_h_tc .eq. .False.) .and. (.not. noL_standard)) then
TC_HF_energy += diag_three_elem_hf TC_HF_three_e_energy = 0.d0
else
TC_HF_three_e_energy = noL_0e
endif
TC_HF_energy += TC_HF_one_e_energy + TC_HF_two_e_energy + TC_HF_three_e_energy
END_PROVIDER END_PROVIDER

80
plugins/local/tc_scf/tcscf_energy_naive.irp.f
View File

@ -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
! ---

189
plugins/local/tc_scf/three_e_energy_bi_ortho.irp.f
View File

@ -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

6
plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f
View File

@ -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