1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-04-28 19:34:46 +02:00

Fix bug in HPC laplacian AO

This commit is contained in:
Anthony Scemama 2024-01-11 14:33:40 +01:00
parent 7fc10a47a1
commit d257e28b92
3 changed files with 20 additions and 13 deletions

View File

@ -114,6 +114,12 @@ int main() {
#+end_src #+end_src
#+NAME:fortran_cutoff
#+begin_src fortran
cutoff = qmckl_get_numprec_precision(context)
cutoff = dlog(2.d0**(cutoff-2))
#+end_src
* Context * Context
** Constant data ** Constant data
@ -3605,7 +3611,7 @@ function qmckl_compute_ao_basis_primitive_gaussian_vgl &
use qmckl_constants use qmckl_constants
use qmckl, only: qmckl_get_numprec_epsilon use qmckl, only: qmckl_get_numprec_precision
implicit none implicit none
integer (qmckl_context), intent(in) , value :: context integer (qmckl_context), intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: prim_num integer (c_int64_t) , intent(in) , value :: prim_num
@ -3624,7 +3630,7 @@ function qmckl_compute_ao_basis_primitive_gaussian_vgl &
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
! Don't compute exponentials when the result will be almost zero. ! Don't compute exponentials when the result will be almost zero.
cutoff = -dlog(qmckl_get_numprec_epsilon(context)) <<fortran_cutoff>>
do inucl=1,nucl_num do inucl=1,nucl_num
! C is zero-based, so shift bounds by one ! C is zero-based, so shift bounds by one
@ -3908,7 +3914,7 @@ function qmckl_compute_ao_basis_shell_gaussian_vgl( &
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
use qmckl_constants use qmckl_constants
use qmckl, only: qmckl_get_numprec_epsilon use qmckl, only: qmckl_get_numprec_precision
implicit none implicit none
integer (qmckl_context), intent(in) , value :: context integer (qmckl_context), intent(in) , value :: context
@ -3936,7 +3942,7 @@ function qmckl_compute_ao_basis_shell_gaussian_vgl( &
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
! Don't compute exponentials when the result will be almost zero. ! Don't compute exponentials when the result will be almost zero.
cutoff = -dlog(qmckl_get_numprec_epsilon(context)) <<fortran_cutoff>>
do ipoint = 1, point_num do ipoint = 1, point_num
@ -4358,7 +4364,7 @@ function test_qmckl_ao_power(context) bind(C)
integer*8 :: i,j integer*8 :: i,j
double precision :: epsilon double precision :: epsilon
epsilon = qmckl_get_numprec_epsilon(context) epsilon = qmckl_get_numprec_precision(context)
n = 100; n = 100;
LDP = 10; LDP = 10;
@ -5247,7 +5253,7 @@ function test_qmckl_ao_polynomial_vgl(context) bind(C)
double precision :: epsilon double precision :: epsilon
integer(qmckl_exit_code) :: test_qmckl_ao_polynomial_vgl integer(qmckl_exit_code) :: test_qmckl_ao_polynomial_vgl
epsilon = qmckl_get_numprec_epsilon(context) epsilon = qmckl_get_numprec_precision(context)
X = (/ 1.1 , 2.2 , 3.3 /) X = (/ 1.1 , 2.2 , 3.3 /)
R = (/ 0.1 , 1.2 , -2.3 /) R = (/ 0.1 , 1.2 , -2.3 /)
@ -5405,7 +5411,7 @@ function qmckl_compute_nucleus_range_gaussian(context, &
bind(C) result(info) bind(C) result(info)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
use qmckl_constants use qmckl_constants
use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_epsilon use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
implicit none implicit none
integer (qmckl_context), intent(in) , value :: context integer (qmckl_context), intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: ao_num integer (c_int64_t) , intent(in) , value :: ao_num
@ -5566,7 +5572,7 @@ function qmckl_compute_ao_value_doc(context, &
bind(C) result(info) bind(C) result(info)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
use qmckl_constants use qmckl_constants
use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_epsilon use qmckl, only: qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
implicit none implicit none
integer (qmckl_context), intent(in) , value :: context integer (qmckl_context), intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: ao_num integer (c_int64_t) , intent(in) , value :: ao_num
@ -5618,7 +5624,7 @@ function qmckl_compute_ao_value_doc(context, &
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
! Don't compute polynomials when the radial part is zero. ! Don't compute polynomials when the radial part is zero.
cutoff = -dlog(qmckl_get_numprec_epsilon(context)) <<fortran_cutoff>>
do ipoint = 1, point_num do ipoint = 1, point_num
e_coord(1) = coord(ipoint,1) e_coord(1) = coord(ipoint,1)
@ -6354,7 +6360,7 @@ function qmckl_compute_ao_vgl_doc(context, &
bind(C) result(info) bind(C) result(info)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
use qmckl_constants use qmckl_constants
use qmckl, only : qmckl_ao_polynomial_vgl, qmckl_get_numprec_epsilon use qmckl, only : qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision
implicit none implicit none
integer (qmckl_context), intent(in) , value :: context integer (qmckl_context), intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: ao_num integer (c_int64_t) , intent(in) , value :: ao_num
@ -6405,7 +6411,7 @@ function qmckl_compute_ao_vgl_doc(context, &
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
! Don't compute polynomials when the radial part is zero. ! Don't compute polynomials when the radial part is zero.
cutoff = -dlog(qmckl_get_numprec_epsilon(context)) <<fortran_cutoff>>
do ipoint = 1, point_num do ipoint = 1, point_num
e_coord(1) = coord(ipoint,1) e_coord(1) = coord(ipoint,1)
@ -6805,7 +6811,7 @@ IVDEP
ao_vgl_2[0] = s2 * f[0]; ao_vgl_2[0] = s2 * f[0];
ao_vgl_3[0] = s3 * f[0]; ao_vgl_3[0] = s3 * f[0];
ao_vgl_4[0] = s4 * f[0]; ao_vgl_4[0] = s4 * f[0];
ao_vgl_5[0] = s5; ao_vgl_5[0] = s5 * f[0];
break; break;
case 3: case 3:
IVDEP IVDEP

View File

@ -119,6 +119,7 @@ int main() {
#include <string.h> #include <string.h>
#include <stdbool.h> #include <stdbool.h>
#include <assert.h> #include <assert.h>
#include <math.h>
#include "qmckl.h" #include "qmckl.h"
#include "qmckl_context_private_type.h" #include "qmckl_context_private_type.h"

View File

@ -368,7 +368,7 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
*** Number of atomic orbitals *** Number of atomic orbitals
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c)
int64_t ao_num = 0LL; int64_t ao_num = 0;
rcio = trexio_read_ao_num_64(file, &ao_num); rcio = trexio_read_ao_num_64(file, &ao_num);
if (rcio != TREXIO_SUCCESS) { if (rcio != TREXIO_SUCCESS) {