From 045ebc94418b39c5f335100cf955ffef9117f641 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 20 Dec 2018 18:43:29 +0100 Subject: [PATCH] added kohn_sham --- .../e_c_e_x_and_pot_general_prov.irp.f | 32 +-- .../e_c_e_x_specific_prov.irp.f | 4 +- .../potentials_ao_specific_prov.irp.f | 188 -------------- .../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/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 | 3 +- src/kohn_sham_range_separated/rs_ks_scf.irp.f | 6 +- 12 files changed, 719 insertions(+), 213 deletions(-) delete mode 100644 src/dft_utils_one_body/potentials_ao_specific_prov.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/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 diff --git a/src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f b/src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f index 4c2ff682..e1a9439e 100644 --- a/src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f +++ b/src/dft_utils_one_body/e_c_e_x_and_pot_general_prov.irp.f @@ -9,11 +9,11 @@ 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 @@ -96,11 +96,11 @@ END_PROVIDER ! 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 @@ -111,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 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 index c766af54..a27eb167 100644 --- 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 @@ -22,8 +22,8 @@ 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) + 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 diff --git a/src/dft_utils_one_body/potentials_ao_specific_prov.irp.f b/src/dft_utils_one_body/potentials_ao_specific_prov.irp.f deleted file mode 100644 index bfd8b0b5..00000000 --- a/src/dft_utils_one_body/potentials_ao_specific_prov.irp.f +++ /dev/null @@ -1,188 +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_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) - 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, 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, 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/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/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 index 73d45583..643ecece 100644 --- a/src/kohn_sham_range_separated/NEED +++ b/src/kohn_sham_range_separated/NEED @@ -1 +1,2 @@ -DFT_Utils_one_body Integrals_erf SCF_Utils Integrals_Bielec +dft_utils_one_body +scf_utils diff --git a/src/kohn_sham_range_separated/rs_ks_scf.irp.f b/src/kohn_sham_range_separated/rs_ks_scf.irp.f index 10dbf6e5..afa5d9bc 100644 --- a/src/kohn_sham_range_separated/rs_ks_scf.irp.f +++ b/src/kohn_sham_range_separated/rs_ks_scf.irp.f @@ -1,6 +1,6 @@ -program scf +program srs_ks_cf BEGIN_DOC -! Produce `Kohn_Sham` MO orbital +! 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 @@ -10,7 +10,7 @@ program scf touch density_for_dft touch read_wf print*, '**************************' - print*, 'mu_erf = ',mu_erf + print*, 'mu_erf_dft = ',mu_erf_dft print*, '**************************' call check_coherence_functional call create_guess