From dbdaeae65a82d44a4d065e88bb05d4e7c74f31a3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 21 Dec 2018 00:55:43 +0100 Subject: [PATCH] Merge with manu (#79) * added becke_numerical_grid/example.irp.f * added exemples for determinants and bitmask * added density_for_dft, dft_keywords and data_energy_and_density * added dft_utils_one_body * added dft_utils_two_body and it compilates * added script_CIPSI_RSH.sh * rm slave_cipsi * rm slave_cipsi * added dft_utils_two_body which works * added scf_utils * added slater_rules_mono_bielec.irp.f * remove integrals_bielec_erf * hatree_fock is ok with scf_utils * added kohn_sham_range_separated * added kohn_sham --- src/ao_basis/aos_transp.irp.f | 60 +++ src/ao_basis/aos_value.irp.f | 261 +++++++++- .../mono_excitations_bielec.irp.f | 136 +++++ src/determinants/psi_energy_mono_elec.irp.f | 1 + .../slater_rules_mono_bielec.irp.f | 361 +++++++++++++ src/dft_utils_one_body/NEED | 2 + src/dft_utils_one_body/README.rst | 8 - ...p.f => e_c_e_x_and_pot_general_prov.irp.f} | 68 ++- .../e_c_e_x_specific_prov.irp.f | 86 ++++ src/dft_utils_one_body/mu_erf_dft.irp.f | 5 + ...gy.irp.f => one_body_psi_energy_dft.irp.f} | 6 +- .../potentials_ao_on_grids.irp.f | 272 ---------- .../short_range_coulomb.irp.f | 57 +-- .../sr_e_c_e_x_specific_prov.irp.f | 86 ++++ .../sr_potentials_ao_specific_prov.irp.f | 188 +++++++ src/dft_utils_one_body/utils.irp.f | 59 ++- src/dft_utils_two_body/INSTALL.cfg | 2 + src/dft_utils_two_body/NEED | 3 + src/dft_utils_two_body/mr_dft_energy.irp.f | 44 ++ .../print_rsdft_variational_energy.irp.f | 16 + src/dft_utils_two_body/psi_energy_erf.irp.f | 81 +++ .../routines_save_integrals_dft.irp.f | 26 + src/dft_utils_two_body/slater_rules_erf.irp.f | 226 ++++++++ src/dft_utils_two_body/u0_w_erf_u0.irp.f | 482 ++++++++++++++++++ .../write_effective_rsdft_hamiltonian.irp.f | 34 ++ src/fci/selection.irp.f | 1 + src/hartree_fock/EZFIO.cfg | 53 +- src/hartree_fock/NEED | 5 +- src/hartree_fock/README.rst | 33 -- ...fock_matrix.irp.f => fock_matrix_hf.irp.f} | 174 +------ src/hartree_fock/hf_density_matrix_ao.irp.f | 41 -- src/hartree_fock/hf_energy.irp.f | 22 + src/hartree_fock/scf.irp.f | 2 +- src/hartree_fock/scf_old.irp.f | 61 +++ src/kohn_sham/NEED | 2 + src/kohn_sham/fock_matrix_ks.irp.f | 244 +++++++++ src/kohn_sham/ks_scf.irp.f | 95 ++++ src/kohn_sham/potential_functional.irp.f | 25 + src/kohn_sham_range_separated/NEED | 2 + .../fock_matrix_rs_ks.irp.f | 294 +++++++++++ .../potential_functional.irp.f | 25 + src/kohn_sham_range_separated/rs_ks_scf.irp.f | 103 ++++ src/mo_basis/track_orb.irp.f | 64 +++ src/mo_one_e_integrals/read_write.irp.f | 2 +- src/mo_two_e_integrals/map_integrals.irp.f | 55 ++ src/nuclei/nuclei.irp.f | 1 - src/scf_utils/EZFIO.cfg | 53 ++ src/scf_utils/NEED | 2 + src/scf_utils/README.rst | 6 + src/scf_utils/damping_scf.irp.f | 146 ++++++ .../diagonalize_fock.irp.f | 0 src/{hartree_fock => scf_utils}/diis.irp.f | 4 +- src/scf_utils/fock_matrix.irp.f | 144 ++++++ src/{hartree_fock => scf_utils}/huckel.irp.f | 0 .../roothaan_hall_scf.irp.f | 15 +- src/scf_utils/scf_density_matrix_ao.irp.f | 41 ++ src/tools/diagonalize_h.irp.f | 16 + 57 files changed, 3638 insertions(+), 663 deletions(-) create mode 100644 src/ao_basis/aos_transp.irp.f create mode 100644 src/determinants/mono_excitations_bielec.irp.f create mode 100644 src/determinants/slater_rules_mono_bielec.irp.f rename src/dft_utils_one_body/{e_c_e_x_and_potentials_mo.irp.f => e_c_e_x_and_pot_general_prov.irp.f} (65%) create mode 100644 src/dft_utils_one_body/e_c_e_x_specific_prov.irp.f create mode 100644 src/dft_utils_one_body/mu_erf_dft.irp.f rename src/dft_utils_one_body/{one_body_psi_energy.irp.f => one_body_psi_energy_dft.irp.f} (81%) delete mode 100644 src/dft_utils_one_body/potentials_ao_on_grids.irp.f create mode 100644 src/dft_utils_one_body/sr_e_c_e_x_specific_prov.irp.f create mode 100644 src/dft_utils_one_body/sr_potentials_ao_specific_prov.irp.f create mode 100644 src/dft_utils_two_body/INSTALL.cfg create mode 100644 src/dft_utils_two_body/NEED create mode 100644 src/dft_utils_two_body/mr_dft_energy.irp.f create mode 100644 src/dft_utils_two_body/print_rsdft_variational_energy.irp.f create mode 100644 src/dft_utils_two_body/psi_energy_erf.irp.f create mode 100644 src/dft_utils_two_body/routines_save_integrals_dft.irp.f create mode 100644 src/dft_utils_two_body/slater_rules_erf.irp.f create mode 100644 src/dft_utils_two_body/u0_w_erf_u0.irp.f create mode 100644 src/dft_utils_two_body/write_effective_rsdft_hamiltonian.irp.f delete mode 100644 src/hartree_fock/README.rst rename src/hartree_fock/{fock_matrix.irp.f => fock_matrix_hf.irp.f} (54%) delete mode 100644 src/hartree_fock/hf_density_matrix_ao.irp.f create mode 100644 src/hartree_fock/hf_energy.irp.f create mode 100644 src/hartree_fock/scf_old.irp.f create mode 100644 src/kohn_sham/NEED create mode 100644 src/kohn_sham/fock_matrix_ks.irp.f create mode 100644 src/kohn_sham/ks_scf.irp.f create mode 100644 src/kohn_sham/potential_functional.irp.f create mode 100644 src/kohn_sham_range_separated/NEED create mode 100644 src/kohn_sham_range_separated/fock_matrix_rs_ks.irp.f create mode 100644 src/kohn_sham_range_separated/potential_functional.irp.f create mode 100644 src/kohn_sham_range_separated/rs_ks_scf.irp.f create mode 100644 src/mo_basis/track_orb.irp.f create mode 100644 src/scf_utils/EZFIO.cfg create mode 100644 src/scf_utils/NEED create mode 100644 src/scf_utils/README.rst create mode 100644 src/scf_utils/damping_scf.irp.f rename src/{hartree_fock => scf_utils}/diagonalize_fock.irp.f (100%) rename src/{hartree_fock => scf_utils}/diis.irp.f (97%) create mode 100644 src/scf_utils/fock_matrix.irp.f rename src/{hartree_fock => scf_utils}/huckel.irp.f (100%) rename src/{hartree_fock => scf_utils}/roothaan_hall_scf.irp.f (96%) create mode 100644 src/scf_utils/scf_density_matrix_ao.irp.f create mode 100644 src/tools/diagonalize_h.irp.f diff --git a/src/ao_basis/aos_transp.irp.f b/src/ao_basis/aos_transp.irp.f new file mode 100644 index 00000000..f88938a2 --- /dev/null +++ b/src/ao_basis/aos_transp.irp.f @@ -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 + diff --git a/src/ao_basis/aos_value.irp.f b/src/ao_basis/aos_value.irp.f index 631065b9..f75f0544 100644 --- a/src/ao_basis/aos_value.irp.f +++ b/src/ao_basis/aos_value.irp.f @@ -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 + + diff --git a/src/determinants/mono_excitations_bielec.irp.f b/src/determinants/mono_excitations_bielec.irp.f new file mode 100644 index 00000000..dd3ff689 --- /dev/null +++ b/src/determinants/mono_excitations_bielec.irp.f @@ -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 + + diff --git a/src/determinants/psi_energy_mono_elec.irp.f b/src/determinants/psi_energy_mono_elec.irp.f index 494865db..64a4d71e 100644 --- a/src/determinants/psi_energy_mono_elec.irp.f +++ b/src/determinants/psi_energy_mono_elec.irp.f @@ -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) diff --git a/src/determinants/slater_rules_mono_bielec.irp.f b/src/determinants/slater_rules_mono_bielec.irp.f new file mode 100644 index 00000000..c9babbe3 --- /dev/null +++ b/src/determinants/slater_rules_mono_bielec.irp.f @@ -0,0 +1,361 @@ + +subroutine i_H_j_mono_spin_bielec(key_i,key_j,Nint,spin,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns 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 mo_bielec_integrals_in_map + + call get_mono_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) + call get_mono_excitation_from_fock_bielec(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij) +end + + +double precision function diag_H_mat_elem_bielec(det_in,Nint) + implicit none + BEGIN_DOC + ! Computes + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: particle(Nint,2) + integer :: i, nexc(2), ispin + integer :: occ_particle(Nint*bit_kind_size,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer(bit_kind) :: det_tmp(Nint,2) + integer :: na, nb + + ASSERT (Nint > 0) + ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num) + + nexc(1) = 0 + nexc(2) = 0 + do i=1,Nint + hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2)) + particle(i,1) = iand(hole(i,1),det_in(i,1)) + particle(i,2) = iand(hole(i,2),det_in(i,2)) + hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) + hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) + nexc(1) = nexc(1) + popcnt(hole(i,1)) + nexc(2) = nexc(2) + popcnt(hole(i,2)) + enddo + + diag_H_mat_elem_bielec = bi_elec_ref_bitmask_energy + if (nexc(1)+nexc(2) == 0) then + return + endif + + !call debug_det(det_in,Nint) + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(particle, occ_particle, tmp, Nint) + ASSERT (tmp(1) == nexc(1)) + ASSERT (tmp(2) == nexc(2)) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(hole, occ_hole, tmp, Nint) + ASSERT (tmp(1) == nexc(1)) + ASSERT (tmp(2) == nexc(2)) + + det_tmp = ref_bitmask + do ispin=1,2 + na = elec_num_tab(ispin) + nb = elec_num_tab(iand(ispin,1)+1) + do i=1,nexc(ispin) + !DIR$ FORCEINLINE + call ac_operator_bielec( occ_particle(i,ispin), ispin, det_tmp, diag_H_mat_elem_bielec, Nint,na,nb) + !DIR$ FORCEINLINE + call a_operator_bielec ( occ_hole (i,ispin), ispin, det_tmp, diag_H_mat_elem_bielec, Nint,na,nb) + enddo + enddo +end + + +subroutine a_operator_bielec(iorb,ispin,key,hjj,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Needed for diag_H_mat_elem + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hjj + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + integer :: tmp(2) + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k > 0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibclr(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + na = na-1 + + ! Same spin + do i=1,na + hjj = hjj - mo_bielec_integral_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + hjj = hjj - mo_bielec_integral_jj(occ(i,other_spin),iorb) + enddo + +end + + +subroutine ac_operator_bielec(iorb,ispin,key,hjj,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Needed for diag_H_mat_elem + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hjj + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + ASSERT (tmp(1) == elec_alpha_num) + ASSERT (tmp(2) == elec_beta_num) + + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k > 0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibset(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + + ! Same spin + do i=1,na + hjj = hjj + mo_bielec_integral_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + hjj = hjj + mo_bielec_integral_jj(occ(i,other_spin),iorb) + enddo + na = na+1 +end + + + +subroutine i_H_j_mono_spin_monoelec(key_i,key_j,Nint,spin,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns 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 + + call get_mono_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) + integer :: m,p + m = exc(1,1) + p = exc(1,2) + hij = phase * mo_mono_elec_integral(m,p) +end + + +double precision function diag_H_mat_elem_monoelec(det_in,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: particle(Nint,2) + integer :: i, nexc(2), ispin + integer :: occ_particle(Nint*bit_kind_size,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer(bit_kind) :: det_tmp(Nint,2) + integer :: na, nb + + ASSERT (Nint > 0) + ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num) + + diag_H_mat_elem_monoelec = 0.d0 + + !call debug_det(det_in,Nint) + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(det_in, occ_particle, tmp, Nint) + do ispin = 1,2 + do i = 1, tmp(ispin) + diag_H_mat_elem_monoelec += mo_mono_elec_integral(occ_particle(i,ispin),occ_particle(i,ispin)) + enddo + enddo + +end + +subroutine i_H_j_monoelec(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + 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 :: degree,m,p + double precision :: diag_H_mat_elem_monoelec,phase + integer :: exc(0:2,2,2) + call get_excitation_degree(key_i,key_j,degree,Nint) + hij = 0.d0 + if(degree>1)then + return + endif + if(degree==0)then + hij = diag_H_mat_elem_monoelec(key_i,N_int) + else + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + endif + hij = phase * mo_mono_elec_integral(m,p) + endif + +end + +subroutine i_H_j_bielec(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + 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 + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase,phase_2 + integer :: n_occ_ab(2) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map big_array_exchange_integrals bi_elec_ref_bitmask_energy + + 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(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(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_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 + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + endif + call get_mono_excitation_from_fock_bielec(key_i,key_j,p,m,spin,phase,hij) + case (0) + double precision :: diag_H_mat_elem_bielec + hij = diag_H_mat_elem_bielec(key_i,Nint) + end select +end + diff --git a/src/dft_utils_one_body/NEED b/src/dft_utils_one_body/NEED index 7a6314d7..3565cb66 100644 --- a/src/dft_utils_one_body/NEED +++ b/src/dft_utils_one_body/NEED @@ -2,3 +2,5 @@ density_for_dft dft_utils_on_grid mo_one_e_integrals mo_two_e_integrals +ao_one_e_integrals +ao_two_e_integrals diff --git a/src/dft_utils_one_body/README.rst b/src/dft_utils_one_body/README.rst index b51169bc..9b292d8e 100644 --- a/src/dft_utils_one_body/README.rst +++ b/src/dft_utils_one_body/README.rst @@ -2,11 +2,3 @@ RSDFT_Utils =========== -Needed Modules -============== -.. Do not edit this section It was auto-generated -.. by the `update_README.py` script. -Documentation -============= -.. Do not edit this section It was auto-generated -.. by the `update_README.py` script. diff --git a/src/dft_utils_one_body/e_c_e_x_and_potentials_mo.irp.f b/src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f similarity index 65% rename from src/dft_utils_one_body/e_c_e_x_and_potentials_mo.irp.f rename to src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f index 54802a84..e1a9439e 100644 --- a/src/dft_utils_one_body/e_c_e_x_and_potentials_mo.irp.f +++ b/src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f @@ -5,15 +5,15 @@ &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 - potential_x_alpha_ao = potential_x_alpha_ao_LDA - potential_x_beta_ao = potential_x_beta_ao_LDA + potential_x_alpha_ao = potential_sr_x_alpha_ao_LDA + potential_x_beta_ao = potential_sr_x_beta_ao_LDA else if(exchange_functional.EQ."short_range_PBE")then - potential_x_alpha_ao = potential_x_alpha_ao_PBE - potential_x_beta_ao = potential_x_beta_ao_PBE + potential_x_alpha_ao = potential_sr_x_alpha_ao_PBE + potential_x_beta_ao = potential_sr_x_beta_ao_PBE else if(exchange_functional.EQ."None")then potential_x_alpha_ao = 0.d0 potential_x_beta_ao = 0.d0 @@ -24,11 +24,11 @@ endif if(trim(correlation_functional)=="short_range_LDA")then - potential_c_alpha_ao = potential_c_alpha_ao_LDA - potential_c_beta_ao = potential_c_beta_ao_LDA + potential_c_alpha_ao = potential_sr_c_alpha_ao_LDA + potential_c_beta_ao = potential_sr_c_beta_ao_LDA else if(correlation_functional.EQ."short_range_PBE")then - potential_c_alpha_ao = potential_c_alpha_ao_PBE - potential_c_beta_ao = potential_c_beta_ao_PBE + potential_c_alpha_ao = potential_sr_c_alpha_ao_PBE + potential_c_beta_ao = potential_sr_c_beta_ao_PBE else if(correlation_functional.EQ."None")then potential_c_alpha_ao = 0.d0 potential_c_beta_ao = 0.d0 @@ -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,12 +92,15 @@ 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 + energy_x = energy_sr_x_LDA + energy_x = energy_sr_x_LDA else if(exchange_functional.EQ."short_range_PBE")then - energy_x = energy_x_PBE - energy_x = energy_x_PBE + energy_x = energy_sr_x_PBE + energy_x = energy_sr_x_PBE else if(exchange_functional.EQ."None")then energy_x = 0.d0 energy_x = 0.d0 @@ -108,11 +111,11 @@ END_PROVIDER endif if(trim(correlation_functional)=="short_range_LDA")then - energy_c = energy_c_LDA - energy_c = energy_c_LDA + energy_c = energy_sr_c_LDA + energy_c = energy_sr_c_LDA else if(correlation_functional.EQ."short_range_PBE")then - energy_c = energy_c_PBE - energy_c = energy_c_PBE + energy_c = energy_sr_c_PBE + energy_c = energy_sr_c_PBE else if(correlation_functional.EQ."None")then energy_c = 0.d0 energy_c = 0.d0 @@ -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 + diff --git a/src/dft_utils_one_body/e_c_e_x_specific_prov.irp.f b/src/dft_utils_one_body/e_c_e_x_specific_prov.irp.f new file mode 100644 index 00000000..a27eb167 --- /dev/null +++ b/src/dft_utils_one_body/e_c_e_x_specific_prov.irp.f @@ -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(rhoa(istate),rhob(istate),e_c,vc_a,vc_b) + call ex_LDA(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 + diff --git a/src/dft_utils_one_body/mu_erf_dft.irp.f b/src/dft_utils_one_body/mu_erf_dft.irp.f new file mode 100644 index 00000000..8999a429 --- /dev/null +++ b/src/dft_utils_one_body/mu_erf_dft.irp.f @@ -0,0 +1,5 @@ +BEGIN_PROVIDER [double precision, mu_erf_dft] + implicit none + mu_erf_dft = mu_erf + +END_PROVIDER diff --git a/src/dft_utils_one_body/one_body_psi_energy.irp.f b/src/dft_utils_one_body/one_body_psi_energy_dft.irp.f similarity index 81% rename from src/dft_utils_one_body/one_body_psi_energy.irp.f rename to src/dft_utils_one_body/one_body_psi_energy_dft.irp.f index d18cdcb6..8aca035b 100644 --- a/src/dft_utils_one_body/one_body_psi_energy.irp.f +++ b/src/dft_utils_one_body/one_body_psi_energy_dft.irp.f @@ -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 diff --git a/src/dft_utils_one_body/potentials_ao_on_grids.irp.f b/src/dft_utils_one_body/potentials_ao_on_grids.irp.f deleted file mode 100644 index b8c8b7ff..00000000 --- a/src/dft_utils_one_body/potentials_ao_on_grids.irp.f +++ /dev/null @@ -1,272 +0,0 @@ - BEGIN_PROVIDER[double precision, aos_vc_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)] -&BEGIN_PROVIDER[double precision, aos_vc_beta_LDA_w, (n_points_final_grid,ao_num,N_states)] -&BEGIN_PROVIDER[double precision, aos_vx_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)] -&BEGIN_PROVIDER[double precision, aos_vx_beta_LDA_w, (n_points_final_grid,ao_num,N_states)] - implicit none - BEGIN_DOC -! aos_vxc_alpha_LDA_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j) - 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)) - 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) - do j =1, ao_num - aos_vc_alpha_LDA_w(i,j,istate) = vc_a * aos_in_r_array(j,i)*weight - aos_vc_beta_LDA_w(i,j,istate) = vc_b * aos_in_r_array(j,i)*weight - aos_vx_alpha_LDA_w(i,j,istate) = vx_a * aos_in_r_array(j,i)*weight - aos_vx_beta_LDA_w(i,j,istate) = vx_b * aos_in_r_array(j,i)*weight - enddo - enddo - enddo - - 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)] -&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_LDA,(ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, potential_c_beta_ao_LDA,(ao_num,ao_num,N_states)] - implicit none - BEGIN_DOC -! short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis - END_DOC - integer :: istate - double precision :: wall_1,wall_2 - call wall_time(wall_1) - do istate = 1, N_states - call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vc_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_c_alpha_ao_LDA(1,1,istate),ao_num) - call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vc_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_c_beta_ao_LDA(1,1,istate),ao_num) - call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vx_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_x_alpha_ao_LDA(1,1,istate),ao_num) - call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_vx_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_x_beta_ao_LDA(1,1,istate),ao_num) - enddo - call wall_time(wall_2) - print*,'time to provide potential_x/c_alpha/beta_ao_LDA = ',wall_2 - wall_1 - - END_PROVIDER - - BEGIN_PROVIDER[double precision, aos_vc_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)] -&BEGIN_PROVIDER[double precision, aos_vc_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)] -&BEGIN_PROVIDER[double precision, aos_vx_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)] -&BEGIN_PROVIDER[double precision, aos_vx_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] -&BEGIN_PROVIDER[double precision, grad_aos_dvc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] -&BEGIN_PROVIDER[double precision, grad_aos_dvc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] -&BEGIN_PROVIDER[double precision, grad_aos_dvx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] -&BEGIN_PROVIDER[double precision, grad_aos_dvx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] - implicit none - BEGIN_DOC -! aos_vxc_alpha_PBE_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j) - 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)) - allocate(contrib_grad_xa(3,N_states),contrib_grad_xb(3,N_states),contrib_grad_ca(3,N_states),contrib_grad_cb(3,N_states)) - 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 ) - vx_rho_a(istate) *= weight - vc_rho_a(istate) *= weight - vx_rho_b(istate) *= weight - vc_rho_b(istate) *= weight - do m= 1,3 - contrib_grad_ca(m,istate) = weight * (2.d0 * vc_grad_rho_a_2(istate) * grad_rho_a(m,istate) + vc_grad_rho_a_b(istate) * grad_rho_b(m,istate)) - contrib_grad_xa(m,istate) = weight * (2.d0 * vx_grad_rho_a_2(istate) * grad_rho_a(m,istate) + vx_grad_rho_a_b(istate) * grad_rho_b(m,istate)) - contrib_grad_cb(m,istate) = weight * (2.d0 * vc_grad_rho_b_2(istate) * grad_rho_b(m,istate) + vc_grad_rho_a_b(istate) * grad_rho_a(m,istate)) - contrib_grad_xb(m,istate) = weight * (2.d0 * vx_grad_rho_b_2(istate) * grad_rho_b(m,istate) + vx_grad_rho_a_b(istate) * grad_rho_a(m,istate)) - enddo - do j = 1, ao_num - aos_vc_alpha_PBE_w(j,i,istate) = vc_rho_a(istate) * aos_in_r_array(j,i) - aos_vc_beta_PBE_w (j,i,istate) = vc_rho_b(istate) * aos_in_r_array(j,i) - aos_vx_alpha_PBE_w(j,i,istate) = vx_rho_a(istate) * aos_in_r_array(j,i) - aos_vx_beta_PBE_w (j,i,istate) = vx_rho_b(istate) * aos_in_r_array(j,i) - do m = 1,3 - aos_dvc_alpha_PBE_w(j,i,m,istate) = contrib_grad_ca(m,istate) * aos_in_r_array(j,i) - aos_dvc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_in_r_array(j,i) - aos_dvx_alpha_PBE_w(j,i,m,istate) = contrib_grad_xa(m,istate) * aos_in_r_array(j,i) - aos_dvx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_in_r_array(j,i) - - grad_aos_dvc_alpha_PBE_w (j,i,m,istate) = contrib_grad_ca(m,istate) * aos_grad_in_r_array(j,i,m) - grad_aos_dvc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_grad_in_r_array(j,i,m) - grad_aos_dvx_alpha_PBE_w (j,i,m,istate) = contrib_grad_xa(m,istate) * aos_grad_in_r_array(j,i,m) - grad_aos_dvx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_grad_in_r_array(j,i,m) - enddo - enddo - 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 - - 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)] -&BEGIN_PROVIDER [double precision, potential_c_beta_ao_PBE,(ao_num,ao_num,N_states)] - implicit none - BEGIN_DOC -! exchange/correlation alpha/beta potentials with the short range PBE functional on the AO basis - END_DOC - integer :: istate, m - double precision :: wall_1,wall_2 - call wall_time(wall_1) - potential_c_alpha_ao_PBE = 0.d0 - potential_x_alpha_ao_PBE = 0.d0 - potential_c_beta_ao_PBE = 0.d0 - potential_x_beta_ao_PBE = 0.d0 - do istate = 1, N_states - ! correlation alpha - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vc_alpha_PBE_w(1,1,istate),size(aos_vc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1)) - ! correlation beta - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vc_beta_PBE_w(1,1,istate),size(aos_vc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1)) - ! exchange alpha - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vx_alpha_PBE_w(1,1,istate),size(aos_vx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1)) - ! exchange beta - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_vx_beta_PBE_w(1,1,istate),size(aos_vx_beta_PBE_w,1), aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate), size(potential_x_beta_ao_PBE,1)) - do m= 1,3 - ! correlation alpha - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvc_alpha_PBE_w(1,1,m,istate),size(aos_dvc_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1)) - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvc_alpha_PBE_w(1,1,m,istate),size(grad_aos_dvc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_alpha_ao_PBE(1,1,istate),size(potential_c_alpha_ao_PBE,1)) - ! correlation beta - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvc_beta_PBE_w(1,1,m,istate),size(aos_dvc_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1)) - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvc_beta_PBE_w(1,1,m,istate),size(grad_aos_dvc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_c_beta_ao_PBE(1,1,istate),size(potential_c_beta_ao_PBE,1)) - ! exchange alpha - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvx_alpha_PBE_w(1,1,m,istate),size(aos_dvx_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1)) - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvx_alpha_PBE_w(1,1,m,istate),size(grad_aos_dvx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_alpha_ao_PBE(1,1,istate),size(potential_x_alpha_ao_PBE,1)) - ! exchange beta - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dvx_beta_PBE_w(1,1,m,istate),size(aos_dvx_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate),size(potential_x_beta_ao_PBE,1)) - call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dvx_beta_PBE_w(1,1,m,istate),size(grad_aos_dvx_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_x_beta_ao_PBE(1,1,istate),size(potential_x_beta_ao_PBE,1)) - enddo - enddo - - call wall_time(wall_2) - -END_PROVIDER diff --git a/src/dft_utils_one_body/short_range_coulomb.irp.f b/src/dft_utils_one_body/short_range_coulomb.irp.f index 703d988a..51863983 100644 --- a/src/dft_utils_one_body/short_range_coulomb.irp.f +++ b/src/dft_utils_one_body/short_range_coulomb.irp.f @@ -37,7 +37,6 @@ 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 @@ -56,63 +55,9 @@ 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) = + \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 + diff --git a/src/dft_utils_one_body/sr_e_c_e_x_specific_prov.irp.f b/src/dft_utils_one_body/sr_e_c_e_x_specific_prov.irp.f new file mode 100644 index 00000000..6c54769b --- /dev/null +++ b/src/dft_utils_one_body/sr_e_c_e_x_specific_prov.irp.f @@ -0,0 +1,86 @@ + + + BEGIN_PROVIDER[double precision, energy_sr_x_LDA, (N_states) ] +&BEGIN_PROVIDER[double precision, energy_sr_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_sr_x_LDA = 0.d0 + energy_sr_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_dft,rhoa(istate),rhob(istate),e_c,vc_a,vc_b) + call ex_LDA_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,vx_a,vx_b) + energy_sr_x_LDA(istate) += weight * e_x + energy_sr_c_LDA(istate) += weight * e_c + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER[double precision, energy_sr_x_PBE, (N_states) ] +&BEGIN_PROVIDER[double precision, energy_sr_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_sr_x_PBE = 0.d0 + energy_sr_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_sr_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_sr_x_PBE += ex * weight + energy_sr_c_PBE += ec * weight + enddo + enddo + + +END_PROVIDER + diff --git a/src/dft_utils_one_body/sr_potentials_ao_specific_prov.irp.f b/src/dft_utils_one_body/sr_potentials_ao_specific_prov.irp.f new file mode 100644 index 00000000..8e937a35 --- /dev/null +++ b/src/dft_utils_one_body/sr_potentials_ao_specific_prov.irp.f @@ -0,0 +1,188 @@ + BEGIN_PROVIDER[double precision, aos_sr_vc_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vc_beta_LDA_w, (n_points_final_grid,ao_num,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vx_alpha_LDA_w, (n_points_final_grid,ao_num,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vx_beta_LDA_w, (n_points_final_grid,ao_num,N_states)] + implicit none + BEGIN_DOC +! aos_sr_vxc_alpha_LDA_w(j,i) = ao_i(r_j) * (sr_v^x_alpha(r_j) + sr_v^c_alpha(r_j)) * W(r_j) + END_DOC + integer :: istate,i,j + double precision :: r(3) + double precision :: mu,weight + double precision :: e_c,sr_vc_a,sr_vc_b,e_x,sr_vx_a,sr_vx_b + double precision, allocatable :: rhoa(:),rhob(:) + allocate(rhoa(N_states), rhob(N_states)) + 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_dft,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b) + call ex_LDA_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b) + do j =1, ao_num + aos_sr_vc_alpha_LDA_w(i,j,istate) = sr_vc_a * aos_in_r_array(j,i)*weight + aos_sr_vc_beta_LDA_w(i,j,istate) = sr_vc_b * aos_in_r_array(j,i)*weight + aos_sr_vx_alpha_LDA_w(i,j,istate) = sr_vx_a * aos_in_r_array(j,i)*weight + aos_sr_vx_beta_LDA_w(i,j,istate) = sr_vx_b * aos_in_r_array(j,i)*weight + enddo + enddo + enddo + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, potential_sr_x_alpha_ao_LDA,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, potential_sr_x_beta_ao_LDA,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, potential_sr_c_alpha_ao_LDA,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, potential_sr_c_beta_ao_LDA,(ao_num,ao_num,N_states)] + implicit none + BEGIN_DOC +! short range exchange/correlation alpha/beta potentials with LDA functional on the AO basis + END_DOC + integer :: istate + double precision :: wall_1,wall_2 + call wall_time(wall_1) + do istate = 1, N_states + call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_sr_vc_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_sr_c_alpha_ao_LDA(1,1,istate),ao_num) + call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_sr_vc_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_sr_c_beta_ao_LDA(1,1,istate),ao_num) + call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_sr_vx_alpha_LDA_w(1,1,istate),n_points_final_grid,0.d0,potential_sr_x_alpha_ao_LDA(1,1,istate),ao_num) + call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0,aos_in_r_array,ao_num,aos_sr_vx_beta_LDA_w(1,1,istate) ,n_points_final_grid,0.d0,potential_sr_x_beta_ao_LDA(1,1,istate),ao_num) + enddo + call wall_time(wall_2) + print*,'time to provide potential_sr_x/c_alpha/beta_ao_LDA = ',wall_2 - wall_1 + + END_PROVIDER + + BEGIN_PROVIDER[double precision, aos_sr_vc_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vc_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vx_alpha_PBE_w , (ao_num,n_points_final_grid,N_states)] !(n_points_final_grid,ao_num,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vx_beta_PBE_w , (ao_num,n_points_final_grid,N_states)]!(n_points_final_grid,ao_num,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] +&BEGIN_PROVIDER[double precision, grad_aos_dsr_vc_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] +&BEGIN_PROVIDER[double precision, grad_aos_dsr_vc_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] +&BEGIN_PROVIDER[double precision, grad_aos_dsr_vx_alpha_PBE_w , (ao_num,n_points_final_grid,3,N_states)] +&BEGIN_PROVIDER[double precision, grad_aos_dsr_vx_beta_PBE_w , (ao_num,n_points_final_grid,3,N_states)] + implicit none + BEGIN_DOC +! aos_vxc_alpha_PBE_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j) + 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 :: sr_vc_rho_a(:), sr_vc_rho_b(:), sr_vx_rho_a(:), sr_vx_rho_b(:) + double precision, allocatable :: sr_vx_grad_rho_a_2(:), sr_vx_grad_rho_b_2(:), sr_vx_grad_rho_a_b(:), sr_vc_grad_rho_a_2(:), sr_vc_grad_rho_b_2(:), sr_vc_grad_rho_a_b(:) + allocate(sr_vc_rho_a(N_states), sr_vc_rho_b(N_states), sr_vx_rho_a(N_states), sr_vx_rho_b(N_states)) + allocate(sr_vx_grad_rho_a_2(N_states), sr_vx_grad_rho_b_2(N_states), sr_vx_grad_rho_a_b(N_states), sr_vc_grad_rho_a_2(N_states), sr_vc_grad_rho_b_2(N_states), sr_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)) + allocate(contrib_grad_xa(3,N_states),contrib_grad_xb(3,N_states),contrib_grad_ca(3,N_states),contrib_grad_cb(3,N_states)) + 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_sr_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + ex,sr_vx_rho_a,sr_vx_rho_b,sr_vx_grad_rho_a_2,sr_vx_grad_rho_b_2,sr_vx_grad_rho_a_b, & ! outputs correlation + ec,sr_vc_rho_a,sr_vc_rho_b,sr_vc_grad_rho_a_2,sr_vc_grad_rho_b_2,sr_vc_grad_rho_a_b ) + sr_vx_rho_a(istate) *= weight + sr_vc_rho_a(istate) *= weight + sr_vx_rho_b(istate) *= weight + sr_vc_rho_b(istate) *= weight + do m= 1,3 + contrib_grad_ca(m,istate) = weight * (2.d0 * sr_vc_grad_rho_a_2(istate) * grad_rho_a(m,istate) + sr_vc_grad_rho_a_b(istate) * grad_rho_b(m,istate)) + contrib_grad_xa(m,istate) = weight * (2.d0 * sr_vx_grad_rho_a_2(istate) * grad_rho_a(m,istate) + sr_vx_grad_rho_a_b(istate) * grad_rho_b(m,istate)) + contrib_grad_cb(m,istate) = weight * (2.d0 * sr_vc_grad_rho_b_2(istate) * grad_rho_b(m,istate) + sr_vc_grad_rho_a_b(istate) * grad_rho_a(m,istate)) + contrib_grad_xb(m,istate) = weight * (2.d0 * sr_vx_grad_rho_b_2(istate) * grad_rho_b(m,istate) + sr_vx_grad_rho_a_b(istate) * grad_rho_a(m,istate)) + enddo + do j = 1, ao_num + aos_sr_vc_alpha_PBE_w(j,i,istate) = sr_vc_rho_a(istate) * aos_in_r_array(j,i) + aos_sr_vc_beta_PBE_w (j,i,istate) = sr_vc_rho_b(istate) * aos_in_r_array(j,i) + aos_sr_vx_alpha_PBE_w(j,i,istate) = sr_vx_rho_a(istate) * aos_in_r_array(j,i) + aos_sr_vx_beta_PBE_w (j,i,istate) = sr_vx_rho_b(istate) * aos_in_r_array(j,i) + do m = 1,3 + aos_dsr_vc_alpha_PBE_w(j,i,m,istate) = contrib_grad_ca(m,istate) * aos_in_r_array(j,i) + aos_dsr_vc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_in_r_array(j,i) + aos_dsr_vx_alpha_PBE_w(j,i,m,istate) = contrib_grad_xa(m,istate) * aos_in_r_array(j,i) + aos_dsr_vx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_in_r_array(j,i) + + grad_aos_dsr_vc_alpha_PBE_w (j,i,m,istate) = contrib_grad_ca(m,istate) * aos_grad_in_r_array(j,i,m) + grad_aos_dsr_vc_beta_PBE_w (j,i,m,istate) = contrib_grad_cb(m,istate) * aos_grad_in_r_array(j,i,m) + grad_aos_dsr_vx_alpha_PBE_w (j,i,m,istate) = contrib_grad_xa(m,istate) * aos_grad_in_r_array(j,i,m) + grad_aos_dsr_vx_beta_PBE_w (j,i,m,istate) = contrib_grad_xb(m,istate) * aos_grad_in_r_array(j,i,m) + enddo + enddo + enddo + enddo + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, potential_sr_x_alpha_ao_PBE,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, potential_sr_x_beta_ao_PBE,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, potential_sr_c_alpha_ao_PBE,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, potential_sr_c_beta_ao_PBE,(ao_num,ao_num,N_states)] + implicit none + BEGIN_DOC +! exchange/correlation alpha/beta potentials with the short range PBE functional on the AO basis + END_DOC + integer :: istate, m + double precision :: wall_1,wall_2 + call wall_time(wall_1) + potential_sr_c_alpha_ao_PBE = 0.d0 + potential_sr_x_alpha_ao_PBE = 0.d0 + potential_sr_c_beta_ao_PBE = 0.d0 + potential_sr_x_beta_ao_PBE = 0.d0 + do istate = 1, N_states + ! correlation alpha + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_sr_vc_alpha_PBE_w(1,1,istate),size(aos_sr_vc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_c_alpha_ao_PBE(1,1,istate),size(potential_sr_c_alpha_ao_PBE,1)) + ! correlation beta + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_sr_vc_beta_PBE_w(1,1,istate),size(aos_sr_vc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_c_beta_ao_PBE(1,1,istate),size(potential_sr_c_beta_ao_PBE,1)) + ! exchange alpha + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_sr_vx_alpha_PBE_w(1,1,istate),size(aos_sr_vx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_x_alpha_ao_PBE(1,1,istate),size(potential_sr_x_alpha_ao_PBE,1)) + ! exchange beta + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_sr_vx_beta_PBE_w(1,1,istate),size(aos_sr_vx_beta_PBE_w,1), aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_x_beta_ao_PBE(1,1,istate), size(potential_sr_x_beta_ao_PBE,1)) + do m= 1,3 + ! correlation alpha + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dsr_vc_alpha_PBE_w(1,1,m,istate),size(aos_dsr_vc_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_sr_c_alpha_ao_PBE(1,1,istate),size(potential_sr_c_alpha_ao_PBE,1)) + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dsr_vc_alpha_PBE_w(1,1,m,istate),size(grad_aos_dsr_vc_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_c_alpha_ao_PBE(1,1,istate),size(potential_sr_c_alpha_ao_PBE,1)) + ! correlation beta + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dsr_vc_beta_PBE_w(1,1,m,istate),size(aos_dsr_vc_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_sr_c_beta_ao_PBE(1,1,istate),size(potential_sr_c_beta_ao_PBE,1)) + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dsr_vc_beta_PBE_w(1,1,m,istate),size(grad_aos_dsr_vc_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_c_beta_ao_PBE(1,1,istate),size(potential_sr_c_beta_ao_PBE,1)) + ! exchange alpha + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dsr_vx_alpha_PBE_w(1,1,m,istate),size(aos_dsr_vx_alpha_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_sr_x_alpha_ao_PBE(1,1,istate),size(potential_sr_x_alpha_ao_PBE,1)) + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dsr_vx_alpha_PBE_w(1,1,m,istate),size(grad_aos_dsr_vx_alpha_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_x_alpha_ao_PBE(1,1,istate),size(potential_sr_x_alpha_ao_PBE,1)) + ! exchange beta + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,aos_dsr_vx_beta_PBE_w(1,1,m,istate),size(aos_dsr_vx_beta_PBE_w,1),aos_grad_in_r_array(1,1,m),size(aos_grad_in_r_array,1),1.d0,potential_sr_x_beta_ao_PBE(1,1,istate),size(potential_sr_x_beta_ao_PBE,1)) + call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0,grad_aos_dsr_vx_beta_PBE_w(1,1,m,istate),size(grad_aos_dsr_vx_beta_PBE_w,1),aos_in_r_array,size(aos_in_r_array,1),1.d0,potential_sr_x_beta_ao_PBE(1,1,istate),size(potential_sr_x_beta_ao_PBE,1)) + enddo + enddo + + call wall_time(wall_2) + +END_PROVIDER diff --git a/src/dft_utils_one_body/utils.irp.f b/src/dft_utils_one_body/utils.irp.f index cfa261be..9775ebb0 100644 --- a/src/dft_utils_one_body/utils.irp.f +++ b/src/dft_utils_one_body/utils.irp.f @@ -1,5 +1,5 @@ -subroutine GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & +subroutine GGA_sr_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) implicit none @@ -10,7 +10,7 @@ subroutine GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho double precision :: r2(3),dr2(3), local_potential,r12,dx2,mu do istate = 1, N_states if(exchange_functional.EQ."short_range_PBE")then - call ex_pbe_sr(mu_erf,rho_a(istate),rho_b(istate),grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),ex(istate),vx_rho_a(istate),vx_rho_b(istate),vx_grad_rho_a_2(istate),vx_grad_rho_b_2(istate),vx_grad_rho_a_b(istate)) + call ex_pbe_sr(mu_erf_dft,rho_a(istate),rho_b(istate),grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),ex(istate),vx_rho_a(istate),vx_rho_b(istate),vx_grad_rho_a_2(istate),vx_grad_rho_b_2(istate),vx_grad_rho_a_b(istate)) else if(exchange_functional.EQ."None")then ex = 0.d0 vx_rho_a = 0.d0 @@ -30,7 +30,60 @@ subroutine GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc) call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco) - call ec_pbe_sr(mu_erf,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec(istate),vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo) + call ec_pbe_sr(mu_erf_dft,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec(istate),vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo) + + call v_rho_oc_to_v_rho_ab(vrhoo,vrhoc,vc_rho_a(istate),vc_rho_b(istate)) + call v_grad_rho_oc_to_v_grad_rho_ab(vsigmaoo,vsigmacc,vsigmaco,vc_grad_rho_a_2(istate),vc_grad_rho_b_2(istate),vc_grad_rho_a_b(istate)) + else if(correlation_functional.EQ."None")then + ec = 0.d0 + vc_rho_a = 0.d0 + vc_rho_b = 0.d0 + vc_grad_rho_a_2 = 0.d0 + vc_grad_rho_a_b = 0.d0 + vc_grad_rho_b_2 = 0.d0 + else + print*, 'Correlation functional required does not exist ...' + print*, 'correlation_functional',correlation_functional + stop + endif + enddo +end + + +subroutine GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & + ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & + ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) + implicit none + double precision, intent(in) :: r(3),rho_a(N_states),rho_b(N_states),grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states) + double precision, intent(out) :: ex(N_states),vx_rho_a(N_states),vx_rho_b(N_states),vx_grad_rho_a_2(N_states),vx_grad_rho_b_2(N_states),vx_grad_rho_a_b(N_states) + double precision, intent(out) :: ec(N_states),vc_rho_a(N_states),vc_rho_b(N_states),vc_grad_rho_a_2(N_states),vc_grad_rho_b_2(N_states),vc_grad_rho_a_b(N_states) + integer :: istate + double precision :: r2(3),dr2(3), local_potential,r12,dx2,mu + double precision :: mu_local + mu_local = 1.d+9 + do istate = 1, N_states + if(exchange_functional.EQ."short_range_PBE")then + call ex_pbe_sr(mu_local,rho_a(istate),rho_b(istate),grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),ex(istate),vx_rho_a(istate),vx_rho_b(istate),vx_grad_rho_a_2(istate),vx_grad_rho_b_2(istate),vx_grad_rho_a_b(istate)) + else if(exchange_functional.EQ."None")then + ex = 0.d0 + vx_rho_a = 0.d0 + vx_rho_b = 0.d0 + vx_grad_rho_a_2 = 0.d0 + vx_grad_rho_a_b = 0.d0 + vx_grad_rho_b_2 = 0.d0 + else + print*, 'Exchange functional required does not exist ...' + print*,'exchange_functional',exchange_functional + stop + endif + + double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo + if(correlation_functional.EQ."short_range_PBE")then + ! convertion from (alpha,beta) formalism to (closed, open) formalism + call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc) + call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco) + + call ec_pbe_sr(mu_local,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec(istate),vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo) call v_rho_oc_to_v_rho_ab(vrhoo,vrhoc,vc_rho_a(istate),vc_rho_b(istate)) call v_grad_rho_oc_to_v_grad_rho_ab(vsigmaoo,vsigmacc,vsigmaco,vc_grad_rho_a_2(istate),vc_grad_rho_b_2(istate),vc_grad_rho_a_b(istate)) diff --git a/src/dft_utils_two_body/INSTALL.cfg b/src/dft_utils_two_body/INSTALL.cfg new file mode 100644 index 00000000..2b126bf9 --- /dev/null +++ b/src/dft_utils_two_body/INSTALL.cfg @@ -0,0 +1,2 @@ +[scripts] +qp_cipsi_rsh diff --git a/src/dft_utils_two_body/NEED b/src/dft_utils_two_body/NEED new file mode 100644 index 00000000..8db4e8fc --- /dev/null +++ b/src/dft_utils_two_body/NEED @@ -0,0 +1,3 @@ +dft_utils_one_body +determinants +davidson_undressed diff --git a/src/dft_utils_two_body/mr_dft_energy.irp.f b/src/dft_utils_two_body/mr_dft_energy.irp.f new file mode 100644 index 00000000..ae85ce8e --- /dev/null +++ b/src/dft_utils_two_body/mr_dft_energy.irp.f @@ -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_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 + diff --git a/src/dft_utils_two_body/print_rsdft_variational_energy.irp.f b/src/dft_utils_two_body/print_rsdft_variational_energy.irp.f new file mode 100644 index 00000000..e00e3b4b --- /dev/null +++ b/src/dft_utils_two_body/print_rsdft_variational_energy.irp.f @@ -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 diff --git a/src/dft_utils_two_body/psi_energy_erf.irp.f b/src/dft_utils_two_body/psi_energy_erf.irp.f new file mode 100644 index 00000000..9c995ec7 --- /dev/null +++ b/src/dft_utils_two_body/psi_energy_erf.irp.f @@ -0,0 +1,81 @@ + +BEGIN_PROVIDER [ double precision, psi_energy_erf, (N_states) ] + use bitmasks + implicit none + BEGIN_DOC + ! Computes e_0 = / + ! + 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 = + 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 = + (1/2) + + 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) + + 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 + diff --git a/src/dft_utils_two_body/routines_save_integrals_dft.irp.f b/src/dft_utils_two_body/routines_save_integrals_dft.irp.f new file mode 100644 index 00000000..6809de86 --- /dev/null +++ b/src/dft_utils_two_body/routines_save_integrals_dft.irp.f @@ -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_erf_bi_elec_integrals_mo + call save_erf_bi_elec_integrals_ao +end + diff --git a/src/dft_utils_two_body/slater_rules_erf.irp.f b/src/dft_utils_two_body/slater_rules_erf.irp.f new file mode 100644 index 00000000..91c15c07 --- /dev/null +++ b/src/dft_utils_two_body/slater_rules_erf.irp.f @@ -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 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 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 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 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 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 + + diff --git a/src/dft_utils_two_body/u0_w_erf_u0.irp.f b/src/dft_utils_two_body/u0_w_erf_u0.irp.f new file mode 100644 index 00000000..daa1501d --- /dev/null +++ b/src/dft_utils_two_body/u0_w_erf_u0.irp.f @@ -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 = / + ! + ! 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 + + diff --git a/src/dft_utils_two_body/write_effective_rsdft_hamiltonian.irp.f b/src/dft_utils_two_body/write_effective_rsdft_hamiltonian.irp.f new file mode 100644 index 00000000..342e83b5 --- /dev/null +++ b/src/dft_utils_two_body/write_effective_rsdft_hamiltonian.irp.f @@ -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 + + diff --git a/src/fci/selection.irp.f b/src/fci/selection.irp.f index 84608f6a..db6b0484 100644 --- a/src/fci/selection.irp.f +++ b/src/fci/selection.irp.f @@ -1375,3 +1375,4 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) enddo end + diff --git a/src/hartree_fock/EZFIO.cfg b/src/hartree_fock/EZFIO.cfg index 22d2447e..35f5c69b 100644 --- a/src/hartree_fock/EZFIO.cfg +++ b/src/hartree_fock/EZFIO.cfg @@ -1,53 +1,6 @@ -[max_dim_diis] -type: integer -doc: Maximum size of the |DIIS| extrapolation procedure -interface: ezfio,provider,ocaml -default: 15 - -[threshold_diis] +[energy] type: Threshold -doc: Threshold on the convergence of the |DIIS| error vector during a Hartree-Fock calculation. If 0. is chosen, the square root of thresh_scf will be used. -interface: ezfio,provider,ocaml +doc: Energy HF +interface: ezfio default: 0. -[thresh_scf] -type: Threshold -doc: Threshold on the convergence of the Hartree Fock energy. -interface: ezfio,provider,ocaml -default: 1.e-10 - -[n_it_scf_max] -type: Strictly_positive_int -doc: Maximum number of |SCF| iterations -interface: ezfio,provider,ocaml -default: 500 - -[level_shift] -type: Positive_float -doc: Initial value of the energy shift on the virtual |MOs| -interface: ezfio,provider,ocaml -default: 0.0 - -[scf_algorithm] -type: character*(32) -doc: Type of |SCF| algorithm used. Possible choices are [ Simple | DIIS] -interface: ezfio,provider,ocaml -default: DIIS - -[mo_guess_type] -type: MO_guess -doc: Initial MO guess. Can be [ Huckel | HCore ] -interface: ezfio,provider,ocaml -default: Huckel - -[energy] -type: double precision -doc: Calculated HF energy -interface: ezfio - -[no_oa_or_av_opt] -type: logical -doc: If |true|, skip the (inactive+core) --> (active) and the (active) --> (virtual) orbital rotations within the |SCF| procedure -interface: ezfio,provider,ocaml -default: False - diff --git a/src/hartree_fock/NEED b/src/hartree_fock/NEED index c45dc538..233ffdb5 100644 --- a/src/hartree_fock/NEED +++ b/src/hartree_fock/NEED @@ -1,4 +1 @@ -ao_one_e_integrals -ao_two_e_integrals -mo_guess -bitmask +scf_utils ao_one_e_integrals ao_two_e_integrals diff --git a/src/hartree_fock/README.rst b/src/hartree_fock/README.rst deleted file mode 100644 index 373a611d..00000000 --- a/src/hartree_fock/README.rst +++ /dev/null @@ -1,33 +0,0 @@ -============ -Hartree-Fock -============ - - -The Hartree-Fock module performs *Restricted* Hartree-Fock calculations (the -spatial part of the |MOs| is common for alpha and beta spinorbitals). - -The Hartree-Fock program does the following: - -#. Compute/Read all the one- and two-electron integrals, and store them in memory -#. Check in the |EZFIO| database if there is a set of |MOs|. If there is, it - will read them as initial guess. Otherwise, it will create a guess. -#. Perform the |SCF| iterations - -At each iteration, the |MOs| are saved in the |EZFIO| database. Hence, if the calculation -crashes for any unexpected reason, the calculation can be restarted by running again -the |SCF| with the same |EZFIO| database. - -The `DIIS`_ algorithm is implemented, as well as the `level-shifting`_ method. -If the |SCF| does not converge, try again with a higher value of :option:`level_shift`. - -To start a calculation from scratch, the simplest way is to remove the -``mo_basis`` directory from the |EZFIO| database, and run the |SCF| again. - - - - -.. _DIIS: https://en.wikipedia.org/w/index.php?title=DIIS -.. _level-shifting: https://doi.org/10.1002/qua.560070407 - - - diff --git a/src/hartree_fock/fock_matrix.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f similarity index 54% rename from src/hartree_fock/fock_matrix.irp.f rename to src/hartree_fock/fock_matrix_hf.irp.f index 7f473f7a..97bbc614 100644 --- a/src/hartree_fock/fock_matrix.irp.f +++ b/src/hartree_fock/fock_matrix_hf.irp.f @@ -1,103 +1,3 @@ - BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num,mo_tot_num) ] -&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo, (mo_tot_num)] - implicit none - BEGIN_DOC - ! Fock matrix on the MO basis. - ! For open shells, the ROHF Fock Matrix is - ! - ! | F-K | F + K/2 | F | - ! |---------------------------------| - ! | F + K/2 | F | F - K/2 | - ! |---------------------------------| - ! | F | F - K/2 | F + K | - ! - ! F = 1/2 (Fa + Fb) - ! - ! K = Fb - Fa - ! - END_DOC - integer :: i,j,n - if (elec_alpha_num == elec_beta_num) then - Fock_matrix_mo = Fock_matrix_mo_alpha - else - - do j=1,elec_beta_num - ! F-K - do i=1,elec_beta_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F+K/2 - do i=elec_beta_num+1,elec_alpha_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F - do i=elec_alpha_num+1, mo_tot_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) - enddo - enddo - - do j=elec_beta_num+1,elec_alpha_num - ! F+K/2 - do i=1,elec_beta_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F - do i=elec_beta_num+1,elec_alpha_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) - enddo - ! F-K/2 - do i=elec_alpha_num+1, mo_tot_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - enddo - - do j=elec_alpha_num+1, mo_tot_num - ! F - do i=1,elec_beta_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) - enddo - ! F-K/2 - do i=elec_beta_num+1,elec_alpha_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F+K - do i=elec_alpha_num+1,mo_tot_num - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) & - + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - enddo - - endif - - do i = 1, mo_tot_num - Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) - enddo -END_PROVIDER - - - - BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ] -&BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num, ao_num) ] - implicit none - BEGIN_DOC - ! Alpha Fock matrix in AO basis set - END_DOC - - integer :: i,j - do j=1,ao_num - do i=1,ao_num - Fock_matrix_ao_alpha(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_alpha(i,j) - Fock_matrix_ao_beta (i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_beta (i,j) - enddo - enddo - -END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_alpha, (ao_num, ao_num) ] &BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_beta , (ao_num, ao_num) ] @@ -123,7 +23,7 @@ END_PROVIDER !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s,i0,j0,k0,l0, & !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp, c0, c1, c2, & !$OMP local_threshold)& - !$OMP SHARED(ao_num,HF_density_matrix_ao_alpha,HF_density_matrix_ao_beta,& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& !$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, & !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) @@ -170,9 +70,9 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - c0 = HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l) - c1 = HF_density_matrix_ao_alpha(k,i) - c2 = HF_density_matrix_ao_beta(k,i) + c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l) + c1 = SCF_density_matrix_ao_alpha(k,i) + c2 = SCF_density_matrix_ao_beta(k,i) if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then cycle endif @@ -209,7 +109,7 @@ END_PROVIDER !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & !$OMP n_elements,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)& - !$OMP SHARED(ao_num,HF_density_matrix_ao_alpha,HF_density_matrix_ao_beta,& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& !$OMP ao_integrals_map, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) @@ -235,12 +135,12 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * values(k1) + integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1) ao_bi_elec_integral_alpha_tmp(i,j) += integral ao_bi_elec_integral_beta_tmp (i,j) += integral integral = values(k1) - ao_bi_elec_integral_alpha_tmp(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral - ao_bi_elec_integral_beta_tmp (l,j) -= HF_density_matrix_ao_beta (k,i) * integral + ao_bi_elec_integral_alpha_tmp(l,j) -= SCF_density_matrix_ao_alpha(k,i) * integral + ao_bi_elec_integral_beta_tmp (l,j) -= SCF_density_matrix_ao_beta (k,i) * integral enddo enddo enddo @@ -258,63 +158,19 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num,mo_tot_num) ] - implicit none - BEGIN_DOC - ! Fock matrix on the MO basis - END_DOC - call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & - Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1)) -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num,mo_tot_num) ] - implicit none - BEGIN_DOC - ! Fock matrix on the MO basis - END_DOC - call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & - Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1)) -END_PROVIDER - -BEGIN_PROVIDER [ double precision, HF_energy ] + BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num, ao_num) ] implicit none BEGIN_DOC - ! Hartree-Fock energy + ! Alpha Fock matrix in AO basis set END_DOC - HF_energy = nuclear_repulsion - + integer :: i,j do j=1,ao_num do i=1,ao_num - HF_energy += 0.5d0 * ( & - (ao_mono_elec_integral(i,j) + Fock_matrix_ao_alpha(i,j) ) * HF_density_matrix_ao_alpha(i,j) +& - (ao_mono_elec_integral(i,j) + Fock_matrix_ao_beta (i,j) ) * HF_density_matrix_ao_beta (i,j) ) + Fock_matrix_ao_alpha(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_alpha(i,j) + Fock_matrix_ao_beta (i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_beta (i,j) enddo enddo - + END_PROVIDER - - -BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] - implicit none - BEGIN_DOC - ! Fock matrix in AO basis set - END_DOC - - if ( (elec_alpha_num == elec_beta_num).and. & - (level_shift == 0.) ) & - then - integer :: i,j - do j=1,ao_num - do i=1,ao_num - Fock_matrix_ao(i,j) = Fock_matrix_ao_alpha(i,j) - enddo - enddo - else - call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & - Fock_matrix_ao,size(Fock_matrix_ao,1)) - endif -END_PROVIDER - - diff --git a/src/hartree_fock/hf_density_matrix_ao.irp.f b/src/hartree_fock/hf_density_matrix_ao.irp.f deleted file mode 100644 index 736d1a97..00000000 --- a/src/hartree_fock/hf_density_matrix_ao.irp.f +++ /dev/null @@ -1,41 +0,0 @@ -BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ] - implicit none - BEGIN_DOC - ! S^{-1}.P_alpha.S^{-1} - END_DOC - - call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & - mo_coef, size(mo_coef,1), & - mo_coef, size(mo_coef,1), 0.d0, & - HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num,ao_num) ] - implicit none - BEGIN_DOC - ! S^{-1}.P_beta.S^{-1} - END_DOC - - call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & - mo_coef, size(mo_coef,1), & - mo_coef, size(mo_coef,1), 0.d0, & - HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1)) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ] - implicit none - BEGIN_DOC - ! S^{-1}.P.S^{-1} where P = C.C^t - END_DOC - ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1)) - if (elec_alpha_num== elec_beta_num) then - HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_alpha - else - ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_beta ,1)) - HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_beta - endif - -END_PROVIDER - diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f new file mode 100644 index 00000000..a4fe52be --- /dev/null +++ b/src/hartree_fock/hf_energy.irp.f @@ -0,0 +1,22 @@ +BEGIN_PROVIDER [double precision, extra_energy_contrib_from_density] + implicit none + extra_energy_contrib_from_density = 0.D0 + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, HF_energy] +&BEGIN_PROVIDER [ double precision, HF_two_electron_energy] +&BEGIN_PROVIDER [ double precision, HF_one_electron_energy] + implicit none + integer :: i,j + HF_energy = nuclear_repulsion + do j=1,ao_num + do i=1,ao_num + HF_two_electron_energy += 0.5d0 * ( ao_bi_elec_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) & + +ao_bi_elec_integral_beta(i,j) * SCF_density_matrix_ao_beta(i,j) ) + HF_one_electron_energy += ao_mono_elec_integral(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + enddo + enddo + HF_energy += HF_two_electron_energy + HF_one_electron_energy +END_PROVIDER + diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index a04fcc31..c8e474c0 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -47,7 +47,7 @@ subroutine run double precision :: EHF integer :: i_it, i, j, k - EHF = HF_energy + EHF = SCF_energy mo_label = "Canonical" diff --git a/src/hartree_fock/scf_old.irp.f b/src/hartree_fock/scf_old.irp.f new file mode 100644 index 00000000..852fa128 --- /dev/null +++ b/src/hartree_fock/scf_old.irp.f @@ -0,0 +1,61 @@ +program scf + BEGIN_DOC +! Produce `Hartree_Fock` MO orbital +! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ +! output: hartree_fock.energy +! optional: mo_basis.mo_coef + END_DOC + call create_guess + call orthonormalize_mos + call run +end + +subroutine create_guess + implicit none + BEGIN_DOC +! Create a MO guess if no MOs are present in the EZFIO directory + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_mo_basis_mo_coef(exists) + if (.not.exists) then + if (mo_guess_type == "HCore") then + mo_coef = ao_ortho_lowdin_coef + TOUCH mo_coef + mo_label = 'Guess' + call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label) + SOFT_TOUCH mo_coef mo_label + else if (mo_guess_type == "Huckel") then + call huckel_guess + else + print *, 'Unrecognized MO guess type : '//mo_guess_type + stop 1 + endif + endif +end + +subroutine run + + BEGIN_DOC +! Run SCF calculation + END_DOC + + use bitmasks + implicit none + + double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem + double precision :: EHF + integer :: i_it, i, j, k + + EHF = SCF_energy + + mo_label = "Canonical" + +! Choose SCF algorithm + + call damping_SCF ! Deprecated routine +! call Roothaan_Hall_SCF + +end + + diff --git a/src/kohn_sham/NEED b/src/kohn_sham/NEED new file mode 100644 index 00000000..83407613 --- /dev/null +++ b/src/kohn_sham/NEED @@ -0,0 +1,2 @@ +dft_utils_one_body +scf_utils diff --git a/src/kohn_sham/fock_matrix_ks.irp.f b/src/kohn_sham/fock_matrix_ks.irp.f new file mode 100644 index 00000000..fe4c29db --- /dev/null +++ b/src/kohn_sham/fock_matrix_ks.irp.f @@ -0,0 +1,244 @@ + + + BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_alpha, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_beta , (ao_num, ao_num) ] + use map_module + implicit none + BEGIN_DOC + ! Alpha Fock matrix in ao basis set + END_DOC + + integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 + integer*8 :: p,q + double precision :: integral, c0, c1, c2 + double precision :: ao_bielec_integral, local_threshold + double precision, allocatable :: ao_bi_elec_integral_alpha_tmp(:,:) + double precision, allocatable :: ao_bi_elec_integral_beta_tmp(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_beta_tmp + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_alpha_tmp + + ao_bi_elec_integral_alpha = 0.d0 + ao_bi_elec_integral_beta = 0.d0 + if (do_direct_integrals) then + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s,i0,j0,k0,l0, & + !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp, c0, c1, c2, & + !$OMP local_threshold)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& + !$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, & + !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) + + allocate(keys(1), values(1)) + allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), & + ao_bi_elec_integral_beta_tmp(ao_num,ao_num)) + ao_bi_elec_integral_alpha_tmp = 0.d0 + ao_bi_elec_integral_beta_tmp = 0.d0 + + q = ao_num*ao_num*ao_num*ao_num + !$OMP DO SCHEDULE(dynamic) + do p=1_8,q + call bielec_integrals_index_reverse(kk,ii,ll,jj,p) + if ( (kk(1)>ao_num).or. & + (ii(1)>ao_num).or. & + (jj(1)>ao_num).or. & + (ll(1)>ao_num) ) then + cycle + endif + k = kk(1) + i = ii(1) + l = ll(1) + j = jj(1) + + if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & + < ao_integrals_threshold) then + cycle + endif + local_threshold = ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j) + if (local_threshold < ao_integrals_threshold) then + cycle + endif + i0 = i + j0 = j + k0 = k + l0 = l + values(1) = 0.d0 + local_threshold = ao_integrals_threshold/local_threshold + do k2=1,8 + if (kk(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l) + c1 = SCF_density_matrix_ao_alpha(k,i) + c2 = SCF_density_matrix_ao_beta(k,i) + if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then + cycle + endif + if (values(1) == 0.d0) then + values(1) = ao_bielec_integral(k0,l0,i0,j0) + endif + integral = c0 * values(1) + ao_bi_elec_integral_alpha_tmp(i,j) += integral + ao_bi_elec_integral_beta_tmp (i,j) += integral + integral = values(1) + ao_bi_elec_integral_alpha_tmp(l,j) -= c1 * integral + ao_bi_elec_integral_beta_tmp (l,j) -= c2 * integral + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_bi_elec_integral_alpha += ao_bi_elec_integral_alpha_tmp + !$OMP END CRITICAL + !$OMP CRITICAL + ao_bi_elec_integral_beta += ao_bi_elec_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp) + !$OMP END PARALLEL + else + PROVIDE ao_bielec_integrals_in_map + + integer(omp_lock_kind) :: lck(ao_num) + integer*8 :: i8 + integer :: ii(8), jj(8), kk(8), ll(8), k2 + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + + + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& + !$OMP ao_integrals_map, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta,HF_exchange) + + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), & + ao_bi_elec_integral_beta_tmp(ao_num,ao_num)) + ao_bi_elec_integral_alpha_tmp = 0.d0 + ao_bi_elec_integral_beta_tmp = 0.d0 + + !$OMP DO SCHEDULE(dynamic,64) + !DIR$ NOVECTOR + do i8=0_8,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + do k1=1,n_elements + call bielec_integrals_index_reverse(kk,ii,ll,jj,keys(k1)) + + do k2=1,8 + if (kk(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1) + ao_bi_elec_integral_alpha_tmp(i,j) += integral + ao_bi_elec_integral_beta_tmp (i,j) += integral + integral = values(k1) + ao_bi_elec_integral_alpha_tmp(l,j) -= HF_exchange * (SCF_density_matrix_ao_alpha(k,i) * integral) + ao_bi_elec_integral_beta_tmp (l,j) -= HF_exchange * (SCF_density_matrix_ao_beta (k,i) * integral) + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_bi_elec_integral_alpha += ao_bi_elec_integral_alpha_tmp + !$OMP END CRITICAL + !$OMP CRITICAL + ao_bi_elec_integral_beta += ao_bi_elec_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp) + !$OMP END PARALLEL + + endif + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Alpha Fock matrix in ao basis set + END_DOC + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao_alpha(i,j) = Fock_matrix_alpha_no_xc_ao(i,j) + ao_potential_alpha_xc(i,j) + Fock_matrix_ao_beta(i,j) = Fock_matrix_beta_no_xc_ao(i,j) + ao_potential_beta_xc(i,j) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_no_xc_ao, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_no_xc_ao, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Mono electronic an Coulomb matrix in ao basis set + END_DOC + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_alpha_no_xc_ao(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_alpha(i,j) + Fock_matrix_beta_no_xc_ao(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_beta (i,j) + enddo + enddo + +END_PROVIDER + + + + + BEGIN_PROVIDER [ double precision, KS_energy] +&BEGIN_PROVIDER [ double precision, two_electron_energy] +&BEGIN_PROVIDER [ double precision, one_electron_energy] +&BEGIN_PROVIDER [ double precision, Fock_matrix_energy] +&BEGIN_PROVIDER [ double precision, trace_potential_xc ] + implicit none + BEGIN_DOC + ! Hartree-Fock energy + END_DOC + + integer :: i,j + double precision :: accu_mono,accu_fock + KS_energy = nuclear_repulsion + one_electron_energy = 0.d0 + two_electron_energy = 0.d0 + Fock_matrix_energy = 0.d0 + trace_potential_xc = 0.d0 + do j=1,ao_num + do i=1,ao_num + Fock_matrix_energy += Fock_matrix_ao_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) + & + Fock_matrix_ao_beta(i,j) * SCF_density_matrix_ao_beta(i,j) + two_electron_energy += 0.5d0 * ( ao_bi_elec_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) & + +ao_bi_elec_integral_beta(i,j) * SCF_density_matrix_ao_beta(i,j) ) + one_electron_energy += ao_mono_elec_integral(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) +! possible bug fix for open-shell +! trace_potential_xc += (ao_potential_alpha_xc(i,j) + ao_potential_beta_xc(i,j) ) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + trace_potential_xc += ao_potential_alpha_xc(i,j) * SCF_density_matrix_ao_alpha(i,j) + ao_potential_beta_xc(i,j) * SCF_density_matrix_ao_beta (i,j) + enddo + enddo + + KS_energy += e_exchange_dft + e_correlation_dft + one_electron_energy + two_electron_energy +END_PROVIDER + +BEGIN_PROVIDER [double precision, extra_energy_contrib_from_density] + implicit none +! possible bug fix for open-shell: +! extra_energy_contrib_from_density = e_exchange_dft + e_correlation_dft - 0.25d0 * trace_potential_xc + extra_energy_contrib_from_density = e_exchange_dft + e_correlation_dft - 0.5d0 * trace_potential_xc +END_PROVIDER + diff --git a/src/kohn_sham/ks_scf.irp.f b/src/kohn_sham/ks_scf.irp.f new file mode 100644 index 00000000..aad5527e --- /dev/null +++ b/src/kohn_sham/ks_scf.irp.f @@ -0,0 +1,95 @@ +program srs_ks_cf + BEGIN_DOC +! Produce `Kohn_Sham` MO orbital +! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ +! output: kohn_sham.energy +! optional: mo_basis.mo_coef + END_DOC + read_wf = .False. + density_for_dft ="WFT" + touch density_for_dft + touch read_wf + print*, '**************************' + print*, 'mu_erf_dft = ',mu_erf_dft + print*, '**************************' + call check_coherence_functional + call create_guess + call orthonormalize_mos + call run +end + +subroutine check_coherence_functional + implicit none + integer :: ifound_x,ifound_c + if(exchange_functional.eq."None")then + ifound_x = 1 + else + ifound_x = index(exchange_functional,"short_range") + endif + + if(correlation_functional.eq."None")then + ifound_c = 1 + else + ifound_c = index(correlation_functional,"short_range") + endif + print*,ifound_x,ifound_c + if(ifound_x .eq.0 .or. ifound_c .eq. 0)then + print*,'YOU ARE USING THE RANGE SEPARATED KS PROGRAM BUT YOUR INPUT KEYWORD FOR ' + print*,'exchange_functional is ',exchange_functional + print*,'correlation_functional is ',correlation_functional + print*,'CHANGE THE exchange_functional and correlation_functional keywords to range separated functionals' + print*,'or switch to the KS_SCF program that uses regular functionals' + stop + endif + +end + + + +subroutine create_guess + implicit none + BEGIN_DOC +! Create a MO guess if no MOs are present in the EZFIO directory + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_mo_basis_mo_coef(exists) + if (.not.exists) then + if (mo_guess_type == "HCore") then + mo_coef = ao_ortho_lowdin_coef + TOUCH mo_coef + mo_label = 'Guess' + call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.) + SOFT_TOUCH mo_coef mo_label + else if (mo_guess_type == "Huckel") then + call huckel_guess + else + print *, 'Unrecognized MO guess type : '//mo_guess_type + stop 1 + endif + endif +end + +subroutine run + + BEGIN_DOC +! Run SCF calculation + END_DOC + + use bitmasks + implicit none + + double precision :: EHF + + EHF = KS_energy + + mo_label = "Canonical" + +! Choose SCF algorithm + +! call damping_SCF ! Deprecated routine + call Roothaan_Hall_SCF + +end + + diff --git a/src/kohn_sham/potential_functional.irp.f b/src/kohn_sham/potential_functional.irp.f new file mode 100644 index 00000000..c0ef5d74 --- /dev/null +++ b/src/kohn_sham/potential_functional.irp.f @@ -0,0 +1,25 @@ + BEGIN_PROVIDER [double precision, ao_potential_alpha_xc, (ao_num, ao_num)] +&BEGIN_PROVIDER [double precision, ao_potential_beta_xc, (ao_num, ao_num)] + implicit none + integer :: i,j,k,l + ao_potential_alpha_xc = 0.d0 + ao_potential_beta_xc = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + ao_potential_alpha_xc(i,j) = potential_c_alpha_ao(i,j,1) + potential_x_alpha_ao(i,j,1) + ao_potential_beta_xc(i,j) = potential_c_beta_ao(i,j,1) + potential_x_beta_ao(i,j,1) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, e_exchange_dft] + implicit none + e_exchange_dft = energy_x(1) + +END_PROVIDER + +BEGIN_PROVIDER [double precision, e_correlation_dft] + implicit none + e_correlation_dft = energy_c(1) + +END_PROVIDER diff --git a/src/kohn_sham_range_separated/NEED b/src/kohn_sham_range_separated/NEED new file mode 100644 index 00000000..643ecece --- /dev/null +++ b/src/kohn_sham_range_separated/NEED @@ -0,0 +1,2 @@ +dft_utils_one_body +scf_utils diff --git a/src/kohn_sham_range_separated/fock_matrix_rs_ks.irp.f b/src/kohn_sham_range_separated/fock_matrix_rs_ks.irp.f new file mode 100644 index 00000000..fdc26dba --- /dev/null +++ b/src/kohn_sham_range_separated/fock_matrix_rs_ks.irp.f @@ -0,0 +1,294 @@ + BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_alpha, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_beta , (ao_num, ao_num) ] + use map_module + implicit none + BEGIN_DOC + ! Alpha Fock matrix in AO basis set + END_DOC + + integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 + integer*8 :: p,q + double precision :: integral, c0, c1, c2 + double precision :: ao_bielec_integral, local_threshold + double precision, allocatable :: ao_bi_elec_integral_alpha_tmp(:,:) + double precision, allocatable :: ao_bi_elec_integral_beta_tmp(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_beta_tmp + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_alpha_tmp + + ao_bi_elec_integral_alpha = 0.d0 + ao_bi_elec_integral_beta = 0.d0 + if (do_direct_integrals) then + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s,i0,j0,k0,l0, & + !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp, c0, c1, c2, & + !$OMP local_threshold)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& + !$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, & + !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) + + allocate(keys(1), values(1)) + allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), & + ao_bi_elec_integral_beta_tmp(ao_num,ao_num)) + ao_bi_elec_integral_alpha_tmp = 0.d0 + ao_bi_elec_integral_beta_tmp = 0.d0 + + q = ao_num*ao_num*ao_num*ao_num + !$OMP DO SCHEDULE(dynamic) + do p=1_8,q + call bielec_integrals_index_reverse(kk,ii,ll,jj,p) + if ( (kk(1)>ao_num).or. & + (ii(1)>ao_num).or. & + (jj(1)>ao_num).or. & + (ll(1)>ao_num) ) then + cycle + endif + k = kk(1) + i = ii(1) + l = ll(1) + j = jj(1) + + if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) & + < ao_integrals_threshold) then + cycle + endif + local_threshold = ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j) + if (local_threshold < ao_integrals_threshold) then + cycle + endif + i0 = i + j0 = j + k0 = k + l0 = l + values(1) = 0.d0 + local_threshold = ao_integrals_threshold/local_threshold + do k2=1,8 + if (kk(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l) + c1 = SCF_density_matrix_ao_alpha(k,i) + c2 = SCF_density_matrix_ao_beta(k,i) + if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then + cycle + endif + if (values(1) == 0.d0) then + values(1) = ao_bielec_integral(k0,l0,i0,j0) + endif + integral = c0 * values(1) + ao_bi_elec_integral_alpha_tmp(i,j) += integral + ao_bi_elec_integral_beta_tmp (i,j) += integral + integral = values(1) + ao_bi_elec_integral_alpha_tmp(l,j) -= c1 * integral + ao_bi_elec_integral_beta_tmp (l,j) -= c2 * integral + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_bi_elec_integral_alpha += ao_bi_elec_integral_alpha_tmp + !$OMP END CRITICAL + !$OMP CRITICAL + ao_bi_elec_integral_beta += ao_bi_elec_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp) + !$OMP END PARALLEL + else + PROVIDE ao_bielec_integrals_in_map + PROVIDE ao_bielec_integrals_erf_in_map + + integer(omp_lock_kind) :: lck(ao_num) + integer*8 :: i8 + integer :: ii(8), jj(8), kk(8), ll(8), k2 + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + integer(cache_map_size_kind) :: n_elements_max_erf, n_elements_erf + integer(key_kind), allocatable :: keys_erf(:) + double precision, allocatable :: values_erf(:) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& + !$OMP ao_integrals_map, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) + + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), & + ao_bi_elec_integral_beta_tmp(ao_num,ao_num)) + ao_bi_elec_integral_alpha_tmp = 0.d0 + ao_bi_elec_integral_beta_tmp = 0.d0 + + !$OMP DO SCHEDULE(dynamic,64) + !DIR$ NOVECTOR + do i8=0_8,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + do k1=1,n_elements + call bielec_integrals_index_reverse(kk,ii,ll,jj,keys(k1)) + + do k2=1,8 + if (kk(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1) + ao_bi_elec_integral_alpha_tmp(i,j) += integral + ao_bi_elec_integral_beta_tmp (i,j) += integral + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_bi_elec_integral_alpha += ao_bi_elec_integral_alpha_tmp + !$OMP END CRITICAL + !$OMP CRITICAL + ao_bi_elec_integral_beta += ao_bi_elec_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp) + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral_erf,ii,jj,kk,ll,i8,keys_erf,values_erf,n_elements_max_erf, & + !$OMP n_elements_erf,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& + !$OMP ao_integrals_erf_map, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) + + + call get_cache_map_n_elements_max(ao_integrals_erf_map,n_elements_max_erf) + allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), & + ao_bi_elec_integral_beta_tmp(ao_num,ao_num)) + allocate(keys_Erf(n_elements_max_erf), values_erf(n_elements_max_erf)) + + ao_bi_elec_integral_alpha_tmp = 0.d0 + ao_bi_elec_integral_beta_tmp = 0.d0 + !$OMP DO SCHEDULE(dynamic,64) + !DIR$ NOVECTOR + do i8=0_8,ao_integrals_erf_map%map_size + n_elements_erf = n_elements_max_erf + call get_cache_map(ao_integrals_erf_map,i8,keys_erf,values_erf,n_elements_erf) + do k1=1,n_elements_erf + call bielec_integrals_index_reverse(kk,ii,ll,jj,keys_erf(k1)) + + do k2=1,8 + if (kk(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + double precision :: integral_erf + integral_erf = values_erf(k1) + ao_bi_elec_integral_alpha_tmp(l,j) -= (SCF_density_matrix_ao_alpha(k,i) * integral_erf) + ao_bi_elec_integral_beta_tmp (l,j) -= (SCF_density_matrix_ao_beta (k,i) * integral_erf) + enddo + enddo + enddo + + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_bi_elec_integral_alpha = ao_bi_elec_integral_alpha + ao_bi_elec_integral_alpha_tmp + !$OMP END CRITICAL + !$OMP CRITICAL + ao_bi_elec_integral_beta = ao_bi_elec_integral_beta + ao_bi_elec_integral_beta_tmp + !$OMP END CRITICAL + deallocate(ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp) + deallocate(keys_erf,values_erf) + !$OMP END PARALLEL + + endif + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Alpha Fock matrix in AO basis set + END_DOC + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao_alpha(i,j) = Fock_matrix_alpha_no_xc_ao(i,j) + ao_potential_alpha_xc(i,j) + Fock_matrix_ao_beta (i,j) = Fock_matrix_beta_no_xc_ao(i,j) + ao_potential_beta_xc(i,j) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_no_xc_ao, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_no_xc_ao, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Mono electronic an Coulomb matrix in AO basis set + END_DOC + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_alpha_no_xc_ao(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_alpha(i,j) + Fock_matrix_beta_no_xc_ao(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_beta (i,j) + enddo + enddo + +END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, RS_KS_energy ] +!BEGIN_PROVIDER [ double precision, SCF_energy ] +&BEGIN_PROVIDER [ double precision, two_electron_energy] +&BEGIN_PROVIDER [ double precision, one_electron_energy] +&BEGIN_PROVIDER [ double precision, Fock_matrix_energy] +&BEGIN_PROVIDER [ double precision, trace_potential_xc ] + implicit none + BEGIN_DOC + ! Range-separated Kohn-Sham energy + END_DOC + RS_KS_energy = nuclear_repulsion + + integer :: i,j + double precision :: accu_mono,accu_fock + one_electron_energy = 0.d0 + two_electron_energy = 0.d0 + Fock_matrix_energy = 0.d0 + trace_potential_xc = 0.d0 + do j=1,ao_num + do i=1,ao_num + Fock_matrix_energy += Fock_matrix_ao_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) + & + Fock_matrix_ao_beta(i,j) * SCF_density_matrix_ao_beta(i,j) + two_electron_energy += 0.5d0 * ( ao_bi_elec_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) & + +ao_bi_elec_integral_beta(i,j) * SCF_density_matrix_ao_beta(i,j) ) + one_electron_energy += ao_mono_elec_integral(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) +! possible bug fix for open-shell +! trace_potential_xc += (ao_potential_alpha_xc(i,j) + ao_potential_beta_xc(i,j) ) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + trace_potential_xc += ao_potential_alpha_xc(i,j) * SCF_density_matrix_ao_alpha(i,j) + ao_potential_beta_xc(i,j) * SCF_density_matrix_ao_beta (i,j) + enddo + enddo + RS_KS_energy += e_exchange_dft + e_correlation_dft + one_electron_energy + two_electron_energy +!SCF_energy = RS_KS_energy +END_PROVIDER + +BEGIN_PROVIDER [double precision, extra_energy_contrib_from_density] + implicit none +! possible bug fix for open-shell: +! extra_energy_contrib_from_density = e_exchange_dft + e_correlation_dft - 0.25d0 * trace_potential_xc + extra_energy_contrib_from_density = e_exchange_dft + e_correlation_dft - 0.5d0 * trace_potential_xc +END_PROVIDER + +!BEGIN_PROVIDER [ double precision, SCF_energy ] +! implicit none +! SCF_energy = RS_KS_energy +!END_PROVIDER diff --git a/src/kohn_sham_range_separated/potential_functional.irp.f b/src/kohn_sham_range_separated/potential_functional.irp.f new file mode 100644 index 00000000..c0ef5d74 --- /dev/null +++ b/src/kohn_sham_range_separated/potential_functional.irp.f @@ -0,0 +1,25 @@ + BEGIN_PROVIDER [double precision, ao_potential_alpha_xc, (ao_num, ao_num)] +&BEGIN_PROVIDER [double precision, ao_potential_beta_xc, (ao_num, ao_num)] + implicit none + integer :: i,j,k,l + ao_potential_alpha_xc = 0.d0 + ao_potential_beta_xc = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + ao_potential_alpha_xc(i,j) = potential_c_alpha_ao(i,j,1) + potential_x_alpha_ao(i,j,1) + ao_potential_beta_xc(i,j) = potential_c_beta_ao(i,j,1) + potential_x_beta_ao(i,j,1) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, e_exchange_dft] + implicit none + e_exchange_dft = energy_x(1) + +END_PROVIDER + +BEGIN_PROVIDER [double precision, e_correlation_dft] + implicit none + e_correlation_dft = energy_c(1) + +END_PROVIDER diff --git a/src/kohn_sham_range_separated/rs_ks_scf.irp.f b/src/kohn_sham_range_separated/rs_ks_scf.irp.f new file mode 100644 index 00000000..afa5d9bc --- /dev/null +++ b/src/kohn_sham_range_separated/rs_ks_scf.irp.f @@ -0,0 +1,103 @@ +program srs_ks_cf + BEGIN_DOC +! Produce `Range_separated_Kohn_Sham` MO orbital +! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ +! output: kohn_sham.energy +! optional: mo_basis.mo_coef + END_DOC + read_wf = .False. + density_for_dft ="WFT" + touch density_for_dft + touch read_wf + print*, '**************************' + print*, 'mu_erf_dft = ',mu_erf_dft + print*, '**************************' + call check_coherence_functional + call create_guess + call orthonormalize_mos + call run +end + +subroutine check_coherence_functional + implicit none + integer :: ifound_x,ifound_c + if(exchange_functional.eq."None")then + ifound_x = 1 + else + ifound_x = index(exchange_functional,"short_range") + endif + + if(correlation_functional.eq."None")then + ifound_c = 1 + else + ifound_c = index(correlation_functional,"short_range") + endif + print*,ifound_x,ifound_c + if(ifound_x .eq.0 .or. ifound_c .eq. 0)then + print*,'YOU ARE USING THE RANGE SEPARATED KS PROGRAM BUT YOUR INPUT KEYWORD FOR ' + print*,'exchange_functional is ',exchange_functional + print*,'correlation_functional is ',correlation_functional + print*,'CHANGE THE exchange_functional and correlation_functional keywords to range separated functionals' + print*,'or switch to the KS_SCF program that uses regular functionals' + stop + endif + +end + + +subroutine create_guess + implicit none + BEGIN_DOC +! Create a MO guess if no MOs are present in the EZFIO directory + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_mo_basis_mo_coef(exists) + if (.not.exists) then + print*,'Creating a guess for the MOs' + print*,'mo_guess_type = ',mo_guess_type + if (mo_guess_type == "HCore") then + mo_coef = ao_ortho_lowdin_coef + TOUCH mo_coef + mo_label = 'Guess' + call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.) + SOFT_TOUCH mo_coef mo_label + else if (mo_guess_type == "Huckel") then + call huckel_guess + else + print *, 'Unrecognized MO guess type : '//mo_guess_type + stop 1 + endif + endif +end + +subroutine run + + BEGIN_DOC +! Run SCF calculation + END_DOC + + use bitmasks + implicit none + + double precision :: EHF + + EHF = RS_KS_energy + + mo_label = "Canonical" + +! Choose SCF algorithm + +! call damping_SCF ! Deprecated routine + call Roothaan_Hall_SCF + + write(*, '(A22,X,F16.10)') 'one_electron_energy = ',one_electron_energy + write(*, '(A22,X,F16.10)') 'two_electron_energy = ',two_electron_energy + write(*, '(A22,X,F16.10)') 'e_exchange_dft = ',e_exchange_dft + write(*, '(A22,X,F16.10)') 'e_correlation_dft = ',e_correlation_dft + write(*, '(A22,X,F16.10)') 'Fock_matrix_energy = ',Fock_matrix_energy + + +end + + diff --git a/src/mo_basis/track_orb.irp.f b/src/mo_basis/track_orb.irp.f new file mode 100644 index 00000000..6eb50e11 --- /dev/null +++ b/src/mo_basis/track_orb.irp.f @@ -0,0 +1,64 @@ +BEGIN_PROVIDER [ double precision, mo_coef_begin_iteration, (ao_num,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Void provider to store the coefficients of the |MO| basis at the beginning of the SCF iteration + ! + ! Usefull to track some orbitals + END_DOC +END_PROVIDER + +subroutine initialize_mo_coef_begin_iteration + implicit none + BEGIN_DOC + ! + ! Initialize :c:data:`mo_coef_begin_iteration` to the current :c:data:`mo_coef` + END_DOC + mo_coef_begin_iteration = mo_coef +end + +subroutine reorder_active_orb + implicit none + BEGIN_DOC +! routines that takes the current :c:data:`mo_coef` and reorder the active orbitals (see :c:data:`list_act` and :c:data:`n_act_orb`) according to the overlap with :c:data:`mo_coef_begin_iteration` + END_DOC + integer :: i,j,iorb + integer :: k,l + double precision, allocatable :: accu(:) + integer, allocatable :: index_active_orb(:),iorder(:) + double precision, allocatable :: mo_coef_tmp(:,:) + allocate(accu(mo_tot_num),index_active_orb(n_act_orb),iorder(mo_tot_num)) + allocate(mo_coef_tmp(ao_num,mo_tot_num)) + + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, mo_tot_num + accu(j) = 0.d0 + iorder(j) = j + do k = 1, ao_num + do l = 1, ao_num + accu(j) += mo_coef_begin_iteration(k,iorb) * mo_coef(l,j) * ao_overlap(k,l) + enddo + enddo + accu(j) = -dabs(accu(j)) + enddo + call dsort(accu,iorder,mo_tot_num) + index_active_orb(i) = iorder(1) + enddo + + double precision :: x + integer :: i1,i2 + print*, 'swapping the active MOs' + do j = 1, n_act_orb + i1 = list_act(j) + i2 = index_active_orb(j) + print*, i1,i2 + do i=1,ao_num + x = mo_coef(i,i1) + mo_coef(i,i1) = mo_coef(i,i2) + mo_coef(i,i2) = x + enddo + enddo +!call loc_cele_routine + + deallocate(accu,index_active_orb, iorder) +end diff --git a/src/mo_one_e_integrals/read_write.irp.f b/src/mo_one_e_integrals/read_write.irp.f index 302ddc77..d9488357 100644 --- a/src/mo_one_e_integrals/read_write.irp.f +++ b/src/mo_one_e_integrals/read_write.irp.f @@ -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 diff --git a/src/mo_two_e_integrals/map_integrals.irp.f b/src/mo_two_e_integrals/map_integrals.irp.f index 2fdb03b9..bb5b2a85 100644 --- a/src/mo_two_e_integrals/map_integrals.irp.f +++ b/src/mo_two_e_integrals/map_integrals.irp.f @@ -239,6 +239,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 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 diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index a343374a..f97c9415 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -135,7 +135,6 @@ END_PROVIDER ASSERT (nucl_dist(ie1,ie2) > 0.d0) enddo enddo - END_PROVIDER BEGIN_PROVIDER [ double precision, nuclear_repulsion ] diff --git a/src/scf_utils/EZFIO.cfg b/src/scf_utils/EZFIO.cfg new file mode 100644 index 00000000..af42dddc --- /dev/null +++ b/src/scf_utils/EZFIO.cfg @@ -0,0 +1,53 @@ +[max_dim_diis] +type: integer +doc: Maximum size of the DIIS extrapolation procedure +interface: ezfio,provider,ocaml +default: 15 + +[threshold_diis] +type: Threshold +doc: Threshold on the convergence of the DIIS error vector during a Hartree-Fock calculation. If 0. is chosen, the square root of thresh_scf will be used. +interface: ezfio,provider,ocaml +default: 0. + +[thresh_scf] +type: Threshold +doc: Threshold on the convergence of the Hartree Fock energy. +interface: ezfio,provider,ocaml +default: 1.e-10 + +[n_it_scf_max] +type: Strictly_positive_int +doc: Maximum number of SCF iterations +interface: ezfio,provider,ocaml +default: 500 + +[level_shift] +type: Positive_float +doc: Energy shift on the virtual MOs to improve SCF convergence +interface: ezfio,provider,ocaml +default: 0.1 + +[scf_algorithm] +type: character*(32) +doc: Type of SCF algorithm used. Possible choices are [ Simple | DIIS] +interface: ezfio,provider,ocaml +default: DIIS + +[mo_guess_type] +type: MO_guess +doc: Initial MO guess. Can be [ Huckel | HCore ] +interface: ezfio,provider,ocaml +default: Huckel + +[energy] +type: double precision +doc: Calculated HF energy +interface: ezfio + +[no_oa_or_av_opt] +type: logical +doc: If true, leave the active orbitals untouched in the SCF procedure +interface: ezfio,provider,ocaml +default: False + diff --git a/src/scf_utils/NEED b/src/scf_utils/NEED new file mode 100644 index 00000000..b89695da --- /dev/null +++ b/src/scf_utils/NEED @@ -0,0 +1,2 @@ +mo_guess +bitmask diff --git a/src/scf_utils/README.rst b/src/scf_utils/README.rst new file mode 100644 index 00000000..a672b7f6 --- /dev/null +++ b/src/scf_utils/README.rst @@ -0,0 +1,6 @@ +========= +scf_utils +========= + + +TODO diff --git a/src/scf_utils/damping_scf.irp.f b/src/scf_utils/damping_scf.irp.f new file mode 100644 index 00000000..6706fb0b --- /dev/null +++ b/src/scf_utils/damping_scf.irp.f @@ -0,0 +1,146 @@ +subroutine damping_SCF + implicit none + double precision :: E + double precision, allocatable :: D_alpha(:,:), D_beta(:,:) + double precision :: E_new + double precision, allocatable :: D_new_alpha(:,:), D_new_beta(:,:), F_new(:,:) + double precision, allocatable :: delta_alpha(:,:), delta_beta(:,:) + double precision :: lambda, E_half, a, b, delta_D, delta_E, E_min + + integer :: i,j,k + logical :: saving + character :: save_char + + allocate( & + D_alpha( ao_num, ao_num ), & + D_beta( ao_num, ao_num ), & + F_new( ao_num, ao_num ), & + D_new_alpha( ao_num, ao_num ), & + D_new_beta( ao_num, ao_num ), & + delta_alpha( ao_num, ao_num ), & + delta_beta( ao_num, ao_num )) + + do j=1,ao_num + do i=1,ao_num + D_alpha(i,j) = SCF_density_matrix_ao_alpha(i,j) + D_beta (i,j) = SCF_density_matrix_ao_beta (i,j) + enddo + enddo + + + call write_time(6) + + write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') & + '====','================','================','================', '====' + write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') & + ' N ', 'Energy ', 'Energy diff ', 'Density diff ', 'Save' + write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') & + '====','================','================','================', '====' + + E = SCF_energy + 1.d0 + E_min = SCF_energy + delta_D = 0.d0 + do k=1,n_it_scf_max + + delta_E = SCF_energy - E + E = SCF_energy + + if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then + exit + endif + + saving = E < E_min + if (saving) then + call save_mos + save_char = 'X' + E_min = E + else + save_char = ' ' + endif + + write(6,'(I4,1X,F16.10, 1X, F16.10, 1X, F16.10, 3X, A )') & + k, E, delta_E, delta_D, save_char + + if(no_oa_or_av_opt)then + call initialize_mo_coef_begin_iteration + endif + D_alpha = SCF_density_matrix_ao_alpha + D_beta = SCF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + if(no_oa_or_av_opt)then + call reorder_active_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH mo_coef + + D_new_alpha = SCF_density_matrix_ao_alpha + D_new_beta = SCF_density_matrix_ao_beta + F_new = Fock_matrix_ao + E_new = SCF_energy + + delta_alpha = D_new_alpha - D_alpha + delta_beta = D_new_beta - D_beta + + lambda = .5d0 + E_half = 0.d0 + do while (E_half > E) + SCF_density_matrix_ao_alpha = D_alpha + lambda * delta_alpha + SCF_density_matrix_ao_beta = D_beta + lambda * delta_beta + TOUCH SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + if(no_oa_or_av_opt)then + call reorder_active_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH mo_coef + E_half = SCF_energy + if ((E_half > E).and.(E_new < E)) then + lambda = 1.d0 + exit + else if ((E_half > E).and.(lambda > 5.d-4)) then + lambda = 0.5d0 * lambda + E_new = E_half + else + exit + endif + enddo + + a = (E_new + E - 2.d0*E_half)*2.d0 + b = -E_new - 3.d0*E + 4.d0*E_half + lambda = -lambda*b/(a+1.d-16) + D_alpha = (1.d0-lambda) * D_alpha + lambda * D_new_alpha + D_beta = (1.d0-lambda) * D_beta + lambda * D_new_beta + delta_E = SCF_energy - E + do j=1,ao_num + do i=1,ao_num + delta_D = delta_D + & + (D_alpha(i,j) - SCF_density_matrix_ao_alpha(i,j))*(D_alpha(i,j) - SCF_density_matrix_ao_alpha(i,j)) + & + (D_beta (i,j) - SCF_density_matrix_ao_beta (i,j))*(D_beta (i,j) - SCF_density_matrix_ao_beta (i,j)) + enddo + enddo + delta_D = dsqrt(delta_D/dble(ao_num)**2) + SCF_density_matrix_ao_alpha = D_alpha + SCF_density_matrix_ao_beta = D_beta + TOUCH SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + if(no_oa_or_av_opt)then + call reorder_active_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH mo_coef + + enddo + write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') '====','================','================','================', '====' + write(6,*) + + if(.not.no_oa_or_av_opt)then + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.) + endif + + call write_double(6, E_min, 'Hartree-Fock energy') + call ezfio_set_hartree_fock_energy(E_min) + + call write_time(6) + + deallocate(D_alpha,D_beta,F_new,D_new_alpha,D_new_beta,delta_alpha,delta_beta) +end diff --git a/src/hartree_fock/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f similarity index 100% rename from src/hartree_fock/diagonalize_fock.irp.f rename to src/scf_utils/diagonalize_fock.irp.f diff --git a/src/hartree_fock/diis.irp.f b/src/scf_utils/diis.irp.f similarity index 97% rename from src/hartree_fock/diis.irp.f rename to src/scf_utils/diis.irp.f index 6ebb5879..dfc43015 100644 --- a/src/hartree_fock/diis.irp.f +++ b/src/scf_utils/diis.irp.f @@ -27,7 +27,7 @@ BEGIN_PROVIDER [double precision, FPS_SPF_Matrix_AO, (AO_num, AO_num)] call dgemm('N','N',AO_num,AO_num,AO_num, & 1.d0, & Fock_Matrix_AO,Size(Fock_Matrix_AO,1), & - HF_Density_Matrix_AO,Size(HF_Density_Matrix_AO,1), & + SCF_Density_Matrix_AO,Size(SCF_Density_Matrix_AO,1), & 0.d0, & scratch,Size(scratch,1)) @@ -45,7 +45,7 @@ BEGIN_PROVIDER [double precision, FPS_SPF_Matrix_AO, (AO_num, AO_num)] call dgemm('N','N',AO_num,AO_num,AO_num, & 1.d0, & AO_Overlap,Size(AO_Overlap,1), & - HF_Density_Matrix_AO,Size(HF_Density_Matrix_AO,1), & + SCF_Density_Matrix_AO,Size(SCF_Density_Matrix_AO,1), & 0.d0, & scratch,Size(scratch,1)) diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f new file mode 100644 index 00000000..63141d04 --- /dev/null +++ b/src/scf_utils/fock_matrix.irp.f @@ -0,0 +1,144 @@ + BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo, (mo_tot_num)] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis. + ! For open shells, the ROHF Fock Matrix is + ! + ! | F-K | F + K/2 | F | + ! |---------------------------------| + ! | F + K/2 | F | F - K/2 | + ! |---------------------------------| + ! | F | F - K/2 | F + K | + ! + ! F = 1/2 (Fa + Fb) + ! + ! K = Fb - Fa + ! + END_DOC + integer :: i,j,n + if (elec_alpha_num == elec_beta_num) then + Fock_matrix_mo = Fock_matrix_mo_alpha + else + + do j=1,elec_beta_num + ! F-K + do i=1,elec_beta_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + enddo + ! F+K/2 + do i=elec_beta_num+1,elec_alpha_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + enddo + ! F + do i=elec_alpha_num+1, mo_tot_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) + enddo + enddo + + do j=elec_beta_num+1,elec_alpha_num + ! F+K/2 + do i=1,elec_beta_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + enddo + ! F + do i=elec_beta_num+1,elec_alpha_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) + enddo + ! F-K/2 + do i=elec_alpha_num+1, mo_tot_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + enddo + enddo + + do j=elec_alpha_num+1, mo_tot_num + ! F + do i=1,elec_beta_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) + enddo + ! F-K/2 + do i=elec_beta_num+1,elec_alpha_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + enddo + ! F+K + do i=elec_alpha_num+1,mo_tot_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) & + + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + enddo + enddo + + endif + + do i = 1, mo_tot_num + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + enddo +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & + Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1)) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & + Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1)) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao(i,j) = Fock_matrix_ao_alpha(i,j) + enddo + enddo + else + call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & + Fock_matrix_ao,size(Fock_matrix_ao,1)) + endif +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, SCF_energy ] + implicit none + BEGIN_DOC + ! Hartree-Fock energy + END_DOC + SCF_energy = nuclear_repulsion + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + SCF_energy += 0.5d0 * ( & + (ao_mono_elec_integral(i,j) + Fock_matrix_ao_alpha(i,j) ) * SCF_density_matrix_ao_alpha(i,j) +& + (ao_mono_elec_integral(i,j) + Fock_matrix_ao_beta (i,j) ) * SCF_density_matrix_ao_beta (i,j) ) + enddo + enddo + SCF_energy += extra_energy_contrib_from_density + +END_PROVIDER + diff --git a/src/hartree_fock/huckel.irp.f b/src/scf_utils/huckel.irp.f similarity index 100% rename from src/hartree_fock/huckel.irp.f rename to src/scf_utils/huckel.irp.f diff --git a/src/hartree_fock/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f similarity index 96% rename from src/hartree_fock/roothaan_hall_scf.irp.f rename to src/scf_utils/roothaan_hall_scf.irp.f index e860496f..76c0b1fe 100644 --- a/src/hartree_fock/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -24,6 +24,7 @@ END_DOC call write_time(6) + print*,'Energy of the guess = ',SCF_energy write(6,'(A4, 1X, A16, 1X, A16, 1X, A16)') & '====','================','================','================' write(6,'(A4, 1X, A16, 1X, A16, 1X, A16)') & @@ -32,8 +33,7 @@ END_DOC '====','================','================','================' ! Initialize energies and density matrices - - energy_SCF_previous = HF_energy + energy_SCF_previous = SCF_energy Delta_energy_SCF = 1.d0 iteration_SCF = 0 dim_DIIS = 0 @@ -88,7 +88,7 @@ END_DOC ! SCF energy - energy_SCF = HF_energy + energy_SCF = SCF_energy Delta_Energy_SCF = energy_SCF - energy_SCF_previous if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then Fock_matrix_AO(1:ao_num,1:ao_num) = Fock_matrix_DIIS (1:ao_num,1:ao_num,index_dim_DIIS) @@ -100,14 +100,14 @@ END_DOC double precision :: level_shift_save level_shift_save = level_shift mo_coef_save(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num) - do while (Delta_Energy_SCF > 0.d0) + do while (Delta_energy_SCF .ge. 0.d0) mo_coef(1:ao_num,1:mo_tot_num) = mo_coef_save TOUCH mo_coef - level_shift = level_shift + 1.0d0 + level_shift = level_shift + 0.1d0 mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_tot_num) TOUCH mo_coef level_shift - Delta_Energy_SCF = HF_energy - energy_SCF_previous - energy_SCF = HF_energy + Delta_Energy_SCF = SCF_energy - energy_SCF_previous + energy_SCF = SCF_energy if (level_shift-level_shift_save > 50.d0) then level_shift = level_shift_save SOFT_TOUCH level_shift @@ -140,6 +140,7 @@ END_DOC if(.not.no_oa_or_av_opt)then call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.) + call save_mos endif call write_double(6, Energy_SCF, 'Hartree-Fock energy') diff --git a/src/scf_utils/scf_density_matrix_ao.irp.f b/src/scf_utils/scf_density_matrix_ao.irp.f new file mode 100644 index 00000000..d085fb92 --- /dev/null +++ b/src/scf_utils/scf_density_matrix_ao.irp.f @@ -0,0 +1,41 @@ +BEGIN_PROVIDER [double precision, SCF_density_matrix_ao_alpha, (ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! S^{-1}.P_alpha.S^{-1} + END_DOC + + call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + SCF_density_matrix_ao_alpha, size(SCF_density_matrix_ao_alpha,1)) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, SCF_density_matrix_ao_beta, (ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! S^{-1}.P_beta.S^{-1} + END_DOC + + call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + SCF_density_matrix_ao_beta, size(SCF_density_matrix_ao_beta,1)) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, SCF_density_matrix_ao, (ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! S^{-1}.P.S^{-1} where P = C.C^t + END_DOC + ASSERT (size(SCF_density_matrix_ao,1) == size(SCF_density_matrix_ao_alpha,1)) + if (elec_alpha_num== elec_beta_num) then + SCF_density_matrix_ao = SCF_density_matrix_ao_alpha + SCF_density_matrix_ao_alpha + else + ASSERT (size(SCF_density_matrix_ao,1) == size(SCF_density_matrix_ao_beta ,1)) + SCF_density_matrix_ao = SCF_density_matrix_ao_alpha + SCF_density_matrix_ao_beta + endif + +END_PROVIDER + diff --git a/src/tools/diagonalize_h.irp.f b/src/tools/diagonalize_h.irp.f new file mode 100644 index 00000000..8f564ca9 --- /dev/null +++ b/src/tools/diagonalize_h.irp.f @@ -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