From a65902aa33da727fa48ce69c81a0fd4caad7d75d Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Tue, 25 Jan 2022 18:20:03 +0100
Subject: [PATCH] fixed bug in Laplacians

---
 src/dft_utils_in_r/ao_in_r.irp.f | 34 +++++++++++++++++++++-----------
 src/dft_utils_in_r/mo_in_r.irp.f |  2 +-
 2 files changed, 24 insertions(+), 12 deletions(-)

diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f
index 6fa6a4c7..38478d21 100644
--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@@ -91,7 +91,19 @@
  enddo
  END_PROVIDER
 
- BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)]
+ BEGIN_PROVIDER [double precision, aos_lapl_in_r_array_transp, (ao_num, n_points_final_grid,3)]
+ implicit none
+ integer :: i,j,m
+ do i = 1, n_points_final_grid
+  do j = 1, ao_num
+   do m = 1, 3
+    aos_lapl_in_r_array_transp(j,i,m) =  aos_lapl_in_r_array(m,j,i)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, aos_lapl_in_r_array, (3,ao_num,n_points_final_grid)]
  implicit none
  BEGIN_DOC
  ! aos_lapl_in_r_array(i,j,k)   = value of the kth component of the laplacian of jth ao on the ith grid point
@@ -100,20 +112,20 @@
  END_DOC
  integer :: i,j,m
  double precision :: aos_array(ao_num), r(3)
- double precision :: aos_grad_array(ao_num,3)
- double precision :: aos_lapl_array(ao_num,3)
+ double precision :: aos_grad_array(3,ao_num)
+ double precision :: aos_lapl_array(3,ao_num)
  !$OMP PARALLEL DO &
  !$OMP DEFAULT (NONE)  &
  !$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) & 
  !$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points)
- do m = 1, 3
-  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)
-   call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
-   do j = 1, ao_num
-    aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m)
+ 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)
+  call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
+  do j = 1, ao_num
+   do m = 1, 3
+    aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j)
    enddo
   enddo
  enddo
diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f
index 0a8b4d52..192cb25a 100644
--- a/src/dft_utils_in_r/mo_in_r.irp.f
+++ b/src/dft_utils_in_r/mo_in_r.irp.f
@@ -138,7 +138,7 @@
  integer :: m
  mos_lapl_in_r_array = 0.d0
  do m=1,3
-  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num)
+  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array_transp(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num)
  enddo
  END_PROVIDER