The qmckl.h header file installed in the ${prefix}/include directory
@@ -385,12 +385,12 @@ Both files are located in the include/ directory.
In a traditional source code, most of the lines of source files of a program @@ -435,8 +435,8 @@ interactively, in the same spirit as Jupyter notebooks.
For a tutorial on literate programming with org-mode, follow this link. @@ -467,8 +467,8 @@ org-mode.
Most of the codes of the TREX CoE are written in Fortran with some scripts in @@ -516,8 +516,8 @@ For more guidelines on using Fortran to generate a C interface, see
The authors should follow the recommendations of the C99
@@ -535,8 +535,8 @@ Compliance can be checked with cppcheck as:
The proposed API should allow the library to: deal with memory transfers @@ -547,8 +547,8 @@ functions (see below).
To avoid namespace collisions, we use qmckl_ as a prefix for all exported
@@ -573,8 +573,8 @@ form is allowed.
In the C language, the number of bits used by the integer types can change @@ -606,15 +606,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
@@ -717,8 +717,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 @@ -727,14 +727,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. @@ -747,27 +747,27 @@ temporary storage, to simplify the use of accelerators.
The high-level functions should be pure, unless the introduction
of non-purity is justified. All the side effects should be made in
-the context variable.
+the context variable.
The number of bits of precision required for a function 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
+functions will use the precision specified in the context
variable.
Reducing the scaling of an algorithm usually implies also reducing @@ -783,7 +783,7 @@ implemented adapted to different problem sizes.
The following arrays are stored in the context: @@ -387,6 +395,88 @@ The following arrays are stored in the context:
type |
++ | Gaussian ('G') or Slater ('S') |
+
shell_num |
++ | Number of shells | +
prim_num |
++ | Total number of primitives | +
nucleus_index |
+[nucl_num] |
+Index of the first shell of each nucleus | +
shell_ang_mom |
+[shell_num] |
+Angular momentum of each shell | +
shell_prim_num |
+[shell_num] |
+Number of primitives in each shell | +
shell_prim_index |
+[shell_num] |
+Address of the first primitive of each shell in the EXPONENT array |
+
shell_factor |
+[shell_num] |
+Normalization factor for each shell | +
exponent |
+[prim_num] |
+Array of exponents | +
coefficient |
+[prim_num] |
+Array of coefficients | +
prim_factor |
+[prim_num] |
+Normalization factors of the primtives | +
+Computed data +
+type |
-- | Gaussian ('G') or Slater ('S') |
-|||
|---|---|---|---|---|---|
shell_num |
-- | Number of shells | -|||
prim_num |
-- | Total number of primitives | -|||
nucleus_index |
+nucleus_prim_index |
[nucl_num] |
-Index of the first shell of each nucleus | +Index of the first primitive for each nucleus | |
shell_ang_mom |
-[shell_num] |
-Angular momentum of each shell | +primitive_vgl |
+[prim_num][5][walk_num][elec_num] |
+Value, gradients, Laplacian of the primitives at electron positions |
shell_prim_num |
-[shell_num] |
-Number of primitives in each shell | -|||
shell_prim_index |
-[shell_num] |
-Address of the first primitive of each shell in the EXPONENT array |
-|||
shell_factor |
-[shell_num] |
-Normalization factor for each shell | -|||
exponent |
-[prim_num] |
-Array of exponents | -|||
coefficient |
-[prim_num] |
-Array of coefficients | -|||
prim_factor |
-[prim_num] |
-Normalization factors of the primtives | +primitive_vgl_date |
+uint64_t |
+Late modification date of Value, gradients, Laplacian of the primitives at electron positions |
| qmcklcontext | -context | -in | -Global state | -
| int64t | -n | -in | -Number of values | -
| double | -X[n] | -in | -Array containing the input values | -
| int32t | -LMAX[n] | -in | -Array containing the maximum power for each value | -
| double | -P[n][ldp] | -out | -Array containing all the powers of X |
-
| int64t | -ldp | -in | -Leading dimension of array P |
-
context is not QMCKL_NULL_CONTEXTn > 0X is allocated with at least \(n \times 8\) bytesLMAX is allocated with at least \(n \times 4\) bytesP is allocated with at least \(n \times \max_i \text{LMAX}_i \times 8\) bytesLDP >= \(\max_i\) LMAX[i]qmckl_exit_code qmckl_ao_power ( - const qmckl_context context, - const int64_t n, - const double* X, - const int32_t* LMAX, - double* const P, - const int64_t ldp ); --
integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info) - use qmckl - implicit none - integer*8 , intent(in) :: context - integer*8 , intent(in) :: n - real*8 , intent(in) :: X(n) - integer , intent(in) :: LMAX(n) - real*8 , intent(out) :: P(ldp,n) - integer*8 , intent(in) :: ldp - - integer*8 :: i,k - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (n <= ldp) then - info = QMCKL_INVALID_ARG_2 - return - endif - - k = MAXVAL(LMAX) - if (LDP < k) then - info = QMCKL_INVALID_ARG_6 - return - endif - - if (k <= 0) then - info = QMCKL_INVALID_ARG_4 - return - endif - - do i=1,n - P(1,i) = X(i) - do k=2,LMAX(i) - P(k,i) = P(k-1,i) * X(i) - end do - end do - -end function qmckl_ao_power_f --
integer(c_int32_t) function test_qmckl_ao_power(context) bind(C) - use qmckl - implicit none - - integer(qmckl_context), intent(in), value :: context - - integer*8 :: n, LDP - integer, allocatable :: LMAX(:) - double precision, allocatable :: X(:), P(:,:) - integer*8 :: i,j - double precision :: epsilon - - epsilon = qmckl_get_numprec_epsilon(context) - - n = 100; - LDP = 10; - - allocate(X(n), P(LDP,n), LMAX(n)) - - do j=1,n - X(j) = -5.d0 + 0.1d0 * dble(j) - LMAX(j) = 1 + int(mod(j, 5),4) - end do - - test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP) - if (test_qmckl_ao_power /= QMCKL_SUCCESS) return - - test_qmckl_ao_power = QMCKL_FAILURE - - do j=1,n - do i=1,LMAX(j) - if ( X(j)**i == 0.d0 ) then - if ( P(i,j) /= 0.d0) return - else - if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return - end if - end do - end do - - test_qmckl_ao_power = QMCKL_SUCCESS - deallocate(X,P,LMAX) -end function test_qmckl_ao_power --
-A polynomial is centered on a nucleus \(\mathbf{R}_i\) -
- --\[ - P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c - \] -
- --The gradients with respect to electron coordinates are -
- -\begin{eqnarray*} -\frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & - = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\ -\frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & - = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\ -\frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & - = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\ -\end{eqnarray*} - --and the Laplacian is -
- -\begin{eqnarray*} -\left( \frac{\partial }{\partial x^2} + - \frac{\partial }{\partial y^2} + - \frac{\partial }{\partial z^2} \right) P_l - \left(\mathbf{r},\mathbf{R}_i \right) & = & - a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\ - && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\ - && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}. -\end{eqnarray*} - -
-qmckl_ao_polynomial_vgl computes the values, gradients and
-Laplacians at a given point in space, of all polynomials with an
-angular momentum up to lmax.
-
| qmcklcontext | -context | -in | -Global state | -
| double | -X[3] | -in | -Array containing the coordinates of the points | -
| double | -R[3] | -in | -Array containing the x,y,z coordinates of the center | -
| int32t | -lmax | -in | -Maximum angular momentum | -
| int64t | -n | -inout | -Number of computed polynomials | -
| int32t | -L[n][ldl] | -out | -Contains a,b,c for all n results |
-
| int64t | -ldl | -in | -Leading dimension of L |
-
| double | -VGL[n][ldv] | -out | -Value, gradients and Laplacian of the polynomials | -
| int64t | -ldv | -in | -Leading dimension of array VGL |
-
context is not QMCKL_NULL_CONTEXTn > 0lmax >= 0ldl >= 3ldv >= 5X is allocated with at least \(3 \times 8\) bytesR is allocated with at least \(3 \times 8\) bytesn >= (lmax+1)(lmax+2)(lmax+3)/6L is allocated with at least \(3 \times n \times 4\) bytesVGL is allocated with at least \(5 \times n \times 8\) bytesn should be equal to (lmax+1)(lmax+2)(lmax+3)/6ll, alphabetical order of the
-string made by a*"x" + b*"y" + c*"z" (in Python notation).
-For example, with a=0, b=2 and c=1 the string is "yyz"qmckl_exit_code qmckl_ao_polynomial_vgl ( - const qmckl_context context, - const double* X, - const double* R, - const int32_t lmax, - int64_t* n, - int32_t* const L, - const int64_t ldl, - double* const VGL, - const int64_t ldv ); --
integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info) - use qmckl - implicit none - integer*8 , intent(in) :: context - real*8 , intent(in) :: X(3), R(3) - integer , intent(in) :: lmax - integer*8 , intent(out) :: n - integer , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6) - integer*8 , intent(in) :: ldl - real*8 , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6) - integer*8 , intent(in) :: ldv - - integer*8 :: i,j - integer :: a,b,c,d - real*8 :: Y(3) - integer :: lmax_array(3) - real*8 :: pows(-2:lmax,3) - integer, external :: qmckl_ao_power_f - double precision :: xy, yz, xz - double precision :: da, db, dc, dd - - info = 0 - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (lmax < 0) then - info = QMCKL_INVALID_ARG_4 - return - endif - - if (ldl < 3) then - info = QMCKL_INVALID_ARG_7 - return - endif - - if (ldv < 5) then - info = QMCKL_INVALID_ARG_9 - return - endif - - - do i=1,3 - Y(i) = X(i) - R(i) - end do - - lmax_array(1:3) = lmax - if (lmax == 0) then - VGL(1,1) = 1.d0 - vgL(2:5,1) = 0.d0 - l(1:3,1) = 0 - n=1 - else if (lmax > 0) then - pows(-2:0,1:3) = 1.d0 - do i=1,lmax - pows(i,1) = pows(i-1,1) * Y(1) - pows(i,2) = pows(i-1,2) * Y(2) - pows(i,3) = pows(i-1,3) * Y(3) - end do - - VGL(1:5,1:4) = 0.d0 - l (1:3,1:4) = 0 - - VGL(1 ,1 ) = 1.d0 - vgl(1:5,2:4) = 0.d0 - - l (1,2) = 1 - vgl(1,2) = pows(1,1) - vgL(2,2) = 1.d0 - - l (2,3) = 1 - vgl(1,3) = pows(1,2) - vgL(3,3) = 1.d0 - - l (3,4) = 1 - vgl(1,4) = pows(1,3) - vgL(4,4) = 1.d0 - - n=4 - endif - - ! l>=2 - dd = 2.d0 - do d=2,lmax - da = dd - do a=d,0,-1 - db = dd-da - do b=d-a,0,-1 - c = d - a - b - dc = dd - da - db - n = n+1 - - l(1,n) = a - l(2,n) = b - l(3,n) = c - - xy = pows(a,1) * pows(b,2) - yz = pows(b,2) * pows(c,3) - xz = pows(a,1) * pows(c,3) - - vgl(1,n) = xy * pows(c,3) - - xy = dc * xy - xz = db * xz - yz = da * yz - - vgl(2,n) = pows(a-1,1) * yz - vgl(3,n) = pows(b-1,2) * xz - vgl(4,n) = pows(c-1,3) * xy - - vgl(5,n) = & - (da-1.d0) * pows(a-2,1) * yz + & - (db-1.d0) * pows(b-2,2) * xz + & - (dc-1.d0) * pows(c-2,3) * xy - - db = db - 1.d0 - end do - da = da - 1.d0 - end do - dd = dd + 1.d0 - end do - - info = QMCKL_SUCCESS - -end function qmckl_ao_polynomial_vgl_f --
integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) - use qmckl - implicit none - - integer(c_int64_t), intent(in), value :: context - - integer :: lmax, d, i - integer, allocatable :: L(:,:) - integer*8 :: n, ldl, ldv, j - double precision :: X(3), R(3), Y(3) - double precision, allocatable :: VGL(:,:) - double precision :: w - double precision :: epsilon - - epsilon = qmckl_get_numprec_epsilon(context) - - X = (/ 1.1 , 2.2 , 3.3 /) - R = (/ 0.1 , 1.2 , -2.3 /) - Y(:) = X(:) - R(:) - - lmax = 4; - ldl = 3; - ldv = 100; - - d = (lmax+1)*(lmax+2)*(lmax+3)/6 - - allocate (L(ldl,d), VGL(ldv,d)) - - test_qmckl_ao_polynomial_vgl = & - qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) - - if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return - if (n /= d) return - - do j=1,n - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - do i=1,3 - if (L(i,j) < 0) return - end do - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (dabs(1.d0 - VGL(1,j) / (& - Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j) & - )) > epsilon ) return - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (L(1,j) < 1) then - if (VGL(2,j) /= 0.d0) return - else - if (dabs(1.d0 - VGL(2,j) / (& - L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) & - )) > epsilon ) return - end if - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (L(2,j) < 1) then - if (VGL(3,j) /= 0.d0) return - else - if (dabs(1.d0 - VGL(3,j) / (& - L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) & - )) > epsilon ) return - end if - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (L(3,j) < 1) then - if (VGL(4,j) /= 0.d0) return - else - if (dabs(1.d0 - VGL(4,j) / (& - L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) & - )) > epsilon ) return - end if - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - w = 0.d0 - if (L(1,j) > 1) then - w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j) - end if - if (L(2,j) > 1) then - w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j) - end if - if (L(3,j) > 1) then - w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2) - end if - if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return - end do - - test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS - - deallocate(L,VGL) -end function test_qmckl_ao_polynomial_vgl --
int test_qmckl_ao_polynomial_vgl(qmckl_context context); -assert(0 == test_qmckl_ao_polynomial_vgl(context)); --
qmckl_ao_gaussian_vgl computes the values, gradients and
Laplacians at a given point of n Gaussian functions centered at
the same point:
@@ -1508,20 +946,892 @@ Requirements
qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double* const primitive_vgl); +
| qmcklcontext | +context | +in | +Global state | +
| int64t | +primnum | +in | +Number of primitives | +
| int64t | +elecnum | +in | +Number of electrons | +
| int64t | +nuclnum | +in | +Number of nuclei | +
| int64t | +walknum | +in | +Number of walkers | +
| int64t | +nucleusprimindex[nuclnum] | +in | +Index of the 1st primitive of each nucleus | +
| double | +eleccoord[walknum][3][elecnum] | +in | +Electron coordinates | +
| double | +nuclcoord[3][elecnum] | +in | +Nuclear coordinates | +
| double | +expo[primnum] | +in | +Exponents of the primitives | +
| double | +primitivevgl[primnum][5][walknum][elecnum] | +out | +Value, gradients and Laplacian of the primitives | +
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, & + prim_num, elec_num, nucl_num, walk_num, & + nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: prim_num + integer*8 , intent(in) :: nucl_num + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: walk_num + integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1) + double precision , intent(in) :: elec_coord(elec_num,3,walk_num) + double precision , intent(in) :: nucl_coord(nucl_num,3) + double precision , intent(in) :: expo(prim_num) + double precision , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num) + + integer*8 :: inucl, iprim, iwalk, ielec + double precision :: x, y, z, two_a, ar2, r2, v + double precision :: r2_cut(elec_num,walk_num) + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (prim_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_6 + return + endif + + do inucl=1,nucl_num + ! C is zero-based, so shift bounds by one + do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1) + do iwalk = 1, walk_num + do ielec = 1, elec_num + x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1) + y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2) + z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3) + + r2 = x*x + y*y + z*z + ar2 = expo(iprim)*r2 + v = dexp(-ar2) + two_a = -2.d0 * expo(iprim) * v + + primitive_vgl(ielec, iwalk, 1, iprim) = v + primitive_vgl(ielec, iwalk, 2, iprim) = two_a * x + primitive_vgl(ielec, iwalk, 3, iprim) = two_a * y + primitive_vgl(ielec, iwalk, 4, iprim) = two_a * z + primitive_vgl(ielec, iwalk, 5, iprim) = two_a * (3.d0 - 2.d0*ar2) + + end do + end do + end do + end do + +end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f ++
#define walk_num chbrclf_walk_num +#define elec_num chbrclf_elec_num +#define prim_num chbrclf_prim_num + +int64_t elec_up_num = chbrclf_elec_up_num; +int64_t elec_dn_num = chbrclf_elec_dn_num; +double* elec_coord = &(chbrclf_elec_coord[0][0][0]); + +rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_electron_walk_num (context, walk_num); +assert (rc == QMCKL_SUCCESS); + +assert(qmckl_electron_provided(context)); + +rc = qmckl_set_electron_coord (context, 'N', elec_coord); +assert(rc == QMCKL_SUCCESS); + + + +double prim_vgl[prim_num][5][walk_num][elec_num]; + +rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +assert( fabs(prim_vgl[7][0][0][26] - ( 1.0501570432064878E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[7][1][0][26] - (-7.5014974095310560E-004)) < 1.e-14 ); +assert( fabs(prim_vgl[7][2][0][26] - (-3.8250692897610380E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[7][3][0][26] - ( 3.4950559194080275E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[7][4][0][26] - ( 2.0392163767356572E-002)) < 1.e-14 ); + +assert( fabs(prim_vgl[39][0][1][15] - ( 1.0825844173157661E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[39][1][1][15] - ( 2.3774237611651531E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[39][2][1][15] - ( 2.1423191526963063E-003)) < 1.e-14 ); +assert( fabs(prim_vgl[39][3][1][15] - ( 4.3312003523048492E-004)) < 1.e-14 ); +assert( fabs(prim_vgl[39][4][1][15] - ( 7.5174404780004771E-003)) < 1.e-14 ); + + + + ++
+The qmckl_ao_power function computes all the powers of the n
+input data up to the given maximum value given in input for each of
+the \(n\) points:
+
+\[ P_{ik} = X_i^k \] +
+ +| qmcklcontext | +context | +in | +Global state | +
| int64t | +n | +in | +Number of values | +
| double | +X[n] | +in | +Array containing the input values | +
| int32t | +LMAX[n] | +in | +Array containing the maximum power for each value | +
| double | +P[n][ldp] | +out | +Array containing all the powers of X |
+
| int64t | +ldp | +in | +Leading dimension of array P |
+
context is not QMCKL_NULL_CONTEXTn > 0X is allocated with at least \(n \times 8\) bytesLMAX is allocated with at least \(n \times 4\) bytesP is allocated with at least \(n \times \max_i \text{LMAX}_i \times 8\) bytesLDP >= \(\max_i\) LMAX[i]qmckl_exit_code qmckl_ao_power ( + const qmckl_context context, + const int64_t n, + const double* X, + const int32_t* LMAX, + double* const P, + const int64_t ldp ); ++
integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info) + use qmckl + implicit none + integer*8 , intent(in) :: context + integer*8 , intent(in) :: n + real*8 , intent(in) :: X(n) + integer , intent(in) :: LMAX(n) + real*8 , intent(out) :: P(ldp,n) + integer*8 , intent(in) :: ldp + + integer*8 :: i,k + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (n <= ldp) then + info = QMCKL_INVALID_ARG_2 + return + endif + + k = MAXVAL(LMAX) + if (LDP < k) then + info = QMCKL_INVALID_ARG_6 + return + endif + + if (k <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + do i=1,n + P(1,i) = X(i) + do k=2,LMAX(i) + P(k,i) = P(k-1,i) * X(i) + end do + end do + +end function qmckl_ao_power_f ++
integer(c_int32_t) function test_qmckl_ao_power(context) bind(C) + use qmckl + implicit none + + integer(qmckl_context), intent(in), value :: context + + integer*8 :: n, LDP + integer, allocatable :: LMAX(:) + double precision, allocatable :: X(:), P(:,:) + integer*8 :: i,j + double precision :: epsilon + + epsilon = qmckl_get_numprec_epsilon(context) + + n = 100; + LDP = 10; + + allocate(X(n), P(LDP,n), LMAX(n)) + + do j=1,n + X(j) = -5.d0 + 0.1d0 * dble(j) + LMAX(j) = 1 + int(mod(j, 5),4) + end do + + test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP) + if (test_qmckl_ao_power /= QMCKL_SUCCESS) return + + test_qmckl_ao_power = QMCKL_FAILURE + + do j=1,n + do i=1,LMAX(j) + if ( X(j)**i == 0.d0 ) then + if ( P(i,j) /= 0.d0) return + else + if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return + end if + end do + end do + + test_qmckl_ao_power = QMCKL_SUCCESS + deallocate(X,P,LMAX) +end function test_qmckl_ao_power ++
+A polynomial is centered on a nucleus \(\mathbf{R}_i\) +
+ ++\[ + P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c + \] +
+ ++The gradients with respect to electron coordinates are +
+ +\begin{eqnarray*} +\frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & + = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\ +\frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & + = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\ +\frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & + = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\ +\end{eqnarray*} + ++and the Laplacian is +
+ +\begin{eqnarray*} +\left( \frac{\partial }{\partial x^2} + + \frac{\partial }{\partial y^2} + + \frac{\partial }{\partial z^2} \right) P_l + \left(\mathbf{r},\mathbf{R}_i \right) & = & + a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\ + && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\ + && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}. +\end{eqnarray*} + +
+qmckl_ao_polynomial_vgl computes the values, gradients and
+Laplacians at a given point in space, of all polynomials with an
+angular momentum up to lmax.
+
| qmcklcontext | +context | +in | +Global state | +
| double | +X[3] | +in | +Array containing the coordinates of the points | +
| double | +R[3] | +in | +Array containing the x,y,z coordinates of the center | +
| int32t | +lmax | +in | +Maximum angular momentum | +
| int64t | +n | +inout | +Number of computed polynomials | +
| int32t | +L[n][ldl] | +out | +Contains a,b,c for all n results |
+
| int64t | +ldl | +in | +Leading dimension of L |
+
| double | +VGL[n][ldv] | +out | +Value, gradients and Laplacian of the polynomials | +
| int64t | +ldv | +in | +Leading dimension of array VGL |
+
context is not QMCKL_NULL_CONTEXTn > 0lmax >= 0ldl >= 3ldv >= 5X is allocated with at least \(3 \times 8\) bytesR is allocated with at least \(3 \times 8\) bytesn >= (lmax+1)(lmax+2)(lmax+3)/6L is allocated with at least \(3 \times n \times 4\) bytesVGL is allocated with at least \(5 \times n \times 8\) bytesn should be equal to (lmax+1)(lmax+2)(lmax+3)/6ll, alphabetical order of the
+string made by a*"x" + b*"y" + c*"z" (in Python notation).
+For example, with a=0, b=2 and c=1 the string is "yyz"qmckl_exit_code qmckl_ao_polynomial_vgl ( + const qmckl_context context, + const double* X, + const double* R, + const int32_t lmax, + int64_t* n, + int32_t* const L, + const int64_t ldl, + double* const VGL, + const int64_t ldv ); ++
integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info) + use qmckl + implicit none + integer*8 , intent(in) :: context + real*8 , intent(in) :: X(3), R(3) + integer , intent(in) :: lmax + integer*8 , intent(out) :: n + integer , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6) + integer*8 , intent(in) :: ldl + real*8 , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6) + integer*8 , intent(in) :: ldv + + integer*8 :: i,j + integer :: a,b,c,d + real*8 :: Y(3) + integer :: lmax_array(3) + real*8 :: pows(-2:lmax,3) + integer, external :: qmckl_ao_power_f + double precision :: xy, yz, xz + double precision :: da, db, dc, dd + + info = 0 + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (lmax < 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (ldl < 3) then + info = QMCKL_INVALID_ARG_7 + return + endif + + if (ldv < 5) then + info = QMCKL_INVALID_ARG_9 + return + endif + + + do i=1,3 + Y(i) = X(i) - R(i) + end do + + lmax_array(1:3) = lmax + if (lmax == 0) then + VGL(1,1) = 1.d0 + vgL(2:5,1) = 0.d0 + l(1:3,1) = 0 + n=1 + else if (lmax > 0) then + pows(-2:0,1:3) = 1.d0 + do i=1,lmax + pows(i,1) = pows(i-1,1) * Y(1) + pows(i,2) = pows(i-1,2) * Y(2) + pows(i,3) = pows(i-1,3) * Y(3) + end do + + VGL(1:5,1:4) = 0.d0 + l (1:3,1:4) = 0 + + VGL(1 ,1 ) = 1.d0 + vgl(1:5,2:4) = 0.d0 + + l (1,2) = 1 + vgl(1,2) = pows(1,1) + vgL(2,2) = 1.d0 + + l (2,3) = 1 + vgl(1,3) = pows(1,2) + vgL(3,3) = 1.d0 + + l (3,4) = 1 + vgl(1,4) = pows(1,3) + vgL(4,4) = 1.d0 + + n=4 + endif + + ! l>=2 + dd = 2.d0 + do d=2,lmax + da = dd + do a=d,0,-1 + db = dd-da + do b=d-a,0,-1 + c = d - a - b + dc = dd - da - db + n = n+1 + + l(1,n) = a + l(2,n) = b + l(3,n) = c + + xy = pows(a,1) * pows(b,2) + yz = pows(b,2) * pows(c,3) + xz = pows(a,1) * pows(c,3) + + vgl(1,n) = xy * pows(c,3) + + xy = dc * xy + xz = db * xz + yz = da * yz + + vgl(2,n) = pows(a-1,1) * yz + vgl(3,n) = pows(b-1,2) * xz + vgl(4,n) = pows(c-1,3) * xy + + vgl(5,n) = & + (da-1.d0) * pows(a-2,1) * yz + & + (db-1.d0) * pows(b-2,2) * xz + & + (dc-1.d0) * pows(c-2,3) * xy + + db = db - 1.d0 + end do + da = da - 1.d0 + end do + dd = dd + 1.d0 + end do + + info = QMCKL_SUCCESS + +end function qmckl_ao_polynomial_vgl_f ++
integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) + use qmckl + implicit none + + integer(c_int64_t), intent(in), value :: context + + integer :: lmax, d, i + integer, allocatable :: L(:,:) + integer*8 :: n, ldl, ldv, j + double precision :: X(3), R(3), Y(3) + double precision, allocatable :: VGL(:,:) + double precision :: w + double precision :: epsilon + + epsilon = qmckl_get_numprec_epsilon(context) + + X = (/ 1.1 , 2.2 , 3.3 /) + R = (/ 0.1 , 1.2 , -2.3 /) + Y(:) = X(:) - R(:) + + lmax = 4; + ldl = 3; + ldv = 100; + + d = (lmax+1)*(lmax+2)*(lmax+3)/6 + + allocate (L(ldl,d), VGL(ldv,d)) + + test_qmckl_ao_polynomial_vgl = & + qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) + + if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return + if (n /= d) return + + do j=1,n + test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE + do i=1,3 + if (L(i,j) < 0) return + end do + test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE + if (dabs(1.d0 - VGL(1,j) / (& + Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j) & + )) > epsilon ) return + + test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE + if (L(1,j) < 1) then + if (VGL(2,j) /= 0.d0) return + else + if (dabs(1.d0 - VGL(2,j) / (& + L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) & + )) > epsilon ) return + end if + + test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE + if (L(2,j) < 1) then + if (VGL(3,j) /= 0.d0) return + else + if (dabs(1.d0 - VGL(3,j) / (& + L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) & + )) > epsilon ) return + end if + + test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE + if (L(3,j) < 1) then + if (VGL(4,j) /= 0.d0) return + else + if (dabs(1.d0 - VGL(4,j) / (& + L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) & + )) > epsilon ) return + end if + + test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE + w = 0.d0 + if (L(1,j) > 1) then + w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j) + end if + if (L(2,j) > 1) then + w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j) + end if + if (L(3,j) > 1) then + w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2) + end if + if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return + end do + + test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS + + deallocate(L,VGL) +end function test_qmckl_ao_polynomial_vgl ++
int test_qmckl_ao_polynomial_vgl(qmckl_context context); +assert(0 == test_qmckl_ao_polynomial_vgl(context)); ++
The context variable is a handle for the state of the library,
@@ -338,7 +338,7 @@ A value of QMCKL_NULL_CONTEXT for the context is equivalent to a
typedef int64_t qmckl_context ; +typedef int64_t qmckl_context ; #define QMCKL_NULL_CONTEXT (qmckl_context) 0
ctx is a qmckl_context_struct* pointer.
The context keeps a ``date'' that allows to check which data needs @@ -367,7 +367,7 @@ coordinates are updated.
When a new element is added to the context, the functions -qmcklcontextcreate, qmcklcontextdestroy and qmcklcontextcopy +qmcklcontextcreate, qmcklcontextdestroy and qmcklcontextcopy should be updated inorder to make deep copies.
@@ -416,8 +416,8 @@ if the context is valid,QMCKL_NULL_CONTEXT otherwise.
To create a new context, qmckl_context_create() should be used.
@@ -460,17 +460,26 @@ To create a new context, qmckl_context_create() should be used.
}
/* Initialize data */
- ctx->tag = VALID_TAG;
+ {
+ ctx->tag = VALID_TAG;
- const qmckl_context context = (const qmckl_context) ctx;
- assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT );
+ const qmckl_context context = (const qmckl_context) ctx;
+ assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT );
- ctx->numprec.precision = QMCKL_DEFAULT_PRECISION;
- ctx->numprec.range = QMCKL_DEFAULT_RANGE;
+ qmckl_exit_code rc;
- ctx->ao_basis.uninitialized = (1 << 12) - 1;
- ctx->nucleus.uninitialized = (1 << 4) - 1;
- ctx->electron.uninitialized = (1 << 3) - 1;
+ ctx->numprec.precision = QMCKL_DEFAULT_PRECISION;
+ ctx->numprec.range = QMCKL_DEFAULT_RANGE;
+
+ rc = qmckl_init_electron(context);
+ assert (rc == QMCKL_SUCCESS);
+
+ rc = qmckl_init_nucleus(context);
+ assert (rc == QMCKL_SUCCESS);
+
+ rc = qmckl_init_ao_basis(context);
+ assert (rc == QMCKL_SUCCESS);
+ }
/* Allocate qmckl_memory_struct */
{
@@ -493,8 +502,8 @@ To create a new context, qmckl_context_create() should be used.
For thread safety, the context may be locked/unlocked. The lock is @@ -539,8 +548,8 @@ number of times the thread has locked it is saved in the
qmckl_context_copy makes a deep copy of a context. It returns
@@ -588,8 +597,8 @@ number of times the thread has locked it is saved in the
The context is destroyed with qmckl_context_destroy, leaving the ancestors untouched.
@@ -643,7 +652,7 @@ It frees the context, and returns the previous context.
qmckl_distance_sqqmckl_distance_sq
qmckl_distance_sq computes the matrix of the squared distances
@@ -391,7 +391,7 @@ between all pairs of points in two sets, one point within each set:
\]
context is not QMCKL_NULL_CONTEXTqmckl_exit_code qmckl_distance_rescaled ( @@ -524,8 +524,8 @@ between all pairs of points in two sets, one point within each set:
integer function qmckl_distance_sq_f(context, transa, transb, m, n, & @@ -660,8 +660,8 @@ between all pairs of points in two sets, one point within each set:
This function is more efficient when A and B are
@@ -671,12 +671,12 @@ transposed.
qmckl_distanceqmckl_distance
qmckl_distance computes the matrix of the distances between all
@@ -694,7 +694,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 ( @@ -827,8 +827,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
integer function qmckl_distance_f(context, transa, transb, m, n, & @@ -995,8 +995,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
This function is more efficient when A and B are transposed.
@@ -1006,12 +1006,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
@@ -1029,7 +1029,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 ( @@ -1170,8 +1170,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, & @@ -1341,8 +1341,8 @@ the leading dimension:[n][3]in C and(3,n)in Fortra
This function is more efficient when A and B are transposed.
@@ -1354,7 +1354,7 @@ This function is more efficient when A and B are trans
The following data stored in the context: @@ -382,123 +382,123 @@ The following data stored in the context:
uninitializedint32_tnumint64_tup_numint64_tdown_numint64_twalk_numint64_trescale_factor_kappa_eedoublerescale_factor_kappa_endoubleprovidedboolelectron is validcoord_newdouble[walk_num][3][num]coord_olddouble[walk_num][3][num]coord_new_dateuint64_tee_distancedouble[walk_num][num][num]ee_distance_dateuint64_ten_distancedouble[walk_num][nucl_num][num]en_distance_dateuint64_tee_distance_rescaleddouble[walk_num][num][num]ee_distance_rescaled_dateuint64_ten_distance_rescaleddouble[walk_num][nucl_num][num]en_distance_rescaled_dateuint64_ttypedef struct qmckl_electron_struct { @@ -522,20 +522,45 @@ The following data stored in the context: int32_t uninitialized; bool provided; } qmckl_electron_struct; +
The uninitialized integer contains one bit set to one for each
-initialization function which has not bee called. It becomes equal
+initialization function which has not been called. It becomes equal
to zero after all initialization functions have been called. The
struct is then initialized and provided == true.
+Some values are initialized by default, and are not concerned by
+this mechanism.
-When all the data relative to electrons have been set, the
-following function returns true.
-
qmckl_exit_code qmckl_init_electron(qmckl_context context); ++
qmckl_exit_code qmckl_init_electron(qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return false; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + ctx->electron.uninitialized = (1 << 2) - 1; + + /* Default values */ + ctx->electron.rescale_factor_kappa_ee = 1.0; + ctx->electron.rescale_factor_kappa_en = 1.0; + + return QMCKL_SUCCESS; +} ++
bool qmckl_electron_provided (const qmckl_context context);
@@ -544,8 +569,8 @@ following function returns true.
Access functions return QMCKL_SUCCESS when the data has been
@@ -557,12 +582,12 @@ contains the requested data. Otherwise, this variable is untouched.
A walker is a set of electron coordinates that are arguments of
@@ -571,12 +596,12 @@ the wave function. walk_num is the number of walkers.
Returns the current electron coordinates. The pointer is assumed @@ -620,8 +645,8 @@ The order of the indices is:
To set the data relative to the electrons in the context, the @@ -631,10 +656,12 @@ both allocated.
qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); -qmckl_exit_code qmckl_set_kappa (qmckl_context context, const double rescale_factor_kappa_ee, const double rescale_factor_kappa_en); -qmckl_exit_code qmckl_set_electron_walk_num (qmckl_context context, const int64_t walk_num); -qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord); +qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); +qmckl_exit_code qmckl_set_electron_walk_num (qmckl_context context, const int64_t walk_num); +qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const double* coord); + +qmckl_exit_code qmckl_set_electron_rescale_factor_ee (qmckl_context context, const double kappa_ee); +qmckl_exit_code qmckl_set_electron_rescale_factor_en (qmckl_context context, const double kappa_en);
/* Reference input data */ @@ -675,14 +702,13 @@ in the context. int64_t elec_num = chbrclf_elec_num; int64_t elec_up_num = chbrclf_elec_up_num; int64_t elec_dn_num = chbrclf_elec_dn_num; -double rescale_factor_kappa_ee = 1.0; // TODO Get rescale_factor_kappa_ee from chbrclf -double rescale_factor_kappa_en = 1.0; // TODO Get rescale_factor_kappa_en from chbrclf +double rescale_factor_kappa_ee = 2.0; +double rescale_factor_kappa_en = 3.0; double* elec_coord = &(chbrclf_elec_coord[0][0][0]); int64_t nucl_num = chbrclf_nucl_num; double* charge = chbrclf_charge; double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -double nucl_rescale_factor_kappa = 1.0; // TODO Change get rescale_factor_kappa from chbrclf example /* --- */ @@ -717,23 +743,27 @@ rc = qmckl_get_electron_num (context, &n); assert(rc == QMCKL_SUCCESS); assert(n == elec_num); -double k_ee; -double k_en; -rc = qmckl_get_kappa_ee (context, &k_ee); -assert(rc == QMCKL_NOT_PROVIDED); - -rc = qmckl_get_kappa_en (context, &k_en); -assert(rc == QMCKL_NOT_PROVIDED); - -rc = qmckl_set_kappa (context, rescale_factor_kappa_ee, rescale_factor_kappa_en); +double k_ee = 0.; +double k_en = 0.; +rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); assert(rc == QMCKL_SUCCESS); -assert(!qmckl_electron_provided(context)); +assert(k_ee == 1.0); -rc = qmckl_get_kappa_ee (context, &k_ee); +rc = qmckl_get_electron_rescale_factor_en (context, &k_en); +assert(rc == QMCKL_SUCCESS); +assert(k_en == 1.0); + +rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); assert(rc == QMCKL_SUCCESS); assert(k_ee == rescale_factor_kappa_ee); -rc = qmckl_get_kappa_en (context, &k_en); +rc = qmckl_get_electron_rescale_factor_en (context, &k_en); assert(rc == QMCKL_SUCCESS); assert(k_en == rescale_factor_kappa_en); @@ -769,8 +799,8 @@ rc = qmckl_get_electron_coord (context, 'N'
The computed data is stored in the context so that it can be reused @@ -783,12 +813,12 @@ current date is stored.
qmckl_exit_code qmckl_get_electron_ee_distance(qmckl_context context, double* const distance); @@ -797,10 +827,10 @@ current date is stored.
typedef struct qmckl_nucleus_struct { @@ -481,6 +481,7 @@ The following data stored in the context: int32_t uninitialized; bool provided; } qmckl_nucleus_struct; +
uninitialized integer contains one bit set to one for each
initialization function which has not been called. It becomes equal
to zero after all initialization functions have been called. The
struct is then initialized and provided == true.
+Some values are initialized by default, and are not concerned by
+this mechanism.
+
+qmckl_exit_code qmckl_init_nucleus(qmckl_context context); ++
qmckl_exit_code qmckl_init_nucleus(qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return false; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + ctx->nucleus.uninitialized = (1 << 3) - 1; + + /* Default values */ + ctx->nucleus.rescale_factor_kappa = 1.0; + + return QMCKL_SUCCESS; +} ++
When all the data relative to nuclei have been set, the following
@@ -508,8 +538,8 @@ function returns true.
To set the data relative to the nuclei in the context, the @@ -519,8 +549,9 @@ following functions need to be called.
qmckl_exit_code qmckl_set_nucleus_num (qmckl_context context, const int64_t num); qmckl_exit_code qmckl_set_nucleus_charge (qmckl_context context, const double* charge); -qmckl_exit_code qmckl_set_nucleus_kappa (qmckl_context context, const double rescale_factor_kappa); qmckl_exit_code qmckl_set_nucleus_coord (qmckl_context context, const char transp, const double* coord); + +qmckl_exit_code qmckl_set_nucleus_rescale_factor (qmckl_context context, const double rescale_factor_kappa);
const int64_t nucl_num = chbrclf_nucl_num; const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -const double nucl_rescale_factor_kappa = 1.0; // TODO Change get rescale_factor_kappa from chbrclf example +const double nucl_rescale_factor_kappa = 2.0; /* --- */ @@ -572,15 +603,15 @@ rc = qmckl_get_nucleus_num (context, &n); assert(n == nucl_num); double k; -rc = qmckl_get_nucleus_kappa (context, &k); -assert(rc == QMCKL_NOT_PROVIDED); - - -rc = qmckl_set_nucleus_kappa (context, nucl_rescale_factor_kappa); +rc = qmckl_get_nucleus_rescale_factor (context, &k); assert(rc == QMCKL_SUCCESS); -assert(!qmckl_nucleus_provided(context)); +assert(k == 1.0); -rc = qmckl_get_nucleus_kappa (context, &k); + +rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_get_nucleus_rescale_factor (context, &k); assert(rc == QMCKL_SUCCESS); assert(k == nucl_rescale_factor_kappa); @@ -628,8 +659,8 @@ rc = qmckl_get_nucleus_charge(context, nucl_charge2);
The computed data is stored in the context so that it can be reused @@ -642,12 +673,12 @@ current date is stored.
qmckl_exit_code qmckl_get_nucleus_nn_distance(qmckl_context context, double* distance); @@ -656,10 +687,10 @@ current date is stored.
context is not QMCKL_NULL_CONTEXTqmckl_exit_code qmckl_transpose ( @@ -456,8 +456,8 @@ Transposes a matrix: \(B_{ji} = A_{ij}\)
integer function qmckl_transpose_f(context, m, n, A, LDA, B, LDB) & @@ -516,7 +516,7 @@ Transposes a matrix: \(B_{ji} = A_{ij}\)