Merge branch 'main' of github.com:TREX-CoE/qmckl into main

src/Makefile

@@ -1,3 +1,8 @@
+COMPILER=GNU
+#COMPILER=INTEL
+#COMPILER=LLVM
+
+ifeq ($(COMPILER),GNU)
 CC=gcc -g
 CFLAGS=-fPIC -fexceptions -Wall -Werror -Wpedantic -Wextra
 
@@ -5,14 +10,28 @@ FC=gfortran -g
 FFLAGS=-fPIC  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
 
 LIBS=-lgfortran -lm
+endif
 
-#CC=icc -xHost
+ifeq ($(COMPILER),INTEL)
-#CFLAGS=-fPIC -g -O2
+CC=icc -xHost
-#
+CFLAGS=-fPIC -g -O2
-#FC=ifort -xHost
+
-#FFLAGS=-fPIC -g -O2
+FC=ifort -xHost
-#
+FFLAGS=-fPIC -g -O2
-#LIBS=-lm -lifcore -lirc
+
+LIBS=-lm -lifcore -lirc
+endif
+
+#TODO
+ifeq ($(COMPILER),LLVM)
+CC=clang
+CFLAGS=-fPIC -g -O2
+
+FC=flang
+FFLAGS=fPIC -g -O2
+
+LIBS=-lm
+endif
 
 
 export CC CFLAGS FC FFLAGS LIBS

src/qmckl.org

@@ -12,6 +12,7 @@
 #define QMCKL_H
 #include <stdlib.h>
 #include <stdint.h>
+#include <math.h>
     #+END_SRC
 
     #+BEGIN_SRC f90 :tangle qmckl_f.f90

src/qmckl_ao.org

@@ -11,7 +11,6 @@
 
 *** Test                                                           :noexport:
       #+BEGIN_SRC C :tangle test_qmckl_ao.c
-#include <math.h>
 #include "qmckl.h"
 #include "munit.h"
 MunitResult test_qmckl_ao() {
@@ -208,7 +207,7 @@ munit_assert_int(0, ==, test_qmckl_ao_power(context));
       | =n=          | output | Number of computed polynomials                       |
       | =L(ldl,n)=   | output | Contains a,b,c for all =n= results                   |
       | =ldl=        | input  | Leading dimension of =L=                             |
-      | =VGL(ldv,5)= | output | Value, gradients and Laplacian of the polynomials    |
+      | =VGL(ldv,n)= | output | Value, gradients and Laplacian of the polynomials    |
       | =ldv=        | input  | Leading dimension of array =VGL=                     |
     
 ***** Requirements
@@ -217,12 +216,25 @@ munit_assert_int(0, ==, test_qmckl_ao_power(context));
       - =n= > 0
       - =lmax= >= 0
       - =ldl= >= 3
-      - =ldv= >= (=lmax=+1)(=lmax=+2)(=lmax=+3)/6
+      - =ldv= >= 5
       - =X= is allocated with at least $3 \times 8$ bytes
       - =R= is allocated with at least $3 \times 8$ bytes
+      - =n= >= =(lmax+1)(lmax+2)(lmax+3)/6=
       - =L= is allocated with at least $3 \times n \times 4$ bytes
-      - =VGL= is allocated with at least $n \times 5 \times 8$ bytes
+      - =VGL= is allocated with at least $5 \times n \times 8$ bytes
-      - On output, =n= should be equal to (=lmax=+1)(=lmax=+2)(=lmax=+3)/6
+      - On output, =n= should be equal to =(lmax+1)(lmax+2)(lmax+3)/6=
+      - On output, the powers are given in the following order (l=a+b+c):
+        - Increase values of =l=
+        - Within a given value of =l=, 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"
+        
+***** Error codes
+
+      | -1 | Null context        |
+      | -2 | Inconsistent =ldl=  |
+      | -3 | Inconsistent =ldv=  |
+      | -4 | Inconsistent =lmax= |
 
 ***** Header
       #+BEGIN_SRC C :tangle qmckl.h
@@ -243,7 +255,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
   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,5)
+  real*8    , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
   integer*8 , intent(in)  :: ldv
 
   integer*8         :: i,j
@@ -267,7 +279,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
      return
   endif
   
-  if (ldv < (lmax+1)*(lmax+2)*(lmax+3)/6) then
+  if (ldv < 5) then
      info = -3
      return
   endif
@@ -281,30 +293,53 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
   do i=1,3
      Y(i) = X(i) - R(i)
   end do
-  pows(-2:-1,1:3) = 0.d0
+
-  pows(0,1:3) = 1.d0
   lmax_array(1:3) = lmax
-  info = qmckl_ao_power_f(context, 1_8, Y(1), (/lmax/), pows(1,1), size(pows,1,kind=8)) 
+  if (lmax == 0) then
-  if (info /= 0) return
+    VGL(1,1) = 1.d0
-  info = qmckl_ao_power_f(context, 1_8, Y(2), (/lmax/), pows(1,2), size(pows,1,kind=8)) 
+    vgL(2:5,1) = 0.d0
-  if (info /= 0) return
+    l(1:3,1) = 0
-  info = qmckl_ao_power_f(context, 1_8, Y(3), (/lmax/), pows(1,3), size(pows,1,kind=8)) 
+    n=1
-  if (info /= 0) return
+  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,1) = 1.d0
-  vgl(1,2:5) = 0.d0
+    vgl(1:5,2:4) = 0.d0
-  l(1:3,1) = 0
+
-  n=1
+    l(1,2) = 1
-  dd = 1.d0
+    vgl(1,2) = pows(1,1)
-  do d=1,lmax
+    vgL(2,2) = 1.d0
-     da = 0.d0
+
-     do a=0,d
+    l(2,3) = 1
-        db = 0.d0
+    vgl(1,3) = pows(1,2)
-        do b=0,d-a
+    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
@@ -313,33 +348,28 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
            yz = pows(b,2) * pows(c,3)
            xz = pows(a,1) * pows(c,3)
            
-           vgl(n,1) = xy * pows(c,3)
+           vgl(1,n) = xy * pows(c,3)
            
            xy = dc * xy
            xz = db * xz
            yz = da * yz
            
-           vgl(n,2) = pows(a-1,1) * yz
+           vgl(2,n) = pows(a-1,1) * yz
-           vgl(n,3) = pows(b-1,2) * xz
+           vgl(3,n) = pows(b-1,2) * xz
-           vgl(n,4) = pows(c-1,3) * xy
+           vgl(4,n) = pows(c-1,3) * xy
            
-           vgl(n,5) = &
+           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
+           db = db - 1.d0
         end do
-        da = da + 1.d0
+        da = da - 1.d0
      end do
      dd = dd + 1.d0
   end do
 
-  if (n /= (lmax+1)*(lmax+2)*(lmax+3)/6) then
-    info = -5
-    return
-  endif
-
   info = 0
 
 end function qmckl_ao_polynomial_vgl_f
@@ -357,7 +387,7 @@ integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, l
   integer (c_int64_t) , intent(out)        :: n
   integer (c_int32_t) , intent(out)        :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
   integer (c_int64_t) , intent(in) , value :: ldl
-  real    (c_double)  , intent(out)        :: VGL(ldv,5)
+  real    (c_double)  , intent(out)        :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
   integer (c_int64_t) , intent(in) , value :: ldv
 
   integer, external :: qmckl_ao_polynomial_vgl_f
@@ -365,6 +395,7 @@ integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, l
 end function qmckl_ao_polynomial_vgl
       #+END_SRC
 
+***** Fortran interface                                            :noexport:
       #+BEGIN_SRC f90 :tangle qmckl_f.f90
   interface
      integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) &
@@ -377,7 +408,7 @@ end function qmckl_ao_polynomial_vgl
        real    (c_double)  , intent(in)         :: X(3), R(3)
        integer (c_int64_t) , intent(out)        :: n
        integer (c_int32_t) , intent(out)        :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
-       real    (c_double)  , intent(out)        :: VGL(ldv,5)
+       real    (c_double)  , intent(out)        :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
      end function qmckl_ao_polynomial_vgl
   end interface
       #+END_SRC
@@ -407,10 +438,10 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
   n = 0;
   ldl = 3;
   ldv = 100;
-
+  
   d = (lmax+1)*(lmax+2)*(lmax+3)/6
 
-  allocate (L(ldl,100), VGL(ldv,5))
+  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)
@@ -426,33 +457,33 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
         if (L(i,j) < 0) return
      end do
      test_qmckl_ao_polynomial_vgl = -12
-     if (dabs(1.d0 - VGL(j,1) / (&
+     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 = -13
      if (L(1,j) < 1) then
-        if (VGL(j,2) /= 0.d0) return
+        if (VGL(2,j) /= 0.d0) return
      else
-        if (dabs(1.d0 - VGL(j,2) / (&
+        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 = -14
      if (L(2,j) < 1) then
-        if (VGL(j,3) /= 0.d0) return
+        if (VGL(3,j) /= 0.d0) return
      else
-        if (dabs(1.d0 - VGL(j,3) / (&
+        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 = -15
      if (L(3,j) < 1) then
-        if (VGL(j,4) /= 0.d0) return
+        if (VGL(4,j) /= 0.d0) return
      else
-        if (dabs(1.d0 - VGL(j,4) / (&
+        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
@@ -468,7 +499,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
      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(j,5) / w) > epsilon ) return
+     if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
   end do
 
   test_qmckl_ao_polynomial_vgl = 0

src/qmckl_context.org

@@ -779,7 +779,7 @@ int qmckl_context_get_range(const qmckl_context context) {
 ***** TODO Tests                                                   :noexport:
 
 **** =qmckl_context_get_epsilon=
-     Returns $\epsilon = 2 / \log_{10} 2^{n-1}$ where =n= is the precision
+     Returns $\epsilon = 2^{1-n}$ where =n= is the precision
      #+BEGIN_SRC C :comments org :tangle qmckl.h
 double qmckl_context_get_epsilon(const qmckl_context context);
      #+END_SRC
@@ -788,7 +788,7 @@ double qmckl_context_get_epsilon(const qmckl_context context);
       #+BEGIN_SRC C :tangle qmckl_context.c
 double qmckl_context_get_epsilon(const qmckl_context context) {
   const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
-  return 1.0 / ((double) ((int64_t) 1 << (ctx->precision-1)));
+  return pow(2.0,(double) 1-ctx->precision);
 }
       #+END_SRC
 

src/qmckl_distance.org

@@ -9,7 +9,6 @@
 
 **** Headers                                                       :noexport:
      #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
-#include <math.h>
 #include "qmckl.h"
 #include "munit.h"
 MunitResult test_qmckl_distance() {