From da8eac81e01e9ee558351195aba1f964ed5fbc0b Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Wed, 1 May 2024 21:52:00 +0200 Subject: [PATCH] TC-SCF CLEANED --- plugins/local/bi_ort_ints/no_dressing.irp.f | 7 +- plugins/local/tc_scf/EZFIO.cfg | 36 + plugins/local/tc_scf/fock_hermit.irp.f | 107 --- plugins/local/tc_scf/fock_tc.irp.f | 40 +- plugins/local/tc_scf/fock_tc_mo_tot.irp.f | 19 +- plugins/local/tc_scf/fock_three_hermit.irp.f | 771 ------------------ .../local/tc_scf/integrals_in_r_stuff.irp.f | 391 --------- plugins/local/tc_scf/jast_schmos_90.irp.f | 318 -------- plugins/local/tc_scf/plot_j_schMos.irp.f | 69 -- plugins/local/tc_scf/print_fit_param.irp.f | 59 -- plugins/local/tc_scf/print_tcscf_energy.irp.f | 55 -- plugins/local/tc_scf/rh_tcscf_simple.irp.f | 129 --- .../local/tc_scf/rotate_tcscf_orbitals.irp.f | 369 --------- .../local/tc_scf/tc_petermann_factor.irp.f | 91 --- plugins/local/tc_scf/tc_scf.irp.f | 25 +- plugins/local/tc_scf/tc_scf_dm.irp.f | 24 +- plugins/local/tc_scf/tc_scf_energy.irp.f | 423 ++++++++++ plugins/local/tc_scf/tcscf_energy_naive.irp.f | 80 -- .../tc_scf/three_e_energy_bi_ortho.irp.f | 189 ----- .../local/tc_scf/write_ao_2e_tc_integ.irp.f | 6 +- 20 files changed, 502 insertions(+), 2706 deletions(-) delete mode 100644 plugins/local/tc_scf/fock_hermit.irp.f delete mode 100644 plugins/local/tc_scf/fock_three_hermit.irp.f delete mode 100644 plugins/local/tc_scf/integrals_in_r_stuff.irp.f delete mode 100644 plugins/local/tc_scf/jast_schmos_90.irp.f delete mode 100644 plugins/local/tc_scf/plot_j_schMos.irp.f delete mode 100644 plugins/local/tc_scf/print_fit_param.irp.f delete mode 100644 plugins/local/tc_scf/print_tcscf_energy.irp.f delete mode 100644 plugins/local/tc_scf/rh_tcscf_simple.irp.f delete mode 100644 plugins/local/tc_scf/rotate_tcscf_orbitals.irp.f delete mode 100644 plugins/local/tc_scf/tc_petermann_factor.irp.f delete mode 100644 plugins/local/tc_scf/tcscf_energy_naive.irp.f delete mode 100644 plugins/local/tc_scf/three_e_energy_bi_ortho.irp.f diff --git a/plugins/local/bi_ort_ints/no_dressing.irp.f b/plugins/local/bi_ort_ints/no_dressing.irp.f index bd225274..721ac0f8 100644 --- a/plugins/local/bi_ort_ints/no_dressing.irp.f +++ b/plugins/local/bi_ort_ints/no_dressing.irp.f @@ -322,6 +322,12 @@ END_PROVIDER BEGIN_PROVIDER [double precision, noL_0e] + BEGIN_DOC + ! + ! < Phi_left | L | Phi_right > + ! + END_DOC + implicit none integer :: i, j, k, ipoint 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_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:) - call wall_time(t0) print*, " Providing noL_0e ..." diff --git a/plugins/local/tc_scf/EZFIO.cfg b/plugins/local/tc_scf/EZFIO.cfg index 510c777c..6820a8b0 100644 --- a/plugins/local/tc_scf/EZFIO.cfg +++ b/plugins/local/tc_scf/EZFIO.cfg @@ -9,3 +9,39 @@ doc: If |true|, tc-scf has converged interface: ezfio,provider,ocaml 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 + diff --git a/plugins/local/tc_scf/fock_hermit.irp.f b/plugins/local/tc_scf/fock_hermit.irp.f deleted file mode 100644 index 5a51b324..00000000 --- a/plugins/local/tc_scf/fock_hermit.irp.f +++ /dev/null @@ -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 - -! --- - diff --git a/plugins/local/tc_scf/fock_tc.irp.f b/plugins/local/tc_scf/fock_tc.irp.f index 508f3cd7..16bb5c87 100644 --- a/plugins/local/tc_scf/fock_tc.irp.f +++ b/plugins/local/tc_scf/fock_tc.irp.f @@ -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, 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) ) - - 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) ) + 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 END_PROVIDER @@ -142,21 +133,12 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] implicit none 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) ) - - 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) ) + 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 END_PROVIDER diff --git a/plugins/local/tc_scf/fock_tc_mo_tot.irp.f b/plugins/local/tc_scf/fock_tc_mo_tot.irp.f index 2df2421e..fd490af6 100644 --- a/plugins/local/tc_scf/fock_tc_mo_tot.irp.f +++ b/plugins/local/tc_scf/fock_tc_mo_tot.irp.f @@ -132,7 +132,7 @@ enddo endif - if(no_oa_or_av_opt)then + if(no_oa_or_av_opt) then do i = 1, n_act_orb iorb = list_act(i) do j = 1, n_inact_orb @@ -153,8 +153,21 @@ enddo endif - if(.not.bi_ortho .and. three_body_h_tc)then - Fock_matrix_tc_mo_tot += fock_3_mat + if(tc_Brillouin_Right) then + + 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 END_PROVIDER diff --git a/plugins/local/tc_scf/fock_three_hermit.irp.f b/plugins/local/tc_scf/fock_three_hermit.irp.f deleted file mode 100644 index 00d47fae..00000000 --- a/plugins/local/tc_scf/fock_three_hermit.irp.f +++ /dev/null @@ -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) ! - 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 - - - - - diff --git a/plugins/local/tc_scf/integrals_in_r_stuff.irp.f b/plugins/local/tc_scf/integrals_in_r_stuff.irp.f deleted file mode 100644 index 3ce85a97..00000000 --- a/plugins/local/tc_scf/integrals_in_r_stuff.irp.f +++ /dev/null @@ -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 - -! --- - diff --git a/plugins/local/tc_scf/jast_schmos_90.irp.f b/plugins/local/tc_scf/jast_schmos_90.irp.f deleted file mode 100644 index 5c5e625f..00000000 --- a/plugins/local/tc_scf/jast_schmos_90.irp.f +++ /dev/null @@ -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 diff --git a/plugins/local/tc_scf/plot_j_schMos.irp.f b/plugins/local/tc_scf/plot_j_schMos.irp.f deleted file mode 100644 index eda0dd25..00000000 --- a/plugins/local/tc_scf/plot_j_schMos.irp.f +++ /dev/null @@ -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 diff --git a/plugins/local/tc_scf/print_fit_param.irp.f b/plugins/local/tc_scf/print_fit_param.irp.f deleted file mode 100644 index e62f0dde..00000000 --- a/plugins/local/tc_scf/print_fit_param.irp.f +++ /dev/null @@ -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 - -! --- - diff --git a/plugins/local/tc_scf/print_tcscf_energy.irp.f b/plugins/local/tc_scf/print_tcscf_energy.irp.f deleted file mode 100644 index 6f9afd9a..00000000 --- a/plugins/local/tc_scf/print_tcscf_energy.irp.f +++ /dev/null @@ -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 - -! --- - diff --git a/plugins/local/tc_scf/rh_tcscf_simple.irp.f b/plugins/local/tc_scf/rh_tcscf_simple.irp.f deleted file mode 100644 index 2c2cf2c2..00000000 --- a/plugins/local/tc_scf/rh_tcscf_simple.irp.f +++ /dev/null @@ -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 - -! --- - diff --git a/plugins/local/tc_scf/rotate_tcscf_orbitals.irp.f b/plugins/local/tc_scf/rotate_tcscf_orbitals.irp.f deleted file mode 100644 index 0f2663e5..00000000 --- a/plugins/local/tc_scf/rotate_tcscf_orbitals.irp.f +++ /dev/null @@ -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 - -! --- diff --git a/plugins/local/tc_scf/tc_petermann_factor.irp.f b/plugins/local/tc_scf/tc_petermann_factor.irp.f deleted file mode 100644 index 14fff898..00000000 --- a/plugins/local/tc_scf/tc_petermann_factor.irp.f +++ /dev/null @@ -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 - -! --- - diff --git a/plugins/local/tc_scf/tc_scf.irp.f b/plugins/local/tc_scf/tc_scf.irp.f index ee8e8dad..f099b90e 100644 --- a/plugins/local/tc_scf/tc_scf.irp.f +++ b/plugins/local/tc_scf/tc_scf.irp.f @@ -10,13 +10,10 @@ program tc_scf integer :: i logical :: good_angles - PROVIDE j1e_type - PROVIDE j2e_type - PROVIDE tcscf_algorithm - print *, ' TC-SCF with:' - print *, ' j1e_type = ', j1e_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' @@ -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_a_grid, 'angular external grid over') - if(tc_integ_type .eq. "numeric") then my_extra_grid_becke = .True. 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') endif - !call create_guess() - !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 + call rh_tcscf_diis() PROVIDE Fock_matrix_tc_diag_mo_tot print*, ' Eigenvalues:' @@ -59,14 +45,11 @@ program tc_scf print*, i, Fock_matrix_tc_diag_mo_tot(i) enddo - ! TODO - ! rotate angles in separate code only if necessary - if(minimize_lr_angles)then + if(minimize_lr_angles) then call minimize_tc_orb_angles() endif call print_energy_and_mos(good_angles) - write(json_unit,json_array_close_fmtx) call json_close diff --git a/plugins/local/tc_scf/tc_scf_dm.irp.f b/plugins/local/tc_scf/tc_scf_dm.irp.f index bf31a4a1..5d25fce2 100644 --- a/plugins/local/tc_scf/tc_scf_dm.irp.f +++ b/plugins/local/tc_scf/tc_scf_dm.irp.f @@ -10,16 +10,8 @@ BEGIN_PROVIDER [double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num) implicit none - if(bi_ortho) then - - 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 + PROVIDE mo_l_coef mo_r_coef + TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta END_PROVIDER @@ -35,16 +27,8 @@ BEGIN_PROVIDER [double precision, TCSCF_density_matrix_ao_alpha, (ao_num, ao_num implicit none - if(bi_ortho) then - - 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 + PROVIDE mo_l_coef mo_r_coef + TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha END_PROVIDER diff --git a/plugins/local/tc_scf/tc_scf_energy.irp.f b/plugins/local/tc_scf/tc_scf_energy.irp.f index 0266c605..c9366195 100644 --- a/plugins/local/tc_scf/tc_scf_energy.irp.f +++ b/plugins/local/tc_scf/tc_scf_energy.irp.f @@ -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 + +! --- + diff --git a/plugins/local/tc_scf/tcscf_energy_naive.irp.f b/plugins/local/tc_scf/tcscf_energy_naive.irp.f deleted file mode 100644 index 82bb8799..00000000 --- a/plugins/local/tc_scf/tcscf_energy_naive.irp.f +++ /dev/null @@ -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 - -! --- - diff --git a/plugins/local/tc_scf/three_e_energy_bi_ortho.irp.f b/plugins/local/tc_scf/three_e_energy_bi_ortho.irp.f deleted file mode 100644 index 0c9ebbd7..00000000 --- a/plugins/local/tc_scf/three_e_energy_bi_ortho.irp.f +++ /dev/null @@ -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 - - diff --git a/plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f b/plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f index 7ce57578..ec5167d1 100644 --- a/plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f +++ b/plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f @@ -4,11 +4,9 @@ program write_ao_2e_tc_integ implicit none - PROVIDE j1e_type - PROVIDE j2e_type - - print *, ' j1e_type = ', j1e_type print *, ' j2e_type = ', j2e_type + print *, ' j1e_type = ', j1e_type + print *, ' env_type = ', env_type my_grid_becke = .True. PROVIDE tc_grid1_a tc_grid1_r