From 226ca23af860dd9622002c22089d79fbf3fd1af5 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 14 Apr 2017 01:16:21 +0200 Subject: [PATCH] DFT integration works for non overlapping densities --- plugins/DFT_Utils/grid_density.irp.f | 27 ++++++++++++------- .../test_integration_3d_density.irp.f | 3 +++ src/Utils/angular_integration.irp.f | 2 +- 3 files changed, 21 insertions(+), 11 deletions(-) diff --git a/plugins/DFT_Utils/grid_density.irp.f b/plugins/DFT_Utils/grid_density.irp.f index 3f4f1183..59500bb5 100644 --- a/plugins/DFT_Utils/grid_density.irp.f +++ b/plugins/DFT_Utils/grid_density.irp.f @@ -5,7 +5,7 @@ BEGIN_PROVIDER [integer, n_points_radial_grid] implicit none - n_points_radial_grid = 100000 + n_points_radial_grid = 10000 END_PROVIDER @@ -24,15 +24,22 @@ END_PROVIDER double precision :: degre_rad degre_rad = 180.d0/pi accu = 0.d0 - do i = 1, n_points_integration_angular_lebedev - 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)) & - * dsin ( degre_rad * phi_angular_integration_lebedev(i)) - angular_quadrature_points(i,2) = dsin ( degre_rad * theta_angular_integration_lebedev(i)) & - * dsin ( degre_rad * phi_angular_integration_lebedev(i)) - angular_quadrature_points(i,3) = dcos ( degre_rad * phi_angular_integration_lebedev(i)) - enddo +!do i = 1, n_points_integration_angular_lebedev +! 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)) & +! * dsin ( degre_rad * phi_angular_integration_lebedev(i)) +! angular_quadrature_points(i,2) = dsin ( degre_rad * theta_angular_integration_lebedev(i)) & +! * dsin ( degre_rad * phi_angular_integration_lebedev(i)) +! angular_quadrature_points(i,3) = dcos ( degre_rad * phi_angular_integration_lebedev(i)) + +!!weights_angular_points(i) = weights_angular_integration_lebedev(i) +!!angular_quadrature_points(i,1) = dcos ( degre_rad * phi_angular_integration_lebedev(i)) & +!! * dsin ( degre_rad * theta_angular_integration_lebedev(i)) +!!angular_quadrature_points(i,2) = dsin ( degre_rad * phi_angular_integration_lebedev(i)) & +!! * dsin ( degre_rad * theta_angular_integration_lebedev(i)) +!!angular_quadrature_points(i,3) = dcos ( degre_rad * theta_angular_integration_lebedev(i)) +!enddo print*,'ANGULAR' print*,'' print*,'accu = ',accu diff --git a/plugins/DFT_Utils/test_integration_3d_density.irp.f b/plugins/DFT_Utils/test_integration_3d_density.irp.f index 9ff15621..3a1428d0 100644 --- a/plugins/DFT_Utils/test_integration_3d_density.irp.f +++ b/plugins/DFT_Utils/test_integration_3d_density.irp.f @@ -11,6 +11,9 @@ subroutine routine integer :: i double precision :: accu(2) accu = 0.d0 + do i = 1, N_det + call debug_det(psi_det(1,1,i),N_int) + enddo 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) diff --git a/src/Utils/angular_integration.irp.f b/src/Utils/angular_integration.irp.f index a0a7207d..7ffbd01b 100644 --- a/src/Utils/angular_integration.irp.f +++ b/src/Utils/angular_integration.irp.f @@ -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= 3 + degree_max_integration_lebedev= 7 END_PROVIDER