mirror of
https://github.com/LCPQ/quantum_package
synced 2024-10-19 22:41:48 +02:00
added kohn_sham
This commit is contained in:
Emmanuel Giner
parent
e2c380a1f3
commit
045ebc9441
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
86
src/dft_utils_one_body/sr_e_c_e_x_specific_prov.irp.f
Normal file
86
src/dft_utils_one_body/sr_e_c_e_x_specific_prov.irp.f
Normal file
@ -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
|
||||
|
||||
188
src/dft_utils_one_body/sr_potentials_ao_specific_prov.irp.f
Normal file
188
src/dft_utils_one_body/sr_potentials_ao_specific_prov.irp.f
Normal file
@ -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
|
||||
@ -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))
|
||||
|
||||
2
src/kohn_sham/NEED
Normal file
2
src/kohn_sham/NEED
Normal file
@ -0,0 +1,2 @@
|
||||
dft_utils_one_body
|
||||
scf_utils
|
||||
244
src/kohn_sham/fock_matrix_ks.irp.f
Normal file
244
src/kohn_sham/fock_matrix_ks.irp.f
Normal file
@ -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
|
||||
|
||||
95
src/kohn_sham/ks_scf.irp.f
Normal file
95
src/kohn_sham/ks_scf.irp.f
Normal file
@ -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
|
||||
|
||||
|
||||
25
src/kohn_sham/potential_functional.irp.f
Normal file
25
src/kohn_sham/potential_functional.irp.f
Normal file
@ -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
|
||||
@ -1 +1,2 @@
|
||||
DFT_Utils_one_body Integrals_erf SCF_Utils Integrals_Bielec
|
||||
dft_utils_one_body
|
||||
scf_utils
|
||||
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user