diff --git a/config/ifort.cfg b/config/ifort.cfg index 6e6dd389..a738a83c 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [OPTION] -MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below CACHE : 1 ; Enable cache_compile.py OPENMP : 1 ; Append OpenMP flags @@ -50,7 +50,7 @@ FCFLAGS : -xSSE4.2 -O2 -ip -ftz # -xSSE2 : Valgrind needs a very simple x86 executable # [DEBUG] -FC : -g -traceback +FC : -g -traceback -fpe0 FCFLAGS : -xSSE2 -C IRPF90_FLAGS : --openmp diff --git a/plugins/DFT_Utils/grid_density.irp.f b/plugins/DFT_Utils/grid_density.irp.f index 7cb93a91..ac3b702f 100644 --- a/plugins/DFT_Utils/grid_density.irp.f +++ b/plugins/DFT_Utils/grid_density.irp.f @@ -1,11 +1,11 @@ BEGIN_PROVIDER [integer, n_points_angular_grid] implicit none - n_points_angular_grid = 18 + n_points_angular_grid = 51 END_PROVIDER BEGIN_PROVIDER [integer, n_points_radial_grid] implicit none - n_points_radial_grid = 10 + n_points_radial_grid = 1000 END_PROVIDER @@ -16,7 +16,30 @@ END_PROVIDER ! weights and grid points for the integration on the angular variables on ! the unit sphere centered on (0,0,0) END_DOC - call cal_quad(n_points_aangular_grid, angular_quadrature_points,weights_angular_points) + call cal_quad(n_points_angular_grid, angular_quadrature_points,weights_angular_points) + +END_PROVIDER + +BEGIN_PROVIDER [integer , m_knowles] + implicit none + BEGIN_DOC +! value of the "m" parameter in the equation (7) of the paper of Knowles (JCP, 104, 1996) + END_DOC + m_knowles = 3 +END_PROVIDER + + BEGIN_PROVIDER [double precision, grid_points_radial, (n_points_radial_grid)] +&BEGIN_PROVIDER [double precision, dr_radial_integral] + + implicit none + BEGIN_DOC +! points in [0,1] to map the radial integral [0,\infty] + END_DOC + dr_radial_integral = 1.d0/dble(n_points_radial_grid-1) + integer :: i + do i = 1, n_points_radial_grid-1 + grid_points_radial(i) = (i-1) * dr_radial_integral + enddo END_PROVIDER @@ -27,46 +50,51 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid implicit none integer :: i,j,k double precision :: dr,x_ref,y_ref,z_ref - dr = 1.d0/dble(n_points_radial_grid-1) + double precision :: knowles_function do i = 1, nucl_num x_ref = nucl_coord(i,1) y_ref = nucl_coord(i,2) z_ref = nucl_coord(i,3) do j = 1, n_points_radial_grid - do k = 1, n_points_angular_grid - grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * dr - grid_points_per_atom(2,k,j,i) = y_ref + angular_quadrature_points(k,2) * dr - grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * dr + 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_angular_grid ! 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 enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (nucl_num,n_points_angular_grid,n_points_radial_grid) ] +BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_angular_grid,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 (j=1,nucl_num) +! 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) - do j = 1, nucl_num - do k = 1, n_points_radial_grid - do l = 1, n_points_angular_grid - r(1) = grid_points_per_atom(1,j,k,l) - r(2) = grid_points_per_atom(2,j,k,l) - r(3) = grid_points_per_atom(3,j,k,l) + ! run over all points in space + do j = 1, nucl_num ! that are referred to each atom + do k = 1, n_points_radial_grid !for each radial grid attached to the "jth" atom + do l = 1, n_points_angular_grid ! 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) accu = 0.d0 - do i = 1, nucl_num - tmp_array(i) = cell_function_becke(r,i) + do i = 1, nucl_num ! For each of these points in space, ou need to evaluate the P_n(r) + ! function defined for each atom "i" by equation (13) and (21) with k == 3 + tmp_array(i) = cell_function_becke(r,i) ! P_n(r) + ! Then you compute the summ the P_n(r) function for each of the "r" points accu += tmp_array(i) enddo accu = 1.d0/accu - do i = 1, nucl_num - weight_functions_at_grid_points(i,j,k,l) = tmp_array(i)*accu - enddo + weight_functions_at_grid_points(l,k,j) = tmp_array(j) * accu enddo enddo enddo @@ -80,22 +108,34 @@ END_PROVIDER integer :: i,j,k,l,m double precision :: contrib double precision :: r(3) - double precision :: aos_array(ao_num) + double precision :: aos_array(ao_num),mos_array(mo_tot_num) do j = 1, nucl_num do k = 1, n_points_radial_grid do l = 1, n_points_angular_grid - r(1) = grid_points_per_atom(1,j,k,l) - r(2) = grid_points_per_atom(2,j,k,l) - r(3) = grid_points_per_atom(3,j,k,l) - call give_all_aos_at_r(r,aos_array) - one_body_dm_mo_alpha_at_grid_points(j,k,l) = 0.d0 - do i = 1, ao_num - do m = 1, ao_num - contrib = aos_array(i) * aos_array(m) - one_body_dm_mo_alpha_at_grid_points(j,k,l) += one_body_dm_ao_alpha(i,m) * contrib - one_body_dm_mo_beta_at_grid_points(j,k,l) += one_body_dm_ao_beta(i,m) * contrib - enddo + 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 + 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) + +! call give_all_aos_at_r(r,aos_array) +! do i = 1, ao_num +! do m = 1, ao_num +! contrib = aos_array(i) * aos_array(m) +! one_body_dm_mo_alpha_at_grid_points(l,k,j) += one_body_dm_ao_alpha(i,m) * contrib +! one_body_dm_mo_beta_at_grid_points(l,k,j) += one_body_dm_ao_beta(i,m) * contrib +! enddo +! enddo + + call give_all_mos_at_r(r,mos_array) + do i = 1, mo_tot_num + do m = 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(i,m) * contrib + one_body_dm_mo_beta_at_grid_points(l,k,j) += one_body_dm_mo_beta(i,m) * contrib + enddo enddo + enddo enddo enddo diff --git a/plugins/DFT_Utils/integration_3d.irp.f b/plugins/DFT_Utils/integration_3d.irp.f index f4088302..e2c6d8d0 100644 --- a/plugins/DFT_Utils/integration_3d.irp.f +++ b/plugins/DFT_Utils/integration_3d.irp.f @@ -5,7 +5,7 @@ double precision function step_function_becke(x) integer :: i,n_max_becke step_function_becke = f_function_becke(x) - n_max_becke = 2 + n_max_becke = 3 do i = 1, n_max_becke step_function_becke = f_function_becke(step_function_becke) enddo diff --git a/plugins/DFT_Utils/integration_radial.irp.f b/plugins/DFT_Utils/integration_radial.irp.f index 59874b8b..3670db14 100644 --- a/plugins/DFT_Utils/integration_radial.irp.f +++ b/plugins/DFT_Utils/integration_radial.irp.f @@ -3,46 +3,108 @@ implicit none double precision :: accu integer :: i,j,k,l - integer :: m_param_knowles - double precision :: dx,x - integer :: n_pt_int_radial + double precision :: x double precision :: integrand(n_points_angular_grid), weights(n_points_angular_grid) double precision :: f_average_angular_alpha,f_average_angular_beta double precision :: derivative_knowles_function,knowles_function - n_pt_int_radial = 10 - dx = 1.d0/dble(n_pt_int_radial-1) - x = 0.d0 - m_param_knowles = 3 - do j = 1, nucl_num - integral_density_alpha_knowles_becke_per_atom(j) = 0.d0 - do i = 1, n_points_radial_grid - ! Angular integration - f_average_angular_alpha = 0.d0 - f_average_angular_beta = 0.d0 - do k = 1, n_points_angular_grid - 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) + + ! Run over all nuclei in order to perform the Voronoi partition + ! 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_angular_grid (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 + do i = 1, n_points_radial_grid + ! 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_angular_grid + 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) + enddo + ! + x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1] + double precision :: contrib_integration +! print*,m_knowles + 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 + integral_density_alpha_knowles_becke_per_atom(j) += contrib_integration *f_average_angular_alpha + integral_density_beta_knowles_becke_per_atom(j) += contrib_integration *f_average_angular_beta + enddo + integral_density_alpha_knowles_becke_per_atom(j) *= dr_radial_integral + integral_density_beta_knowles_becke_per_atom(j) *= dr_radial_integral enddo - integral_density_alpha_knowles_becke_per_atom(j) += derivative_knowles_function(alpha,m_param_knowles,x) & - *knowles_function(alpha,m_param_knowles,x)**2 & - *f_average_angular_alpha - x += dx - enddo - integral_density_alpha_knowles_becke_per_atom(j) *= dx - enddo END_PROVIDER double precision function knowles_function(alpha,m,x) implicit none + BEGIN_DOC +! function proposed by Knowles (JCP, 104, 1996) for distributing the radial points : +! the Log "m" function ( equation (7) in the paper ) + END_DOC double precision, intent(in) :: alpha,x integer, intent(in) :: m knowles_function = -alpha * dlog(1.d0-x**m) +!knowles_function = 1.d0 end double precision function derivative_knowles_function(alpha,m,x) implicit none + BEGIN_DOC +! derivative of the function proposed by Knowles (JCP, 104, 1996) for distributing the radial points + END_DOC double precision, intent(in) :: alpha,x integer, intent(in) :: m - derivative_knowles_function = m x**(m-1) / (alpha * dlog(1.d0-x**m)) + derivative_knowles_function = alpha * dble(m) * x**(m-1) / (1.d0 - x**m) end + + BEGIN_PROVIDER [double precision, alpha_knowles, (100)] + implicit none + integer :: i + BEGIN_DOC +! recommended values for the alpha parameters according to the paper of Knowles (JCP, 104, 1996) +! as a function of the nuclear charge + END_DOC + + ! H-He + alpha_knowles(1) = 5.d0 + alpha_knowles(2) = 5.d0 + + ! Li-Be + alpha_knowles(3) = 7.d0 + alpha_knowles(4) = 7.d0 + + ! B-Ne + do i = 5, 10 + alpha_knowles(i) = 5.d0 + enddo + + ! Na-Mg + do i = 11, 12 + alpha_knowles(i) = 7.d0 + enddo + + ! Al-Ar + do i = 13, 18 + alpha_knowles(i) = 5.d0 + enddo + + ! K-Ca + do i = 19, 20 + alpha_knowles(i) = 7.d0 + enddo + + ! Sc-Zn + do i = 21, 30 + alpha_knowles(i) = 5.d0 + enddo + + ! Ga-Kr + do i = 31, 36 + alpha_knowles(i) = 7.d0 + enddo + + END_PROVIDER diff --git a/plugins/DFT_Utils/test_integration_3d_density.irp.f b/plugins/DFT_Utils/test_integration_3d_density.irp.f new file mode 100644 index 00000000..93ce58f4 --- /dev/null +++ b/plugins/DFT_Utils/test_integration_3d_density.irp.f @@ -0,0 +1,24 @@ +program pouet + print*,'coucou' + read_wf = .True. + touch read_wf + print*,'m_knowles = ',m_knowles + call routine + +end +subroutine routine + implicit none + integer :: i + double precision :: accu(2) + accu = 0.d0 + do i = 1, nucl_num + accu(1) += integral_density_alpha_knowles_becke_per_atom(i) + accu(2) += integral_density_beta_knowles_becke_per_atom(i) + enddo + print*,'accu(1) = ',accu(1) + print*,'Nalpha = ',elec_alpha_num + print*,'accu(2) = ',accu(2) + print*,'Nalpha = ',elec_beta_num + + +end diff --git a/plugins/Properties/eval_ab_two_body_dm.irp.f b/plugins/Properties/eval_ab_two_body_dm.irp.f deleted file mode 100644 index 595b0087..00000000 --- a/plugins/Properties/eval_ab_two_body_dm.irp.f +++ /dev/null @@ -1,80 +0,0 @@ -program two_bod_ab_dm - implicit none - read_wf = .True. - touch read_wf - call routine - -end - -subroutine routine - implicit none - integer :: i,j,k,l - double precision :: dx1,dx2 - double precision :: interval_1,interval_2 - integer :: nx1,nx2 - double precision :: r1(3),r2(3) - double precision :: xmin_1,xmax_1,xmin_2,xmax_2 - double precision :: compute_extra_diag_two_body_dm_ab,two_bod_extra_diag - double precision :: compute_diag_two_body_dm_ab,two_bod_diag - - double precision,allocatable :: mos_array_r1(:),mos_array_r2(:) - double precision :: test_diag, test_extra_diag - double precision, allocatable :: x_array(:),y_array(:),z_diag_array(:),z_extra_diag_array(:),z_total_array(:) - allocate(mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num)) - - - - ! He triplet S - ! r1 = x - ! r2 = z - nx1 = 100 - nx2 = 100 - allocate(x_array(nx1),y_array(nx2),z_diag_array(nx1),z_extra_diag_array(nx1),z_total_array(nx1)) - xmin_1 = -2.d0 - xmax_1 = 2.d0 - xmin_2 = -2.d0 - xmax_2 = 2.d0 - interval_1 = xmax_1 - xmin_1 - interval_2 = xmax_2 - xmin_2 - dx1 = interval_1/dble(nx1) - dx2 = interval_2/dble(nx2) - r1 = 0.d0 - r2 = 0.d0 - - double precision :: x_tmp,y_tmp - x_tmp = xmin_1 - do i = 1, nx1 - x_array(i) = x_tmp - write(i_unit_x_two_body_dm_ab,'(10000(F16.10,X))')x_array(i) - x_tmp += dx1 - enddo - x_tmp = xmin_2 - do i = 1, nx1 - y_array(i) = x_tmp - write(i_unit_y_two_body_dm_ab,'(10000(F16.10,X))')x_array(i) - x_tmp += dx2 - enddo - - - ! initialization - r1(1) = xmin_1 - do i = 1, nx1 - r2(3) = xmin_2 - do j = 1, nx2 - two_bod_extra_diag = compute_extra_diag_two_body_dm_ab(r1,r2) - two_bod_diag = compute_diag_two_body_dm_ab(r1,r2) - z_diag_array(j) = two_bod_diag - z_extra_diag_array(j) = two_bod_extra_diag - z_total_array(j) = two_bod_extra_diag + two_bod_diag -! write(i_unit_two_body_dm_ab,'(100(F16.10,X))')r1(1),r2(3),two_bod_diag,two_bod_extra_diag,two_bod_diag+two_bod_extra_diag - r2(3) += dx2 - enddo - write(i_unit_z_two_body_diag_dm_ab,'(10000(F16.10,X))')z_diag_array(:) - write(i_unit_z_two_body_extra_diag_dm_ab,'(10000(F16.10,X))')z_extra_diag_array(:) - write(i_unit_z_two_body_total_dm_ab,'(10000(F16.10,X))')z_total_array(:) - r1(1) += dx1 - enddo - - deallocate(mos_array_r1,mos_array_r2) - -end diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f index c89b1125..f88d6ea3 100644 --- a/src/Determinants/two_body_dm_map.irp.f +++ b/src/Determinants/two_body_dm_map.irp.f @@ -242,7 +242,24 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb if(degree>2)cycle call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! if(i==3.or.j==3)then +! print*,'i,j = ',i,j +! call debug_det(psi_det(1,1,i),N_int) +! call debug_det(psi_det(1,1,j),N_int) +! print*,degree,s1,s2 +! print*,h1,p1,h2,p2 +! print*,phase +! pause +! endif contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase +! print*,'coucou' +! print*,'i,j = ',i,j +! print*,'contrib = ',contrib +! print*,h1,p1,h2,p2 +! print*,'s1,s2',s1,s2 +! call debug_det(psi_det(1,1,i),N_int) +! call debug_det(psi_det(1,1,j),N_int) +! pause if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************ if(s1==s2)cycle ! Only the alpha/beta two body density matrix ! * c_I * c_J @@ -288,7 +305,7 @@ end double precision function compute_extra_diag_two_body_dm_ab(r1,r2) implicit none - double precision :: r1(3),r2(3) + double precision, intent(in) :: r1(3),r2(3) integer :: i,j,k,l double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num) double precision :: contrib @@ -298,16 +315,16 @@ double precision function compute_extra_diag_two_body_dm_ab(r1,r2) call give_all_act_mos_at_r(r1,mos_array_r1) call give_all_act_mos_at_r(r2,mos_array_r2) do l = 1, n_act_orb ! p2 - contrib = mos_array_r2(l) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle do k = 1, n_act_orb ! h2 -! contrib *= mos_array_r2(k) -! if(dabs(contrib*mos_array_r2(k)).lt.threshld_two_bod_dm)cycle do j = 1, n_act_orb ! p1 -! contrib *= mos_array_r1(j) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle do i = 1,n_act_orb ! h1 double precision :: contrib_tmp +! print*,'i,j',i,j +! print*,mos_array_r1(i) , mos_array_r1(j) +! print*,'k,l',k,l +! print*,mos_array_r2(k) * mos_array_r2(l) +! print*,'gama = ',two_body_dm_ab_big_array(i,j,k,l) +! pause contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l) compute_extra_diag_two_body_dm_ab += two_body_dm_ab_big_array(i,j,k,l) * contrib_tmp enddo