LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-10-19 22:41:48 +02:00
Code Issues Releases Wiki Activity

added dft_utils_two_body which works

Browse Source
This commit is contained in:
Emmanuel Giner 2018-12-20 16:07:06 +01:00
parent c1a6002d8c
commit acf5bd4d52
28 changed files with 1576 additions and 176 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
TODO
View File

@ -5,6 +5,7 @@
* dft_utils_one_body
* dm_for_dft
* integrals_erf
* Faire un README pour les scripts workflow
* Molden format : http://cheminf.cmbi.ru.nl/molden/molden_format.html
* Virer tous les modules qui sont dans plugins
* Permettre aux utilisateurs de facilement deposer des plugins dans plugins via une commande

2
quantum_package.rc.default
View File

@ -20,7 +20,7 @@ else
fi
}
export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")
export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)
export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml:"${QP_ROOT}"/scripts/workflow)
export LD_LIBRARY_PATH=$(qp_prepend_export "LD_LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)
export LIBRARY_PATH=$(qp_prepend_export "LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)
export C_INCLUDE_PATH=$(qp_prepend_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)

60
src/ao_basis/aos_transp.irp.f Normal file
View File

@ -0,0 +1,60 @@
BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)]
implicit none
BEGIN_DOC
! List of AOs attached on each atom
END_DOC
integer :: i
integer, allocatable :: nucl_tmp(:)
allocate(nucl_tmp(nucl_num))
nucl_tmp = 0
Nucl_Aos = 0
do i = 1, ao_num
nucl_tmp(ao_nucl(i))+=1
Nucl_Aos_transposed(nucl_tmp(ao_nucl(i)),ao_nucl(i)) = i
enddo
deallocate(nucl_tmp)
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_expo_ordered_transp_per_nucl, (ao_prim_num_max,N_AOs_max,nucl_num) ]
implicit none
integer :: i,j,k,l
do i = 1, nucl_num
do j = 1,Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i)
do l = 1, ao_prim_num(k)
ao_expo_ordered_transp_per_nucl(l,j,i) = ao_expo_ordered_transp(l,k)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_power_ordered_transp_per_nucl, (3,N_AOs_max,nucl_num) ]
implicit none
integer :: i,j,k,l
do i = 1, nucl_num
do j = 1,Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i)
do l = 1, 3
ao_power_ordered_transp_per_nucl(l,j,i) = ao_power(k,l)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered_transp_per_nucl, (ao_prim_num_max,N_AOs_max,nucl_num) ]
implicit none
integer :: i,j,k,l
do i = 1, nucl_num
do j = 1,Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i)
do l = 1, ao_prim_num(k)
ao_coef_normalized_ordered_transp_per_nucl(l,j,i) = ao_coef_normalized_ordered_transp(l,k)
enddo
enddo
enddo
END_PROVIDER

261
src/ao_basis/aos_value.irp.f
View File

@ -1,7 +1,7 @@
double precision function ao_value(i,r)
implicit none
BEGIN_DOC
! Returns the value of the i-th |AO| at point r
! return the value of the ith ao at point r
END_DOC
double precision, intent(in) :: r(3)
integer, intent(in) :: i
@ -10,10 +10,10 @@ double precision function ao_value(i,r)
double precision :: center_ao(3)
double precision :: beta
integer :: power_ao(3)
double precision :: accu,dx,dy,dz,r2
num_ao = ao_nucl(i)
power_ao(1:3)= ao_power(i,1:3)
center_ao(1:3) = nucl_coord(num_ao,1:3)
double precision :: accu,dx,dy,dz,r2
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
@ -31,10 +31,44 @@ double precision function ao_value(i,r)
end
subroutine give_all_aos_at_r(r,aos_array)
double precision function primitive_value(i,j,r)
implicit none
BEGIN_DOC
! return the value of the jth primitive of ith ao at point r WITHOUT THE COEF
END_DOC
double precision, intent(in) :: r(3)
integer, intent(in) :: i,j
integer :: m,num_ao
double precision :: center_ao(3)
double precision :: beta
integer :: power_ao(3)
double precision :: accu,dx,dy,dz,r2
num_ao = ao_nucl(i)
power_ao(1:3)= ao_power(i,1:3)
center_ao(1:3) = nucl_coord(num_ao,1:3)
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
r2 = dx*dx + dy*dy + dz*dz
dx = dx**power_ao(1)
dy = dy**power_ao(2)
dz = dz**power_ao(3)
accu = 0.d0
m=j
beta = ao_expo_ordered_transp(m,i)
accu += dexp(-beta*r2)
primitive_value = accu * dx * dy * dz
end
subroutine give_all_aos_at_r_old(r,aos_array)
implicit none
BEGIN_dOC
! Gives the values of |AOs| at a given point r
! gives the values of aos at a given point r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num)
@ -43,5 +77,222 @@ subroutine give_all_aos_at_r(r,aos_array)
do i = 1, ao_num
aos_array(i) = ao_value(i,r)
enddo
end
subroutine give_all_aos_at_r(r,aos_array)
implicit none
BEGIN_dOC
! input : r == r(1) = x and so on
! aos_array(i) = aos(i) evaluated in r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num)
integer :: power_ao(3)
integer :: i,j,k,l,m
double precision :: dx,dy,dz,r2
double precision :: dx2,dy2,dz2
double precision :: center_ao(3)
double precision :: beta
do i = 1, nucl_num
center_ao(1:3) = nucl_coord(i,1:3)
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
r2 = dx*dx + dy*dy + dz*dz
do j = 1,Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
aos_array(k) = 0.d0
power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
dx2 = dx**power_ao(1)
dy2 = dy**power_ao(2)
dz2 = dz**power_ao(3)
do l = 1,ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i)
aos_array(k)+= ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
enddo
aos_array(k) = aos_array(k) * dx2 * dy2 * dz2
enddo
enddo
end
subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
implicit none
BEGIN_DOC
! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
! output : aos_array(i) = ao(i) evaluated at r
! : aos_grad_array(1,i) = gradient X of the ao(i) evaluated at r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num)
double precision, intent(out) :: aos_grad_array(3,ao_num)
integer :: power_ao(3)
integer :: i,j,k,l,m
double precision :: dx,dy,dz,r2
double precision :: dx2,dy2,dz2
double precision :: dx1,dy1,dz1
double precision :: center_ao(3)
double precision :: beta,accu_1,accu_2,contrib
do i = 1, nucl_num
center_ao(1:3) = nucl_coord(i,1:3)
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
r2 = dx*dx + dy*dy + dz*dz
do j = 1,Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
aos_array(k) = 0.d0
aos_grad_array(1,k) = 0.d0
aos_grad_array(2,k) = 0.d0
aos_grad_array(3,k) = 0.d0
power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
dx2 = dx**power_ao(1)
dy2 = dy**power_ao(2)
dz2 = dz**power_ao(3)
if(power_ao(1) .ne. 0)then
dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
else
dx1 = 0.d0
endif
if(power_ao(2) .ne. 0)then
dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
else
dy1 = 0.d0
endif
if(power_ao(3) .ne. 0)then
dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
else
dz1 = 0.d0
endif
accu_1 = 0.d0
accu_2 = 0.d0
do l = 1,ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i)
contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
accu_1 += contrib
accu_2 += contrib * beta
enddo
aos_array(k) = accu_1 * dx2 * dy2 * dz2
aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2- 2.d0 * dx2 * dx * dy2 * dz2 * accu_2
aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2- 2.d0 * dx2 * dy2 * dy * dz2 * accu_2
aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1- 2.d0 * dx2 * dy2 * dz2 * dz * accu_2
enddo
enddo
end
subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
implicit none
BEGIN_DOC
! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
! output : aos_array(i) = ao(i) evaluated at r
! : aos_grad_array(1,i) = gradient X of the ao(i) evaluated at r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num)
double precision, intent(out) :: aos_grad_array(ao_num,3)
double precision, intent(out) :: aos_lapl_array(ao_num,3)
integer :: power_ao(3)
integer :: i,j,k,l,m
double precision :: dx,dy,dz,r2
double precision :: dx2,dy2,dz2
double precision :: dx1,dy1,dz1
double precision :: dx3,dy3,dz3
double precision :: dx4,dy4,dz4
double precision :: dx5,dy5,dz5
double precision :: center_ao(3)
double precision :: beta,accu_1,accu_2,accu_3,contrib
do i = 1, nucl_num
center_ao(1:3) = nucl_coord(i,1:3)
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
r2 = dx*dx + dy*dy + dz*dz
do j = 1,Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
aos_array(k) = 0.d0
aos_grad_array(k,1) = 0.d0
aos_grad_array(k,2) = 0.d0
aos_grad_array(k,3) = 0.d0
aos_lapl_array(k,1) = 0.d0
aos_lapl_array(k,2) = 0.d0
aos_lapl_array(k,3) = 0.d0
power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i)
dx2 = dx**power_ao(1)
dy2 = dy**power_ao(2)
dz2 = dz**power_ao(3)
if(power_ao(1) .ne. 0)then
dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
else
dx1 = 0.d0
endif
! For the Laplacian
if(power_ao(1) .ge. 2)then
dx3 = dble(power_ao(1)) * dble((power_ao(1)-1)) * dx**(power_ao(1)-2)
else
dx3 = 0.d0
endif
dx4 = dble((2 * power_ao(1) + 1)) * dx**(power_ao(1))
dx5 = dx**(power_ao(1)+2)
if(power_ao(2) .ne. 0)then
dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
else
dy1 = 0.d0
endif
! For the Laplacian
if(power_ao(2) .ge. 2)then
dy3 = dble(power_ao(2)) * dble((power_ao(2)-1)) * dy**(power_ao(2)-2)
else
dy3 = 0.d0
endif
dy4 = dble((2 * power_ao(2) + 1)) * dy**(power_ao(2))
dy5 = dy**(power_ao(2)+2)
if(power_ao(3) .ne. 0)then
dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
else
dz1 = 0.d0
endif
! For the Laplacian
if(power_ao(3) .ge. 2)then
dz3 = dble(power_ao(3)) * dble((power_ao(3)-1)) * dz**(power_ao(3)-2)
else
dz3 = 0.d0
endif
dz4 = dble((2 * power_ao(3) + 1)) * dz**(power_ao(3))
dz5 = dz**(power_ao(3)+2)
accu_1 = 0.d0
accu_2 = 0.d0
accu_3 = 0.d0
do l = 1,ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i)
contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
accu_1 += contrib
accu_2 += contrib * beta
accu_3 += contrib * beta**2
enddo
aos_array(k) = accu_1 * dx2 * dy2 * dz2
aos_grad_array(k,1) = accu_1 * dx1 * dy2 * dz2- 2.d0 * dx2 * dx * dy2 * dz2 * accu_2
aos_grad_array(k,2) = accu_1 * dx2 * dy1 * dz2- 2.d0 * dx2 * dy2 * dy * dz2 * accu_2
aos_grad_array(k,3) = accu_1 * dx2 * dy2 * dz1- 2.d0 * dx2 * dy2 * dz2 * dz * accu_2
aos_lapl_array(k,1) = accu_1 * dx3 * dy2 * dz2- 2.d0 * dx4 * dy2 * dz2* accu_2 +4.d0 * dx5 *dy2 * dz2* accu_3
aos_lapl_array(k,2) = accu_1 * dx2 * dy3 * dz2- 2.d0 * dx2 * dy4 * dz2* accu_2 +4.d0 * dx2 *dy5 * dz2* accu_3
aos_lapl_array(k,3) = accu_1 * dx2 * dy2 * dz3- 2.d0 * dx2 * dy2 * dz4* accu_2 +4.d0 * dx2 *dy2 * dz5* accu_3
enddo
enddo
end

136
src/determinants/mono_excitations_bielec.irp.f Normal file
View File

@ -0,0 +1,136 @@
use bitmasks
subroutine get_mono_excitation_from_fock_bielec(det_1,det_2,h,p,spin,phase,hij)
use bitmasks
implicit none
integer,intent(in) :: h,p,spin
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: det_1(N_int,2), det_2(N_int,2)
double precision, intent(out) :: hij
integer(bit_kind) :: differences(N_int,2)
integer(bit_kind) :: hole(N_int,2)
integer(bit_kind) :: partcl(N_int,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_partcl(N_int*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i
do i = 1, N_int
differences(i,1) = xor(det_1(i,1),ref_closed_shell_bitmask(i,1))
differences(i,2) = xor(det_1(i,2),ref_closed_shell_bitmask(i,2))
hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
partcl(i,1) = iand(differences(i,1),det_1(i,1))
partcl(i,2) = iand(differences(i,2),det_1(i,2))
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
hij = fock_operator_bielec_closed_shell_ref_bitmask(h,p)
! holes :: direct terms
do i0 = 1, n_occ_ab_hole(1)
i = occ_hole(i0,1)
hij -= big_array_coulomb_integrals(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
do i0 = 1, n_occ_ab_hole(2)
i = occ_hole(i0,2)
hij -= big_array_coulomb_integrals(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
! holes :: exchange terms
do i0 = 1, n_occ_ab_hole(spin)
i = occ_hole(i0,spin)
hij += big_array_exchange_integrals(i,h,p) ! get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map)
enddo
! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,1)
hij += big_array_coulomb_integrals(i,h,p)!get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
do i0 = 1, n_occ_ab_partcl(2)
i = occ_partcl(i0,2)
hij += big_array_coulomb_integrals(i,h,p) !get_mo_bielec_integral_schwartz(h,i,p,i,mo_integrals_map)
enddo
! particles :: exchange terms
do i0 = 1, n_occ_ab_partcl(spin)
i = occ_partcl(i0,spin)
hij -= big_array_exchange_integrals(i,h,p)!get_mo_bielec_integral_schwartz(h,i,i,p,mo_integrals_map)
enddo
hij = hij * phase
end
BEGIN_PROVIDER [double precision, fock_operator_bielec_closed_shell_ref_bitmask, (mo_tot_num, mo_tot_num) ]
implicit none
integer :: i0,j0,i,j,k0,k
integer :: n_occ_ab(2)
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab_virt(2)
integer :: occ_virt(N_int*bit_kind_size,2)
integer(bit_kind) :: key_test(N_int)
integer(bit_kind) :: key_virt(N_int,2)
call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int)
do i = 1, N_int
key_virt(i,1) = full_ijkl_bitmask(i)
key_virt(i,2) = full_ijkl_bitmask(i)
key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1))
key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2))
enddo
double precision :: array_coulomb(mo_tot_num),array_exchange(mo_tot_num)
call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int)
! docc ---> virt mono excitations
do i0 = 1, n_occ_ab(1)
i=occ(i0,1)
do j0 = 1, n_occ_ab_virt(1)
j = occ_virt(j0,1)
call get_mo_bielec_integrals_coulomb_ii(i,j,mo_tot_num,array_coulomb,mo_integrals_map)
call get_mo_bielec_integrals_exch_ii(i,j,mo_tot_num,array_exchange,mo_integrals_map)
double precision :: accu
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * array_coulomb(k) - array_exchange(k)
enddo
fock_operator_bielec_closed_shell_ref_bitmask(i,j) = accu
fock_operator_bielec_closed_shell_ref_bitmask(j,i) = accu
enddo
enddo
! virt ---> virt mono excitations
do i0 = 1, n_occ_ab_virt(1)
i=occ_virt(i0,1)
do j0 = 1, n_occ_ab_virt(1)
j = occ_virt(j0,1)
call get_mo_bielec_integrals_coulomb_ii(i,j,mo_tot_num,array_coulomb,mo_integrals_map)
call get_mo_bielec_integrals_exch_ii(i,j,mo_tot_num,array_exchange,mo_integrals_map)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * array_coulomb(k) - array_exchange(k)
enddo
fock_operator_bielec_closed_shell_ref_bitmask(i,j) = accu
fock_operator_bielec_closed_shell_ref_bitmask(j,i) = accu
enddo
enddo
! docc ---> docc mono excitations
do i0 = 1, n_occ_ab(1)
i=occ(i0,1)
do j0 = 1, n_occ_ab(1)
j = occ(j0,1)
call get_mo_bielec_integrals_coulomb_ii(i,j,mo_tot_num,array_coulomb,mo_integrals_map)
call get_mo_bielec_integrals_exch_ii(i,j,mo_tot_num,array_exchange,mo_integrals_map)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * array_coulomb(k) - array_exchange(k)
enddo
fock_operator_bielec_closed_shell_ref_bitmask(i,j) = accu
fock_operator_bielec_closed_shell_ref_bitmask(j,i) = accu
enddo
enddo
END_PROVIDER

1
src/determinants/psi_energy_mono_elec.irp.f
View File

@ -16,6 +16,7 @@
enddo
enddo
double precision :: accu
accu = 0.d0
do i = 1, N_states
do j = 1, mo_tot_num
accu += one_body_dm_mo_alpha(j,j,i) + one_body_dm_mo_beta(j,j,i)

36
src/dft_utils_one_body/e_c_e_x_and_potentials_mo.irp.f → src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f
View File

@ -5,7 +5,7 @@
&BEGIN_PROVIDER [double precision, potential_c_beta_ao,(ao_num,ao_num,N_states)]
implicit none
BEGIN_DOC
! alpha/beta exchange/correlation potentials on the AO basis
! general providers for the alpha/beta exchange/correlation potentials on the AO basis
END_DOC
if(trim(exchange_functional)=="short_range_LDA")then
@ -51,7 +51,7 @@ END_PROVIDER
&BEGIN_PROVIDER [double precision, potential_c_beta_mo,(mo_tot_num,mo_tot_num,N_states)]
implicit none
BEGIN_DOC
! alpha/beta exchange/correlation potentials on the MO basis
! general providers for the alpha/beta exchange/correlation potentials on the MO basis
END_DOC
integer :: istate
do istate = 1, N_states
@ -92,6 +92,9 @@ END_PROVIDER
BEGIN_PROVIDER [double precision, energy_x, (N_states)]
&BEGIN_PROVIDER [double precision, energy_c, (N_states)]
implicit none
BEGIN_DOC
! correlation and exchange energies general providers.
END_DOC
if(trim(exchange_functional)=="short_range_LDA")then
energy_x = energy_x_LDA
energy_x = energy_x_LDA
@ -123,3 +126,32 @@ END_PROVIDER
endif
END_PROVIDER
BEGIN_PROVIDER [double precision, Trace_v_xc, (N_states)]
&BEGIN_PROVIDER [double precision, Trace_v_H, (N_states)]
&BEGIN_PROVIDER [double precision, Trace_v_Hxc, (N_states)]
implicit none
integer :: i,j,istate
double precision :: dm
BEGIN_DOC
! Trace_v_xc = \sum_{i,j} (rho_{ij}_\alpha v^{xc}_{ij}^\alpha + rho_{ij}_\beta v^{xc}_{ij}^\beta)
! Trace_v_Hxc = \sum_{i,j} v^{H}_{ij} (rho_{ij}_\alpha + rho_{ij}_\beta)
! Trace_v_Hxc = \sum_{i,j} rho_{ij} v^{Hxc}_{ij}
END_DOC
do istate = 1, N_states
Trace_v_xc(istate) = 0.d0
Trace_v_H(istate) = 0.d0
do i = 1, mo_tot_num
do j = 1, mo_tot_num
Trace_v_xc(istate) += (potential_x_alpha_mo(j,i,istate) + potential_c_alpha_mo(j,i,istate)) * one_body_dm_mo_alpha_for_dft(j,i,istate)
Trace_v_xc(istate) += (potential_x_beta_mo(j,i,istate) + potential_c_beta_mo(j,i,istate) ) * one_body_dm_mo_beta_for_dft(j,i,istate)
dm = one_body_dm_mo_alpha_for_dft(j,i,istate) + one_body_dm_mo_beta_for_dft(j,i,istate)
Trace_v_H(istate) += dm * short_range_Hartree_operator(j,i,istate)
enddo
enddo
Trace_v_Hxc(istate) = Trace_v_xc(istate) + Trace_v_H(istate)
enddo
END_PROVIDER

86
src/dft_utils_one_body/e_c_e_x_specific_prov.irp.f Normal file
View File

@ -0,0 +1,86 @@
BEGIN_PROVIDER[double precision, energy_x_LDA, (N_states) ]
&BEGIN_PROVIDER[double precision, energy_c_LDA, (N_states) ]
implicit none
BEGIN_DOC
! exchange/correlation energy with the short range LDA functional
END_DOC
integer :: istate,i,j
double precision :: r(3)
double precision :: mu,weight
double precision :: e_c,vc_a,vc_b,e_x,vx_a,vx_b
double precision, allocatable :: rhoa(:),rhob(:)
allocate(rhoa(N_states), rhob(N_states))
energy_x_LDA = 0.d0
energy_c_LDA = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rhoa(istate) = one_body_dm_alpha_at_r(i,istate)
rhob(istate) = one_body_dm_beta_at_r(i,istate)
call ec_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
call ex_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
energy_x_LDA(istate) += weight * e_x
energy_c_LDA(istate) += weight * e_c
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, energy_x_PBE, (N_states) ]
&BEGIN_PROVIDER[double precision, energy_c_PBE, (N_states) ]
implicit none
BEGIN_DOC
! exchange/correlation energy with the short range PBE functional
END_DOC
integer :: istate,i,j,m
double precision :: r(3)
double precision :: mu,weight
double precision, allocatable :: ex(:), ec(:)
double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
energy_x_PBE = 0.d0
energy_c_PBE = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rho_a(istate) = one_body_dm_alpha_and_grad_at_r(4,i,istate)
rho_b(istate) = one_body_dm_beta_and_grad_at_r(4,i,istate)
grad_rho_a(1:3,istate) = one_body_dm_alpha_and_grad_at_r(1:3,i,istate)
grad_rho_b(1:3,istate) = one_body_dm_beta_and_grad_at_r(1:3,i,istate)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
enddo
! inputs
call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
energy_x_PBE += ex * weight
energy_c_PBE += ec * weight
enddo
enddo
END_PROVIDER

6
src/dft_utils_one_body/one_body_psi_energy.irp.f → src/dft_utils_one_body/one_body_psi_energy_dft.irp.f
View File

@ -17,13 +17,15 @@
enddo
enddo
enddo
accu = 0.d0
do i = 1, N_states
do j = 1, mo_tot_num
accu += one_body_dm_mo_alpha_for_dft(j,j,i) + one_body_dm_mo_beta_for_dft(j,j,i)
enddo
accu = (elec_alpha_num + elec_beta_num ) / accu
psi_energy_h_core(i) = psi_dft_energy_h_core(i) * accu
psi_dft_energy_kinetic(i) = psi_dft_energy_kinetic(i) * accu
psi_dft_energy_nuclear_elec(i) = psi_dft_energy_nuclear_elec(i) * accu
psi_dft_energy_h_core(i) = psi_dft_energy_nuclear_elec(i) + psi_dft_energy_kinetic(i)
enddo
END_PROVIDER

84
src/dft_utils_one_body/potentials_ao_on_grids.irp.f → src/dft_utils_one_body/potentials_ao_specific_prov.irp.f
View File

@ -33,37 +33,6 @@
END_PROVIDER
BEGIN_PROVIDER[double precision, energy_x_LDA, (N_states) ]
&BEGIN_PROVIDER[double precision, energy_c_LDA, (N_states) ]
implicit none
BEGIN_DOC
! exchange/correlation energy with the short range LDA functional
END_DOC
integer :: istate,i,j
double precision :: r(3)
double precision :: mu,weight
double precision :: e_c,vc_a,vc_b,e_x,vx_a,vx_b
double precision, allocatable :: rhoa(:),rhob(:)
allocate(rhoa(N_states), rhob(N_states))
energy_x_LDA = 0.d0
energy_c_LDA = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rhoa(istate) = one_body_dm_alpha_at_r(i,istate)
rhob(istate) = one_body_dm_beta_at_r(i,istate)
call ec_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_c,vc_a,vc_b)
call ex_LDA_sr(mu_erf,rhoa(istate),rhob(istate),e_x,vx_a,vx_b)
energy_x_LDA(istate) += weight * e_x
energy_c_LDA(istate) += weight * e_c
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, potential_x_alpha_ao_LDA,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_x_beta_ao_LDA,(ao_num,ao_num,N_states)]
@ -174,59 +143,6 @@
END_PROVIDER
BEGIN_PROVIDER[double precision, energy_x_PBE, (N_states) ]
&BEGIN_PROVIDER[double precision, energy_c_PBE, (N_states) ]
implicit none
BEGIN_DOC
! exchange/correlation energy with the short range PBE functional
END_DOC
integer :: istate,i,j,m
double precision :: r(3)
double precision :: mu,weight
double precision, allocatable :: ex(:), ec(:)
double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
energy_x_PBE = 0.d0
energy_c_PBE = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight=final_weight_functions_at_final_grid_points(i)
rho_a(istate) = one_body_dm_alpha_and_grad_at_r(4,i,istate)
rho_b(istate) = one_body_dm_beta_and_grad_at_r(4,i,istate)
grad_rho_a(1:3,istate) = one_body_dm_alpha_and_grad_at_r(1:3,i,istate)
grad_rho_b(1:3,istate) = one_body_dm_beta_and_grad_at_r(1:3,i,istate)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
enddo
! inputs
call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange
ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation
ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
energy_x_PBE += ex * weight
energy_c_PBE += ec * weight
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, potential_x_alpha_ao_PBE,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_x_beta_ao_PBE,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_PBE,(ao_num,ao_num,N_states)]

16
src/dft_utils_one_body/shifted_potential.irp.f
View File

@ -1,16 +0,0 @@
BEGIN_PROVIDER [double precision, shifting_constant, (N_states)]
implicit none
BEGIN_DOC
! shifting_constant = (E_{Hxc} - <\Psi | V_{Hxc} | \Psi>) / N_elec
! constant to add to the potential in order to obtain the variational energy as
! the eigenvalue of the effective long-range Hamiltonian
! (see original paper of Levy PRL 113, 113002 (2014), equation (17) )
END_DOC
integer :: istate
do istate = 1, N_states
shifting_constant(istate) = energy_x(istate) + energy_c(istate) + short_range_Hartree(istate) - Trace_v_Hxc(istate)
enddo
shifting_constant = shifting_constant / dble(elec_num)
END_PROVIDER

59
src/dft_utils_one_body/short_range_coulomb.irp.f
View File

@ -37,15 +37,15 @@ END_PROVIDER
BEGIN_PROVIDER [double precision, effective_one_e_potential, (mo_tot_num, mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, effective_one_e_potential_without_kin, (mo_tot_num, mo_tot_num,N_states)]
&BEGIN_PROVIDER [double precision, shifted_effective_one_e_potential_without_kin, (mo_tot_num, mo_tot_num,N_states)]
implicit none
integer :: i,j,istate
effective_one_e_potential = 0.d0
BEGIN_DOC
! effective_one_e_potential(i,j) = <i| v_{H}^{sr} |j> + <i| h_{core} |j> + <i|v_{xc} |j>
!
! Taking the expectation value does not provide any energy
!
! but effective_one_e_potential(i,j) is the potential coupling DFT and WFT part to be used in any WFT calculation
! shifted_effective_one_e_potential_without_kin = effective_one_e_potential_without_kin + shifting_constant on the diagonal
END_DOC
do istate = 1, N_states
do i = 1, mo_tot_num
@ -56,63 +56,8 @@ END_PROVIDER
effective_one_e_potential_without_kin(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_nucl_elec_integral(i,j) &
+ 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate) &
+ potential_x_beta_mo(i,j,istate) + potential_c_beta_mo(i,j,istate) )
shifted_effective_one_e_potential_without_kin(j,i,istate) = effective_one_e_potential_without_kin(j,i,istate)
enddo
enddo
do i = 1, mo_tot_num
shifted_effective_one_e_potential_without_kin(i,i,istate) += shifting_constant(istate)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, Fock_matrix_expectation_value]
implicit none
call get_average(effective_one_e_potential,one_body_dm_average_mo_for_dft,Fock_matrix_expectation_value)
END_PROVIDER
BEGIN_PROVIDER [double precision, Trace_v_xc, (N_states)]
&BEGIN_PROVIDER [double precision, Trace_v_H, (N_states)]
&BEGIN_PROVIDER [double precision, Trace_v_Hxc, (N_states)]
implicit none
integer :: i,j,istate
double precision :: dm
BEGIN_DOC
! Trace_v_xc = \sum_{i,j} (rho_{ij}_\alpha v^{xc}_{ij}^\alpha + rho_{ij}_\beta v^{xc}_{ij}^\beta)
! Trace_v_Hxc = \sum_{i,j} v^{H}_{ij} (rho_{ij}_\alpha + rho_{ij}_\beta)
! Trace_v_Hxc = \sum_{i,j} rho_{ij} v^{Hxc}_{ij}
END_DOC
do istate = 1, N_states
Trace_v_xc(istate) = 0.d0
Trace_v_H(istate) = 0.d0
do i = 1, mo_tot_num
do j = 1, mo_tot_num
Trace_v_xc(istate) += (potential_x_alpha_mo(j,i,istate) + potential_c_alpha_mo(j,i,istate)) * one_body_dm_mo_alpha_for_dft(j,i,istate)
Trace_v_xc(istate) += (potential_x_beta_mo(j,i,istate) + potential_c_beta_mo(j,i,istate) ) * one_body_dm_mo_beta_for_dft(j,i,istate)
dm = one_body_dm_mo_alpha_for_dft(j,i,istate) + one_body_dm_mo_beta_for_dft(j,i,istate)
Trace_v_H(istate) += dm * short_range_Hartree_operator(j,i,istate)
enddo
enddo
Trace_v_Hxc(istate) = Trace_v_xc(istate) + Trace_v_H(istate)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, DFT_one_e_energy_potential, (mo_tot_num, mo_tot_num,N_states)]
implicit none
integer :: i,j,istate
BEGIN_DOC
! one_e_energy_potential(i,j) = <i|h_{core}|j> + \int dr i(r)j(r) \int r' \rho(r') W_{ee}^{sr}
! If one take the expectation value over Psi, one gets the total one body energy
END_DOC
do istate = 1, N_states
do i = 1, mo_tot_num
do j = 1, mo_tot_num
DFT_one_e_energy_potential(j,i,istate) = mo_nucl_elec_integral(j,i) + mo_kinetic_integral(j,i) + short_range_Hartree_operator(j,i,istate) * 0.5d0
enddo
enddo
enddo
END_PROVIDER

2
src/dft_utils_two_body/INSTALL.cfg Normal file
View File

@ -0,0 +1,2 @@
[scripts]
qp_cipsi_rsh

4
src/dft_utils_two_body/NEED Normal file
View File

@ -0,0 +1,4 @@
dft_utils_one_body
integrals_bielec_erf
determinants
davidson_undressed

44
src/dft_utils_two_body/mr_dft_energy.irp.f Normal file
View File

@ -0,0 +1,44 @@
BEGIN_PROVIDER [double precision, electronic_energy_mr_dft, (N_states)]
implicit none
BEGIN_DOC
! Energy for the multi determinantal DFT calculation
END_DOC
print*,'You are using a variational method which uses the wave function stored in the EZFIO folder'
electronic_energy_mr_dft = total_range_separated_electronic_energy
END_PROVIDER
subroutine print_variational_energy_dft
implicit none
print*,'/////////////////////////'
print*, '****************************************'
print*,'///////////////////'
print*, ' Regular range separated DFT energy '
write(*, '(A22,X,F32.10)') 'mu_erf = ',mu_erf
write(*, '(A22,X,F16.10)') 'TOTAL ENERGY = ',electronic_energy_mr_dft+nuclear_repulsion
print*, ''
print*, 'Component of the energy ....'
print*, ''
write(*, '(A22,X,F16.10)') 'nuclear_repulsion = ',nuclear_repulsion
write(*, '(A22,X,F16.10)') 'psi_energy_erf = ',psi_energy_erf
write(*, '(A22,X,F16.10)') 'psi_dft_energy_h_cor= ',psi_dft_energy_h_core
write(*, '(A22,X,F16.10)') 'short_range_Hartree = ',short_range_Hartree
write(*, '(A22,X,F16.10)') 'two_elec_energy = ',two_elec_energy_dft
write(*, '(A22,X,F16.10)') 'energy_x = ',energy_x
write(*, '(A22,X,F16.10)') 'energy_c = ',energy_c
write(*, '(A22,X,F16.10)') 'E_xc = ',energy_x + energy_c
write(*, '(A22,X,F16.10)') 'E_Hxc = ',energy_x + energy_c + short_range_Hartree
print*, ''
print*, '****************************************'
print*, ''
write(*, '(A22,X,F16.10)') 'Approx eigenvalue = ',electronic_energy_mr_dft+nuclear_repulsion + Trace_v_Hxc - (short_range_Hartree + energy_x + energy_c)
write(*, '(A22,X,F16.10)') 'Trace_v_xc = ',Trace_v_xc
write(*, '(A22,X,F16.10)') 'Trace_v_Hxc = ',Trace_v_Hxc
write(*, '(A22,X,F16.10)') '<Psi| H | Psi> = ',psi_energy
write(*, '(A22,X,F16.10)') 'psi_energy_bielec = ',psi_energy_bielec
write(*, '(A22,X,F16.10)') 'psi_energy_h_core = ',psi_energy_h_core
end

16
src/dft_utils_two_body/print_rsdft_variational_energy.irp.f Normal file
View File

@ -0,0 +1,16 @@
program DFT_Utils_two_body_main
implicit none
read_wf = .true.
touch read_wf
disk_access_mo_one_integrals = "None"
touch disk_access_mo_one_integrals
disk_access_mo_integrals = "None"
touch disk_access_mo_integrals
disk_access_ao_integrals = "None"
touch disk_access_ao_integrals
density_for_dft = "WFT"
touch density_for_dft
call print_variational_energy_dft
end

81
src/dft_utils_two_body/psi_energy_erf.irp.f Normal file
View File

@ -0,0 +1,81 @@
BEGIN_PROVIDER [ double precision, psi_energy_erf, (N_states) ]
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <Psi|W_{ee}^{lr}|Psi>/<Psi|Psi>
!
END_DOC
integer :: i
call u_0_H_u_0_erf(psi_energy_erf,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size)
do i=N_det+1,N_states
psi_energy_erf(i) = 0.d0
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_energy_h_core_and_sr_hartree, (N_states) ]
implicit none
BEGIN_DOC
! psi_energy_h_core = <Psi| h_{core} + v_{H}^{sr}|Psi>
END_DOC
psi_energy_h_core_and_sr_hartree = psi_energy_h_core + short_range_Hartree
END_PROVIDER
BEGIN_PROVIDER [ double precision, total_range_separated_electronic_energy, (N_states) ]
implicit none
BEGIN_DOC
! Total_range_separated_electronic_energy = <Psi| h_{core} |Psi> + (1/2) <Psi| v_{H}^{sr} |Psi> + <i|W_{ee}^{lr}|i> + E_{x} + E_{c}
END_DOC
total_range_separated_electronic_energy = psi_energy_h_core + short_range_Hartree + psi_energy_erf + energy_x + energy_c
END_PROVIDER
BEGIN_PROVIDER [ double precision, two_elec_energy_dft, (N_states) ]
implicit none
BEGIN_DOC
! two_elec_energy_dft = (1/2) <Psi| v_{H}^{sr} |Psi> + <Psi|W_{ee}^{lr}|Psi>
END_DOC
two_elec_energy_dft = short_range_Hartree + psi_energy_erf
END_PROVIDER
BEGIN_PROVIDER [ double precision, ref_bitmask_energy_erf ]
&BEGIN_PROVIDER [ double precision, bi_elec_ref_bitmask_energy_erf ]
use bitmasks
implicit none
BEGIN_DOC
! Energy with the LONG RANGE INTERACTION of the reference bitmask used in Slater rules
END_DOC
integer :: occ(N_int*bit_kind_size,2)
integer :: i,j
call bitstring_to_list(ref_bitmask(1,1), occ(1,1), i, N_int)
call bitstring_to_list(ref_bitmask(1,2), occ(1,2), i, N_int)
ref_bitmask_energy_erf = 0.d0
bi_elec_ref_bitmask_energy_erf = 0.d0
do j= 1, elec_alpha_num
do i = j+1, elec_alpha_num
bi_elec_ref_bitmask_energy_erf += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1))
ref_bitmask_energy_erf += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1))
enddo
enddo
do j= 1, elec_beta_num
do i = j+1, elec_beta_num
bi_elec_ref_bitmask_energy_erf += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2))
ref_bitmask_energy_erf += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2))
enddo
do i= 1, elec_alpha_num
bi_elec_ref_bitmask_energy_erf += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2))
ref_bitmask_energy_erf += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2))
enddo
enddo
END_PROVIDER

0
scripts/script_CIPSI_RSH.sh → src/dft_utils_two_body/qp_cipsi_rsh
View File

26
src/dft_utils_two_body/routines_save_integrals_dft.irp.f Normal file
View File

@ -0,0 +1,26 @@
subroutine save_one_e_effective_potential
implicit none
BEGIN_DOC
! used to save the effective_one_e_potential into the one-body integrals in the ezfio folder
! this effective_one_e_potential is computed with the current density
! and will couple the WFT with DFT for the next regular WFT calculation
END_DOC
call ezfio_set_mo_one_e_integrals_integral_nuclear(effective_one_e_potential_without_kin)
call ezfio_set_mo_one_e_integrals_integral_kinetic(mo_kinetic_integral)
print *, 'Effective DFT potential is written on disk on the mo_ne_integral integrals'
call ezfio_set_mo_one_e_integrals_disk_access_mo_one_integrals("Read")
end
subroutine write_all_integrals_for_mrdft
implicit none
BEGIN_DOC
! saves all integrals needed for RS-DFT-MRCI calculation: one-body effective potential and two-elec erf integrals
END_DOC
call save_one_e_effective_potential
call save_bielec_ints_erf_mo_into_ints_mo
call save_bielec_ints_erf_mo_into_ints_ao
end

226
src/dft_utils_two_body/slater_rules_erf.irp.f Normal file
View File

@ -0,0 +1,226 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! THIS FILE CONTAINS EVERYTHING YOU NEED TO COMPUTE THE LONG RANGE PART OF THE INTERACTION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine i_H_j_erf(key_i,key_j,Nint,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|W_{ee}^{lr}|j> where i and j are determinants
! and the W_{ee}^{lr} is the long range two-body interaction
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij
integer :: exc(0:2,2,2)
integer :: degree
double precision :: get_mo_bielec_integral_erf
integer :: m,n,p,q
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem_erf, phase,phase_2
integer :: n_occ_ab(2)
PROVIDE mo_bielec_integrals_erf_in_map mo_integrals_erf_map big_array_exchange_integrals_erf
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
hij = 0.d0
!DIR$ FORCEINLINE
call get_excitation_degree(key_i,key_j,degree,Nint)
integer :: spin
select case (degree)
case (2)
call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
if(exc(1,1,1) == exc(1,2,2) )then
hij = phase * big_array_exchange_integrals_erf(exc(1,1,1),exc(1,1,2),exc(1,2,1))
else if (exc(1,2,1) ==exc(1,1,2))then
hij = phase * big_array_exchange_integrals_erf(exc(1,2,1),exc(1,1,1),exc(1,2,2))
else
hij = phase*get_mo_bielec_integral_erf( &
exc(1,1,1), &
exc(1,1,2), &
exc(1,2,1), &
exc(1,2,2) ,mo_integrals_erf_map)
endif
else if (exc(0,1,1) == 2) then
! Double alpha
hij = phase*(get_mo_bielec_integral_erf( &
exc(1,1,1), &
exc(2,1,1), &
exc(1,2,1), &
exc(2,2,1) ,mo_integrals_erf_map) - &
get_mo_bielec_integral_erf( &
exc(1,1,1), &
exc(2,1,1), &
exc(2,2,1), &
exc(1,2,1) ,mo_integrals_erf_map) )
else if (exc(0,1,2) == 2) then
! Double beta
hij = phase*(get_mo_bielec_integral_erf( &
exc(1,1,2), &
exc(2,1,2), &
exc(1,2,2), &
exc(2,2,2) ,mo_integrals_erf_map) - &
get_mo_bielec_integral_erf( &
exc(1,1,2), &
exc(2,1,2), &
exc(2,2,2), &
exc(1,2,2) ,mo_integrals_erf_map) )
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
if (exc(0,1,1) == 1) then
! Mono alpha
m = exc(1,1,1)
p = exc(1,2,1)
spin = 1
do i = 1, n_occ_ab(1)
hij += -big_array_exchange_integrals_erf(occ(i,1),m,p) + big_array_coulomb_integrals_erf(occ(i,1),m,p)
enddo
do i = 1, n_occ_ab(2)
hij += big_array_coulomb_integrals_erf(occ(i,2),m,p)
enddo
else
! Mono beta
m = exc(1,1,2)
p = exc(1,2,2)
spin = 2
do i = 1, n_occ_ab(2)
hij += -big_array_exchange_integrals_erf(occ(i,2),m,p) + big_array_coulomb_integrals_erf(occ(i,2),m,p)
enddo
do i = 1, n_occ_ab(1)
hij += big_array_coulomb_integrals_erf(occ(i,1),m,p)
enddo
endif
hij = hij * phase
case (0)
hij = diag_H_mat_elem_erf(key_i,Nint)
end select
end
double precision function diag_H_mat_elem_erf(key_i,Nint)
BEGIN_DOC
! returns <i|W_{ee}^{lr}|i> where |i> is a determinant and
! W_{ee}^{lr} is the two body long-range interaction
END_DOC
implicit none
integer(bit_kind), intent(in) :: key_i(N_int,2)
integer, intent(in) :: Nint
integer :: i,j
integer :: occ(Nint*bit_kind_size,2)
integer :: n_occ_ab(2)
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
diag_H_mat_elem_erf = 0.d0
! alpha - alpha
do i = 1, n_occ_ab(1)
do j = i+1, n_occ_ab(1)
diag_H_mat_elem_erf += mo_bielec_integral_erf_jj_anti(occ(i,1),occ(j,1))
enddo
enddo
! beta - beta
do i = 1, n_occ_ab(2)
do j = i+1, n_occ_ab(2)
diag_H_mat_elem_erf += mo_bielec_integral_erf_jj_anti(occ(i,2),occ(j,2))
enddo
enddo
! alpha - beta
do i = 1, n_occ_ab(1)
do j = 1, n_occ_ab(2)
diag_H_mat_elem_erf += mo_bielec_integral_erf_jj(occ(i,1),occ(j,2))
enddo
enddo
end
subroutine i_H_j_mono_spin_erf(key_i,key_j,Nint,spin,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|H|j> where i and j are determinants differing by a single excitation
END_DOC
integer, intent(in) :: Nint, spin
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij
integer :: exc(0:2,2)
double precision :: phase
PROVIDE big_array_exchange_integrals_erf mo_bielec_integrals_erf_in_map
call i_H_j_erf(key_i,key_j,Nint,hij)
end
subroutine i_H_j_double_spin_erf(key_i,key_j,Nint,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|H|j> where i and j are determinants differing by a same-spin double excitation
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint), key_j(Nint)
double precision, intent(out) :: hij
integer :: exc(0:2,2)
double precision :: phase
double precision, external :: get_mo_bielec_integral_erf
PROVIDE big_array_exchange_integrals_erf mo_bielec_integrals_erf_in_map
call get_double_excitation_spin(key_i,key_j,exc,phase,Nint)
hij = phase*(get_mo_bielec_integral_erf( &
exc(1,1), &
exc(2,1), &
exc(1,2), &
exc(2,2), mo_integrals_erf_map) - &
get_mo_bielec_integral_erf( &
exc(1,1), &
exc(2,1), &
exc(2,2), &
exc(1,2), mo_integrals_erf_map) )
end
subroutine i_H_j_double_alpha_beta_erf(key_i,key_j,Nint,hij)
use bitmasks
implicit none
BEGIN_DOC
! Returns <i|H|j> where i and j are determinants differing by an opposite-spin double excitation
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hij
integer :: exc(0:2,2,2)
double precision :: phase, phase2
double precision, external :: get_mo_bielec_integral_erf
PROVIDE big_array_exchange_integrals_erf mo_bielec_integrals_erf_in_map
call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint)
call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint)
phase = phase*phase2
if (exc(1,1,1) == exc(1,2,2)) then
hij = phase * big_array_exchange_integrals_erf(exc(1,1,1),exc(1,1,2),exc(1,2,1))
else if (exc(1,2,1) == exc(1,1,2)) then
hij = phase * big_array_exchange_integrals_erf(exc(1,2,1),exc(1,1,1),exc(1,2,2))
else
hij = phase*get_mo_bielec_integral_erf( &
exc(1,1,1), &
exc(1,1,2), &
exc(1,2,1), &
exc(1,2,2) ,mo_integrals_erf_map)
endif
end

482
src/dft_utils_two_body/u0_w_erf_u0.irp.f Normal file
View File

@ -0,0 +1,482 @@
subroutine u_0_H_u_0_erf(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|W_{ee}^{lr}|u_0>/<u_0|u_0>
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint, N_st, sze
double precision, intent(out) :: e_0(N_st)
double precision, intent(inout) :: u_0(sze,N_st)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision, allocatable :: v_0(:,:), s_0(:,:), u_1(:,:)
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
integer :: i,j
allocate (v_0(sze,N_st),s_0(sze,N_st))
call H_S2_u_0_erf_nstates_openmp(v_0,s_0,u_0,N_st,sze)
double precision :: norm
do i=1,N_st
norm = u_dot_u(u_0(1,i),n)
if (norm /= 0.d0) then
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
else
e_0(i) = 0.d0
endif
enddo
deallocate (s_0, v_0)
end
subroutine H_S2_u_0_erf_nstates_openmp(v_0,s_0,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(inout) :: v_0(sze,N_st), s_0(sze,N_st), u_0(sze,N_st)
integer :: k
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det),v_t(N_st,N_det),s_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
v_t = 0.d0
s_t = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_st)
call H_S2_u_0_erf_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,1,N_det,0,1)
deallocate(u_t)
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_st, N_det)
call dtranspose( &
s_t, &
size(s_t, 1), &
s_0, &
size(s_0, 1), &
N_st, N_det)
deallocate(v_t,s_t)
do k=1,N_st
call dset_order(v_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
call dset_order(s_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine H_S2_u_0_erf_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0> and s_0 = S^2 |u_0>
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
PROVIDE ref_bitmask_energy_erf N_int short_range_Hartree
select case (N_int)
case (1)
call H_S2_u_0_erf_nstates_openmp_work_1(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case (2)
call H_S2_u_0_erf_nstates_openmp_work_2(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case (3)
call H_S2_u_0_erf_nstates_openmp_work_3(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case (4)
call H_S2_u_0_erf_nstates_openmp_work_4(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
case default
call H_S2_u_0_erf_nstates_openmp_work_N_int(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine H_S2_u_0_erf_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_t = H|u_t> and s_t = S^2 |u_t>
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
double precision :: hij, sij
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b, m_a, m_b
integer :: istate
integer :: krow, kcol, krow_b, kcol_b
integer :: lrow, lcol
integer :: mrow, mcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev, nmax
integer*8 :: k8
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson
!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det, &
!$OMP psi_bilinear_matrix_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique, &
!$OMP n_det_alpha_unique, n_det_beta_unique, N_int, &
!$OMP psi_bilinear_matrix_transp_rows, &
!$OMP psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_transp_order, N_st, &
!$OMP psi_bilinear_matrix_order_transp_reverse, &
!$OMP psi_bilinear_matrix_columns_loc, &
!$OMP psi_bilinear_matrix_transp_rows_loc, &
!$OMP istart, iend, istep, irp_here, v_t, s_t, &
!$OMP ishift, idx0, u_t, maxab) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, &
!$OMP lcol, lrow, l_a, l_b, &
!$OMP buffer, doubles, n_doubles, &
!$OMP tmp_det2, hij, sij, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, &
!$OMP n_singles_b, k8)
! Alpha/Beta double excitations
! =============================
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
kcol_prev=-1
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
endif
kcol_prev = kcol
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call i_H_j_double_alpha_beta_erf(tmp_det,tmp_det2,$N_int,hij)
call get_s2(tmp_det,tmp_det2,$N_int,sij)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha excitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
call i_H_j_mono_spin_erf( tmp_det, tmp_det2, $N_int, 1, hij)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! single => sij = 0
enddo
enddo
! Compute Hij for all alpha doubles
! ----------------------------------
do i=1,n_doubles
l_a = doubles(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
call i_H_j_double_spin_erf( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! same spin => sij = 0
enddo
enddo
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
call i_H_j_mono_spin_erf( tmp_det, tmp_det2, $N_int, 2, hij)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! single => sij = 0
enddo
enddo
! Compute Hij for all beta doubles
! ----------------------------------
do i=1,n_doubles
l_b = doubles(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
call i_H_j_double_spin_erf( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! same spin => sij = 0
enddo
enddo
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
double precision, external :: diag_H_mat_elem_erf, diag_S_mat_elem
hij = diag_H_mat_elem_erf(tmp_det,$N_int)
sij = diag_S_mat_elem(tmp_det,$N_int)
do l=1,N_st
v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,k_a)
s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,k_a)
enddo
end do
!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx)
!$OMP END PARALLEL
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE

34
src/dft_utils_two_body/write_effective_rsdft_hamiltonian.irp.f Normal file
View File

@ -0,0 +1,34 @@
program write_effective_RSDFT_hamiltonian
implicit none
BEGIN_DOC
! This programs writes the effective RS-DFT Hamiltonian into the EZFIO folder.
! The next programs that will run unto the EZFIO folder will, by default, have the one- and two-body integrals loaded from the EZFIO data.
END_DOC
read_wf = .true.
touch read_wf
disk_access_mo_one_integrals = "None"
touch disk_access_mo_one_integrals
disk_access_mo_integrals = "None"
touch disk_access_mo_integrals
disk_access_ao_integrals = "None"
touch disk_access_ao_integrals
call routines_write_int
call routines_compute_energy
end
subroutine routines_write_int
implicit none
call write_all_integrals_for_mrdft
density_for_dft = "WFT"
touch density_for_dft
end
subroutine routines_compute_energy
implicit none
call print_variational_energy_dft
call ezfio_set_data_energy_and_density_data_one_body_alpha_dm_mo(one_body_dm_mo_alpha)
call ezfio_set_data_energy_and_density_data_one_body_beta_dm_mo(one_body_dm_mo_beta)
end

55
src/integrals_bielec/map_integrals.irp.f
View File

@ -533,6 +533,61 @@ subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map)
deallocate(pairs,hash,iorder,tmp_val)
end
subroutine get_mo_bielec_integrals_i1j1(k,l,sze,out_array,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple integrals <ik|jl> in the MO basis, all
! i(1)j(1) 1/r12 k(2)l(2)
! i, j for k,l fixed.
END_DOC
integer, intent(in) :: k,l, sze
double precision, intent(out) :: out_array(sze,sze)
type(map_type), intent(inout) :: map
integer :: i,j,kk,ll,m
integer(key_kind),allocatable :: hash(:)
integer ,allocatable :: pairs(:,:), iorder(:)
real(integral_kind), allocatable :: tmp_val(:)
PROVIDE mo_bielec_integrals_in_map
allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), &
tmp_val(sze*sze))
kk=0
out_array = 0.d0
do j=1,sze
do i=1,sze
kk += 1
!DIR$ FORCEINLINE
call bielec_integrals_index(i,k,j,l,hash(kk))
pairs(1,kk) = i
pairs(2,kk) = j
iorder(kk) = kk
enddo
enddo
logical :: integral_is_in_map
if (key_kind == 8) then
call i8radix_sort(hash,iorder,kk,-1)
else if (key_kind == 4) then
call iradix_sort(hash,iorder,kk,-1)
else if (key_kind == 2) then
call i2radix_sort(hash,iorder,kk,-1)
endif
call map_get_many(mo_integrals_map, hash, tmp_val, kk)
do ll=1,kk
m = iorder(ll)
i=pairs(1,m)
j=pairs(2,m)
out_array(i,j) = tmp_val(ll)
enddo
deallocate(pairs,hash,iorder,tmp_val)
end
subroutine get_mo_bielec_integrals_coulomb_ii(k,l,sze,out_val,map)
use map_module
implicit none

8
src/integrals_bielec_erf/routines_save_integrals_erf.irp.f
View File

@ -1,4 +1,4 @@
subroutine save_erf_bi_elec_integrals_mo
subroutine save_bielec_ints_erf_mo
implicit none
integer :: i,j,k,l
PROVIDE mo_bielec_integrals_erf_in_map
@ -7,7 +7,7 @@ subroutine save_erf_bi_elec_integrals_mo
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
end
subroutine save_erf_bi_elec_integrals_ao
subroutine save_bielec_ints_erf_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map
@ -16,7 +16,7 @@ subroutine save_erf_bi_elec_integrals_ao
call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
end
subroutine save_erf_bielec_ints_mo_into_ints_mo
subroutine save_bielec_ints_erf_mo_into_ints_mo
implicit none
integer :: i,j,k,l
PROVIDE mo_bielec_integrals_erf_in_map
@ -25,7 +25,7 @@ subroutine save_erf_bielec_ints_mo_into_ints_mo
call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
end
subroutine save_erf_bielec_ints_mo_into_ints_ao
subroutine save_bielec_ints_erf_mo_into_ints_ao
implicit none
integer :: i,j,k,l
PROVIDE ao_bielec_integrals_erf_in_map

4
src/integrals_bielec_erf/write_integrals_erf.irp.f
View File

@ -13,8 +13,8 @@ end
subroutine routine
implicit none
call save_erf_bi_elec_integrals_ao
call save_erf_bi_elec_integrals_mo
call save_bielec_ints_erf_ao
call save_bielec_ints_erf_mo
end

4
src/integrals_bielec_erf/write_ints_erf_to_regular_ints.irp.f
View File

@ -14,8 +14,8 @@ end
subroutine routine
implicit none
call save_erf_bielec_ints_mo_into_ints_ao
call save_erf_bielec_ints_mo_into_ints_mo
call save_bielec_ints_erf_mo_into_ints_ao
call save_bielec_ints_erf_mo_into_ints_mo
end

2
src/mo_one_e_integrals/read_write.irp.f
View File

@ -19,7 +19,7 @@
write_mo_one_integrals = .False.
else
print *, 'bielec_integrals/disk_access_mo_integrals has a wrong type'
print *, 'mo_one_e_integrals/disk_access_mo_one_integrals has a wrong type'
stop 1
endif

16
src/tools/diagonalize_h.irp.f Normal file
View File

@ -0,0 +1,16 @@
program diagonalize_h
implicit none
BEGIN_DOC
! program that extracts the N_states lowest states of the Hamiltonian within the set of Slater determinants stored in the EZFIO folder
END_DOC
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
call diagonalize_CI
print*,'N_det = ',N_det
call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
end