From 51cf96a506a87c807b88f705e68f6522142d256c Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 7 Apr 2020 11:01:24 +0200 Subject: [PATCH] added mu_of_r --- src/mu_of_r/EZFIO.cfg | 18 ++ src/mu_of_r/NEED | 1 + src/mu_of_r/README.rst | 4 + src/mu_of_r/basis_def.irp.f | 106 +++++++++ src/mu_of_r/example.irp.f | 202 +++++++++++++++++ src/mu_of_r/f_hf_utils.irp.f | 138 ++++++++++++ src/mu_of_r/f_psi_i_a_v_utils.irp.f | 313 +++++++++++++++++++++++++++ src/mu_of_r/f_psi_old.irp.f | 39 ++++ src/mu_of_r/f_psi_utils.irp.f | 151 +++++++++++++ src/mu_of_r/f_val_general.irp.f | 90 ++++++++ src/mu_of_r/mu_of_r_conditions.irp.f | 150 +++++++++++++ src/mu_of_r/test_proj_op.irp.f | 21 ++ 12 files changed, 1233 insertions(+) create mode 100644 src/mu_of_r/EZFIO.cfg create mode 100644 src/mu_of_r/NEED create mode 100644 src/mu_of_r/README.rst create mode 100644 src/mu_of_r/basis_def.irp.f create mode 100644 src/mu_of_r/example.irp.f create mode 100644 src/mu_of_r/f_hf_utils.irp.f create mode 100644 src/mu_of_r/f_psi_i_a_v_utils.irp.f create mode 100644 src/mu_of_r/f_psi_old.irp.f create mode 100644 src/mu_of_r/f_psi_utils.irp.f create mode 100644 src/mu_of_r/f_val_general.irp.f create mode 100644 src/mu_of_r/mu_of_r_conditions.irp.f create mode 100644 src/mu_of_r/test_proj_op.irp.f diff --git a/src/mu_of_r/EZFIO.cfg b/src/mu_of_r/EZFIO.cfg new file mode 100644 index 00000000..5677b3ab --- /dev/null +++ b/src/mu_of_r/EZFIO.cfg @@ -0,0 +1,18 @@ +[mu_of_r_disk] +type: double precision +doc: array of the values of mu(r) +interface: ezfio +size: (becke_numerical_grid.n_points_final_grid,determinants.n_states) + +[mu_of_r_potential] +type: character*(32) +doc: type of potential for the mu(r) interaction: can be [ hf| cas_ful | cas_truncated] +interface: ezfio, provider, ocaml +default: hf + +[io_mu_of_r] +type: Disk_access +doc: Read/Write the mu(r) from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + diff --git a/src/mu_of_r/NEED b/src/mu_of_r/NEED new file mode 100644 index 00000000..22bb3ebe --- /dev/null +++ b/src/mu_of_r/NEED @@ -0,0 +1 @@ +cas_based_on_top diff --git a/src/mu_of_r/README.rst b/src/mu_of_r/README.rst new file mode 100644 index 00000000..b5234e7e --- /dev/null +++ b/src/mu_of_r/README.rst @@ -0,0 +1,4 @@ +================== +mu_of_r_definition +================== + diff --git a/src/mu_of_r/basis_def.irp.f b/src/mu_of_r/basis_def.irp.f new file mode 100644 index 00000000..4da27cb0 --- /dev/null +++ b/src/mu_of_r/basis_def.irp.f @@ -0,0 +1,106 @@ + + + BEGIN_PROVIDER [integer, n_occ_val_orb_for_hf,(2)] +&BEGIN_PROVIDER [integer, n_max_occ_val_orb_for_hf] + implicit none + BEGIN_DOC + ! Number of OCCUPIED VALENCE ORBITALS for each spin to build the f_{HF}(r_1,r_2) function + ! + ! This is typically elec_alpha_num - n_core_orb for alpha electrons and elec_beta_num - n_core_orb for beta electrons + ! + ! This determines the size of the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 + END_DOC + integer :: i + n_occ_val_orb_for_hf = 0 + ! You browse the ALPHA ELECTRONS and check if its not a CORE ORBITAL + do i = 1, elec_alpha_num + if( trim(mo_class(i))=="Inactive" & + .or. trim(mo_class(i))=="Active" & + .or. trim(mo_class(i))=="Virtual" )then + n_occ_val_orb_for_hf(1) +=1 + endif + enddo + + ! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL + do i = 1, elec_beta_num + if( trim(mo_class(i))=="Inactive" & + .or. trim(mo_class(i))=="Active" & + .or. trim(mo_class(i))=="Virtual" )then + n_occ_val_orb_for_hf(2) +=1 + endif + enddo + n_max_occ_val_orb_for_hf = maxval(n_occ_val_orb_for_hf) + +END_PROVIDER + +BEGIN_PROVIDER [integer, list_valence_orb_for_hf, (n_max_occ_val_orb_for_hf,2)] + implicit none + BEGIN_DOC + ! List of OCCUPIED valence orbitals for each spin to build the f_{HF}(r_1,r_2) function + ! + ! This corresponds to ALL OCCUPIED orbitals in the HF wave function, except those defined as "core" + ! + ! This determines the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 + END_DOC + integer :: i,j + j = 0 + ! You browse the ALPHA ELECTRONS and check if its not a CORE ORBITAL + do i = 1, elec_alpha_num + if( trim(mo_class(i))=="Inactive" & + .or. trim(mo_class(i))=="Active" & + .or. trim(mo_class(i))=="Virtual" )then + j +=1 + list_valence_orb_for_hf(j,1) = i + endif + enddo + + j = 0 + ! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL + do i = 1, elec_beta_num + if( trim(mo_class(i))=="Inactive" & + .or. trim(mo_class(i))=="Active" & + .or. trim(mo_class(i))=="Virtual" )then + j +=1 + list_valence_orb_for_hf(j,2) = i + endif + enddo + +END_PROVIDER + +BEGIN_PROVIDER [integer, n_basis_orb] + implicit none + BEGIN_DOC + ! Defines the number of orbitals you will use to explore the basis set + ! + ! This determines the size of the space \mathcal{B} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 + ! + ! It corresponds to all MOs except those defined as "deleted" + END_DOC + n_basis_orb = n_all_but_del_orb +END_PROVIDER + +BEGIN_PROVIDER [integer, list_basis, (n_basis_orb)] + implicit none + BEGIN_DOC + ! Defines the set of orbitals you will use to explore the basis set + ! + ! This determines the space \mathcal{B} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 + ! + ! It corresponds to all MOs except those defined as "deleted" + END_DOC + integer :: i + do i = 1, n_all_but_del_orb + list_basis(i) = list_all_but_del_orb(i) + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_final_grid)] + implicit none + integer :: ipoint,i,ii + do ipoint = 1, n_points_final_grid + do i = 1, n_basis_orb + ii = list_basis(i) + basis_mos_in_r_array(i,ipoint) = mos_in_r_array(ii,ipoint) + enddo + enddo +END_PROVIDER diff --git a/src/mu_of_r/example.irp.f b/src/mu_of_r/example.irp.f new file mode 100644 index 00000000..b32f4f4f --- /dev/null +++ b/src/mu_of_r/example.irp.f @@ -0,0 +1,202 @@ + +subroutine test_f_HF_valence_ab + implicit none + BEGIN_DOC +! routine to test the function f_HF(r1,r2) +! +! the integral over r1,r2 should be equal to the alpha/beta interaction of HF determinant + END_DOC + integer :: ipoint,i,j,i_i,j_j,jpoint + double precision :: accu_val,accu_ful, weight1,weight2, r1(3),integral_psi_val,integral_psi,r2(3),two_bod + accu_2 = 0.d0 + ! You compute the coulomb repulsion between alpha-beta electrons for HF + do i = 1, n_occ_val_orb_for_hf(1) + i_i = list_valence_orb_for_hf(i,1) + do j = 1, n_occ_val_orb_for_hf(2) + j_j = list_valence_orb_for_hf(j,2) + accu_2 += mo_two_e_integrals_jj(j_j,i_i) + enddo + enddo + print*,'' + print*,'' + print*,'' + print*,'**************************' + print*,'**************************' + print*,'Routine to test the f_HF(r1,r2) function' + print*,'**************************' + print*,'' + print*,'' + print*,'' + print*,'**************************' + print*,' = ',accu_2 + print*,'**************************' + + print*,'semi analytical form ' + accu_val = 0.d0 + ! You integrate on r2 the analytical integral over r1 of f_HF(r1,r2) + do ipoint = 1, n_points_final_grid + weight1 =final_weight_at_r_vector(ipoint) + r2(1) = final_grid_points(1,ipoint) + r2(2) = final_grid_points(2,ipoint) + r2(3) = final_grid_points(3,ipoint) + call integral_f_HF_valence_ab(r2,integral_psi_val) + accu_val += integral_psi_val * weight1 + enddo + print*,'**************************' + ! Should give you the alpha-beta repulsion of HF, excluding core contributions, + print*,'int dr1 dr2 f_HF(r1,r2) = ',accu_val + double precision :: accu_2 + + + print*,'pure numerical form (might take quite some time as it grows as N_g^2 * N_e^2 * N_b^2 ...)' + ! You integrate brut force on r1 and r2 + accu_val = 0.d0 + do jpoint = 1, n_points_final_grid + weight1 =final_weight_at_r_vector(jpoint) + r1(1) = final_grid_points(1,jpoint) + r1(2) = final_grid_points(2,jpoint) + r1(3) = final_grid_points(3,jpoint) + do ipoint = 1, n_points_final_grid + weight2 =final_weight_at_r_vector(ipoint) + r2(1) = final_grid_points(1,ipoint) + r2(2) = final_grid_points(2,ipoint) + r2(3) = final_grid_points(3,ipoint) + call f_HF_valence_ab(r1,r2,integral_psi_val,two_bod) + accu_val += integral_psi_val * weight1 * weight2 + enddo + enddo + print*,'int dr1 dr2 f_HF(r1,r2) = ',accu_val + + + print*,'**************************' + print*,'**************************' + print*,'**************************' + accu_val = 0.d0 + r1 = 0.d0 + r1(1) = 0.5d0 + print*,'r1 = ',r1 + ! You compute the integral over r2 of f_HF(r1,r2) + call integral_f_HF_valence_ab(r1,integral_psi) + do ipoint = 1, n_points_final_grid + weight1 =final_weight_at_r_vector(ipoint) + r2(1) = final_grid_points(1,ipoint) + r2(2) = final_grid_points(2,ipoint) + r2(3) = final_grid_points(3,ipoint) + call f_HF_valence_ab(r1,r2,integral_psi_val,two_bod) + accu_val += integral_psi_val * weight1 + enddo + print*,'int dr2 f_HF(r1,r2) = ',integral_psi + print*,'analytical form = ',accu_val + print*,'**************************' +end + + + +subroutine test_f_ii_valence_ab + implicit none + BEGIN_DOC +! routine to test the function f_ii(r1,r2) +! +! it should be the same that f_HF(r1,r2) only for inactive orbitals + END_DOC + integer :: ipoint + double precision :: accu_f, accu_n2, weight, r1(3),r2(3) + double precision :: accu_f_on_top + double precision :: f_HF_val_ab,two_bod_dens_hf,f_ii_val_ab,two_bod_dens_ii + accu_f = 0.d0 + accu_n2 = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + r2 = r1 + call f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens_hf) + call give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens_ii) + accu_f += dabs(f_HF_val_ab - f_ii_val_ab) * weight + accu_n2+= dabs(two_bod_dens_hf - two_bod_dens_ii) * weight + accu_f_on_top += dabs(two_bod_dens_hf) * weight + enddo + print*,'**************************' + print*,'' + print*,'accu_f = ',accu_f + print*,'accu_n2 = ',accu_n2 + print*,'' + print*,'accu_f_on_top = ',accu_f_on_top +end + + +subroutine test_f_ia_valence_ab + implicit none + BEGIN_DOC +! routine to test the function f_ii(r1,r2), f_ia(r1,r2) and f_aa(r1,r2) + END_DOC + integer :: ipoint,istate + double precision :: accu_f, accu_n2, weight, r1(3),r2(3) + double precision :: accu_f_on_top + double precision :: f_ref,f_comp,on_top_ref,on_top_comp + double precision :: f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa + double precision :: accu + accu_f = 0.d0 + accu_n2 = 0.d0 + accu = 0.d0 + istate = 1 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + r1(1) = final_grid_points(1,ipoint) + r1(2) = final_grid_points(2,ipoint) + r1(3) = final_grid_points(3,ipoint) + r2 = r1 + call give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens_ii) + call give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens_ia,istate) + call give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens_aa,istate) + f_ref = f_psi_cas_ab_old(ipoint,istate) + f_comp = f_ii_val_ab + f_ia_val_ab + f_aa_val_ab + on_top_ref = total_cas_on_top_density(ipoint,istate) + on_top_comp= two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa + accu_f += dabs(f_ref - f_comp) * weight + accu_n2+= dabs(on_top_ref - on_top_comp) * weight + accu += f_ref * weight + enddo + print*,'**************************' + print*,'' + print*,'accu_f = ',accu_f + print*,'accu_n2 = ',accu_n2 + print*,'' + print*,'accu = ',accu + +end + +subroutine test_f_ii_ia_aa_valence_ab + implicit none + BEGIN_DOC +! routine to test the function f_Psi(r1,r2) based on core/inactive/active orbitals + END_DOC + integer :: ipoint,istate + double precision :: accu_f, accu_n2, weight, r1(3),r2(3) + double precision :: accu_f_on_top + double precision :: f_ref,f_comp,on_top_ref,on_top_comp + double precision :: f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa + double precision :: accu + accu_f = 0.d0 + accu_n2 = 0.d0 + accu = 0.d0 + istate = 1 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + f_ref = f_psi_cas_ab(ipoint,istate) + f_comp = f_psi_cas_ab_old(ipoint,istate) + on_top_ref = total_cas_on_top_density(ipoint,istate) + on_top_comp= on_top_cas_mu_r(ipoint,istate) + accu_f += dabs(f_ref - f_comp) * weight + accu_n2+= dabs(on_top_ref - on_top_comp) * weight + accu += f_ref * weight + enddo + print*,'**************************' + print*,'' + print*,'accu_f = ',accu_f + print*,'accu_n2 = ',accu_n2 + print*,'' + print*,'accu = ',accu + +end diff --git a/src/mu_of_r/f_hf_utils.irp.f b/src/mu_of_r/f_hf_utils.irp.f new file mode 100644 index 00000000..b89dda18 --- /dev/null +++ b/src/mu_of_r/f_hf_utils.irp.f @@ -0,0 +1,138 @@ +BEGIN_PROVIDER [double precision, two_e_int_hf_f, (n_basis_orb,n_basis_orb,n_max_occ_val_orb_for_hf,n_max_occ_val_orb_for_hf)] + implicit none + BEGIN_DOC +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{HF}(r_1,r_2) +! +! two_e_int_hf_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" + END_DOC + integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n + double precision :: get_two_e_integral + do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1 + m = list_valence_orb_for_hf(orb_m,1) + do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2 + n = list_valence_orb_for_hf(orb_n,1) + do orb_i = 1, n_basis_orb ! electron 1 + i = list_basis(orb_i) + do orb_j = 1, n_basis_orb ! electron 2 + j = list_basis(orb_j) + ! 2 1 2 1 + two_e_int_hf_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) + enddo + enddo + enddo + enddo +END_PROVIDER + +subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens) + implicit none + BEGIN_DOC +! f_HF_val_ab(r1,r2) = function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) +! +! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931 2937) +! +! two_bod_dens = on-top pair density of the HF wave function +! +! < HF | wee_{\alpha\beta} | HF > = \int (r1,r2) f_HF_ab(r1,r2) excluding all contributions from "core" "electrons" + END_DOC + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out):: f_HF_val_ab,two_bod_dens + integer :: i,j,m,n,i_m,i_n + integer :: i_i,i_j + double precision :: mo_two_e_integral + double precision, allocatable :: mos_array_r1(:) + double precision, allocatable :: mos_array_r2(:) + double precision, allocatable :: mos_array_valence_r1(:),mos_array_valence_r2(:) + double precision, allocatable :: mos_array_valence_hf_r1(:),mos_array_valence_hf_r2(:) + double precision :: get_two_e_integral + allocate(mos_array_valence_r1(n_basis_orb) , mos_array_valence_r2(n_basis_orb), mos_array_r1(mo_num), mos_array_r2(mo_num)) + allocate(mos_array_valence_hf_r1(n_occ_val_orb_for_hf(1)) , mos_array_valence_hf_r2(n_occ_val_orb_for_hf(2)) ) + ! You get all orbitals in r_1 and r_2 + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + ! You extract the occupied ALPHA/BETA orbitals belonging to the space \mathcal{A} + do i_m = 1, n_occ_val_orb_for_hf(1) + mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1)) + enddo + do i_m = 1, n_occ_val_orb_for_hf(2) + mos_array_valence_hf_r2(i_m) = mos_array_r2(list_valence_orb_for_hf(i_m,2)) + enddo + + ! You extract the orbitals belonging to the space \mathcal{B} + do i_m = 1, n_basis_orb + mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m)) + mos_array_valence_r2(i_m) = mos_array_r2(list_basis(i_m)) + enddo + + + f_HF_val_ab = 0.d0 + two_bod_dens = 0.d0 + ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space + do m = 1, n_occ_val_orb_for_hf(1)! electron 1 + ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space + do n = 1, n_occ_val_orb_for_hf(2)! electron 2 + ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2) + two_bod_dens += mos_array_valence_hf_r1(m) * mos_array_valence_hf_r1(m) * mos_array_valence_hf_r2(n) * mos_array_valence_hf_r2(n) + ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space + do i = 1, n_basis_orb + do j = 1, n_basis_orb + ! 2 1 2 1 + f_HF_val_ab += two_e_int_hf_f(j,i,n,m) & + * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) & + * mos_array_valence_r2(j) * mos_array_valence_hf_r2(n) + enddo + enddo + enddo + enddo +end + + +subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab) + implicit none + BEGIN_DOC +! in_f_HF_val_ab(r_1) = \int dr_2 f_{\Psi^B}(r_1,r_2) +! +! where f_{\Psi^B}(r_1,r_2) is defined by Eq. (22) of J. Chem. Phys. 149, 194301 (2018) +! +! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931 2937) +! +! Such function can be used to test if the f_HF_val_ab(r_1,r_2) is correctly built. +! +! < HF | wee_{\alpha\beta} | HF > = \int (r1) int_f_HF_val_ab(r_1) + END_DOC + double precision, intent(in) :: r1(3) + double precision, intent(out):: int_f_HF_val_ab + integer :: i,j,m,n,i_m,i_n + integer :: i_i,i_j + double precision :: mo_two_e_integral + double precision :: mos_array_r1(mo_num) + double precision, allocatable :: mos_array_valence_r1(:) + double precision, allocatable :: mos_array_valence_hf_r1(:) + double precision :: get_two_e_integral + call give_all_mos_at_r(r1,mos_array_r1) + allocate(mos_array_valence_r1( n_basis_orb )) + allocate(mos_array_valence_hf_r1( n_occ_val_orb_for_hf(1) ) ) + do i_m = 1, n_occ_val_orb_for_hf(1) + mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1)) + enddo + do i_m = 1, n_basis_orb + mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m)) + enddo + + int_f_HF_val_ab = 0.d0 + ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space + do m = 1, n_occ_val_orb_for_hf(1)! electron 1 + ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space + do n = 1, n_occ_val_orb_for_hf(2)! electron 2 + ! You browse all ORBITALS in the \mathacal{B} space + do i = 1, n_basis_orb + ! due to integration in real-space and the use of orthonormal MOs, a Kronecker delta_jn shoes up + j = n + ! 2 1 2 1 + int_f_HF_val_ab += two_e_int_hf_f(j,i,n,m) & + * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) + enddo + enddo + enddo +end diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f new file mode 100644 index 00000000..21f330d6 --- /dev/null +++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f @@ -0,0 +1,313 @@ +subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens) + implicit none + BEGIN_DOC +! contribution from purely inactive orbitals to f_{\Psi^B}(r_1,r_2) for a CAS wave function + END_DOC + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: f_ii_val_ab,two_bod_dens + integer :: i,j,m,n,i_m,i_n + integer :: i_i,i_j + double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) + double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision :: get_two_e_integral + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + ! You extract the inactive orbitals + allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) + do i_m = 1, n_inact_orb + mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m)) + enddo + do i_m = 1, n_inact_orb + mos_array_inact_r2(i_m) = mos_array_r2(list_inact(i_m)) + enddo + + ! You extract the orbitals belonging to the space \mathcal{B} + allocate(mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) + do i_m = 1, n_basis_orb + mos_array_basis_r1(i_m) = mos_array_r1(list_basis(i_m)) + mos_array_basis_r2(i_m) = mos_array_r2(list_basis(i_m)) + enddo + + f_ii_val_ab = 0.d0 + two_bod_dens = 0.d0 + ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space + do m = 1, n_inact_orb ! electron 1 + ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space + do n = 1, n_inact_orb ! electron 2 + ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2) + two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) + ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space + do i = 1, n_basis_orb + do j = 1, n_basis_orb + ! 2 1 2 1 + f_ii_val_ab += two_e_int_ii_f(j,i,n,m) * mos_array_inact_r1(m) * mos_array_basis_r1(i) & + * mos_array_inact_r2(n) * mos_array_basis_r2(j) + enddo + enddo + enddo + enddo +end + + +subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) + BEGIN_DOC +! contribution from inactive and active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function + END_DOC + implicit none + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: f_ia_val_ab,two_bod_dens + integer :: i,orb_i,a,orb_a,n,m,b + double precision :: rho + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) + double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) + double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) + double precision, allocatable :: integrals_array(:,:),rho_tilde(:,:),v_tilde(:,:) + + f_ia_val_ab = 0.d0 + two_bod_dens = 0.d0 + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + + ! You extract the inactive orbitals + allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) + do i = 1, n_inact_orb + mos_array_inact_r1(i) = mos_array_r1(list_inact(i)) + enddo + do i= 1, n_inact_orb + mos_array_inact_r2(i) = mos_array_r2(list_inact(i)) + enddo + + ! You extract the active orbitals + allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) + do i= 1, n_act_orb + mos_array_act_r1(i) = mos_array_r1(list_act(i)) + enddo + do i= 1, n_act_orb + mos_array_act_r2(i) = mos_array_r2(list_act(i)) + enddo + + ! You extract the orbitals belonging to the space \mathcal{B} + allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) + do i= 1, n_basis_orb + mos_array_basis_r1(i) = mos_array_r1(list_basis(i)) + enddo + do i= 1, n_basis_orb + mos_array_basis_r2(i) = mos_array_r2(list_basis(i)) + enddo + + ! Contracted density : intermediate quantity + ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2) + allocate(rho_tilde(n_inact_orb,n_act_orb)) + two_bod_dens = 0.d0 + do a = 1, n_act_orb + do i = 1, n_inact_orb + rho_tilde(i,a) = 0.d0 + do b = 1, n_act_orb + rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) + two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho + rho_tilde(i,a) += rho * mos_array_inact_r1(i) * mos_array_act_r2(b) + enddo + enddo + enddo + + ! Contracted two-e integrals : intermediate quantity + ! v_tilde(i,a) = \sum_{m,n} phi_m(1) * phi_n(2) < i a | m n > + allocate( v_tilde(n_act_orb,n_act_orb) ) + allocate( integrals_array(mo_num,mo_num) ) + v_tilde = 0.d0 + do a = 1, n_act_orb + orb_a = list_act(a) + do i = 1, n_inact_orb + v_tilde(i,a) = 0.d0 + orb_i = list_inact(i) +! call get_mo_two_e_integrals_ij(orb_i,orb_a,mo_num,integrals_array,mo_integrals_map) + do m = 1, n_basis_orb + do n = 1, n_basis_orb +! v_tilde(i,a) += integrals_array(n,m) * mos_array_basis_r2(n) * mos_array_basis_r1(m) + v_tilde(i,a) += two_e_int_ia_f(n,m,i,a) * mos_array_basis_r2(n) * mos_array_basis_r1(m) + enddo + enddo + enddo + enddo + + do a = 1, n_act_orb + do i = 1, n_inact_orb + f_ia_val_ab += v_tilde(i,a) * rho_tilde(i,a) + enddo + enddo +end + + +subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) + BEGIN_DOC +! contribution from purely active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function + END_DOC + implicit none + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: f_aa_val_ab,two_bod_dens + integer :: i,orb_i,a,orb_a,n,m,b,c,d + double precision :: rho + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) + double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) + double precision, allocatable :: integrals_array(:,:),rho_tilde(:,:),v_tilde(:,:) + + f_aa_val_ab = 0.d0 + two_bod_dens = 0.d0 + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + + ! You extract the active orbitals + allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) + do i= 1, n_act_orb + mos_array_act_r1(i) = mos_array_r1(list_act(i)) + enddo + do i= 1, n_act_orb + mos_array_act_r2(i) = mos_array_r2(list_act(i)) + enddo + + ! You extract the orbitals belonging to the space \mathcal{B} + allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) + do i= 1, n_basis_orb + mos_array_basis_r1(i) = mos_array_r1(list_basis(i)) + enddo + do i= 1, n_basis_orb + mos_array_basis_r2(i) = mos_array_r2(list_basis(i)) + enddo + + ! Contracted density : intermediate quantity + ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2) + allocate(rho_tilde(n_act_orb,n_act_orb)) + two_bod_dens = 0.d0 + rho_tilde = 0.d0 + do a = 1, n_act_orb ! 1 + do b = 1, n_act_orb ! 2 + do c = 1, n_act_orb ! 1 + do d = 1, n_act_orb ! 2 + rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate) + rho_tilde(b,a) += rho + two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) + enddo + enddo + enddo + enddo + + ! Contracted two-e integrals : intermediate quantity + ! v_tilde(i,a) = \sum_{m,n} phi_m(1) * phi_n(2) < i a | m n > + allocate( v_tilde(n_act_orb,n_act_orb) ) + v_tilde = 0.d0 + do a = 1, n_act_orb + do b = 1, n_act_orb + v_tilde(b,a) = 0.d0 + do m = 1, n_basis_orb + do n = 1, n_basis_orb + v_tilde(b,a) += two_e_int_aa_f(n,m,b,a) * mos_array_basis_r2(n) * mos_array_basis_r1(m) + enddo + enddo + enddo + enddo + + do a = 1, n_act_orb + do b = 1, n_act_orb + f_aa_val_ab += v_tilde(b,a) * rho_tilde(b,a) + enddo + enddo +end + +BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{ii}(r_1,r_2) +! +! two_e_int_aa_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" + END_DOC + integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n + double precision :: integrals_array(mo_num,mo_num),get_two_e_integral + do orb_m = 1, n_act_orb ! electron 1 + m = list_act(orb_m) + do orb_n = 1, n_act_orb ! electron 2 + n = list_act(orb_n) + call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) + do orb_i = 1, n_basis_orb ! electron 1 + i = list_basis(orb_i) + do orb_j = 1, n_basis_orb ! electron 2 + j = list_basis(orb_j) + ! 2 1 2 1 + two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) +! two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) + enddo + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, two_e_int_ia_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_act_orb)] + implicit none + BEGIN_DOC +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{ia}(r_1,r_2) +! +! two_e_int_aa_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" + END_DOC + integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n + double precision :: integrals_array(mo_num,mo_num),get_two_e_integral + do orb_m = 1, n_act_orb ! electron 1 + m = list_act(orb_m) + do orb_n = 1, n_inact_orb ! electron 2 + n = list_inact(orb_n) + call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) + do orb_i = 1, n_basis_orb ! electron 1 + i = list_basis(orb_i) + do orb_j = 1, n_basis_orb ! electron 2 + j = list_basis(orb_j) + ! 2 1 2 1 +! two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) + two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) + enddo + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_inact_orb)] + implicit none + BEGIN_DOC +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{ii}(r_1,r_2) +! +! two_e_int_ii_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" + END_DOC + integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n + double precision :: get_two_e_integral,integrals_array(mo_num,mo_num) + do orb_m = 1, n_inact_orb ! electron 1 + m = list_inact(orb_m) + do orb_n = 1, n_inact_orb ! electron 2 + n = list_inact(orb_n) + call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) + do orb_i = 1, n_basis_orb ! electron 1 + i = list_basis(orb_i) + do orb_j = 1, n_basis_orb ! electron 2 + j = list_basis(orb_j) + ! 2 1 2 1 +! two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) + two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) + enddo + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/mu_of_r/f_psi_old.irp.f b/src/mu_of_r/f_psi_old.irp.f new file mode 100644 index 00000000..643b0608 --- /dev/null +++ b/src/mu_of_r/f_psi_old.irp.f @@ -0,0 +1,39 @@ + + +BEGIN_PROVIDER [double precision, f_psi_cas_ab_old, (n_points_final_grid,N_states)] + implicit none + BEGIN_DOC +! +! Function f_{\Psi^B}(r,r) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) on each point of the grid and for all states +! +! Assumes that the wave function in psi_det is developped within an active space defined +! + END_DOC + integer :: ipoint,k,l,istate + double precision :: wall0,wall1 + print*,'Providing f_psi_cas_ab_old ..... ' + provide full_occ_2_rdm_ab_mo + call wall_time(wall0) + provide core_inact_act_V_kl_contracted full_occ_2_rdm_cntrctd + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,k,l,istate) & + !$OMP SHARED (n_core_inact_act_orb, n_points_final_grid, full_occ_2_rdm_cntrctd, core_inact_act_V_kl_contracted, f_psi_cas_ab_old,N_states) + !$OMP DO + do istate = 1, N_states + do ipoint = 1, n_points_final_grid + f_psi_cas_ab_old(ipoint,istate) = 0.d0 + do l = 1, n_core_inact_act_orb ! 2 + do k = 1, n_core_inact_act_orb ! 1 + f_psi_cas_ab_old(ipoint,istate) += core_inact_act_V_kl_contracted(k,l,ipoint) * full_occ_2_rdm_cntrctd(k,l,ipoint,istate) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + call wall_time(wall1) + print*,'Time to provide f_psi_cas_ab_old = ',wall1 - wall0 +END_PROVIDER + + diff --git a/src/mu_of_r/f_psi_utils.irp.f b/src/mu_of_r/f_psi_utils.irp.f new file mode 100644 index 00000000..bdd76f18 --- /dev/null +++ b/src/mu_of_r/f_psi_utils.irp.f @@ -0,0 +1,151 @@ +BEGIN_PROVIDER [double precision, core_inact_act_V_kl_contracted, (n_core_inact_act_orb,n_core_inact_act_orb,n_points_final_grid)] + implicit none + BEGIN_DOC +! core_inact_act_V_kl_contracted(k,l,ipoint) = \sum_{ij} V_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint) +! +! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) +! + END_DOC + integer :: ipoint,k,l + do k = 1, n_core_inact_act_orb + do l = 1, n_core_inact_act_orb + do ipoint = 1, n_points_final_grid + core_inact_act_V_kl_contracted(k,l,ipoint) = full_occ_v_kl_cntrctd(ipoint,k,l) + enddo + enddo + enddo + free full_occ_v_kl_cntrctd + +END_PROVIDER + +BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd, (n_core_inact_act_orb,n_core_inact_act_orb,n_points_final_grid,N_states)] + implicit none + BEGIN_DOC +! full_occ_2_rdm_cntrctd(k,l,ipoint,istate) = \sum_{ij} \Gamma_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint) +! +! where \Gamma_{ij}^{kl}(istate) = +! +! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) +! + END_DOC + integer :: ipoint,k,l,istate + do istate = 1, N_states + do k = 1, n_core_inact_act_orb + do l = 1, n_core_inact_act_orb + do ipoint = 1, n_points_final_grid + full_occ_2_rdm_cntrctd(k,l,ipoint,istate) = full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) + enddo + enddo + enddo + enddo + free full_occ_2_rdm_cntrctd_trans +END_PROVIDER + + + + +BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd_trans, (n_points_final_grid,n_core_inact_act_orb,n_core_inact_act_orb,N_states)] + implicit none + BEGIN_DOC +! full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) = \sum_{ij} \Gamma_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint) +! +! where \Gamma_{ij}^{kl}(istate) = +! +! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) +! + END_DOC + integer :: i,j,k,l,istate + integer :: ipoint + double precision, allocatable :: mos_array_r(:),r(:) + provide full_occ_2_rdm_ab_mo + double precision :: wall0,wall1 + print*,'Providing full_occ_2_rdm_cntrctd_trans ..... ' + call wall_time(wall0) + full_occ_2_rdm_cntrctd_trans = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,k,l,i,j,istate) & + !$OMP SHARED (n_core_inact_act_orb, n_points_final_grid, full_occ_2_rdm_cntrctd_trans, final_grid_points,full_occ_2_rdm_ab_mo,core_inact_act_mos_in_r_array,N_states ) + !$OMP DO + do istate = 1, N_states + do l = 1, n_core_inact_act_orb ! 2 + do k = 1, n_core_inact_act_orb ! 1 + do ipoint = 1, n_points_final_grid + do j = 1, n_core_inact_act_orb + do i = 1, n_core_inact_act_orb + ! 1 2 1 2 + full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint) + enddo + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + call wall_time(wall1) + print*,'Time to provide full_occ_2_rdm_cntrctd_trans = ',wall1 - wall0 + +END_PROVIDER + + + +BEGIN_PROVIDER [double precision, full_occ_v_kl_cntrctd, (n_points_final_grid,n_core_inact_act_orb,n_core_inact_act_orb)] + implicit none + BEGIN_DOC +! full_occ_v_kl_cntrctd(ipoint,k,l) = \sum_{ij} V_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint) +! +! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) +! + END_DOC + integer :: i,j,k,l,kk,ll,ii,jj + integer :: ipoint + double precision, allocatable :: integrals_array(:,:), mos_array_r(:),r(:), integrals_basis(:,:) + ! just not to mess with parallelization + allocate(integrals_array(mo_num,mo_num)) + k = 1 + l = 1 + call get_mo_two_e_integrals_ij(k,l,mo_num,integrals_array,mo_integrals_map) + provide basis_mos_in_r_array + deallocate(integrals_array) + double precision :: wall0,wall1 + call wall_time(wall0) + + full_occ_v_kl_cntrctd = 0.d0 + print*,'Providing full_occ_v_kl_cntrctd ..... ' + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,kk,ll,k,l,i,j,ii,jj,integrals_array,integrals_basis) & + !$OMP SHARED (mo_num, n_points_final_grid, n_basis_orb, list_basis, full_occ_v_kl_cntrctd, mo_integrals_map,final_grid_points,basis_mos_in_r_array, n_core_inact_act_orb, list_core_inact_act) + allocate(integrals_array(mo_num,mo_num), integrals_basis(n_basis_orb,n_basis_orb)) + !$OMP DO + do l = 1, n_core_inact_act_orb! 2 + ll = list_core_inact_act(l) + do k = 1, n_core_inact_act_orb ! 1 + kk = list_core_inact_act(k) + call get_mo_two_e_integrals_ij(kk,ll,mo_num,integrals_array,mo_integrals_map) + do j = 1, n_basis_orb + jj = list_basis(j) + do i = 1, n_basis_orb + ii = list_basis(i) + integrals_basis(i,j) = integrals_array(ii,jj) + enddo + enddo + do ipoint = 1, n_points_final_grid + do j = 1, n_basis_orb ! condition on mo_num in order to ensure the correct CBS limit + do i = 1, n_basis_orb ! + !1 2 1 2 + full_occ_v_kl_cntrctd(ipoint,k,l) += integrals_basis(i,j) * basis_mos_in_r_array(j,ipoint) * basis_mos_in_r_array(i,ipoint) + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + deallocate(integrals_array,integrals_basis) + !$OMP END PARALLEL + + call wall_time(wall1) + print*,'Time to provide full_occ_v_kl_cntrctd = ',wall1 - wall0 +END_PROVIDER + diff --git a/src/mu_of_r/f_val_general.irp.f b/src/mu_of_r/f_val_general.irp.f new file mode 100644 index 00000000..7939dbf0 --- /dev/null +++ b/src/mu_of_r/f_val_general.irp.f @@ -0,0 +1,90 @@ + BEGIN_PROVIDER [double precision, f_psi_cas_ab, (n_points_final_grid,N_states)] +&BEGIN_PROVIDER [double precision, on_top_cas_mu_r, (n_points_final_grid,N_states)] + implicit none + BEGIN_DOC +! +! Function f_{\Psi^B}(r,r) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) on each point of the grid and for all states and for a CAS wave function +! +! Assumes that the wave function in psi_det is developped within an active space defined +! + END_DOC + integer :: ipoint,istate + double precision :: wall0,wall1,r(3) + double precision :: f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa + double precision :: accu + accu = 0.d0 + r = 0.d0 + istate = 1 + ! To initialize parallelization + call give_f_ii_val_ab(r,r,f_ii_val_ab,two_bod_dens_ii) + call give_f_ia_val_ab(r,r,f_ia_val_ab,two_bod_dens_ia,istate) + call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate) + provide final_grid_points act_2_rdm_ab_mo + + print*,'Providing f_psi_cas_ab..... ' + + + call wall_time(wall0) + do istate = 1, N_states + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,r,f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa) & + !$OMP SHARED (n_points_final_grid,f_psi_cas_ab,on_top_cas_mu_r,final_grid_points,istate) + !$OMP DO + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + ! inactive-inactive part of f_psi(r1,r2) + call give_f_ii_val_ab(r,r,f_ii_val_ab,two_bod_dens_ii) + ! inactive-active part of f_psi(r1,r2) + call give_f_ia_val_ab(r,r,f_ia_val_ab,two_bod_dens_ia,istate) + ! active-active part of f_psi(r1,r2) + call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate) + f_psi_cas_ab(ipoint,istate) = f_ii_val_ab + f_ia_val_ab + f_aa_val_ab + on_top_cas_mu_r(ipoint,istate) = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa + enddo + !$OMP END DO + !$OMP END PARALLEL + enddo + call wall_time(wall1) + print*,'Time to provide f_psi_cas_ab = ',wall1 - wall0 + print*,'accu = ',accu + +END_PROVIDER + + BEGIN_PROVIDER [double precision, f_psi_hf_ab, (n_points_final_grid)] +&BEGIN_PROVIDER [double precision, on_top_hf_mu_r, (n_points_final_grid)] + implicit none + BEGIN_DOC +! +! Function f_{\Psi^B}(r,r) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) on each point of the grid for a HF wave function +! + END_DOC + integer :: ipoint + double precision :: wall0,wall1,r(3),f_HF_val_ab,two_bod_dens + f_psi_hf_ab = 0.d0 + r = 0.d0 + ! To initialize parallelization + call f_HF_valence_ab(r,r,f_HF_val_ab,two_bod_dens) + + call wall_time(wall0) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,r,f_HF_val_ab,two_bod_dens) & + !$OMP SHARED (n_points_final_grid,f_psi_hf_ab,on_top_hf_mu_r,final_grid_points) + !$OMP DO + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + call f_HF_valence_ab(r,r,f_HF_val_ab,two_bod_dens) + f_psi_hf_ab(ipoint) = f_HF_val_ab + on_top_hf_mu_r(ipoint) = two_bod_dens + enddo + !$OMP END DO + !$OMP END PARALLEL + call wall_time(wall1) + print*,'Time to provide f_psi_hf_ab = ',wall1 - wall0 + +END_PROVIDER diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f new file mode 100644 index 00000000..05a2a4d7 --- /dev/null +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -0,0 +1,150 @@ + + BEGIN_PROVIDER [double precision, mu_of_r_prov, (n_points_final_grid,N_states) ] + implicit none + BEGIN_DOC + ! general variable for mu(r) + ! + ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the two-body density matrix are excluded + END_DOC + integer :: ipoint,istate + double precision :: wall0,wall1 + print*,'providing mu_of_r ...' +! PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + call wall_time(wall0) + + if (read_mu_of_r) then + print*,'Reading mu(r) from disk ...' + call ezfio_get_mu_of_r_mu_of_r_disk(mu_of_r_prov) + return + endif + + do istate = 1, N_states + do ipoint = 1, n_points_final_grid + if(mu_of_r_potential.EQ."hf")then + mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) + else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then + mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) + else + print*,'you requested the following mu_of_r_potential' + print*,mu_of_r_potential + print*,'which does not correspond to any of the options for such keyword' + stop + endif + enddo + enddo + + if (write_mu_of_r) then + print*,'Writing mu(r) on disk ...' + call ezfio_set_mu_of_r_io_mu_of_r('Read') + call ezfio_set_mu_of_r_mu_of_r_disk(mu_of_r_prov) + endif + + call wall_time(wall1) + print*,'Time to provide mu_of_r = ',wall1-wall0 + END_PROVIDER + + + BEGIN_PROVIDER [double precision, mu_of_r_hf, (n_points_final_grid) ] + implicit none + BEGIN_DOC + ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO) + ! + ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the two-body density matrix are excluded + END_DOC + integer :: ipoint + double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + print*,'providing mu_of_r_hf ...' + call wall_time(wall0) + sqpi = dsqrt(dacos(-1.d0)) + provide f_psi_hf_ab + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & + !$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_psi_hf_ab,on_top_hf_mu_r,sqpi) + do ipoint = 1, n_points_final_grid + f_hf = f_psi_hf_ab(ipoint) + on_top = on_top_hf_mu_r(ipoint) + if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then + w_hf = 1.d+10 + else + w_hf = f_hf / on_top + endif + mu_of_r_hf(ipoint) = w_hf * sqpi * 0.5d0 + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide mu_of_r_hf = ',wall1-wall0 + END_PROVIDER + + BEGIN_PROVIDER [double precision, mu_of_r_psi_cas, (n_points_final_grid,N_states) ] + implicit none + BEGIN_DOC + ! mu(r) computed with a wave function developped in an active space + ! + ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the one- and two-body density matrix are excluded + END_DOC + integer :: ipoint,istate + double precision :: wall0,wall1,f_psi,on_top,w_psi,sqpi + print*,'providing mu_of_r_psi_cas ...' + call wall_time(wall0) + sqpi = dsqrt(dacos(-1.d0)) + + provide f_psi_cas_ab + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & + !$OMP SHARED (n_points_final_grid,mu_of_r_psi_cas,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) + do istate = 1, N_states + do ipoint = 1, n_points_final_grid + f_psi = f_psi_cas_ab(ipoint,istate) + on_top = on_top_cas_mu_r(ipoint,istate) + if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then + w_psi = 1.d+10 + else + w_psi = f_psi / on_top + endif + mu_of_r_psi_cas(ipoint,istate) = w_psi * sqpi * 0.5d0 + enddo + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0 + END_PROVIDER + + + BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)] + implicit none + BEGIN_DOC + ! average value of mu(r) weighted with the total one-e density and divised by the number of electrons + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the one- and two-body density matrix are excluded + END_DOC + integer :: ipoint,istate + double precision :: weight,density + mu_average_prov = 0.d0 + do istate = 1, N_states + do ipoint = 1, n_points_final_grid + weight =final_weight_at_r_vector(ipoint) + density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & + + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) + if(mu_of_r_prov(ipoint,istate).gt.1.d+09)cycle + mu_average_prov(istate) += mu_of_r_prov(ipoint,istate) * weight * density + enddo + mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate) + enddo + END_PROVIDER diff --git a/src/mu_of_r/test_proj_op.irp.f b/src/mu_of_r/test_proj_op.irp.f new file mode 100644 index 00000000..1d46da5e --- /dev/null +++ b/src/mu_of_r/test_proj_op.irp.f @@ -0,0 +1,21 @@ +program projected_operators + implicit none + BEGIN_DOC +! TODO + END_DOC + read_wf = .True. + touch read_wf + ! You specify that you want to avoid any contribution from + ! orbitals coming from core + no_core_density = .True. + touch no_core_density + mu_of_r_potential = "cas_ful" + touch mu_of_r_potential + print*,'Using Valence Only functions' +! call test_f_HF_valence_ab +! call routine_full_mos +! call test_f_ii_valence_ab + call test_f_ia_valence_ab + call test_f_ii_ia_aa_valence_ab +end +