From 2467214b3a0741bd65f7304ac054adabbb3449e0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 5 Nov 2020 00:46:19 +0100 Subject: [PATCH] Introduced Gaussian basis functions --- src/.gitignore | 1 + src/Makefile | 10 +- src/README.org | 5 + src/qmckl_ao.org | 317 +++++++++++++++++++++++++++++++++++++----- src/qmckl_context.org | 193 +++++++++++++++++++------ 5 files changed, 442 insertions(+), 84 deletions(-) diff --git a/src/.gitignore b/src/.gitignore index 90eb50e..0304f4e 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -2,6 +2,7 @@ *.c *.f90 *.h +*.fh *.html *~ *.so diff --git a/src/Makefile b/src/Makefile index 6cf6be8..66a009a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,8 +1,8 @@ -CC=gcc -CFLAGS=-fPIC -fexceptions -Wall -Werror -Wpedantic -Wextra -g +CC=gcc -g +CFLAGS=-fPIC -fexceptions -Wall -Werror -Wpedantic -Wextra -FC=gfortran -FFLAGS=-fPIC -g -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 +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 @@ -33,7 +33,7 @@ doc:$(ORG_SOURCE_FILES) ./create_doc.sh README.org $(ORG_SOURCE_FILES) clean: - rm -f qmckl.h test_qmckl_* test_qmckl.c qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h Makefile.generated libqmckl.so *.html *.fh + rm -f qmckl.h test_qmckl_* test_qmckl.c test_qmckl qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h Makefile.generated libqmckl.so *.html *.fh Makefile.generated: $(ORG_SOURCE_FILES) Makefile create_makefile.sh ./create_makefile.sh $(ORG_SOURCE_FILES) diff --git a/src/README.org b/src/README.org index a895586..3601733 100644 --- a/src/README.org +++ b/src/README.org @@ -80,6 +80,9 @@ rm ${nb}.md Fortran interface files should also be written in a file with a =.fh= extension. + For more guidelines on using Fortran to generate a C interface, see + [[http://fortranwiki.org/fortran/show/Generating+C+Interfaces][this link]] + ** Coding style # TODO: decide on a coding style @@ -124,6 +127,8 @@ rm ${nb}.md 32-bit architectures) - ASCII strings are represented as a pointers to a character arrays and terminated by a zero character (C convention). + + Complex numbers can be represented by an array of 2 floats. # TODO : Link to repositories for bindings To facilitate the use in other languages than C, we provide some diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org index b14f440..5c39a62 100644 --- a/src/qmckl_ao.org +++ b/src/qmckl_ao.org @@ -162,6 +162,7 @@ end function qmckl_ao_powers integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C) use, intrinsic :: iso_c_binding implicit none + include 'qmckl_context.fh' include 'qmckl_ao.fh' integer(c_int64_t), intent(in), value :: context @@ -170,7 +171,10 @@ integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C) integer, allocatable :: LMAX(:) double precision, allocatable :: X(:), P(:,:) integer*8 :: i,j - + double precision :: epsilon + + epsilon = qmckl_context_get_epsilon(context) + n = 100; LDP = 10; @@ -178,7 +182,7 @@ integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C) do j=1,n X(j) = -5.d0 + 0.1d0 * dble(j) - LMAX(j) = 1 + int(mod(j, 9),4) + LMAX(j) = 1 + int(mod(j, 5),4) end do test_qmckl_ao_powers = qmckl_ao_powers(context, n, X, LMAX, P, LDP) @@ -188,7 +192,11 @@ integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C) do j=1,n do i=1,LMAX(j) - if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > 1.d-14 ) return + 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 @@ -202,6 +210,7 @@ int test_qmckl_ao_powers(qmckl_context context); munit_assert_int(0, ==, test_qmckl_ao_powers(context)); #+END_SRC + ** =qmckl_ao_polynomial_vgl= Computes the values, gradients and Laplacians at a given point of @@ -216,21 +225,21 @@ munit_assert_int(0, ==, test_qmckl_ao_powers(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,n)= | output | Value, gradients and Laplacian of the polynomials | + | =VGL(ldv,5)= | output | Value, gradients and Laplacian of the polynomials | | =ldv= | input | Leading dimension of array =VGL= | *** Requirements - =context= is not 0 - =n= > 0 + - =lmax= >= 0 + - =ldl= >= 3 + - =ldv= >= (=lmax=+1)(=lmax=+2)(=lmax=+3)/6 - =X= is allocated with at least $3 \times 8$ bytes - =R= is allocated with at least $3 \times 8$ bytes - - =lmax= >= 0 - - On output, =n= should be equal to (=lmax=+1)(=lmax=+2)(=lmax=+3)/6 - =L= is allocated with at least $3 \times n \times 4$ bytes - - =ldl= >= 3 - - =VGL= is allocated with at least $5 \times n \times 8$ bytes - - =ldv= >= 5 + - =VGL= is allocated with at least $n \times 5 \times 8$ bytes + - On output, =n= should be equal to (=lmax=+1)(=lmax=+2)(=lmax=+3)/6 *** Header #+BEGIN_SRC C :comments link :tangle qmckl_ao.h @@ -238,7 +247,7 @@ qmckl_exit_code qmckl_ao_polynomial_vgl(const qmckl_context context, const double *X, const double *R, const int32_t lmax, const int64_t *n, const int32_t *L, const int64_t ldl, - const double *VGL, const int64_t ldv); + const double *VGL, const int64_t ldv); #+END_SRC *** Source @@ -251,7 +260,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,(lmax+1)*(lmax+2)*(lmax+3)/6) + real*8 , intent(out) :: VGL(ldv,5) integer*8 , intent(in) :: ldv integer*8 :: i,j @@ -270,18 +279,21 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, return endif - n = (lmax+1)*(lmax+2)*(lmax+3)/6 - if (ldl < 3) then info = -2 return endif - if (ldv < 5) then + if (ldv < (lmax+1)*(lmax+2)*(lmax+3)/6) then info = -3 return endif + if (lmax <= 0) then + info = -4 + return + endif + do i=1,3 Y(i) = X(i) - R(i) @@ -297,10 +309,10 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, if (info /= 0) return + vgl(1,1) = 1.d0 + vgl(1,2:5) = 0.d0 + l(1:3,1) = 0 n=1 - vgl(1:5,1:n) = 0.d0 - l(1:3,n) = 0 - vgl(1,n) = 1.d0 dd = 1.d0 do d=1,lmax da = 0.d0 @@ -318,17 +330,17 @@ 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(1,n) = xy * pows(c,3) + vgl(n,1) = 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(n,2) = pows(a-1,1) * yz + vgl(n,3) = pows(b-1,2) * xz + vgl(n,4) = pows(c-1,3) * xy - vgl(5,n) = & + vgl(n,5) = & (da-1.d0) * pows(a-2,1) * yz + & (db-1.d0) * pows(b-2,2) * xz + & (dc-1.d0) * pows(c-2,3) * xy @@ -340,6 +352,13 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, 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 #+END_SRC @@ -355,7 +374,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,(lmax+1)*(lmax+2)*(lmax+3)/6) + real (c_double) , intent(out) :: VGL(ldv,5) integer (c_int64_t) , intent(in) , value :: ldv integer, external :: qmckl_ao_polynomial_vgl_f @@ -375,7 +394,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,(lmax+1)*(lmax+2)*(lmax+3)/6) + real (c_double) , intent(out) :: VGL(ldv,5) end function qmckl_ao_polynomial_vgl end interface #+END_SRC @@ -384,6 +403,7 @@ end function qmckl_ao_polynomial_vgl integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) use, intrinsic :: iso_c_binding implicit none + include 'qmckl_context.fh' include 'qmckl_ao.fh' integer(c_int64_t), intent(in), value :: context @@ -394,6 +414,9 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) double precision :: X(3), R(3), Y(3) double precision, allocatable :: VGL(:,:) double precision :: w + double precision :: epsilon + + epsilon = qmckl_context_get_epsilon(context) X = (/ 1.1 , 2.2 , 3.3 /) R = (/ 0.1 , 1.2 , -2.3 /) @@ -406,7 +429,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) d = (lmax+1)*(lmax+2)*(lmax+3)/6 - allocate (L(ldl,100), VGL(ldv,100)) + allocate (L(ldl,100), VGL(ldv,5)) test_qmckl_ao_polynomial_vgl = & qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) @@ -417,37 +440,43 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) if (n /= d) return do j=1,n + test_qmckl_ao_polynomial_vgl = -11 do i=1,3 if (L(i,j) < 0) return end do - if (dabs(1.d0 - VGL(1,j) / (& + test_qmckl_ao_polynomial_vgl = -12 + if (dabs(1.d0 - VGL(j,1) / (& Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j) & - )) > 1.d-14 ) return + )) > epsilon ) return + test_qmckl_ao_polynomial_vgl = -13 if (L(1,j) < 1) then - if (VGL(2,j) /= 0.d0) return + if (VGL(j,2) /= 0.d0) return else - if (dabs(1.d0 - VGL(2,j) / (& + if (dabs(1.d0 - VGL(j,2) / (& L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) & - )) > 1.d-14 ) return + )) > epsilon ) return end if + test_qmckl_ao_polynomial_vgl = -14 if (L(2,j) < 1) then - if (VGL(3,j) /= 0.d0) return + if (VGL(j,3) /= 0.d0) return else - if (dabs(1.d0 - VGL(3,j) / (& + if (dabs(1.d0 - VGL(j,3) / (& L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) & - )) > 1.d-14 ) return + )) > epsilon ) return end if + test_qmckl_ao_polynomial_vgl = -15 if (L(3,j) < 1) then - if (VGL(4,j) /= 0.d0) return + if (VGL(j,4) /= 0.d0) return else - if (dabs(1.d0 - VGL(4,j) / (& + if (dabs(1.d0 - VGL(j,4) / (& L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) & - )) > 1.d-14 ) return + )) > epsilon ) return end if + test_qmckl_ao_polynomial_vgl = -16 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) @@ -458,7 +487,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(5,j) / w) > 1.d-14 ) return + if (dabs(1.d0 - VGL(j,5) / w) > epsilon ) return end do test_qmckl_ao_polynomial_vgl = 0 @@ -474,8 +503,220 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context)); #+END_SRC -* TODO Gaussian basis functions +* Gaussian basis functions + +** =qmckl_ao_gaussians_vgl= + + Computes the values, gradients and Laplacians at a given point of + =n= Gaussian functions centered at the same point: + + \[ v_i = exp(-a_i |X-R|^2) \] + \[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \] + \[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \] + \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] + \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \] + +*** Arguments + + | =context= | input | Global state | + | =X(3)= | input | Array containing the coordinates of the points | + | =R(3)= | input | Array containing the x,y,z coordinates of the center | + | =n= | input | Number of computed gaussians | + | =A(n)= | input | Exponents of the Gaussians | + | =VGL(ldv,5)= | output | Value, gradients and Laplacian of the Gaussians | + | =ldv= | input | Leading dimension of array =VGL= | + +*** Requirements + + - =context= is not 0 + - =n= > 0 + - =ldv= >= 5 + - =A(i)= > 0 for all =i= + - =X= is allocated with at least $3 \times 8$ bytes + - =R= is allocated with at least $3 \times 8$ bytes + - =A= is allocated with at least $n \times 8$ bytes + - =VGL= is allocated with at least $n \times 5 \times 8$ bytes + +*** Header + #+BEGIN_SRC C :comments link :tangle qmckl_ao.h +qmckl_exit_code qmckl_ao_gaussians_vgl(const qmckl_context context, + const double *X, const double *R, + const int64_t *n, const int64_t *A, + const double *VGL, const int64_t ldv); + #+END_SRC + +*** Source + #+BEGIN_SRC f90 :comments link :tangle qmckl_ao.f90 +integer function qmckl_ao_gaussians_vgl_f(context, X, R, n, A, VGL, ldv) result(info) + implicit none + integer*8 , intent(in) :: context + real*8 , intent(in) :: X(3), R(3) + integer*8 , intent(in) :: n + real*8 , intent(in) :: A(n) + real*8 , intent(out) :: VGL(ldv,5) + integer*8 , intent(in) :: ldv + + integer*8 :: i,j + real*8 :: Y(3), r2, t, u, v + + info = 0 + + if (context == 0_8) then + info = -1 + return + endif + + if (n <= 0) then + info = -2 + return + endif + + if (ldv < n) then + info = -3 + return + endif + + + do i=1,3 + Y(i) = X(i) - R(i) + end do + r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3) + + do i=1,n + VGL(i,1) = dexp(-A(i) * r2) + end do + + do i=1,n + VGL(i,5) = A(i) * VGL(i,1) + end do + + t = -2.d0 * ( X(1) - R(1) ) + u = -2.d0 * ( X(2) - R(2) ) + v = -2.d0 * ( X(3) - R(3) ) + + do i=1,n + VGL(i,2) = t * VGL(i,5) + VGL(i,3) = u * VGL(i,5) + VGL(i,4) = v * VGL(i,5) + end do + + t = 4.d0 * r2 + do i=1,n + VGL(i,5) = (t * A(i) - 6.d0) * VGL(i,5) + end do + +end function qmckl_ao_gaussians_vgl_f + #+END_SRC + +*** C interface :noexport: + #+BEGIN_SRC f90 :comments link :tangle qmckl_ao.f90 +integer(c_int32_t) function qmckl_ao_gaussians_vgl(context, X, R, n, A, VGL, ldv) & + bind(C) result(info) + use, intrinsic :: iso_c_binding + implicit none + integer (c_int64_t) , intent(in) , value :: context + real (c_double) , intent(in) :: X(3), R(3) + integer (c_int64_t) , intent(in) , value :: n + real (c_double) , intent(in) :: A(n) + real (c_double) , intent(out) :: VGL(ldv,5) + integer (c_int64_t) , intent(in) , value :: ldv + + integer, external :: qmckl_ao_gaussians_vgl_f + info = qmckl_ao_gaussians_vgl_f(context, X, R, n, A, VGL, ldv) +end function qmckl_ao_gaussians_vgl + #+END_SRC + + #+BEGIN_SRC f90 :comments link :tangle qmckl_ao.fh + interface + integer(c_int32_t) function qmckl_ao_gaussians_vgl(context, X, R, n, A, VGL, ldv) & + bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: ldv + integer (c_int64_t) , intent(in) , value :: n + real (c_double) , intent(in) :: X(3), R(3), A(n) + real (c_double) , intent(out) :: VGL(ldv,5) + end function qmckl_ao_gaussians_vgl + end interface + #+END_SRC +*** Test :noexport: + #+BEGIN_SRC f90 :comments link :tangle test_qmckl_ao_f.f90 +integer(c_int32_t) function test_qmckl_ao_gaussians_vgl(context) bind(C) + use, intrinsic :: iso_c_binding + implicit none + include 'qmckl_context.fh' + include 'qmckl_ao.fh' + + integer(c_int64_t), intent(in), value :: context + + integer*8 :: n, ldv, j, i + double precision :: X(3), R(3), Y(3), r2 + double precision, allocatable :: VGL(:,:), A(:) + double precision :: epsilon + + epsilon = qmckl_context_get_epsilon(context) + + X = (/ 1.1 , 2.2 , 3.3 /) + R = (/ 0.1 , 1.2 , -2.3 /) + Y(:) = X(:) - R(:) + r2 = Y(1)**2 + Y(2)**2 + Y(3)**2 + + n = 10; + ldv = 100; + + allocate (A(n), VGL(ldv,5)) + do i=1,n + A(i) = 0.0013 * dble(ishft(1,i)) + end do + + + test_qmckl_ao_gaussians_vgl = & + qmckl_ao_gaussians_vgl(context, X, R, n, A, VGL, ldv) + if (test_qmckl_ao_gaussians_vgl /= 0) return + + test_qmckl_ao_gaussians_vgl = -1 + + do i=1,n + test_qmckl_ao_gaussians_vgl = -11 + if (dabs(1.d0 - VGL(i,1) / (& + dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussians_vgl = -12 + if (dabs(1.d0 - VGL(i,2) / (& + -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussians_vgl = -13 + if (dabs(1.d0 - VGL(i,3) / (& + -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussians_vgl = -14 + if (dabs(1.d0 - VGL(i,4) / (& + -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) & + )) > epsilon ) return + + test_qmckl_ao_gaussians_vgl = -15 + if (dabs(1.d0 - VGL(i,5) / (& + A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) & + )) > epsilon ) return + end do + + test_qmckl_ao_gaussians_vgl = 0 + + deallocate(VGL) +end function test_qmckl_ao_gaussians_vgl + #+END_SRC + + #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c +int test_qmckl_ao_gaussians_vgl(qmckl_context context); +munit_assert_int(0, ==, test_qmckl_ao_gaussians_vgl(context)); + #+END_SRC + #+END_SRC + + * TODO Slater basis functions * End of files :noexport: diff --git a/src/qmckl_context.org b/src/qmckl_context.org index ad35918..635da3e 100644 --- a/src/qmckl_context.org +++ b/src/qmckl_context.org @@ -83,11 +83,11 @@ qmckl_context qmckl_context_check(const qmckl_context context) ; *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_context qmckl_context_check(const qmckl_context context) { - qmckl_context_struct * ctx; if (context == (qmckl_context) 0) return (qmckl_context) 0; - ctx = (qmckl_context_struct*) context; + const qmckl_context_struct * ctx = (qmckl_context_struct*) context; + if (ctx->tag != VALID_TAG) return (qmckl_context) 0; return context; @@ -109,9 +109,8 @@ qmckl_context qmckl_context_create(); #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_context qmckl_context_create() { - qmckl_context_struct* context; - - context = (qmckl_context_struct*) qmckl_malloc ((qmckl_context) 0, sizeof(qmckl_context_struct)); + qmckl_context_struct* context = + (qmckl_context_struct*) qmckl_malloc ((qmckl_context) 0, sizeof(qmckl_context_struct)); if (context == NULL) { return (qmckl_context) 0; } @@ -125,6 +124,15 @@ qmckl_context qmckl_context_create() { } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int64_t) function qmckl_context_create() bind(C) + use, intrinsic :: iso_c_binding + end function qmckl_context_create + end interface + #+END_SRC + *** Test :noexport: #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c context = qmckl_context_create(); @@ -149,23 +157,20 @@ qmckl_context qmckl_context_copy(const qmckl_context context); #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_context qmckl_context_copy(const qmckl_context context) { - qmckl_context_struct* old_context; - qmckl_context_struct* new_context; - qmckl_context checked_context; - - checked_context = qmckl_context_check(context); + const qmckl_context checked_context = qmckl_context_check(context); if (checked_context == (qmckl_context) 0) { return (qmckl_context) 0; } - new_context = (qmckl_context_struct*) qmckl_malloc (context, sizeof(qmckl_context_struct)); + qmckl_context_struct* old_context = (qmckl_context_struct*) checked_context; + + qmckl_context_struct* new_context = + (qmckl_context_struct*) qmckl_malloc (context, sizeof(qmckl_context_struct)); if (new_context == NULL) { return (qmckl_context) 0; } - old_context = (qmckl_context_struct*) checked_context; - new_context->prev = old_context; new_context->precision = old_context->precision; new_context->range = old_context->range; @@ -176,6 +181,16 @@ qmckl_context qmckl_context_copy(const qmckl_context context) { #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int64_t) function qmckl_context_copy(context) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + end function qmckl_context_copy + end interface + #+END_SRC + *** Test :noexport: #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c new_context = qmckl_context_copy(context); @@ -200,19 +215,26 @@ qmckl_context qmckl_context_previous(const qmckl_context context); #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_context qmckl_context_previous(const qmckl_context context) { - qmckl_context checked_context; - qmckl_context_struct* ctx; - - checked_context = qmckl_context_check(context); + const qmckl_context checked_context = qmckl_context_check(context); if (checked_context == (qmckl_context) 0) { return (qmckl_context) 0; } - ctx = (qmckl_context_struct*) checked_context; + const qmckl_context_struct* ctx = (qmckl_context_struct*) checked_context; return qmckl_context_check((qmckl_context) ctx->prev); } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int64_t) function qmckl_context_previous(context) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + end function qmckl_context_previous + end interface + #+END_SRC + *** Test :noexport: #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c munit_assert_int64(qmckl_context_previous(new_context), !=, (qmckl_context) 0); @@ -236,15 +258,12 @@ qmckl_exit_code qmckl_context_destroy(qmckl_context context); *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c -qmckl_exit_code qmckl_context_destroy(qmckl_context context) { +qmckl_exit_code qmckl_context_destroy(const qmckl_context context) { - qmckl_context_struct* ctx; - qmckl_context checked_context; - - checked_context = qmckl_context_check(context); + const qmckl_context checked_context = qmckl_context_check(context); if (checked_context == (qmckl_context) 0) return QMCKL_FAILURE; - ctx = (qmckl_context_struct*) context; + qmckl_context_struct* ctx = (qmckl_context_struct*) context; if (ctx == NULL) return QMCKL_FAILURE; ctx->tag = INVALID_TAG; @@ -253,6 +272,16 @@ qmckl_exit_code qmckl_context_destroy(qmckl_context context) { } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int32_t) function qmckl_context_destroy(context) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + end function qmckl_context_destroy + end interface + #+END_SRC + *** Test :noexport: #+BEGIN_SRC C :comments link :tangle test_qmckl_context.c munit_assert_int64(qmckl_context_check(new_context), ==, new_context); @@ -283,12 +312,11 @@ qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, cons *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision) { - qmckl_context_struct* ctx; if (precision < 2) return QMCKL_FAILURE; if (precision > 53) return QMCKL_FAILURE; - ctx = (qmckl_context_struct*) context; + qmckl_context_struct* ctx = (qmckl_context_struct*) context; if (ctx == NULL) return QMCKL_FAILURE; ctx->precision = precision; @@ -296,6 +324,17 @@ qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, cons } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int32_t) function qmckl_context_update_precision(context, precision) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + integer (c_int32_t), intent(in), value :: precision + end function qmckl_context_update_precision + end interface + #+END_SRC + *** TODO Tests :noexport: ** =qmckl_context_update_range= *** Header @@ -306,12 +345,11 @@ qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const in *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range) { - qmckl_context_struct* ctx; if (range < 2) return QMCKL_FAILURE; if (range > 11) return QMCKL_FAILURE; - ctx = (qmckl_context_struct*) context; + qmckl_context_struct* ctx = (qmckl_context_struct*) context; if (ctx == NULL) return QMCKL_FAILURE; ctx->range = range; @@ -319,6 +357,17 @@ qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const in } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int32_t) function qmckl_context_update_range(context, range) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + integer (c_int32_t), intent(in), value :: range + end function qmckl_context_update_range + end interface + #+END_SRC + *** TODO Tests :noexport: ** =qmckl_context_set_precision= *** Header @@ -329,9 +378,7 @@ qmckl_context qmckl_context_set_precision(const qmckl_context context, const int *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision) { - qmckl_context new_context; - - new_context = qmckl_context_copy(context); + qmckl_context new_context = qmckl_context_copy(context); if (new_context == 0) return 0; if (qmckl_context_update_precision(context, precision) == QMCKL_FAILURE) return 0; @@ -340,6 +387,17 @@ qmckl_context qmckl_context_set_precision(const qmckl_context context, const int } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int32_t) function qmckl_context_set_precision(context, precision) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + integer (c_int32_t), intent(in), value :: precision + end function qmckl_context_set_precision + end interface + #+END_SRC + *** TODO Tests :noexport: ** =qmckl_context_set_range= *** Header @@ -350,9 +408,7 @@ qmckl_context qmckl_context_set_range(const qmckl_context context, const int ran *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c qmckl_context qmckl_context_set_range(const qmckl_context context, const int range) { - qmckl_context new_context; - - new_context = qmckl_context_copy(context); + qmckl_context new_context = qmckl_context_copy(context); if (new_context == 0) return 0; if (qmckl_context_update_range(context, range) == QMCKL_FAILURE) return 0; @@ -361,42 +417,97 @@ qmckl_context qmckl_context_set_range(const qmckl_context context, const int ran } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int32_t) function qmckl_context_set_range(context, range) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + integer (c_int32_t), intent(in), value :: range + end function qmckl_context_set_range + end interface + #+END_SRC + *** TODO Tests :noexport: ** =qmckl_context_get_precision= *** Header #+BEGIN_SRC C :comments link :tangle qmckl_context.h -int qmckl_context_get_precision(const qmckl_context context); +int32_t qmckl_context_get_precision(const qmckl_context context); #+END_SRC *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c int qmckl_context_get_precision(const qmckl_context context) { - qmckl_context_struct* ctx; - ctx = (qmckl_context_struct*) context; + const qmckl_context_struct* ctx = (qmckl_context_struct*) context; return ctx->precision; } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int32_t) function qmckl_context_get_precision(context) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + end function qmckl_context_get_precision + end interface + #+END_SRC + *** TODO Tests :noexport: ** =qmckl_context_get_range= *** Header #+BEGIN_SRC C :comments link :tangle qmckl_context.h -int qmckl_context_get_range(const qmckl_context context); +int32_t qmckl_context_get_range(const qmckl_context context); #+END_SRC *** Source #+BEGIN_SRC C :comments link :tangle qmckl_context.c int qmckl_context_get_range(const qmckl_context context) { - qmckl_context_struct* ctx; - ctx = (qmckl_context_struct*) context; + const qmckl_context_struct* ctx = (qmckl_context_struct*) context; return ctx->range; } #+END_SRC +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + integer (c_int32_t) function qmckl_context_get_range(context) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + end function qmckl_context_get_range + end interface + #+END_SRC + +*** TODO Tests :noexport: + +** =qmckl_context_get_epsilon= + Returns $\epsilon = 2 / \log_{10} 2^{n-1}$ where =n= is the precision +*** Header + #+BEGIN_SRC C :comments link :tangle qmckl_context.h +double qmckl_context_get_epsilon(const qmckl_context context); + #+END_SRC + +*** Source + #+BEGIN_SRC C :comments link :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))); +} + #+END_SRC + +*** Fortran interface + #+BEGIN_SRC f90 :comments link :tangle qmckl_context.fh + interface + real (c_double) function qmckl_context_get_epsilon(context) bind(C) + use, intrinsic :: iso_c_binding + integer (c_int64_t), intent(in), value :: context + end function qmckl_context_get_epsilon + end interface + #+END_SRC + *** TODO Tests :noexport: - * Info about the molecular system ** TODO =qmckl_context_set_nucl_coord=