diff --git a/src/dft_utils_one_e/effective_pot.irp.f b/src/dft_utils_one_e/effective_pot.irp.f index 5a61de15..27f4841e 100644 --- a/src/dft_utils_one_e/effective_pot.irp.f +++ b/src/dft_utils_one_e/effective_pot.irp.f @@ -10,10 +10,9 @@ ! ! Taking the expectation value does not provide any energy, but ! -! effective_one_e_potential(i,j) is the potential coupling DFT and WFT part to -! -! be used in any WFT calculation. +! effective_one_e_potential(i,j) is the potential coupling DFT and WFT parts ! +! and it is used in any RS-DFT based calculations END_DOC do istate = 1, N_states do j = 1, mo_num diff --git a/src/dft_utils_one_e/mu_erf_dft.irp.f b/src/dft_utils_one_e/mu_erf_dft.irp.f index 3a3a2f28..53effcb6 100644 --- a/src/dft_utils_one_e/mu_erf_dft.irp.f +++ b/src/dft_utils_one_e/mu_erf_dft.irp.f @@ -1,7 +1,9 @@ BEGIN_PROVIDER [double precision, mu_erf_dft] implicit none BEGIN_DOC -! range separation parameter used in RS-DFT. It is set to mu_erf in order to be consistent with the two electrons integrals erf +! range separation parameter used in RS-DFT. +! +! It is set to mu_erf in order to be consistent with the module "ao_two_e_erf_ints" END_DOC mu_erf_dft = mu_erf diff --git a/src/functionals/pbe.irp.f b/src/functionals/pbe.irp.f index 48b0661d..df32cce2 100644 --- a/src/functionals/pbe.irp.f +++ b/src/functionals/pbe.irp.f @@ -1,124 +1,79 @@ - BEGIN_PROVIDER[double precision, energy_x_pbe, (N_states) ] +&BEGIN_PROVIDER[double precision, energy_c_pbe, (N_states) ] implicit none BEGIN_DOC + ! exchange / correlation energies with the short-range version Perdew-Burke-Ernzerhof GGA functional + ! + ! defined in Chem. Phys.329, 276 (2006) + END_DOC + BEGIN_DOC ! exchange/correlation energy with the short range pbe functional END_DOC integer :: istate,i,j,m 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)) + double precision :: ex, ec + double precision :: rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),grad_rho_a_2,grad_rho_b_2,grad_rho_a_b + double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b + double precision :: 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(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 - do istate = 1, N_states - do i = 1, n_points_final_grid - weight = final_weight_at_r_vector(i) - rho_a(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) - rho_b(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - grad_rho_a(1:3,istate) = one_e_dm_and_grad_alpha_in_r(1:3,i,istate) - grad_rho_b(1:3,istate) = one_e_dm_and_grad_beta_in_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(0.d0,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 - enddo - enddo - - -END_PROVIDER - -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 :: 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_c_pbe = 0.d0 + mu = 0.d0 do istate = 1, N_states do i = 1, n_points_final_grid weight = final_weight_at_r_vector(i) - rho_a(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) - rho_b(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - grad_rho_a(1:3,istate) = one_e_dm_and_grad_alpha_in_r(1:3,i,istate) - grad_rho_b(1:3,istate) = one_e_dm_and_grad_beta_in_r(1:3,i,istate) + rho_a = one_e_dm_and_grad_alpha_in_r(4,i,istate) + rho_b = one_e_dm_and_grad_beta_in_r(4,i,istate) + grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,i,istate) + grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_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) + grad_rho_a_2 += grad_rho_a(m) * grad_rho_a(m) + grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m) + grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo ! inputs - call GGA_sr_type_functionals(0.d0,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu,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_c_pbe += ec * weight + energy_x_pbe(istate) += ex * weight + energy_c_pbe(istate) += 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 potential for alpha / beta electrons with the Perdew-Burke-Ernzerhof GGA functional + ! exchange / correlation potential for alpha / beta electrons with the short-range version Perdew-Burke-Ernzerhof GGA functional + ! + ! defined in Chem. Phys.329, 276 (2006) END_DOC integer :: i,j,istate do istate = 1, n_states do i = 1, ao_num do j = 1, ao_num - potential_x_alpha_ao_pbe(j,i,istate) = pot_scal_x_alpha_ao_pbe(j,i,istate) + pot_grad_x_alpha_ao_pbe(j,i,istate) + pot_grad_x_alpha_ao_pbe(i,j,istate) - potential_x_beta_ao_pbe(j,i,istate) = pot_scal_x_beta_ao_pbe(j,i,istate) + pot_grad_x_beta_ao_pbe(j,i,istate) + pot_grad_x_beta_ao_pbe(i,j,istate) + potential_x_alpha_ao_pbe(j,i,istate) = pot_sr_scal_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(i,j,istate) + potential_x_beta_ao_pbe(j,i,istate) = pot_sr_scal_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(i,j,istate) - potential_c_alpha_ao_pbe(j,i,istate) = pot_scal_c_alpha_ao_pbe(j,i,istate) + pot_grad_c_alpha_ao_pbe(j,i,istate) + pot_grad_c_alpha_ao_pbe(i,j,istate) - potential_c_beta_ao_pbe(j,i,istate) = pot_scal_c_beta_ao_pbe(j,i,istate) + pot_grad_c_beta_ao_pbe(j,i,istate) + pot_grad_c_beta_ao_pbe(i,j,istate) + potential_c_alpha_ao_pbe(j,i,istate) = pot_sr_scal_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(i,j,istate) + potential_c_beta_ao_pbe(j,i,istate) = pot_sr_scal_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(i,j,istate) enddo enddo enddo END_PROVIDER - - BEGIN_PROVIDER [double precision, potential_xc_alpha_ao_pbe,(ao_num,ao_num,N_states)] &BEGIN_PROVIDER [double precision, potential_xc_beta_ao_pbe,(ao_num,ao_num,N_states)] implicit none @@ -129,8 +84,8 @@ END_PROVIDER do istate = 1, n_states do i = 1, ao_num do j = 1, ao_num - potential_xc_alpha_ao_pbe(j,i,istate) = pot_scal_xc_alpha_ao_pbe(j,i,istate) + pot_grad_xc_alpha_ao_pbe(j,i,istate) + pot_grad_xc_alpha_ao_pbe(i,j,istate) - potential_xc_beta_ao_pbe(j,i,istate) = pot_scal_xc_beta_ao_pbe(j,i,istate) + pot_grad_xc_beta_ao_pbe(j,i,istate) + pot_grad_xc_beta_ao_pbe(i,j,istate) + potential_xc_alpha_ao_pbe(j,i,istate) = pot_sr_scal_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(i,j,istate) + potential_xc_beta_ao_pbe(j,i,istate) = pot_sr_scal_xc_beta_ao_pbe(j,i,istate) + pot_sr_grad_xc_beta_ao_pbe(j,i,istate) + pot_sr_grad_xc_beta_ao_pbe(i,j,istate) enddo enddo enddo @@ -138,78 +93,76 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER[double precision, aos_vc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_vc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_vx_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_vx_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvx_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvx_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] + + BEGIN_PROVIDER[double precision, aos_sr_vc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vx_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vx_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vx_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vx_beta_pbe_w , (ao_num,n_points_final_grid,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) +! intermediates to compute the sr_pbe potentials +! +! aos_sr_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 :: 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)) - - aos_dvc_alpha_pbe_w = 0.d0 - aos_dvc_beta_pbe_w = 0.d0 - aos_dvx_alpha_pbe_w = 0.d0 - aos_dvx_beta_pbe_w = 0.d0 - + double precision :: ex, ec + double precision :: rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),grad_rho_a_2,grad_rho_b_2,grad_rho_a_b + double precision :: contrib_grad_xa(3),contrib_grad_xb(3),contrib_grad_ca(3),contrib_grad_cb(3) + double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b + double precision :: 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 + aos_dsr_vc_alpha_pbe_w= 0.d0 + aos_dsr_vc_beta_pbe_w = 0.d0 + aos_dsr_vx_alpha_pbe_w= 0.d0 + aos_dsr_vx_beta_pbe_w = 0.d0 + mu = 0.d0 do istate = 1, N_states do i = 1, n_points_final_grid weight = final_weight_at_r_vector(i) - rho_a(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) - rho_b(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - grad_rho_a(1:3,istate) = one_e_dm_and_grad_alpha_in_r(1:3,i,istate) - grad_rho_b(1:3,istate) = one_e_dm_and_grad_beta_in_r(1:3,i,istate) + + rho_a = one_e_dm_and_grad_alpha_in_r(4,i,istate) + rho_b = one_e_dm_and_grad_beta_in_r(4,i,istate) + grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,i,istate) + grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_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) + grad_rho_a_2 += grad_rho_a(m) * grad_rho_a(m) + grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m) + grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo ! inputs - call GGA_sr_type_functionals(0.d0,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu,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 + vx_rho_a *= weight + vc_rho_a *= weight + vx_rho_b *= weight + vc_rho_b *= 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)) + contrib_grad_ca(m) = weight * (2.d0 * vc_grad_rho_a_2 * grad_rho_a(m) + vc_grad_rho_a_b * grad_rho_b(m) ) + contrib_grad_xa(m) = weight * (2.d0 * vx_grad_rho_a_2 * grad_rho_a(m) + vx_grad_rho_a_b * grad_rho_b(m) ) + contrib_grad_cb(m) = weight * (2.d0 * vc_grad_rho_b_2 * grad_rho_b(m) + vc_grad_rho_a_b * grad_rho_a(m) ) + contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 * grad_rho_b(m) + vx_grad_rho_a_b * grad_rho_a(m) ) 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) + aos_sr_vc_alpha_pbe_w(j,i,istate) = vc_rho_a * aos_in_r_array(j,i) + aos_sr_vc_beta_pbe_w (j,i,istate) = vc_rho_b * aos_in_r_array(j,i) + aos_sr_vx_alpha_pbe_w(j,i,istate) = vx_rho_a * aos_in_r_array(j,i) + aos_sr_vx_beta_pbe_w (j,i,istate) = vx_rho_b * aos_in_r_array(j,i) enddo do j = 1, ao_num do m = 1,3 - aos_dvc_alpha_pbe_w(j,i,istate) += contrib_grad_ca(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i) - aos_dvc_beta_pbe_w (j,i,istate) += contrib_grad_cb(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i) - aos_dvx_alpha_pbe_w(j,i,istate) += contrib_grad_xa(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i) - aos_dvx_beta_pbe_w (j,i,istate) += contrib_grad_xb(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i) + aos_dsr_vc_alpha_pbe_w(j,i,istate) += contrib_grad_ca(m) * aos_grad_in_r_array_transp_xyz(m,j,i) + aos_dsr_vc_beta_pbe_w (j,i,istate) += contrib_grad_cb(m) * aos_grad_in_r_array_transp_xyz(m,j,i) + aos_dsr_vx_alpha_pbe_w(j,i,istate) += contrib_grad_xa(m) * aos_grad_in_r_array_transp_xyz(m,j,i) + aos_dsr_vx_beta_pbe_w (j,i,istate) += contrib_grad_xb(m) * aos_grad_in_r_array_transp_xyz(m,j,i) enddo enddo enddo @@ -218,42 +171,44 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER [double precision, pot_scal_x_alpha_ao_pbe, (ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_scal_c_alpha_ao_pbe, (ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_scal_x_beta_ao_pbe, (ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_scal_c_beta_ao_pbe, (ao_num,ao_num,N_states)] + BEGIN_PROVIDER [double precision, pot_sr_scal_x_alpha_ao_pbe, (ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_scal_c_alpha_ao_pbe, (ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_scal_x_beta_ao_pbe, (ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_scal_c_beta_ao_pbe, (ao_num,ao_num,N_states)] implicit none +! intermediates to compute the sr_pbe potentials +! integer :: istate BEGIN_DOC ! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential END_DOC - pot_scal_c_alpha_ao_pbe = 0.d0 - pot_scal_x_alpha_ao_pbe = 0.d0 - pot_scal_c_beta_ao_pbe = 0.d0 - pot_scal_x_beta_ao_pbe = 0.d0 + pot_sr_scal_c_alpha_ao_pbe = 0.d0 + pot_sr_scal_x_alpha_ao_pbe = 0.d0 + pot_sr_scal_c_beta_ao_pbe = 0.d0 + pot_sr_scal_x_beta_ao_pbe = 0.d0 double precision :: wall_1,wall_2 call wall_time(wall_1) 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_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, & - pot_scal_c_alpha_ao_pbe(1,1,istate),size(pot_scal_c_alpha_ao_pbe,1)) + pot_sr_scal_c_alpha_ao_pbe(1,1,istate),size(pot_sr_scal_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_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, & - pot_scal_c_beta_ao_pbe(1,1,istate),size(pot_scal_c_beta_ao_pbe,1)) + pot_sr_scal_c_beta_ao_pbe(1,1,istate),size(pot_sr_scal_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_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, & - pot_scal_x_alpha_ao_pbe(1,1,istate),size(pot_scal_x_alpha_ao_pbe,1)) + pot_sr_scal_x_alpha_ao_pbe(1,1,istate),size(pot_sr_scal_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_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, & - pot_scal_x_beta_ao_pbe(1,1,istate), size(pot_scal_x_beta_ao_pbe,1)) + pot_sr_scal_x_beta_ao_pbe(1,1,istate), size(pot_sr_scal_x_beta_ao_pbe,1)) enddo call wall_time(wall_2) @@ -261,10 +216,10 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER [double precision, pot_grad_x_alpha_ao_pbe,(ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_grad_x_beta_ao_pbe,(ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_grad_c_alpha_ao_pbe,(ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_grad_c_beta_ao_pbe,(ao_num,ao_num,N_states)] + BEGIN_PROVIDER [double precision, pot_sr_grad_x_alpha_ao_pbe,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_grad_x_beta_ao_pbe,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_grad_c_alpha_ao_pbe,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_grad_c_beta_ao_pbe,(ao_num,ao_num,N_states)] implicit none BEGIN_DOC ! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the gradienst of the density and orbitals @@ -272,31 +227,31 @@ END_PROVIDER integer :: istate double precision :: wall_1,wall_2 call wall_time(wall_1) - pot_grad_c_alpha_ao_pbe = 0.d0 - pot_grad_x_alpha_ao_pbe = 0.d0 - pot_grad_c_beta_ao_pbe = 0.d0 - pot_grad_x_beta_ao_pbe = 0.d0 + pot_sr_grad_c_alpha_ao_pbe = 0.d0 + pot_sr_grad_x_alpha_ao_pbe = 0.d0 + pot_sr_grad_c_beta_ao_pbe = 0.d0 + pot_sr_grad_x_beta_ao_pbe = 0.d0 do istate = 1, N_states ! correlation alpha call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_dvc_alpha_pbe_w(1,1,istate),size(aos_dvc_alpha_pbe_w,1), & + aos_dsr_vc_alpha_pbe_w(1,1,istate),size(aos_dsr_vc_alpha_pbe_w,1), & aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, & - pot_grad_c_alpha_ao_pbe(1,1,istate),size(pot_grad_c_alpha_ao_pbe,1)) + pot_sr_grad_c_alpha_ao_pbe(1,1,istate),size(pot_sr_grad_c_alpha_ao_pbe,1)) ! correlation beta call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_dvc_beta_pbe_w(1,1,istate),size(aos_dvc_beta_pbe_w,1), & + aos_dsr_vc_beta_pbe_w(1,1,istate),size(aos_dsr_vc_beta_pbe_w,1), & aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, & - pot_grad_c_beta_ao_pbe(1,1,istate),size(pot_grad_c_beta_ao_pbe,1)) + pot_sr_grad_c_beta_ao_pbe(1,1,istate),size(pot_sr_grad_c_beta_ao_pbe,1)) ! exchange alpha call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_dvx_alpha_pbe_w(1,1,istate),size(aos_dvx_alpha_pbe_w,1), & + aos_dsr_vx_alpha_pbe_w(1,1,istate),size(aos_dsr_vx_alpha_pbe_w,1), & aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, & - pot_grad_x_alpha_ao_pbe(1,1,istate),size(pot_grad_x_alpha_ao_pbe,1)) + pot_sr_grad_x_alpha_ao_pbe(1,1,istate),size(pot_sr_grad_x_alpha_ao_pbe,1)) ! exchange beta call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_dvx_beta_pbe_w(1,1,istate),size(aos_dvx_beta_pbe_w,1), & + aos_dsr_vx_beta_pbe_w(1,1,istate),size(aos_dsr_vx_beta_pbe_w,1), & aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, & - pot_grad_x_beta_ao_pbe(1,1,istate),size(pot_grad_x_beta_ao_pbe,1)) + pot_sr_grad_x_beta_ao_pbe(1,1,istate),size(pot_sr_grad_x_beta_ao_pbe,1)) enddo call wall_time(wall_2) @@ -304,13 +259,13 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER[double precision, aos_vxc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_vxc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvxc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] -&BEGIN_PROVIDER[double precision, aos_dvxc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] + BEGIN_PROVIDER[double precision, aos_sr_vxc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_sr_vxc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vxc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)] +&BEGIN_PROVIDER[double precision, aos_dsr_vxc_beta_pbe_w , (ao_num,n_points_final_grid,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) +! aos_sr_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 :: mu,weight @@ -320,8 +275,9 @@ END_PROVIDER double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b double precision :: 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 - aos_dvxc_alpha_pbe_w = 0.d0 - aos_dvxc_beta_pbe_w = 0.d0 + mu = 0.d0 + aos_dsr_vxc_alpha_pbe_w = 0.d0 + aos_dsr_vxc_beta_pbe_w = 0.d0 do istate = 1, N_states do i = 1, n_points_final_grid @@ -339,28 +295,28 @@ END_PROVIDER grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo - ! call exc_sr_pbe - call GGA_sr_type_functionals(0.d0,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! inputs - ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs exchange - ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) ! outputs correlation + ! inputs + call GGA_sr_type_functionals(mu,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 *= weight vc_rho_a *= weight vx_rho_b *= weight vc_rho_b *= weight do m= 1,3 - contrib_grad_ca(m) = weight * (2.d0 * vc_grad_rho_a_2 * grad_rho_a(m) + vc_grad_rho_a_b * grad_rho_b(m)) - contrib_grad_xa(m) = weight * (2.d0 * vx_grad_rho_a_2 * grad_rho_a(m) + vx_grad_rho_a_b * grad_rho_b(m)) - contrib_grad_cb(m) = weight * (2.d0 * vc_grad_rho_b_2 * grad_rho_b(m) + vc_grad_rho_a_b * grad_rho_a(m)) - contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 * grad_rho_b(m) + vx_grad_rho_a_b * grad_rho_a(m)) + contrib_grad_ca(m) = weight * (2.d0 * vc_grad_rho_a_2 * grad_rho_a(m) + vc_grad_rho_a_b * grad_rho_b(m) ) + contrib_grad_xa(m) = weight * (2.d0 * vx_grad_rho_a_2 * grad_rho_a(m) + vx_grad_rho_a_b * grad_rho_b(m) ) + contrib_grad_cb(m) = weight * (2.d0 * vc_grad_rho_b_2 * grad_rho_b(m) + vc_grad_rho_a_b * grad_rho_a(m) ) + contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 * grad_rho_b(m) + vx_grad_rho_a_b * grad_rho_a(m) ) enddo do j = 1, ao_num - aos_vxc_alpha_pbe_w(j,i,istate) = ( vc_rho_a + vx_rho_a ) * aos_in_r_array(j,i) - aos_vxc_beta_pbe_w (j,i,istate) = ( vc_rho_b + vx_rho_b ) * aos_in_r_array(j,i) + aos_sr_vxc_alpha_pbe_w(j,i,istate) = ( vc_rho_a + vx_rho_a ) * aos_in_r_array(j,i) + aos_sr_vxc_beta_pbe_w (j,i,istate) = ( vc_rho_b + vx_rho_b ) * aos_in_r_array(j,i) enddo do j = 1, ao_num do m = 1,3 - aos_dvxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i) - aos_dvxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i) + aos_dsr_vxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i) + aos_dsr_vxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i) enddo enddo enddo @@ -369,36 +325,36 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER [double precision, pot_scal_xc_alpha_ao_pbe, (ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_scal_xc_beta_ao_pbe, (ao_num,ao_num,N_states)] + BEGIN_PROVIDER [double precision, pot_sr_scal_xc_alpha_ao_pbe, (ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_scal_xc_beta_ao_pbe, (ao_num,ao_num,N_states)] implicit none integer :: istate BEGIN_DOC ! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential END_DOC - pot_scal_xc_alpha_ao_pbe = 0.d0 - pot_scal_xc_beta_ao_pbe = 0.d0 + pot_sr_scal_xc_alpha_ao_pbe = 0.d0 + pot_sr_scal_xc_beta_ao_pbe = 0.d0 double precision :: wall_1,wall_2 call wall_time(wall_1) do istate = 1, N_states ! exchange - correlation alpha call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_vxc_alpha_pbe_w(1,1,istate),size(aos_vxc_alpha_pbe_w,1), & + aos_sr_vxc_alpha_pbe_w(1,1,istate),size(aos_sr_vxc_alpha_pbe_w,1), & aos_in_r_array,size(aos_in_r_array,1),1.d0, & - pot_scal_xc_alpha_ao_pbe(1,1,istate),size(pot_scal_xc_alpha_ao_pbe,1)) + pot_sr_scal_xc_alpha_ao_pbe(1,1,istate),size(pot_sr_scal_xc_alpha_ao_pbe,1)) ! exchange - correlation beta call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_vxc_beta_pbe_w(1,1,istate),size(aos_vxc_beta_pbe_w,1), & + aos_sr_vxc_beta_pbe_w(1,1,istate),size(aos_sr_vxc_beta_pbe_w,1), & aos_in_r_array,size(aos_in_r_array,1),1.d0, & - pot_scal_xc_beta_ao_pbe(1,1,istate),size(pot_scal_xc_beta_ao_pbe,1)) + pot_sr_scal_xc_beta_ao_pbe(1,1,istate),size(pot_sr_scal_xc_beta_ao_pbe,1)) enddo call wall_time(wall_2) END_PROVIDER - BEGIN_PROVIDER [double precision, pot_grad_xc_alpha_ao_pbe,(ao_num,ao_num,N_states)] -&BEGIN_PROVIDER [double precision, pot_grad_xc_beta_ao_pbe,(ao_num,ao_num,N_states)] + BEGIN_PROVIDER [double precision, pot_sr_grad_xc_alpha_ao_pbe,(ao_num,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, pot_sr_grad_xc_beta_ao_pbe,(ao_num,ao_num,N_states)] implicit none BEGIN_DOC ! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the gradienst of the density and orbitals @@ -406,21 +362,22 @@ END_PROVIDER integer :: istate double precision :: wall_1,wall_2 call wall_time(wall_1) - pot_grad_xc_alpha_ao_pbe = 0.d0 - pot_grad_xc_beta_ao_pbe = 0.d0 + pot_sr_grad_xc_alpha_ao_pbe = 0.d0 + pot_sr_grad_xc_beta_ao_pbe = 0.d0 do istate = 1, N_states - ! correlation alpha + ! exchange - correlation alpha call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_dvxc_alpha_pbe_w(1,1,istate),size(aos_dvxc_alpha_pbe_w,1), & + aos_dsr_vxc_alpha_pbe_w(1,1,istate),size(aos_dsr_vxc_alpha_pbe_w,1), & aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, & - pot_grad_xc_alpha_ao_pbe(1,1,istate),size(pot_grad_xc_alpha_ao_pbe,1)) - ! correlation beta + pot_sr_grad_xc_alpha_ao_pbe(1,1,istate),size(pot_sr_grad_xc_alpha_ao_pbe,1)) + ! exchange - correlation beta call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, & - aos_dvxc_beta_pbe_w(1,1,istate),size(aos_dvxc_beta_pbe_w,1), & + aos_dsr_vxc_beta_pbe_w(1,1,istate),size(aos_dsr_vxc_beta_pbe_w,1), & aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, & - pot_grad_xc_beta_ao_pbe(1,1,istate),size(pot_grad_xc_beta_ao_pbe,1)) + pot_sr_grad_xc_beta_ao_pbe(1,1,istate),size(pot_sr_grad_xc_beta_ao_pbe,1)) enddo call wall_time(wall_2) END_PROVIDER + diff --git a/src/mo_one_e_ints/EZFIO.cfg b/src/mo_one_e_ints/EZFIO.cfg index 1cb77b6e..0f31b16a 100644 --- a/src/mo_one_e_ints/EZFIO.cfg +++ b/src/mo_one_e_ints/EZFIO.cfg @@ -24,7 +24,6 @@ interface: ezfio,provider,ocaml default: None - [mo_integrals_pseudo] type: double precision doc: Pseudopotential integrals in |MO| basis set