mirror of
https://github.com/LCPQ/quantum_package
synced 2025-01-03 10:05:57 +01:00
DFT begins to work with lda
This commit is contained in:
Emmanuel Giner
parent
226ca23af8
commit
5b8175e818
6951
plugins/DFT_Utils/angular.f
Normal file
6951
plugins/DFT_Utils/angular.f
Normal file
File diff suppressed because it is too large
Load Diff
25
plugins/DFT_Utils/functional.irp.f
Normal file
25
plugins/DFT_Utils/functional.irp.f
Normal file
@ -0,0 +1,25 @@
|
||||
double precision function ex_lda(rho)
|
||||
include 'constants.include.F'
|
||||
implicit none
|
||||
double precision, intent(in) :: rho
|
||||
ex_lda = cst_lda * rho**(c_4_3)
|
||||
|
||||
end
|
||||
|
||||
BEGIN_PROVIDER [double precision, lda_exchange, (N_states)]
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
double precision :: ex_lda
|
||||
do l = 1, N_states
|
||||
lda_exchange(l) = 0.d0
|
||||
do j = 1, nucl_num
|
||||
do i = 1, n_points_radial_grid
|
||||
do k = 1, n_points_integration_angular
|
||||
lda_exchange(l) += final_weight_functions_at_grid_points(k,i,j) * &
|
||||
(ex_lda(one_body_dm_mo_alpha_at_grid_points(k,i,j,l)) + ex_lda(one_body_dm_mo_beta_at_grid_points(k,i,j,l)))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
@ -1,7 +1,7 @@
|
||||
!BEGIN_PROVIDER [integer, n_points_integration_angular_lebedev]
|
||||
!implicit none
|
||||
!n_points_integration_angular_lebedev = 50
|
||||
!END_PROVIDER
|
||||
BEGIN_PROVIDER [integer, n_points_integration_angular]
|
||||
implicit none
|
||||
n_points_integration_angular = 110
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, n_points_radial_grid]
|
||||
implicit none
|
||||
@ -9,22 +9,33 @@ BEGIN_PROVIDER [integer, n_points_radial_grid]
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular_lebedev,3) ]
|
||||
&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular_lebedev)]
|
||||
BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_integration_angular,3) ]
|
||||
&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_integration_angular)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! weights and grid points for the integration on the angular variables on
|
||||
! the unit sphere centered on (0,0,0)
|
||||
! According to the LEBEDEV scheme
|
||||
END_DOC
|
||||
call cal_quad(n_points_integration_angular_lebedev, angular_quadrature_points,weights_angular_points)
|
||||
angular_quadrature_points = 0.d0
|
||||
weights_angular_points = 0.d0
|
||||
!call cal_quad(n_points_integration_angular, angular_quadrature_points,weights_angular_points)
|
||||
include 'constants.include.F'
|
||||
integer :: i
|
||||
integer :: i,n
|
||||
double precision :: accu
|
||||
double precision :: degre_rad
|
||||
degre_rad = 180.d0/pi
|
||||
degre_rad = pi/180.d0
|
||||
accu = 0.d0
|
||||
!do i = 1, n_points_integration_angular_lebedev
|
||||
double precision :: x(n_points_integration_angular),y(n_points_integration_angular),z(n_points_integration_angular),w(n_points_integration_angular)
|
||||
call LD0110(X,Y,Z,W,N)
|
||||
do i = 1, n_points_integration_angular
|
||||
angular_quadrature_points(i,1) = x(i)
|
||||
angular_quadrature_points(i,2) = y(i)
|
||||
angular_quadrature_points(i,3) = z(i)
|
||||
weights_angular_points(i) = w(i) * 4.d0 * pi
|
||||
accu += w(i)
|
||||
enddo
|
||||
!do i = 1, n_points_integration_angular
|
||||
! accu += weights_angular_integration_lebedev(i)
|
||||
! weights_angular_points(i) = weights_angular_integration_lebedev(i) * 4.d0 * pi
|
||||
! angular_quadrature_points(i,1) = dcos ( degre_rad * theta_angular_integration_lebedev(i)) &
|
||||
@ -70,7 +81,7 @@ END_PROVIDER
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num)]
|
||||
BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular,n_points_radial_grid,nucl_num)]
|
||||
BEGIN_DOC
|
||||
! points for integration over space
|
||||
END_DOC
|
||||
@ -86,7 +97,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
|
||||
double precision :: x,r
|
||||
x = grid_points_radial(j) ! x value for the mapping of the [0, +\infty] to [0,1]
|
||||
r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) ! value of the radial coordinate for the integration
|
||||
do k = 1, n_points_integration_angular_lebedev ! explicit values of the grid points centered around each atom
|
||||
do k = 1, n_points_integration_angular ! explicit values of the grid points centered around each atom
|
||||
grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * r
|
||||
grid_points_per_atom(2,k,j,i) = y_ref + angular_quadrature_points(k,2) * r
|
||||
grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * r
|
||||
@ -95,7 +106,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
|
||||
BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
|
||||
BEGIN_DOC
|
||||
! Weight function at grid points : w_n(r) according to the equation (22) of Becke original paper (JCP, 88, 1988)
|
||||
! the "n" discrete variable represents the nucleis which in this array is represented by the last dimension
|
||||
@ -109,7 +120,7 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_int
|
||||
! run over all points in space
|
||||
do j = 1, nucl_num ! that are referred to each atom
|
||||
do k = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
|
||||
do l = 1, n_points_integration_angular_lebedev ! for each angular point attached to the "jth" atom
|
||||
do l = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom
|
||||
r(1) = grid_points_per_atom(1,l,k,j)
|
||||
r(2) = grid_points_per_atom(2,l,k,j)
|
||||
r(3) = grid_points_per_atom(3,l,k,j)
|
||||
@ -129,18 +140,47 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_int
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_integration_angular_lebedev,n_points_radial_grid,nucl_num) ]
|
||||
BEGIN_PROVIDER [double precision, final_weight_functions_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
|
||||
BEGIN_DOC
|
||||
! Weight function at grid points : w_n(r) according to the equation (22) of Becke original paper (JCP, 88, 1988)
|
||||
! the "n" discrete variable represents the nucleis which in this array is represented by the last dimension
|
||||
! and the points are labelled by the other dimensions
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j,k,l,m
|
||||
double precision :: r(3)
|
||||
double precision :: accu,cell_function_becke
|
||||
double precision :: tmp_array(nucl_num)
|
||||
double precision :: contrib_integration,x
|
||||
double precision :: derivative_knowles_function,knowles_function
|
||||
! run over all points in space
|
||||
do j = 1, nucl_num ! that are referred to each atom
|
||||
do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
|
||||
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
|
||||
do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom
|
||||
contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x) &
|
||||
*knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
|
||||
final_weight_functions_at_grid_points(k,i,j) = weights_angular_points(k) * weight_functions_at_grid_points(k,i,j) * contrib_integration * dr_radial_integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
|
||||
&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_integration_angular,n_points_radial_grid,nucl_num,N_states) ]
|
||||
implicit none
|
||||
integer :: i,j,k,l,m,i_state
|
||||
double precision :: contrib
|
||||
double precision :: r(3)
|
||||
double precision :: aos_array(ao_num),mos_array(mo_tot_num)
|
||||
do i_state = 1, N_states
|
||||
do j = 1, nucl_num
|
||||
do k = 1, n_points_radial_grid
|
||||
do l = 1, n_points_integration_angular_lebedev
|
||||
one_body_dm_mo_alpha_at_grid_points(l,k,j) = 0.d0
|
||||
one_body_dm_mo_beta_at_grid_points(l,k,j) = 0.d0
|
||||
do l = 1, n_points_integration_angular
|
||||
one_body_dm_mo_alpha_at_grid_points(l,k,j,i_state) = 0.d0
|
||||
one_body_dm_mo_beta_at_grid_points(l,k,j,i_state) = 0.d0
|
||||
r(1) = grid_points_per_atom(1,l,k,j)
|
||||
r(2) = grid_points_per_atom(2,l,k,j)
|
||||
r(3) = grid_points_per_atom(3,l,k,j)
|
||||
@ -149,14 +189,15 @@ END_PROVIDER
|
||||
do m = 1, mo_tot_num
|
||||
do i = 1, mo_tot_num
|
||||
contrib = mos_array(i) * mos_array(m)
|
||||
one_body_dm_mo_alpha_at_grid_points(l,k,j) += one_body_dm_mo_alpha_average(i,m) * contrib
|
||||
one_body_dm_mo_beta_at_grid_points(l,k,j) += one_body_dm_mo_beta_average(i,m) * contrib
|
||||
one_body_dm_mo_alpha_at_grid_points(l,k,j,i_state) += one_body_dm_mo_alpha(i,m,i_state) * contrib
|
||||
one_body_dm_mo_beta_at_grid_points(l,k,j,i_state) += one_body_dm_mo_beta(i,m,i_state) * contrib
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -5,11 +5,10 @@ double precision function step_function_becke(x)
|
||||
integer :: i,n_max_becke
|
||||
|
||||
step_function_becke = f_function_becke(x)
|
||||
do i = 1, 3
|
||||
do i = 1,5
|
||||
step_function_becke = f_function_becke(step_function_becke)
|
||||
enddo
|
||||
step_function_becke = 0.5d0*(1.d0 - step_function_becke)
|
||||
! step_function_becke = (1.d0 - step_function_becke)
|
||||
end
|
||||
|
||||
double precision function f_function_becke(x)
|
||||
|
||||
@ -4,7 +4,7 @@
|
||||
double precision :: accu
|
||||
integer :: i,j,k,l
|
||||
double precision :: x
|
||||
double precision :: integrand(n_points_integration_angular_lebedev), weights(n_points_integration_angular_lebedev)
|
||||
double precision :: integrand(n_points_integration_angular), weights(n_points_integration_angular)
|
||||
double precision :: f_average_angular_alpha,f_average_angular_beta
|
||||
double precision :: derivative_knowles_function,knowles_function
|
||||
|
||||
@ -12,7 +12,7 @@
|
||||
! according ot equation (6) of the paper of Becke (JCP, (88), 1988)
|
||||
! Here the m index is referred to the w_m(r) weight functions of equation (22)
|
||||
! Run over all points of integrations : there are
|
||||
! n_points_radial_grid (i) * n_points_integration_angular_lebedev (k)
|
||||
! n_points_radial_grid (i) * n_points_integration_angular (k)
|
||||
do j = 1, nucl_num
|
||||
integral_density_alpha_knowles_becke_per_atom(j) = 0.d0
|
||||
integral_density_beta_knowles_becke_per_atom(j) = 0.d0
|
||||
@ -20,9 +20,9 @@
|
||||
! Angular integration over the solid angle Omega for a FIXED angular coordinate "r"
|
||||
f_average_angular_alpha = 0.d0
|
||||
f_average_angular_beta = 0.d0
|
||||
do k = 1, n_points_integration_angular_lebedev
|
||||
f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j) * weight_functions_at_grid_points(k,i,j)
|
||||
f_average_angular_beta += weights_angular_points(k) * one_body_dm_mo_beta_at_grid_points(k,i,j) * weight_functions_at_grid_points(k,i,j)
|
||||
do k = 1, n_points_integration_angular
|
||||
f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j,1) * weight_functions_at_grid_points(k,i,j)
|
||||
f_average_angular_beta += weights_angular_points(k) * one_body_dm_mo_beta_at_grid_points(k,i,j,1) * weight_functions_at_grid_points(k,i,j)
|
||||
enddo
|
||||
!
|
||||
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
|
||||
|
||||
@ -4,8 +4,47 @@ program pouet
|
||||
touch read_wf
|
||||
print*,'m_knowles = ',m_knowles
|
||||
call routine
|
||||
call routine3
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine routine3
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
double precision :: accu
|
||||
accu = 0.d0
|
||||
do j = 1, nucl_num ! that are referred to each atom
|
||||
do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
|
||||
do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom
|
||||
accu += final_weight_functions_at_grid_points(k,i,j) * one_body_dm_mo_alpha_at_grid_points(k,i,j,1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*, accu
|
||||
print*, 'lda_exchange',lda_exchange
|
||||
|
||||
end
|
||||
subroutine routine2
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
double precision :: x,y,z
|
||||
double precision :: r
|
||||
double precision :: accu
|
||||
accu = 0.d0
|
||||
r = 1.d0
|
||||
do k = 1, n_points_integration_angular
|
||||
x = angular_quadrature_points(k,1) * r
|
||||
y = angular_quadrature_points(k,2) * r
|
||||
z = angular_quadrature_points(k,3) * r
|
||||
accu += weights_angular_points(k) * (x**2 + y**2 + z**2)
|
||||
enddo
|
||||
print*, accu
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine routine
|
||||
implicit none
|
||||
integer :: i
|
||||
|
||||
@ -4,7 +4,7 @@ BEGIN_PROVIDER [integer, degree_max_integration_lebedev]
|
||||
! needed for the angular integration according to LEBEDEV formulae
|
||||
END_DOC
|
||||
implicit none
|
||||
degree_max_integration_lebedev= 7
|
||||
degree_max_integration_lebedev= 13
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -644,14 +644,14 @@ END_PROVIDER
|
||||
weights_angular_integration_lebedev(16) = 0.016604069565742d0
|
||||
weights_angular_integration_lebedev(17) = 0.016604069565742d0
|
||||
weights_angular_integration_lebedev(18) = 0.016604069565742d0
|
||||
weights_angular_integration_lebedev(19) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(20) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(21) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(22) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(23) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(24) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(25) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(26) = -0.029586038961039d0
|
||||
weights_angular_integration_lebedev(19) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(20) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(21) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(22) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(23) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(24) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(25) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(26) = 0.029586038961039d0
|
||||
weights_angular_integration_lebedev(27) = 0.026576207082159d0
|
||||
weights_angular_integration_lebedev(28) = 0.026576207082159d0
|
||||
weights_angular_integration_lebedev(29) = 0.026576207082159d0
|
||||
|
||||
@ -10,3 +10,7 @@ double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0)
|
||||
double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0))
|
||||
double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
|
||||
double precision, parameter :: thresh = 1.d-15
|
||||
double precision, parameter :: cx_lda = -0.73855876638202234d0
|
||||
double precision, parameter :: c_2_4_3 = 2.5198420997897464d0
|
||||
double precision, parameter :: cst_lda = -0.93052573634909996d0
|
||||
double precision, parameter :: c_4_3 = 1.3333333333333333d0
|
||||
|
||||
Loading…
Reference in New Issue
Block a user