From 96c17686b447cca1dd0a5e97446fc76a6e87280f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 28 May 2019 18:49:21 +0200 Subject: [PATCH] fixed bugs with dummy atom and becke grid --- src/becke_numerical_grid/atomic_number.irp.f | 9 +++++++++ src/becke_numerical_grid/grid_becke.irp.f | 6 +++--- src/becke_numerical_grid/grid_becke_vector.irp.f | 3 ++- src/becke_numerical_grid/step_function_becke.irp.f | 5 ----- src/nuclei/atomic_radii.irp.f | 4 ++-- 5 files changed, 16 insertions(+), 11 deletions(-) create mode 100644 src/becke_numerical_grid/atomic_number.irp.f diff --git a/src/becke_numerical_grid/atomic_number.irp.f b/src/becke_numerical_grid/atomic_number.irp.f new file mode 100644 index 00000000..eea1fad7 --- /dev/null +++ b/src/becke_numerical_grid/atomic_number.irp.f @@ -0,0 +1,9 @@ +BEGIN_PROVIDER [ integer, grid_atomic_number, (nucl_num) ] + implicit none + BEGIN_DOC + ! Atomic number used to adjust the grid + END_DOC + grid_atomic_number(:) = max(1,int(nucl_charge(:))) + +END_PROVIDER + diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f index 38d4053f..e72f6460 100644 --- a/src/becke_numerical_grid/grid_becke.irp.f +++ b/src/becke_numerical_grid/grid_becke.irp.f @@ -146,7 +146,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_ x = grid_points_radial(j) ! value of the radial coordinate for the integration - r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) + r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x) ! explicit values of the grid points centered around each atom do k = 1, n_points_integration_angular @@ -232,8 +232,8 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul 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 + contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)& + *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2 final_weight_at_r(k,i,j) = weights_angular_points(k) * weight_at_r(k,i,j) * contrib_integration * dr_radial_integral if(isnan(final_weight_at_r(k,i,j)))then print*,'isnan(final_weight_at_r(k,i,j))' diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f index 9856bcbd..3c2a6b91 100644 --- a/src/becke_numerical_grid/grid_becke_vector.irp.f +++ b/src/becke_numerical_grid/grid_becke_vector.irp.f @@ -1,5 +1,6 @@ BEGIN_PROVIDER [integer, n_points_final_grid] + implicit none BEGIN_DOC ! Number of points which are non zero END_DOC @@ -8,7 +9,7 @@ BEGIN_PROVIDER [integer, n_points_final_grid] do j = 1, nucl_num do i = 1, n_points_radial_grid -1 do k = 1, n_points_integration_angular - if(dabs(final_weight_at_r(k,i,j)) < tresh_grid)then + if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then cycle endif n_points_final_grid += 1 diff --git a/src/becke_numerical_grid/step_function_becke.irp.f b/src/becke_numerical_grid/step_function_becke.irp.f index b6335c3d..2905c6c0 100644 --- a/src/becke_numerical_grid/step_function_becke.irp.f +++ b/src/becke_numerical_grid/step_function_becke.irp.f @@ -31,10 +31,6 @@ double precision function cell_function_becke(r,atom_number) double precision :: mu_ij,nu_ij double precision :: distance_i,distance_j,step_function_becke integer :: j - if(int(nucl_charge(atom_number))==0)then - cell_function_becke = 0.d0 - return - endif distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number)) distance_i += (r(2) - nucl_coord_transp(2,atom_number) ) * (r(2) - nucl_coord_transp(2,atom_number)) distance_i += (r(3) - nucl_coord_transp(3,atom_number) ) * (r(3) - nucl_coord_transp(3,atom_number)) @@ -42,7 +38,6 @@ double precision function cell_function_becke(r,atom_number) cell_function_becke = 1.d0 do j = 1, nucl_num if(j==atom_number)cycle - if(int(nucl_charge(j))==0)cycle distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j)) distance_j+= (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j)) distance_j+= (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j)) diff --git a/src/nuclei/atomic_radii.irp.f b/src/nuclei/atomic_radii.irp.f index 82487b9d..439b5cec 100644 --- a/src/nuclei/atomic_radii.irp.f +++ b/src/nuclei/atomic_radii.irp.f @@ -66,7 +66,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)] implicit none integer :: i do i = 1, nucl_num - slater_bragg_radii_per_atom(i) = slater_bragg_radii(int(nucl_charge(i))) + slater_bragg_radii_per_atom(i) = slater_bragg_radii(max(1,int(nucl_charge(i)))) enddo END_PROVIDER @@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)] implicit none integer :: i do i = 1, nucl_num - slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(int(nucl_charge(i))) + slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(max(1,int(nucl_charge(i)))) enddo END_PROVIDER