From b14325fef482bdf6cb471b40edf8fa46f2aeac65 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 6 May 2024 18:21:58 +0200 Subject: [PATCH 1/2] Introducing qmckl --- plugins/local/non_h_ints_mu/NEED | 1 + plugins/local/non_h_ints_mu/deb_aos.irp.f | 49 ++++++++-- plugins/local/non_h_ints_mu/qmckl.irp.f | 104 ++++++++++++++++++++++ 3 files changed, 148 insertions(+), 6 deletions(-) diff --git a/plugins/local/non_h_ints_mu/NEED b/plugins/local/non_h_ints_mu/NEED index 48c1c24b..5ca1d543 100644 --- a/plugins/local/non_h_ints_mu/NEED +++ b/plugins/local/non_h_ints_mu/NEED @@ -3,3 +3,4 @@ hamiltonian jastrow ao_tc_eff_map bi_ortho_mos +trexio diff --git a/plugins/local/non_h_ints_mu/deb_aos.irp.f b/plugins/local/non_h_ints_mu/deb_aos.irp.f index c9bc9c9a..86d011fb 100644 --- a/plugins/local/non_h_ints_mu/deb_aos.irp.f +++ b/plugins/local/non_h_ints_mu/deb_aos.irp.f @@ -34,21 +34,58 @@ subroutine print_aos() PROVIDE final_grid_points aos_in_r_array aos_grad_in_r_array aos_lapl_in_r_array - do ipoint = 1, n_points_final_grid - r(:) = final_grid_points(:,ipoint) - print*, r - enddo +! do ipoint = 1, n_points_final_grid +! r(:) = final_grid_points(:,ipoint) +! print*, r +! enddo +double precision :: accu_vgl(5) +double precision :: accu_vgl_nrm(5) do ipoint = 1, n_points_final_grid - r(:) = final_grid_points(:,ipoint) do i = 1, ao_num ao_val = aos_in_r_array (i,ipoint) ao_der(:) = aos_grad_in_r_array(i,ipoint,:) ao_lap = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint) - write(*, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap + write(111, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap enddo enddo + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + ao_val = aos_in_r_array_qmckl (i,ipoint) + ao_der(:) = aos_grad_in_r_array_qmckl(i,ipoint,:) + ao_lap = aos_lapl_in_r_array_qmckl(i,ipoint) + write(222, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap + enddo + enddo + + accu_vgl = 0.d0 + accu_vgl_nrm = 0.d0 + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + ao_val = aos_in_r_array (i,ipoint) + ao_der(:) = aos_grad_in_r_array(i,ipoint,:) + ao_lap = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint) + accu_vgl_nrm(1) += dabs(ao_val) + accu_vgl_nrm(2) += dabs(ao_der(1)) + accu_vgl_nrm(3) += dabs(ao_der(2)) + accu_vgl_nrm(4) += dabs(ao_der(3)) + accu_vgl_nrm(5) += dabs(ao_lap) + + ao_val -= aos_in_r_array_qmckl (i,ipoint) + ao_der(:) -= aos_grad_in_r_array_qmckl(i,ipoint,:) + ao_lap -= aos_lapl_in_r_array_qmckl(i,ipoint) + accu_vgl(1) += dabs(ao_val) + accu_vgl(2) += dabs(ao_der(1)) + accu_vgl(3) += dabs(ao_der(2)) + accu_vgl(4) += dabs(ao_der(3)) + accu_vgl(5) += dabs(ao_lap) + enddo + + enddo + accu_vgl(:) *= 1.d0 / accu_vgl_nrm(:) + print *, accu_vgl + return end diff --git a/plugins/local/non_h_ints_mu/qmckl.irp.f b/plugins/local/non_h_ints_mu/qmckl.irp.f index 1df80457..4d419e24 100644 --- a/plugins/local/non_h_ints_mu/qmckl.irp.f +++ b/plugins/local/non_h_ints_mu/qmckl.irp.f @@ -75,3 +75,107 @@ BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ] endif END_PROVIDER + + + BEGIN_PROVIDER [ double precision, aos_in_r_array_qmckl, (ao_num,n_points_final_grid)] +&BEGIN_PROVIDER [ double precision, aos_grad_in_r_array_qmckl, (ao_num,n_points_final_grid,3)] +&BEGIN_PROVIDER [ double precision, aos_lapl_in_r_array_qmckl, (ao_num, n_points_final_grid)] + implicit none + BEGIN_DOC + ! AOS computed with qmckl + END_DOC + use qmckl + + integer*8 :: qmckl_ctx + integer(qmckl_exit_code) :: rc + + qmckl_ctx = qmckl_context_create() + + rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename))) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in read_trexio' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in set_electron_point' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + double precision, allocatable :: vgl(:,:,:) + allocate( vgl(ao_num,5,n_points_final_grid)) + rc = qmckl_get_ao_basis_ao_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*ao_num*5_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in get_ao_vgl' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + integer :: i,k + do k=1,n_points_final_grid + do i=1,ao_num + aos_in_r_array_qmckl(i,k) = vgl(i,1,k) + aos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k) + aos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k) + aos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k) + aos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k) + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, mos_in_r_array_qmckl, (mo_num,n_points_final_grid)] +&BEGIN_PROVIDER [ double precision, mos_grad_in_r_array_qmckl, (mo_num,n_points_final_grid,3)] +&BEGIN_PROVIDER [ double precision, mos_lapl_in_r_array_qmckl, (mo_num, n_points_final_grid)] + implicit none + BEGIN_DOC + ! moS computed with qmckl + END_DOC + use qmckl + + integer*8 :: qmckl_ctx + integer(qmckl_exit_code) :: rc + + qmckl_ctx = qmckl_context_create() + + rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename))) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in read_trexio' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in set_electron_point' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + double precision, allocatable :: vgl(:,:,:) + allocate( vgl(mo_num,5,n_points_final_grid)) + rc = qmckl_get_mo_basis_mo_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*mo_num*5_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in get_mo_vgl' + rc = qmckl_check(qmckl_ctx, rc) + stop -1 + endif + + integer :: i,k + do k=1,n_points_final_grid + do i=1,mo_num + mos_in_r_array_qmckl(i,k) = vgl(i,1,k) + mos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k) + mos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k) + mos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k) + mos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k) + enddo + enddo + +END_PROVIDER + + From 2a8b9e544b8c9f47ce55dd8f0c4e7df5b0a67ea1 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 7 May 2024 01:56:14 +0200 Subject: [PATCH 2/2] working on aos debug --- plugins/local/non_h_ints_mu/deb_aos.irp.f | 16 +- src/ao_basis/aos_in_r.irp.f | 508 ++++++++++++---------- 2 files changed, 285 insertions(+), 239 deletions(-) diff --git a/plugins/local/non_h_ints_mu/deb_aos.irp.f b/plugins/local/non_h_ints_mu/deb_aos.irp.f index c9bc9c9a..a84e1b91 100644 --- a/plugins/local/non_h_ints_mu/deb_aos.irp.f +++ b/plugins/local/non_h_ints_mu/deb_aos.irp.f @@ -31,12 +31,14 @@ subroutine print_aos() integer :: i, ipoint double precision :: r(3) double precision :: ao_val, ao_der(3), ao_lap + double precision :: mo_val, mo_der(3), mo_lap PROVIDE final_grid_points aos_in_r_array aos_grad_in_r_array aos_lapl_in_r_array + write(1000, *) n_points_final_grid do ipoint = 1, n_points_final_grid r(:) = final_grid_points(:,ipoint) - print*, r + write(1000, '(3(f15.7, 3X))') r enddo do ipoint = 1, n_points_final_grid @@ -45,7 +47,17 @@ subroutine print_aos() ao_val = aos_in_r_array (i,ipoint) ao_der(:) = aos_grad_in_r_array(i,ipoint,:) ao_lap = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint) - write(*, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap + write(1010, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap + enddo + enddo + + do ipoint = 1, n_points_final_grid + r(:) = final_grid_points(:,ipoint) + do i = 1, mo_num + mo_val = mos_in_r_array (i,ipoint) + mo_der(:) = mos_grad_in_r_array(i,ipoint,:) + mo_lap = mos_lapl_in_r_array(i,ipoint,1) + mos_lapl_in_r_array(i,ipoint,2) + mos_lapl_in_r_array(i,ipoint,3) + write(2010, '(5(f15.7, 3X))') mo_val, mo_der, mo_lap enddo enddo diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f index 1b1595a3..053c86a2 100644 --- a/src/ao_basis/aos_in_r.irp.f +++ b/src/ao_basis/aos_in_r.irp.f @@ -1,67 +1,76 @@ -double precision function ao_value(i,r) - implicit none - BEGIN_DOC -! Returns the value of the i-th ao at point $\textbf{r}$ - END_DOC - double precision, intent(in) :: r(3) - integer, intent(in) :: i - integer :: m,num_ao - double precision :: center_ao(3) - double precision :: beta - integer :: power_ao(3) - double precision :: accu,dx,dy,dz,r2 - num_ao = ao_nucl(i) - power_ao(1:3)= ao_power(i,1:3) - center_ao(1:3) = nucl_coord(num_ao,1:3) - dx = (r(1) - center_ao(1)) - dy = (r(2) - center_ao(2)) - dz = (r(3) - center_ao(3)) - r2 = dx*dx + dy*dy + dz*dz - dx = dx**power_ao(1) - dy = dy**power_ao(2) - dz = dz**power_ao(3) +! --- - accu = 0.d0 - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) - enddo - ao_value = accu * dx * dy * dz +double precision function ao_value(i, r) + + BEGIN_DOC + ! Returns the value of the i-th ao at point $\textbf{r}$ + END_DOC + + implicit none + integer, intent(in) :: i + double precision, intent(in) :: r(3) + + integer :: m, num_ao + integer :: power_ao(3) + double precision :: center_ao(3) + double precision :: beta + double precision :: accu, dx, dy, dz, r2 + + num_ao = ao_nucl(i) + power_ao(1:3) = ao_power(i,1:3) + center_ao(1:3) = nucl_coord(num_ao,1:3) + dx = r(1) - center_ao(1) + dy = r(2) - center_ao(2) + dz = r(3) - center_ao(3) + r2 = dx*dx + dy*dy + dz*dz + dx = dx**power_ao(1) + dy = dy**power_ao(2) + dz = dz**power_ao(3) + + accu = 0.d0 + do m = 1, ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) + enddo + ao_value = accu * dx * dy * dz end -double precision function primitive_value(i,j,r) - implicit none - BEGIN_DOC -! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r} -! **without the coefficient** - END_DOC - double precision, intent(in) :: r(3) - integer, intent(in) :: i,j +double precision function primitive_value(i, j, r) - integer :: m,num_ao - double precision :: center_ao(3) - double precision :: beta - integer :: power_ao(3) - double precision :: accu,dx,dy,dz,r2 - num_ao = ao_nucl(i) - power_ao(1:3)= ao_power(i,1:3) - center_ao(1:3) = nucl_coord(num_ao,1:3) - dx = (r(1) - center_ao(1)) - dy = (r(2) - center_ao(2)) - dz = (r(3) - center_ao(3)) - r2 = dx*dx + dy*dy + dz*dz - dx = dx**power_ao(1) - dy = dy**power_ao(2) - dz = dz**power_ao(3) + BEGIN_DOC + ! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r} + ! **without the coefficient** + END_DOC - accu = 0.d0 - m=j - beta = ao_expo_ordered_transp(m,i) - accu += dexp(-beta*r2) - primitive_value = accu * dx * dy * dz + implicit none + integer, intent(in) :: i, j + double precision, intent(in) :: r(3) + + integer :: m, num_ao + integer :: power_ao(3) + double precision :: center_ao(3) + double precision :: beta + double precision :: accu, dx, dy, dz, r2 + + num_ao = ao_nucl(i) + power_ao(1:3)= ao_power(i,1:3) + center_ao(1:3) = nucl_coord(num_ao,1:3) + dx = r(1) - center_ao(1) + dy = r(2) - center_ao(2) + dz = r(3) - center_ao(3) + r2 = dx*dx + dy*dy + dz*dz + dx = dx**power_ao(1) + dy = dy**power_ao(2) + dz = dz**power_ao(3) + + accu = 0.d0 + m = j + beta = ao_expo_ordered_transp(m,i) + accu += dexp(-beta*r2) + primitive_value = accu * dx * dy * dz end @@ -104,9 +113,9 @@ subroutine give_all_aos_at_r(r, tmp_array) dz2 = dz**p_ao(3) tmp_array(k) = 0.d0 - do l = 1,ao_prim_num(k) + do l = 1, ao_prim_num(k) beta = ao_expo_ordered_transp_per_nucl(l,j,i) - if(dabs(beta*r2).gt.40.d0) cycle + if(beta*r2.gt.50.d0) cycle tmp_array(k) += ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) enddo @@ -120,207 +129,232 @@ end ! --- -subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) - implicit none - BEGIN_DOC -! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z -! -! output : -! -! * aos_array(i) = ao(i) evaluated at ro -! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$ -! - END_DOC - double precision, intent(in) :: r(3) - double precision, intent(out) :: aos_array(ao_num) - double precision, intent(out) :: aos_grad_array(3,ao_num) +subroutine give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array) - integer :: power_ao(3) - integer :: i,j,k,l,m - double precision :: dx,dy,dz,r2 - double precision :: dx2,dy2,dz2 - double precision :: dx1,dy1,dz1 - double precision :: center_ao(3) - double precision :: beta,accu_1,accu_2,contrib - do i = 1, nucl_num - center_ao(1:3) = nucl_coord(i,1:3) - dx = (r(1) - center_ao(1)) - dy = (r(2) - center_ao(2)) - dz = (r(3) - center_ao(3)) - r2 = dx*dx + dy*dy + dz*dz - do j = 1,Nucl_N_Aos(i) - k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format - aos_array(k) = 0.d0 - aos_grad_array(1,k) = 0.d0 - aos_grad_array(2,k) = 0.d0 - aos_grad_array(3,k) = 0.d0 - power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i) - dx2 = dx**power_ao(1) - dy2 = dy**power_ao(2) - dz2 = dz**power_ao(3) - if(power_ao(1) .ne. 0)then - dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1) - else - dx1 = 0.d0 - endif - if(power_ao(2) .ne. 0)then - dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1) - else - dy1 = 0.d0 - endif - if(power_ao(3) .ne. 0)then - dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1) - else - dz1 = 0.d0 - endif - accu_1 = 0.d0 - accu_2 = 0.d0 - do l = 1,ao_prim_num(k) - beta = ao_expo_ordered_transp_per_nucl(l,j,i) - contrib = 0.d0 - if(beta*r2.gt.50.d0)cycle - contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) - accu_1 += contrib - accu_2 += contrib * beta - enddo - aos_array(k) = accu_1 * dx2 * dy2 * dz2 - aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2- 2.d0 * dx2 * dx * dy2 * dz2 * accu_2 - aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2- 2.d0 * dx2 * dy2 * dy * dz2 * accu_2 - aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1- 2.d0 * dx2 * dy2 * dz2 * dz * accu_2 + BEGIN_DOC + ! + ! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z + ! + ! output : + ! + ! * aos_array(i) = ao(i) evaluated at ro + ! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$ + ! + END_DOC + + implicit none + double precision, intent(in) :: r(3) + double precision, intent(out) :: aos_array(ao_num) + double precision, intent(out) :: aos_grad_array(3,ao_num) + + integer :: power_ao(3) + integer :: i, j, k, l, m + double precision :: dx, dy, dz, r2 + double precision :: dx1, dy1, dz1 + double precision :: dx2, dy2, dz2 + double precision :: center_ao(3) + double precision :: beta, accu_1, accu_2, contrib + + do i = 1, nucl_num + + center_ao(1:3) = nucl_coord(i,1:3) + + dx = r(1) - center_ao(1) + dy = r(2) - center_ao(2) + dz = r(3) - center_ao(3) + r2 = dx*dx + dy*dy + dz*dz + + do j = 1, Nucl_N_Aos(i) + + k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format + + aos_array(k) = 0.d0 + aos_grad_array(1,k) = 0.d0 + aos_grad_array(2,k) = 0.d0 + aos_grad_array(3,k) = 0.d0 + + power_ao(1:3) = ao_power_ordered_transp_per_nucl(1:3,j,i) + dx2 = dx**power_ao(1) + dy2 = dy**power_ao(2) + dz2 = dz**power_ao(3) + + dx1 = 0.d0 + if(power_ao(1) .ne. 0) then + dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1) + endif + + dy1 = 0.d0 + if(power_ao(2) .ne. 0) then + dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1) + endif + + dz1 = 0.d0 + if(power_ao(3) .ne. 0) then + dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1) + endif + + accu_1 = 0.d0 + accu_2 = 0.d0 + do l = 1, ao_prim_num(k) + beta = ao_expo_ordered_transp_per_nucl(l,j,i) + if(beta*r2.gt.50.d0) cycle + contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) + accu_1 += contrib + accu_2 += contrib * beta + enddo + + aos_array(k) = accu_1 * dx2 * dy2 * dz2 + aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx * dy2 * dz2 * accu_2 + aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy * dz2 * accu_2 + aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz * accu_2 + enddo enddo - enddo + end +! --- -subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) - implicit none - BEGIN_DOC -! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z -! -! output : -! -! * aos_array(i) = ao(i) evaluated at $\textbf{r}$ -! * aos_grad_array(1,i) = $\nabla_x$ of the ao(i) evaluated at $\textbf{r}$ - END_DOC - double precision, intent(in) :: r(3) - double precision, intent(out) :: aos_array(ao_num) - double precision, intent(out) :: aos_grad_array(3,ao_num) - double precision, intent(out) :: aos_lapl_array(3,ao_num) +subroutine give_all_aos_and_grad_and_lapl_at_r(r, aos_array, aos_grad_array, aos_lapl_array) - integer :: power_ao(3) - integer :: i,j,k,l,m - double precision :: dx,dy,dz,r2 - double precision :: dx2,dy2,dz2 - double precision :: dx1,dy1,dz1 - double precision :: dx3,dy3,dz3 - double precision :: dx4,dy4,dz4 - double precision :: dx5,dy5,dz5 - double precision :: center_ao(3) - double precision :: beta,accu_1,accu_2,accu_3,contrib - do i = 1, nucl_num - center_ao(1:3) = nucl_coord(i,1:3) - dx = (r(1) - center_ao(1)) - dy = (r(2) - center_ao(2)) - dz = (r(3) - center_ao(3)) - r2 = dx*dx + dy*dy + dz*dz - do j = 1,Nucl_N_Aos(i) - k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format - aos_array(k) = 0.d0 - aos_grad_array(1,k) = 0.d0 - aos_grad_array(2,k) = 0.d0 - aos_grad_array(3,k) = 0.d0 + BEGIN_DOC + ! + ! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z + ! + ! output : + ! + ! * aos_array(i) = ao(i) evaluated at $\textbf{r}$ + ! * aos_grad_array(1,i) = $\nabla_x$ of the ao(i) evaluated at $\textbf{r}$ + ! + END_DOC - aos_lapl_array(1,k) = 0.d0 - aos_lapl_array(2,k) = 0.d0 - aos_lapl_array(3,k) = 0.d0 + implicit none + double precision, intent(in) :: r(3) + double precision, intent(out) :: aos_array(ao_num) + double precision, intent(out) :: aos_grad_array(3,ao_num) + double precision, intent(out) :: aos_lapl_array(3,ao_num) - power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i) - dx2 = dx**power_ao(1) - dy2 = dy**power_ao(2) - dz2 = dz**power_ao(3) - if(power_ao(1) .ne. 0)then - dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1) - else - dx1 = 0.d0 - endif - ! For the Laplacian - if(power_ao(1) .ge. 2)then - dx3 = dble(power_ao(1)) * dble((power_ao(1)-1)) * dx**(power_ao(1)-2) - else - dx3 = 0.d0 - endif - if(power_ao(1) .ge. 1)then - dx4 = dble((2 * power_ao(1) + 1)) * dx**(power_ao(1)) - else - dx4 = dble((power_ao(1) + 1)) * dx**(power_ao(1)) - endif + integer :: power_ao(3) + integer :: i, j, k, l, m + double precision :: dx, dy, dz, r2 + double precision :: dx1, dy1, dz1 + double precision :: dx2, dy2, dz2 + double precision :: dx3, dy3, dz3 + double precision :: dx4, dy4, dz4 + double precision :: dx5, dy5, dz5 + double precision :: center_ao(3) + double precision :: beta, accu_1, accu_2, accu_3, contrib - dx5 = dx**(power_ao(1)+2) + do i = 1, nucl_num - if(power_ao(2) .ne. 0)then - dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1) - else - dy1 = 0.d0 - endif - ! For the Laplacian - if(power_ao(2) .ge. 2)then - dy3 = dble(power_ao(2)) * dble((power_ao(2)-1)) * dy**(power_ao(2)-2) - else - dy3 = 0.d0 - endif + center_ao(1:3) = nucl_coord(i,1:3) - if(power_ao(2) .ge. 1)then - dy4 = dble((2 * power_ao(2) + 1)) * dy**(power_ao(2)) - else - dy4 = dble((power_ao(2) + 1)) * dy**(power_ao(2)) - endif + dx = r(1) - center_ao(1) + dy = r(2) - center_ao(2) + dz = r(3) - center_ao(3) + r2 = dx*dx + dy*dy + dz*dz + + do j = 1, Nucl_N_Aos(i) - dy5 = dy**(power_ao(2)+2) + k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format + aos_array(k) = 0.d0 + aos_grad_array(1,k) = 0.d0 + aos_grad_array(2,k) = 0.d0 + aos_grad_array(3,k) = 0.d0 + aos_lapl_array(1,k) = 0.d0 + aos_lapl_array(2,k) = 0.d0 + aos_lapl_array(3,k) = 0.d0 + + power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i) + dx2 = dx**power_ao(1) + dy2 = dy**power_ao(2) + dz2 = dz**power_ao(3) - if(power_ao(3) .ne. 0)then - dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1) - else - dz1 = 0.d0 - endif - ! For the Laplacian - if(power_ao(3) .ge. 2)then - dz3 = dble(power_ao(3)) * dble((power_ao(3)-1)) * dz**(power_ao(3)-2) - else - dz3 = 0.d0 - endif + ! --- - if(power_ao(3) .ge. 1)then - dz4 = dble((2 * power_ao(3) + 1)) * dz**(power_ao(3)) - else - dz4 = dble((power_ao(3) + 1)) * dz**(power_ao(3)) - endif + dx1 = 0.d0 + if(power_ao(1) .ne. 0) then + dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1) + endif - dz5 = dz**(power_ao(3)+2) + dx3 = 0.d0 + if(power_ao(1) .ge. 2) then + dx3 = dble(power_ao(1)) * dble((power_ao(1)-1)) * dx**(power_ao(1)-2) + endif + if(power_ao(1) .ge. 1) then + dx4 = dble((2 * power_ao(1) + 1)) * dx**(power_ao(1)) + else + dx4 = dble((power_ao(1) + 1)) * dx**(power_ao(1)) + endif + + dx5 = dx**(power_ao(1)+2) + + ! --- + + dy1 = 0.d0 + if(power_ao(2) .ne. 0) then + dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1) + endif - accu_1 = 0.d0 - accu_2 = 0.d0 - accu_3 = 0.d0 - do l = 1,ao_prim_num(k) - beta = ao_expo_ordered_transp_per_nucl(l,j,i) - contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) - accu_1 += contrib - accu_2 += contrib * beta - accu_3 += contrib * beta**2 - enddo - aos_array(k) = accu_1 * dx2 * dy2 * dz2 + dy3 = 0.d0 + if(power_ao(2) .ge. 2) then + dy3 = dble(power_ao(2)) * dble((power_ao(2)-1)) * dy**(power_ao(2)-2) + endif + + if(power_ao(2) .ge. 1) then + dy4 = dble((2 * power_ao(2) + 1)) * dy**(power_ao(2)) + else + dy4 = dble((power_ao(2) + 1)) * dy**(power_ao(2)) + endif + + dy5 = dy**(power_ao(2)+2) - aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2- 2.d0 * dx2 * dx * dy2 * dz2 * accu_2 - aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2- 2.d0 * dx2 * dy2 * dy * dz2 * accu_2 - aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1- 2.d0 * dx2 * dy2 * dz2 * dz * accu_2 + ! --- + + dz1 = 0.d0 + if(power_ao(3) .ne. 0) then + dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1) + endif - aos_lapl_array(1,k) = accu_1 * dx3 * dy2 * dz2- 2.d0 * dx4 * dy2 * dz2* accu_2 +4.d0 * dx5 *dy2 * dz2* accu_3 - aos_lapl_array(2,k) = accu_1 * dx2 * dy3 * dz2- 2.d0 * dx2 * dy4 * dz2* accu_2 +4.d0 * dx2 *dy5 * dz2* accu_3 - aos_lapl_array(3,k) = accu_1 * dx2 * dy2 * dz3- 2.d0 * dx2 * dy2 * dz4* accu_2 +4.d0 * dx2 *dy2 * dz5* accu_3 + dz3 = 0.d0 + if(power_ao(3) .ge. 2) then + dz3 = dble(power_ao(3)) * dble((power_ao(3)-1)) * dz**(power_ao(3)-2) + endif + + if(power_ao(3) .ge. 1) then + dz4 = dble((2 * power_ao(3) + 1)) * dz**(power_ao(3)) + else + dz4 = dble((power_ao(3) + 1)) * dz**(power_ao(3)) + endif + + dz5 = dz**(power_ao(3)+2) + + ! --- + + accu_1 = 0.d0 + accu_2 = 0.d0 + accu_3 = 0.d0 + do l = 1,ao_prim_num(k) + beta = ao_expo_ordered_transp_per_nucl(l,j,i) + if(beta*r2.gt.50.d0) cycle + contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) + accu_1 += contrib + accu_2 += contrib * beta + accu_3 += contrib * beta**2 + enddo + aos_array(k) = accu_1 * dx2 * dy2 * dz2 + aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx * dy2 * dz2 * accu_2 + aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy * dz2 * accu_2 + aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz * accu_2 + aos_lapl_array(1,k) = accu_1 * dx3 * dy2 * dz2 - 2.d0 * dx4 * dy2 * dz2 * accu_2 + 4.d0 * dx5 * dy2 * dz2 * accu_3 + aos_lapl_array(2,k) = accu_1 * dx2 * dy3 * dz2 - 2.d0 * dx2 * dy4 * dz2 * accu_2 + 4.d0 * dx2 * dy5 * dz2 * accu_3 + aos_lapl_array(3,k) = accu_1 * dx2 * dy2 * dz3 - 2.d0 * dx2 * dy2 * dz4 * accu_2 + 4.d0 * dx2 * dy2 * dz5 * accu_3 + enddo enddo - enddo + end +! ---