The latest version fo QMCkl can be downloaded @@ -371,8 +371,8 @@ The latest version fo QMCkl can be downloaded
QMCkl is built with GNU Autotools, so the usual
@@ -387,8 +387,8 @@ options are defined using CFLAGS and FCFLAGS.
To compile from the source repository, additional dependencies are @@ -409,8 +409,8 @@ to be executed first.
The qmckl.h header file installed in the ${prefix}/include directory
@@ -439,12 +439,12 @@ Both files are located in the include/ directory.
In a traditional source code, most of the lines of source files of a program @@ -494,8 +494,8 @@ tarball contains the generated source code.
For a tutorial on literate programming with org-mode, follow this link. @@ -526,8 +526,8 @@ org-mode.
Most of the codes of the TREX CoE are written in Fortran with some @@ -591,8 +591,8 @@ For more guidelines on using Fortran to generate a C interface, see
The authors should follow the recommendations of the C99 @@ -612,8 +612,8 @@ make cppcheck ; cat cppcheck.out
The proposed API should allow the library to: deal with memory transfers @@ -624,8 +624,8 @@ functions (see below).
To avoid namespace collisions, we use qmckl_ as a prefix for all exported
@@ -646,8 +646,8 @@ form is allowed.
In the C language, the number of bits used by the integer types can change @@ -679,15 +679,15 @@ bindings in other languages in other repositories.
Global variables should be avoided in the library, because it is
possible that one single program needs to use multiple instances
of the library. To solve this problem we propose to use a pointer
to a context variable, built by the library with the
-qmckl_context_create function. The =context= contains the global
+qmckl_context_create function. The =context= contains the global
state of the library, and is used as the first argument of many
QMCkl functions.
A single qmckl.h header to be distributed by the library
@@ -790,8 +790,8 @@ and the types definitions should be written in the *fh_type.f90 fil
Low-level functions are very simple functions which are leaves of @@ -800,14 +800,14 @@ the function call tree (they don't call any other QMCkl function).
These functions are pure, and unaware of the QMCkl
-context. They are not allowed to allocate/deallocate memory, and
+context. They are not allowed to allocate/deallocate memory, and
if they need temporary memory it should be provided in input.
High-level functions are at the top of the function call tree. @@ -819,8 +819,8 @@ temporary storage, to simplify the use of accelerators.
The minimal number of bits of precision required for a function
@@ -828,7 +828,7 @@ should be given as an input of low-level computational
functions. This input will be used to define the values of the
different thresholds that might be used to avoid computing
unnecessary noise. High-level functions will use the precision
-specified in the context variable.
+specified in the context variable.
@@ -896,8 +896,8 @@ following points :
Reducing the scaling of an algorithm usually implies also reducing @@ -913,7 +913,7 @@ implemented adapted to different problem sizes.
The atomic basis set is defined as a list of shells. Each shell \(s\) is @@ -432,19 +438,19 @@ gradients and Laplacian of the atomic basis functions.
The following arrays are stored in the context, and need to be set when initializing the library:
-primitive_vgl |
-double[5][point_num][prim_num] |
+double[point_num][5][prim_num] |
Value, gradients, Laplacian of the primitives at current positions | |
shell_vgl |
-double[5][point_num][shell_num] |
+double[point_num][5][shell_num] |
Value, gradients, Laplacian of the primitives at current positions | |
ao_vgl |
-double[5][point_num][ao_num] |
+double[point_num][5][ao_num] |
Value, gradients, Laplacian of the primitives at current positions |
primitive_vgl |
-double[5][point_num][prim_num] |
+double[point_num][5][prim_num] |
out | Value, gradients and Laplacian of the primitives |
shell_vgl |
-double[5][point_num][shell_num] |
+double[point_num][5][shell_num] |
out | Value, gradients and Laplacian of the shells |
shell_vgl |
-double[5][point_num][shell_num] |
+double[point_num][5][shell_num] |
in | Value, gradients and Laplacian of the shells |
ao_vgl |
-double[5][point_num][ao_num] |
+double[point_num][5][ao_num] |
out | Value, gradients and Laplacian of the AOs |
integer function qmckl_compute_ao_vgl_f(context, & +integer function qmckl_compute_ao_vgl_doc_f(context, & ao_num, shell_num, point_num, nucl_num, & coord, nucl_coord, nucleus_index, nucleus_shell_num, & nucleus_range, nucleus_max_ang_mom, shell_ang_mom, & @@ -2667,12 +2821,12 @@ For example, with a=0, b=2 and c=1 the string is "yyz" integer , intent(in) :: nucleus_max_ang_mom(nucl_num) integer , intent(in) :: shell_ang_mom(shell_num) double precision , intent(in) :: ao_factor(ao_num) - double precision , intent(in) :: shell_vgl(shell_num,point_num,5) - double precision , intent(out) :: ao_vgl(ao_num,point_num,5) + double precision , intent(in) :: shell_vgl(shell_num,5,point_num) + double precision , intent(out) :: ao_vgl(ao_num,5,point_num) double precision :: e_coord(3), n_coord(3) integer*8 :: n_poly - integer :: l, il, k, m, n + integer :: l, il, k integer*8 :: ipoint, inucl, ishell integer*8 :: ishell_start, ishell_end integer :: lstart(0:20) @@ -2682,103 +2836,239 @@ For example, with a=0, b=2 and c=1 the string is "yyz" double precision, allocatable :: poly_vgl(:,:) integer , allocatable :: powers(:,:) - integer , allocatable :: kil(:), knucl(:), kshell(:) - allocate(poly_vgl(8,ao_num), powers(8,ao_num)) - allocate(kil(ao_num), kshell(ao_num), knucl(nucl_num+1)) + allocate(poly_vgl(5,ao_num), powers(3,ao_num)) ! Pre-computed data - - k=1 - do inucl=1,nucl_num - knucl(inucl) = k - ishell_start = nucleus_index(inucl) + 1 - ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) - do ishell = ishell_start, ishell_end - l = shell_ang_mom(ishell) - m = l*(l+1)*(l+2)/6 +1 - n = (l+1)*(l+2)*(l+3)/6 - do il = m, n - kil(k) = il - kshell(k) = ishell - k = k+1 - end do - end do + do l=0,20 + lstart(l) = l*(l+1)*(l+2)/6 +1 end do - knucl(nucl_num+1) = ao_num+1 - info = QMCKL_SUCCESS ! Don't compute polynomials when the radial part is zero. - ! TODO : Use numerical precision here cutoff = -dlog(1.d-12) do ipoint = 1, point_num e_coord(1) = coord(ipoint,1) e_coord(2) = coord(ipoint,2) e_coord(3) = coord(ipoint,3) - - ! Express the radial part in the AO basis - + k=1 do inucl=1,nucl_num - n_coord(1) = nucl_coord(inucl,1) - n_coord(2) = nucl_coord(inucl,2) - n_coord(3) = nucl_coord(inucl,3) + n_coord(1) = nucl_coord(inucl,1) + n_coord(2) = nucl_coord(inucl,2) + n_coord(3) = nucl_coord(inucl,3) ! Test if the point is in the range of the nucleus x = e_coord(1) - n_coord(1) y = e_coord(2) - n_coord(2) z = e_coord(3) - n_coord(3) - r2 = x*x + y*y + z*z + r2 = x*x + z*z + z*z if (r2 > cutoff*nucleus_range(inucl)) then - do k = knucl(inucl), knucl(inucl+1)-1 - ao_vgl(k,ipoint,1) = 0.d0 - ao_vgl(k,ipoint,2) = 0.d0 - ao_vgl(k,ipoint,3) = 0.d0 - ao_vgl(k,ipoint,4) = 0.d0 - ao_vgl(k,ipoint,5) = 0.d0 - end do cycle end if - ! Compute polynomials + ! Compute polynomials info = qmckl_ao_polynomial_vgl_f(context, e_coord, n_coord, & - nucleus_max_ang_mom(inucl), n_poly, powers, 8_8, poly_vgl, 8_8) + nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, & + poly_vgl, 5_8) - do k = knucl(inucl), knucl(inucl+1)-1 - y = shell_vgl(kshell(k),ipoint,1) * ao_factor(k) + ! Loop over shells + ishell_start = nucleus_index(inucl) + 1 + ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) + do ishell = ishell_start, ishell_end + l = shell_ang_mom(ishell) + do il = lstart(l), lstart(l+1)-1 + ! Value + ao_vgl(k,1,ipoint) = & + poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k) - ao_vgl(k,ipoint,1) = y * poly_vgl(1,kil(k)) - ao_vgl(k,ipoint,2) = y * poly_vgl(2,kil(k)) - ao_vgl(k,ipoint,3) = y * poly_vgl(3,kil(k)) - ao_vgl(k,ipoint,4) = y * poly_vgl(4,kil(k)) - ao_vgl(k,ipoint,5) = y * poly_vgl(5,kil(k)) + ! Grad_x + ao_vgl(k,2,ipoint) = ( & + poly_vgl(2,il) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(1,il) * shell_vgl(ishell,2,ipoint) & + ) * ao_factor(k) - x = poly_vgl(1,kil(k)) * ao_factor(k) + ! Grad_y + ao_vgl(k,3,ipoint) = ( & + poly_vgl(3,il) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(1,il) * shell_vgl(ishell,3,ipoint) & + ) * ao_factor(k) - ao_vgl(k,ipoint,2) = ao_vgl(k,ipoint,2) + x * shell_vgl(kshell(k),ipoint,2) - ao_vgl(k,ipoint,3) = ao_vgl(k,ipoint,3) + x * shell_vgl(kshell(k),ipoint,3) - ao_vgl(k,ipoint,4) = ao_vgl(k,ipoint,4) + x * shell_vgl(kshell(k),ipoint,4) - ao_vgl(k,ipoint,5) = ao_vgl(k,ipoint,5) + x * shell_vgl(kshell(k),ipoint,5) + ! Grad_z + ao_vgl(k,4,ipoint) = ( & + poly_vgl(4,il) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(1,il) * shell_vgl(ishell,4,ipoint) & + ) * ao_factor(k) - ao_vgl(k,ipoint,5) = ao_vgl(k,ipoint,5) + & - (ao_factor(k) + ao_factor(k)) * (& - poly_vgl(2,kil(k)) * shell_vgl(kshell(k),ipoint,2) + & - poly_vgl(3,kil(k)) * shell_vgl(kshell(k),ipoint,3) + & - poly_vgl(4,kil(k)) * shell_vgl(kshell(k),ipoint,4) ) + ! Lapl_z + ao_vgl(k,5,ipoint) = ( & + poly_vgl(5,il) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(1,il) * shell_vgl(ishell,5,ipoint) + & + 2.d0 * ( & + poly_vgl(2,il) * shell_vgl(ishell,2,ipoint) + & + poly_vgl(3,il) * shell_vgl(ishell,3,ipoint) + & + poly_vgl(4,il) * shell_vgl(ishell,4,ipoint) ) & + ) * ao_factor(k) + + k = k+1 + end do end do - end do end do - deallocate(poly_vgl, powers, kshell, kil, knucl) -end function qmckl_compute_ao_vgl_f + deallocate(poly_vgl, powers) +end function qmckl_compute_ao_vgl_doc_f
integer function qmckl_compute_ao_vgl_hpc_f(context, & + ao_num, shell_num, point_num, nucl_num, & + coord, nucl_coord, nucleus_index, nucleus_shell_num, & + nucleus_range, nucleus_max_ang_mom, shell_ang_mom, & + ao_factor, shell_vgl, ao_vgl) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: ao_num + integer*8 , intent(in) :: shell_num + integer*8 , intent(in) :: point_num + integer*8 , intent(in) :: nucl_num + double precision , intent(in) :: coord(point_num,3) + double precision , intent(in) :: nucl_coord(nucl_num,3) + integer*8 , intent(in) :: nucleus_index(nucl_num) + integer*8 , intent(in) :: nucleus_shell_num(nucl_num) + double precision , intent(in) :: nucleus_range(nucl_num) + integer , intent(in) :: nucleus_max_ang_mom(nucl_num) + integer , intent(in) :: shell_ang_mom(shell_num) + double precision , intent(in) :: ao_factor(ao_num) + double precision , intent(in) :: shell_vgl(shell_num,5,point_num) + double precision , intent(out) :: ao_vgl(ao_num,5,point_num) + + double precision :: e_coord(3), n_coord(3) + integer*8 :: n_poly + integer :: l, il, k + integer*8 :: ipoint, inucl, ishell + integer*8 :: ishell_start, ishell_end + integer :: lstart(0:20) + double precision :: x, y, z, r2 + double precision :: cutoff + integer, external :: qmckl_ao_polynomial_transp_vgl_f + + double precision, allocatable :: poly_vgl(:,:) + integer , allocatable :: powers(:,:) + + allocate(poly_vgl(ao_num,5), powers(3,ao_num)) + + ! Pre-computed data + do l=0,20 + lstart(l) = l*(l+1)*(l+2)/6 +1 + end do + + info = QMCKL_SUCCESS + + ! Don't compute polynomials when the radial part is zero. + cutoff = -dlog(1.d-12) + + do ipoint = 1, point_num + e_coord(1) = coord(ipoint,1) + e_coord(2) = coord(ipoint,2) + e_coord(3) = coord(ipoint,3) + k=1 + do inucl=1,nucl_num + n_coord(1) = nucl_coord(inucl,1) + n_coord(2) = nucl_coord(inucl,2) + n_coord(3) = nucl_coord(inucl,3) + + ! Test if the point is in the range of the nucleus + x = e_coord(1) - n_coord(1) + y = e_coord(2) - n_coord(2) + z = e_coord(3) - n_coord(3) + + r2 = x*x + z*z + z*z + + if (r2 > cutoff*nucleus_range(inucl)) then + cycle + end if + + ! Compute polynomials + info = qmckl_ao_polynomial_transp_vgl_f(context, e_coord, n_coord, & + nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, & + poly_vgl, int(ao_num,8)) + + ! Loop over shells + ishell_start = nucleus_index(inucl) + 1 + ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) + do ishell = ishell_start, ishell_end + l = shell_ang_mom(ishell) + if (shell_vgl(ishell,1,ipoint) /= 0.d0) then + do il = lstart(l), lstart(l+1)-1 + ! Value + ao_vgl(k,1,ipoint) = & + poly_vgl(il,1) * shell_vgl(ishell,1,ipoint) * ao_factor(k) + + ! Grad_x + ao_vgl(k,2,ipoint) = ( & + poly_vgl(il,2) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(il,1) * shell_vgl(ishell,2,ipoint) & + ) * ao_factor(k) + + ! Grad_y + ao_vgl(k,3,ipoint) = ( & + poly_vgl(il,3) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(il,1) * shell_vgl(ishell,3,ipoint) & + ) * ao_factor(k) + + ! Grad_z + ao_vgl(k,4,ipoint) = ( & + poly_vgl(il,4) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(il,1) * shell_vgl(ishell,4,ipoint) & + ) * ao_factor(k) + + ! Lapl_z + ao_vgl(k,5,ipoint) = ( & + poly_vgl(il,5) * shell_vgl(ishell,1,ipoint) + & + poly_vgl(il,1) * shell_vgl(ishell,5,ipoint) + & + 2.d0 * ( & + poly_vgl(il,2) * shell_vgl(ishell,2,ipoint) + & + poly_vgl(il,3) * shell_vgl(ishell,3,ipoint) + & + poly_vgl(il,4) * shell_vgl(ishell,4,ipoint) ) & + ) * ao_factor(k) + k = k+1 + end do + else + do il = lstart(l), lstart(l+1)-1 + ao_vgl(k,1,ipoint) = 0.d0 + ao_vgl(k,2,ipoint) = 0.d0 + ao_vgl(k,3,ipoint) = 0.d0 + ao_vgl(k,4,ipoint) = 0.d0 + ao_vgl(k,5,ipoint) = 0.d0 + k = k+1 + end do + end if + end do + end do + end do + + deallocate(poly_vgl, powers) +end function qmckl_compute_ao_vgl_hpc_f ++
qmckl_exit_code qmckl_compute_ao_vgl ( const qmckl_context context, @@ -2801,9 +3091,10 @@ For example, with a=0, b=2 and c=1 the string is "yyz"
typedef struct qmckl_determinant_struct { @@ -598,8 +598,8 @@ this mechanism.
When all the data for the slater determinants have been provided, the following
@@ -613,8 +613,8 @@ function returns true.
To set the basis set, all the following functions need to be @@ -638,24 +638,24 @@ computed to accelerate the calculations.
qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double* const det_vgl_alpha); @@ -665,14 +665,14 @@ computed to accelerate the calculations.
context is not QMCKL_NULL_CONTEXT[n][3] in C and (3,n) in Fortra
qmckl_exit_code qmckl_distance ( @@ -834,8 +834,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
integer function qmckl_distance_f(context, transa, transb, m, n, & @@ -1002,8 +1002,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
This function is more efficient when A and B are transposed.
@@ -1013,12 +1013,12 @@ This function is more efficient when A and B are trans
qmckl_distance_rescaledqmckl_distance_rescaled
qmckl_distance_rescaled computes the matrix of the rescaled distances between all
@@ -1036,7 +1036,7 @@ If the input array is normal ('N'), the xyz coordinates are in
the leading dimension: [n][3] in C and (3,n) in Fortran.
context is not QMCKL_NULL_CONTEXT[n][3] in C and (3,n) in Fortra
qmckl_exit_code qmckl_distance_rescaled ( @@ -1185,8 +1185,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, & @@ -1356,8 +1356,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
This function is more efficient when A and B are transposed.
@@ -1366,12 +1366,12 @@ This function is more efficient when A and B are trans
qmckl_distance_rescaled_deriv_eqmckl_distance_rescaled_deriv_e
qmckl_distance_rescaled_deriv_e computes the matrix of the gradient and laplacian of the
@@ -1438,7 +1438,7 @@ If the input array is normal ('N'), the xyz coordinates are in
the leading dimension: [n][3] in C and (3,n) in Fortran.