From a41c67b94fe46cea7b457fb4da07d534f9488adf Mon Sep 17 00:00:00 2001 From: v1j4y Date: Tue, 7 Sep 2021 16:36:26 +0200 Subject: [PATCH 01/23] Working on qmckl_dgemm. --- org/qmckl_ao.org | 2 + org/qmckl_blas.org | 287 ++++++++++++++++++++++++++++++++++++++++++ org/table_of_contents | 1 + tools/lib.org | 2 + 4 files changed, 292 insertions(+) create mode 100644 org/qmckl_blas.org diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index f78350c..011c015 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -135,6 +135,7 @@ int main() { For H_2 with the following basis set, + #+NAME: basis #+BEGIN_EXAMPLE HYDROGEN S 5 @@ -157,6 +158,7 @@ D 1 we have: + #+NAME: params #+BEGIN_EXAMPLE type = 'G' shell_num = 12 diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org new file mode 100644 index 0000000..0869a93 --- /dev/null +++ b/org/qmckl_blas.org @@ -0,0 +1,287 @@ +#+TITLE: BLAS functions +#+SETUPFILE: ../tools/theme.setup +#+INCLUDE: ../tools/lib.org + +* Headers :noexport: + #+begin_src elisp :noexport :results none +(org-babel-lob-ingest "../tools/lib.org") +#+end_src + + #+begin_src c :comments link :tangle (eval c_test) :noweb yes +#include "qmckl.h" +#include "assert.h" +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif +int main() { + qmckl_context context; + context = qmckl_context_create(); + + #+end_src + +* Matrix operations + +** ~qmckl_dgemm~ + + Matrix multiply l$C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function. + + #+NAME: qmckl_dgemm_args + | qmckl_context | context | in | Global state | + | bool | TransA | in | Number of rows of the input matrix | + | bool | TransB | in | Number of rows of the input matrix | + | int64_t | m | in | Number of rows of the input matrix | + | int64_t | n | in | Number of columns of the input matrix | + | int64_t | k | in | Number of columns of the input matrix | + | double | alpha | in | Number of columns of the input matrix | + | double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ | + | int64_t | lda | in | Leading dimension of array ~A~ | + | double | B[][ldb] | in | Array containing the $n \times m$ matrix $B$ | + | int64_t | ldb | in | Leading dimension of array ~B~ | + | double | beta | in | Array containing the $n \times m$ matrix $B$ | + | double | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ | + | int64_t | ldc | in | Leading dimension of array ~B~ | + +*** Requirements + + - ~context~ is not ~QMCKL_NULL_CONTEXT~ + - ~m > 0~ + - ~n > 0~ + - ~k > 0~ + - ~lda >= m~ + - ~ldb >= n~ + - ~ldc >= n~ + - ~A~ is allocated with at least $m \times k \times 8$ bytes + - ~B~ is allocated with at least $k \times n \times 8$ bytes + - ~C~ is allocated with at least $m \times n \times 8$ bytes + +*** C header + + #+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") + + #+RESULTS: + #+BEGIN_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_dgemm ( + const qmckl_context context, + const bool TransA, + const bool TransB, + const int64_t m, + const int64_t n, + const int64_t k, + const double alpha, + const double* A, + const int64_t lda, + const double* B, + const int64_t ldb, + const double beta, + double* const C, + const int64_t ldc ); + #+END_src + + +*** Source + #+begin_src f90 :tangle (eval f) +integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + logical*8 , intent(in) :: TransA, TransN + integer*8 , intent(in) :: m, n, k + integer*8 , intent(in) :: lda + real*8 , intent(in) :: A(lda,*) + integer*8 , intent(in) :: ldb + real*8 , intent(in) :: B(ldb,*) + integer*8 , intent(in) :: ldc + real*8 , intent(out) :: C(ldc,*) + + integer*8 :: i,j + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (m <= 0_8) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (n <= 0_8) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (k <= 0_8) then + info = QMCKL_INVALID_ARG_6 + return + endif + + if (LDA < m) then + info = QMCKL_INVALID_ARG_9 + return + endif + + if (LDB < n) then + info = QMCKL_INVALID_ARG_10 + return + endif + + if (LDB < n) then + info = QMCKL_INVALID_ARG_13 + return + endif + + C = matmul(A,B) +end function qmckl_dgemm_f + #+end_src + +*** C interface :noexport: + + #+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") + + #+RESULTS: + #+BEGIN_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_dgemm & + (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + logical*8 , intent(in) , value :: TransA + logical*8 , intent(in) , value :: TransB + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + integer (c_int64_t) , intent(in) , value :: k + real (c_double ) , intent(in) , value :: alpha + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) , value :: beta + real (c_double ) , intent(out) :: C(ldc,*) + integer (c_int64_t) , intent(in) , value :: ldc + + integer(c_int32_t), external :: qmckl_dgemm_f + info = qmckl_dgemm_f & + (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) + + end function qmckl_dgemm + #+END_src + + + #+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_dgemm & + (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beat, C, ldc) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + logical*8 (c_bool) , intent(in) , value :: TransA + logical*8 (c_bool) , intent(in) , value :: TransB + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + integer (c_int64_t) , intent(in) , value :: k + real (c_double ) , intent(in) :: alpha + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) :: beta + real (c_double ) , intent(out) :: C(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + + end function qmckl_dgemm + end interface + #+end_src + +*** Test :noexport: + #+begin_src f90 :tangle (eval f_test) +integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) + use qmckl + implicit none + integer(qmckl_context), intent(in), value :: context + + double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:) + integer*8 :: m, n, k, LDA, LDB, LDC + integer*8 :: i,j + logical*8 :: TransA, TransB + double precision :: x + + TransA = .False. + TransB = .False. + m = 5 + k = 4 + n = 6 + LDA = m + LDB = k + LDC = m + + allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n)) + + A = 0.d0 + B = 0.d0 + C = 0.d0 + D = 0.d0 + do j=1,m + do i=1,k + A(i,j) = -10.d0 + dble(i+j) + end do + end do + + do j=1,k + do i=1,n + B(i,j) = -10.d0 + dble(i+j) + end do + end do + + test_qmckl_dgemm = qmckl_dgemm(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) + + if (test_qmckl_dgemm /= QMCKL_SUCCESS) return + + test_qmckl_dgemm = QMCKL_FAILURE + + x = 0.d0 + do j=1,m + do i=l,n + do l=1,k + D(i,j) += D(i,j) + (A(i,k)*B(k,j)) + end do + x = x + (D(i,j) - C(i,j))**2 + end do + end do + + if (dabs(x) <= 1.d-15) then + test_qmckl_dgemm = QMCKL_SUCCESS + endif + + deallocate(A,B) +end function test_qmckl_dgemm + #+end_src + + #+begin_src c :comments link :tangle (eval c_test) +qmckl_exit_code test_qmckl_dgemm(qmckl_context context); +assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); + #+end_src + +* End of files :noexport: + + #+begin_src c :comments link :tangle (eval c_test) + assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); + return 0; +} + + #+end_src + + +# -*- mode: org -*- +# vim: syntax=c diff --git a/org/table_of_contents b/org/table_of_contents index 50a6159..01bb7e5 100644 --- a/org/table_of_contents +++ b/org/table_of_contents @@ -9,4 +9,5 @@ qmckl_ao.org qmckl_jastrow.org qmckl_distance.org qmckl_utils.org +qmckl_blas.org qmckl_tests.org diff --git a/tools/lib.org b/tools/lib.org index d8e642a..3fe9fe2 100644 --- a/tools/lib.org +++ b/tools/lib.org @@ -40,6 +40,7 @@ f_of_c_d = { '' : '' , 'qmckl_context' : 'integer (c_int64_t)' , 'qmckl_exit_code' : 'integer (c_int32_t)' + , 'bool' : 'logical*8' , 'int32_t' : 'integer (c_int32_t)' , 'int64_t' : 'integer (c_int64_t)' , 'float' : 'real (c_float )' @@ -58,6 +59,7 @@ ctypeid_d = { '' : '' , 'real' : 'real(c_float)' , 'real*8' : 'real(c_double)' , 'character' : 'character(c_char)' + , 'character' : 'character(c_char)' } #+END_SRC From 42639bcbff7c1e62befbee2bea2a66f14ec13563 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Tue, 7 Sep 2021 16:43:25 +0200 Subject: [PATCH 02/23] Added a TODO. #32 --- org/qmckl_blas.org | 1 + 1 file changed, 1 insertion(+) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 0869a93..72e9cb5 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -24,6 +24,7 @@ int main() { ** ~qmckl_dgemm~ Matrix multiply l$C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function. + TODO: Add description about the external library dependence. #+NAME: qmckl_dgemm_args | qmckl_context | context | in | Global state | From cf3550b6b7c2b1a7a6bd06bafc1e1ac898e84c1a Mon Sep 17 00:00:00 2001 From: v1j4y Date: Wed, 15 Sep 2021 11:55:45 +0200 Subject: [PATCH 03/23] Fixed some bugs, now compiles. --- org/qmckl_blas.org | 44 +++++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 72e9cb5..fd7e34e 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -41,7 +41,7 @@ int main() { | double | beta | in | Array containing the $n \times m$ matrix $B$ | | double | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ | | int64_t | ldc | in | Leading dimension of array ~B~ | - + *** Requirements - ~context~ is not ~QMCKL_NULL_CONTEXT~ @@ -86,14 +86,15 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, use qmckl implicit none integer(qmckl_context) , intent(in) :: context - logical*8 , intent(in) :: TransA, TransN + logical*8 , intent(in) :: TransA, TransB integer*8 , intent(in) :: m, n, k + real*8 , intent(in) :: alpha, beta integer*8 , intent(in) :: lda - real*8 , intent(in) :: A(lda,*) + real*8 , intent(in) :: A(m,n) integer*8 , intent(in) :: ldb - real*8 , intent(in) :: B(ldb,*) + real*8 , intent(in) :: B(n,k) integer*8 , intent(in) :: ldc - real*8 , intent(out) :: C(ldc,*) + real*8 , intent(out) :: C(m,n) integer*8 :: i,j @@ -158,13 +159,13 @@ end function qmckl_dgemm_f integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: k real (c_double ) , intent(in) , value :: alpha - real (c_double ) , intent(in) :: A(lda,*) integer (c_int64_t) , intent(in) , value :: lda - real (c_double ) , intent(in) :: B(ldb,*) + real (c_double ) , intent(in) :: A(lda,*) integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) :: B(ldb,*) real (c_double ) , intent(in) , value :: beta - real (c_double ) , intent(out) :: C(ldc,*) integer (c_int64_t) , intent(in) , value :: ldc + real (c_double ) , intent(out) :: C(ldc,*) integer(c_int32_t), external :: qmckl_dgemm_f info = qmckl_dgemm_f & @@ -180,26 +181,26 @@ end function qmckl_dgemm_f #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(c_int32_t) function qmckl_dgemm & - (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beat, C, ldc) & + (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) & bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context - logical*8 (c_bool) , intent(in) , value :: TransA - logical*8 (c_bool) , intent(in) , value :: TransB + logical*8 , intent(in) , value :: TransA + logical*8 , intent(in) , value :: TransB integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: k real (c_double ) , intent(in) :: alpha - real (c_double ) , intent(in) :: A(lda,*) integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: ldb real (c_double ) , intent(in) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: ldb real (c_double ) , intent(in) :: beta - real (c_double ) , intent(out) :: C(ldb,*) - integer (c_int64_t) , intent(in) , value :: ldb + integer (c_int64_t) , intent(in) , value :: ldc + real (c_double ) , intent(out) :: C(ldc,*) end function qmckl_dgemm end interface @@ -214,15 +215,15 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:) integer*8 :: m, n, k, LDA, LDB, LDC - integer*8 :: i,j + integer*8 :: i,j,l logical*8 :: TransA, TransB - double precision :: x + double precision :: x, alpha, beta TransA = .False. TransB = .False. - m = 5 - k = 4 - n = 6 + m = 5_8 + k = 4_8 + n = 6_8 LDA = m LDB = k LDC = m @@ -255,11 +256,12 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) do j=1,m do i=l,n do l=1,k - D(i,j) += D(i,j) + (A(i,k)*B(k,j)) + D(i,j) = D(i,j) + A(i,k)*B(k,j) end do x = x + (D(i,j) - C(i,j))**2 end do end do + print *,"DABS(X)= ",dabs(x) if (dabs(x) <= 1.d-15) then test_qmckl_dgemm = QMCKL_SUCCESS From eaede28a73043c2df7d296fa5053ed56f54e265b Mon Sep 17 00:00:00 2001 From: v1j4y Date: Wed, 15 Sep 2021 16:21:42 +0200 Subject: [PATCH 04/23] Tests pass for qmckl_dgemm. #32. --- org/qmckl_blas.org | 47 +++++++++++++++++++++++----------------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index fd7e34e..b05cac5 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -90,13 +90,13 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, integer*8 , intent(in) :: m, n, k real*8 , intent(in) :: alpha, beta integer*8 , intent(in) :: lda - real*8 , intent(in) :: A(m,n) + real*8 , intent(in) :: A(m,k) integer*8 , intent(in) :: ldb - real*8 , intent(in) :: B(n,k) + real*8 , intent(in) :: B(k,n) integer*8 , intent(in) :: ldc real*8 , intent(out) :: C(m,n) - integer*8 :: i,j + integer*8 :: i,j,l info = QMCKL_SUCCESS @@ -120,17 +120,17 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, return endif - if (LDA < m) then + if (LDA .ne. m) then info = QMCKL_INVALID_ARG_9 return endif - if (LDB < n) then + if (LDB .ne. k) then info = QMCKL_INVALID_ARG_10 return endif - if (LDB < n) then + if (LDC .ne. m) then info = QMCKL_INVALID_ARG_13 return endif @@ -153,8 +153,8 @@ end function qmckl_dgemm_f implicit none integer (c_int64_t) , intent(in) , value :: context - logical*8 , intent(in) , value :: TransA - logical*8 , intent(in) , value :: TransB + logical , intent(in) , value :: TransA + logical , intent(in) , value :: TransB integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: k @@ -188,17 +188,17 @@ end function qmckl_dgemm_f implicit none integer (c_int64_t) , intent(in) , value :: context - logical*8 , intent(in) , value :: TransA - logical*8 , intent(in) , value :: TransB + logical , intent(in) , value :: TransA + logical , intent(in) , value :: TransB integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: k - real (c_double ) , intent(in) :: alpha + real (c_double ) , intent(in) , value :: alpha integer (c_int64_t) , intent(in) , value :: lda real (c_double ) , intent(in) :: A(lda,*) integer (c_int64_t) , intent(in) , value :: ldb real (c_double ) , intent(in) :: B(ldb,*) - real (c_double ) , intent(in) :: beta + real (c_double ) , intent(in) , value :: beta integer (c_int64_t) , intent(in) , value :: ldc real (c_double ) , intent(out) :: C(ldc,*) @@ -216,11 +216,11 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:) integer*8 :: m, n, k, LDA, LDB, LDC integer*8 :: i,j,l - logical*8 :: TransA, TransB + logical :: TransA, TransB double precision :: x, alpha, beta TransA = .False. - TransB = .False. + TransB = .False. m = 5_8 k = 4_8 n = 6_8 @@ -234,14 +234,16 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) B = 0.d0 C = 0.d0 D = 0.d0 - do j=1,m - do i=1,k + alpha = 1.0d0 + beta = 0.0d0 + do j=1,k + do i=1,m A(i,j) = -10.d0 + dble(i+j) end do end do - do j=1,k - do i=1,n + do j=1,n + do i=1,k B(i,j) = -10.d0 + dble(i+j) end do end do @@ -253,21 +255,20 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) test_qmckl_dgemm = QMCKL_FAILURE x = 0.d0 - do j=1,m - do i=l,n + do j=1,n + do i=1,m do l=1,k - D(i,j) = D(i,j) + A(i,k)*B(k,j) + D(i,j) = D(i,j) + A(i,l)*B(l,j) end do x = x + (D(i,j) - C(i,j))**2 end do end do - print *,"DABS(X)= ",dabs(x) if (dabs(x) <= 1.d-15) then test_qmckl_dgemm = QMCKL_SUCCESS endif - deallocate(A,B) + deallocate(A,B,C,D) end function test_qmckl_dgemm #+end_src From 3fa60a3009bf3e19b2380caa7808a20395f87ebf Mon Sep 17 00:00:00 2001 From: v1j4y Date: Wed, 15 Sep 2021 18:30:41 +0200 Subject: [PATCH 05/23] Implemented MOs. Needs testing. #32 --- org/qmckl_context.org | 3 + org/qmckl_mo.org | 671 ++++++++++++++++++++++++++++++++++++++++++ org/table_of_contents | 1 + 3 files changed, 675 insertions(+) create mode 100644 org/qmckl_mo.org diff --git a/org/qmckl_context.org b/org/qmckl_context.org index d2ca467..bd1a0ce 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -31,10 +31,12 @@ int main() { #include "qmckl_nucleus_private_type.h" #include "qmckl_electron_private_type.h" #include "qmckl_ao_private_type.h" +#include "qmckl_mo_private_type.h" #include "qmckl_jastrow_private_type.h" #include "qmckl_nucleus_private_func.h" #include "qmckl_electron_private_func.h" #include "qmckl_ao_private_func.h" +#include "qmckl_mo_private_func.h" #+end_src #+begin_src c :tangle (eval c) @@ -119,6 +121,7 @@ typedef struct qmckl_context_struct { qmckl_nucleus_struct nucleus; qmckl_electron_struct electron; qmckl_ao_basis_struct ao_basis; + qmckl_mo_basis_struct mo_basis; qmckl_jastrow_struct jastrow; /* To be implemented: diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org new file mode 100644 index 0000000..c6a6b6b --- /dev/null +++ b/org/qmckl_mo.org @@ -0,0 +1,671 @@ +#+TITLE: Molecular Orbitals +#+SETUPFILE: ../tools/theme.setup +#+INCLUDE: ../tools/lib.org + +The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO +coefficient matrix \[C\]. Using these coefficients (e.g. from Hartree Fock SCF method) +the MOs are defined as follows: + +\[ +\phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r}) +\] + + +In this section we demonstrate how to use the QMCkl specific DGEMM +function to calculate the MOs. + + +* Headers :noexport: + #+begin_src elisp :noexport :results none +(org-babel-lob-ingest "../tools/lib.org") + #+end_src + + + #+begin_src c :tangle (eval h_private_type) +#ifndef QMCKL_MO_HPT +#define QMCKL_MO_HPT + +#include + #+end_src + + #+begin_src c :tangle (eval c_test) :noweb yes +#include "qmckl.h" +#include "assert.h" +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include "chbrclf.h" +#include "qmckl_ao_private_func.h" +#include "qmckl_mo_private_func.h" + +int main() { + qmckl_context context; + context = qmckl_context_create(); + + qmckl_exit_code rc; + #+end_src + + #+begin_src c :tangle (eval c) +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_STDINT_H +#include +#elif HAVE_INTTYPES_H +#include +#endif + +#include +#include +#include +#include + +#include "qmckl.h" +#include "qmckl_context_private_type.h" +#include "qmckl_memory_private_type.h" +#include "qmckl_memory_private_func.h" +#include "qmckl_ao_private_type.h" +#include "qmckl_ao_private_func.h" +#include "qmckl_mo_private_type.h" +#include "qmckl_mo_private_func.h" + #+end_src + + +* Context + + The following arrays are stored in the context: + + |---------------------+--------------------+--------------------------------------------------------------| + | ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) | + | ~nucleus_index~ | ~[nucl_num]~ | Index of the first shell of each nucleus | + | ~nucleus_shell_num~ | ~[nucl_num]~ | Number of shells per nucleus | + | ~ao_num~ | | Number of AOs | + | ~ao_cartesian~ | | If true, use polynomials. Otherwise, use spherical harmonics | + | ~ao_factor~ | ~[ao_num]~ | Normalization factor of the AO | + | ~ao_shell~ | ~[ao_num]~ | For each AO, specify to which shell it belongs | + | ~coefficient~ | ~[mo_num, ao_num]~ | Orbital coefficients | + + Computed data: + + |---------------+-----------------------------------+----------------------------------------------------------------------------------------| + |---------------+-----------------------------------+----------------------------------------------------------------------------------------| + | ~mo_vgl~ | ~[5][walk_num][elec_num][mo_num]~ | Value, gradients, Laplacian of the MOs at electron positions | + | ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions | + |---------------+-----------------------------------+----------------------------------------------------------------------------------------| + +** Data structure + + #+begin_src c :comments org :tangle (eval h_private_type) +typedef struct qmckl_mo_basis_struct { + int64_t mo_num; + double * coefficient; + + double * mo_vgl; + int64_t mo_vgl_date; + + int32_t uninitialized; + bool provided; + char type; +} qmckl_mo_basis_struct; + #+end_src + + The ~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. + +** Access functions + + #+begin_src c :comments org :tangle (eval h_private_func) :exports none +char qmckl_get_mo_basis_type (const qmckl_context context); +int64_t qmckl_get_mo_basis_mo_num (const qmckl_context context); +double* qmckl_get_mo_basis_coefficient (const qmckl_context context); + #+end_src + + When all the data for the AOs have been provided, the following + function returns ~true~. + + #+begin_src c :comments org :tangle (eval h_func) +bool qmckl_mo_basis_provided (const qmckl_context context); + #+end_src + + #+NAME:post + #+begin_src c :exports none +if ( (ctx->mo_basis.uninitialized & mask) != 0) { + return NULL; +} + #+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +char qmckl_get_mo_basis_type (const qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1; + + if ( (ctx->mo_basis.uninitialized & mask) != 0) { + return (char) 0; + } + + assert (ctx->mo_basis.type != (char) 0); + return ctx->mo_basis.type; +} + +int64_t qmckl_get_mo_basis_mo_num (const qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (int64_t) 0; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + int32_t mask = 1 << 1; + + if ( (ctx->mo_basis.uninitialized & mask) != 0) { + return (int64_t) 0; + } + + assert (ctx->mo_basis.mo_num > (int64_t) 0); + return ctx->mo_basis.mo_num; +} + + #+end_src + +** Initialization functions + + To set the basis set, all the following functions need to be + called. + + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code qmckl_set_mo_basis_type (qmckl_context context, const char t); +qmckl_exit_code qmckl_set_mo_basis_mo_num (qmckl_context context, const int64_t mo_num); +qmckl_exit_code qmckl_set_mo_basis_coefficient (qmckl_context context, const double * coefficient); + #+end_src + + #+NAME:pre2 + #+begin_src c :exports none +if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + +qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + #+end_src + + #+NAME:post2 + #+begin_src c :exports none +ctx->mo_basis.uninitialized &= ~mask; +ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0); +if (ctx->mo_basis.provided) { + qmckl_exit_code rc_ = qmckl_finalize_basis(context); + if (rc_ != QMCKL_SUCCESS) return rc_; + } + +return QMCKL_SUCCESS; + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_set_mo_basis_type(qmckl_context context, const char t) { + <> + + if (t != 'G' && t != 'S') { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_mo_basis_type", + NULL); + } + + int32_t mask = 1; + ctx->mo_basis.type = t; + + <> +} + +qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t mo_num) { + <> + + if (mo_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_mo_basis_mo_num", + "mo_num <= 0"); + } + + int32_t mask = 1 << 1; + ctx->mo_basis.mo_num = mo_num; + + <> +} + +qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const double* coefficient) { + <> + + int32_t mask = 1 << 2; + + if (ctx->mo_basis.coefficient != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_mo_basis_coefficient", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double); + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if (new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_mo_basis_coefficient", + NULL); + } + + memcpy(new_array, coefficient, mem_info.size); + + ctx->mo_basis.coefficient = new_array; + + <> +} + + #+end_src + + When the basis set is completely entered, other data structures are + computed to accelerate the calculations. + +* Computation + +** Computation of MOs + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ao_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = 5 * ctx->electron.num * ctx->mo_basis.mo_num * ctx->electron.walk_num; + memcpy(mo_vgl, ctx->mo_basis.mo_vgl, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_mo_basis_vgl (context, mo_vgl) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: mo_vgl(*) + end function + end interface + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_mo_basis_vgl(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_mo_basis_vgl(qmckl_context context) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->mo_basis.mo_vgl_date) { + + /* Allocate array */ + if (ctx->mo_basis.mo_vgl == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 5 * ctx->electron.num * ctx->mo_basis.mo_num * + ctx->electron.walk_num * sizeof(double); + double* mo_vgl = (double*) qmckl_malloc(context, mem_info); + + if (mo_vgl == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_mo_basis_vgl", + NULL); + } + ctx->mo_basis.mo_vgl = mo_vgl; + } + + qmckl_exit_code rc; + if (ctx->mo_basis.type == 'G') { + rc = qmckl_compute_mo_basis_gaussian_vgl(context, + ctx->ao_basis.ao_num, + ctx->mo_basis.mo_num, + ctx->electron.num, + ctx->electron.walk_num, + ctx->mo_basis.coefficient, + ctx->ao_basis.ao_vgl, + ctx->mo_basis.mo_vgl); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_mo_basis_vgl", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->mo_basis.mo_vgl_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_mo_basis_gaussian_vgl + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_mo_basis_gaussian_vgl_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~ao_num~ | in | Number of AOs | + | ~int64_t~ | ~mo_num~ | in | Number of MOs | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~double~ | ~coef_normalized[mo_num][ao_num]~ | in | AO to MO transformation matrix | + | ~double~ | ~ao_vgl[5][walk_num][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the AOs | + | ~double~ | ~mo_vgl[5][walk_num][elec_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & + ao_num, mo_num, elec_num, walk_num, & + coef_normalized, ao_vgl, mo_vgl) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: ao_num, mo_num + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: ao_vgl(ao_num,elec_num,walk_num,5) + double precision , intent(in) :: coef_normalized(mo_num,ao_num) + double precision , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5) + logical :: TransA, TransB + double precision :: alpha, beta + integer :: info_qmckl_dgemm_value + integer :: info_qmckl_dgemm_Gx + integer :: info_qmckl_dgemm_Gy + integer :: info_qmckl_dgemm_Gz + integer :: info_qmckl_dgemm_lap + integer*8 :: M, N, K, LDA, LDB, LDC + + integer*8 :: inucl, iprim, iwalk, ielec, ishell + double precision :: x, y, z, two_a, ar2, r2, v, cutoff + TransA = .False. + TransB = .False. + alpha = 1.0d0 + beta = 0.0d0 + + info = QMCKL_SUCCESS + + ! Don't compute exponentials when the result will be almost zero. + ! TODO : Use numerical precision here + cutoff = -dlog(1.d-15) + M = mo_num * 1_8 + N = ao_num * 1_8 + K = 1_8 + + do iwalk = 1, walk_num + do ielec = 1, elec_num + ! Value + info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & + coef_normalized,size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,1),size(mo_vgl,1) * 1_8) + ! Grad_x + info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & + coef_normalized,size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,2),size(mo_vgl,1) * 1_8) + ! Grad_y + info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & + coef_normalized,size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,3),size(mo_vgl,1) * 1_8) + ! Grad_z + info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & + coef_normalized,size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,4),size(mo_vgl,1) * 1_8) + ! Lapl_z + info_qmckl_dgemm_lap = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & + coef_normalized,size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,5),size(mo_vgl,1) * 1_8) + end do + end do + + if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. & + info_qmckl_dgemm_Gx .eq. QMCKL_SUCCESS .and. & + info_qmckl_dgemm_Gy .eq. QMCKL_SUCCESS .and. & + info_qmckl_dgemm_Gz .eq. QMCKL_SUCCESS .and. & + info_qmckl_dgemm_lap .eq. QMCKL_SUCCESS ) then + info = QMCKL_SUCCESS + else + info = QMCKL_FAILURE + end if + +end function qmckl_compute_mo_basis_gaussian_vgl_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl")) + + #+RESULTS: + #+BEGIN_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_mo_basis_gaussian_vgl ( + const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t elec_num, + const int64_t walk_num, + const double* coef_normalized, + const double* ao_vgl, + double* const mo_vgl ); + #+END_src + + + #+CALL: generate_c_interface(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl")) + + #+RESULTS: + #+BEGIN_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_mo_basis_gaussian_vgl & + (context, ao_num, mo_num, elec_num, walk_num, coef_normalized, ao_vgl, mo_vgl) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: ao_num + integer (c_int64_t) , intent(in) , value :: mo_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: coef_normalized(ao_num,mo_num) + real (c_double ) , intent(in) :: ao_vgl(mo_num,elec_num,walk_num,5) + real (c_double ) , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5) + + integer(c_int32_t), external :: qmckl_compute_mo_basis_gaussian_vgl_f + info = qmckl_compute_mo_basis_gaussian_vgl_f & + (context, ao_num, mo_num, elec_num, walk_num, coef_normalized, ao_vgl, mo_vgl) + + end function qmckl_compute_mo_basis_gaussian_vgl + #+END_src + + + #+begin_src python :results output :exports none +import numpy as np + +def f(a,x,y): + return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] ) + +def df(a,x,y,n): + h0 = 1.e-6 + if n == 1: h = np.array([h0,0.,0.]) + elif n == 2: h = np.array([0.,h0,0.]) + elif n == 3: h = np.array([0.,0.,h0]) + return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0) + +def d2f(a,x,y,n): + h0 = 1.e-6 + if n == 1: h = np.array([h0,0.,0.]) + elif n == 2: h = np.array([0.,h0,0.]) + elif n == 3: h = np.array([0.,0.,h0]) + return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2 + +def lf(a,x,y): + return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3) + +elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] ) +elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] ) +nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) +nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) + +#double prim_vgl[prim_num][5][walk_num][elec_num]; +x = elec_26_w1 ; y = nucl_1 +a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ), + ( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ), + ( 2.808000E+02, -4.540000E-03 * 4.8888635917437597e+01 ), + ( 7.927000E+01, -1.813300E-02 * 1.8933972232608955e+01 ), + ( 2.559000E+01, -5.576000E-02 * 8.1089160941724145e+00 ), + ( 8.997000E+00, -1.268950E-01 * 3.7024003863155635e+00 ), + ( 3.319000E+00, -1.703520E-01 * 1.7525302846177560e+00 ), + ( 9.059000E-01, 1.403820E-01 * 6.6179013183966806e-01 ), + ( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ), + ( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )] + +print ( "[1][0][0][26] : %25.15e"% f(a,x,y)) +print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1)) +print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2)) +print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3)) +print ( "[1][4][0][26] : %25.15e"% lf(a,x,y)) + +x = elec_15_w2 ; y = nucl_2 +a = [(3.387000E+01, 6.068000E-03 *1.0006253235944540e+01), + (5.095000E+00, 4.530800E-02 *2.4169531573445120e+00), + (1.159000E+00, 2.028220E-01 *7.9610924849766440e-01), + (3.258000E-01, 5.039030E-01 *3.0734305383061117e-01), + (1.027000E-01, 3.834210E-01 *1.2929684417481876e-01)] + +print ( "[0][1][15][14] : %25.15e"% f(a,x,y)) +print ( "[1][1][15][14] : %25.15e"% df(a,x,y,1)) +print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2)) +print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3)) +print ( "[4][1][15][14] : %25.15e"% lf(a,x,y)) + + #+end_src + + +*** Test + + #+begin_src c :tangle (eval c_test) :exports none + + #+end_src + + +* End of files :noexport: + + #+begin_src c :tangle (eval h_private_type) +#endif + #+end_src + +*** Test + #+begin_src c :tangle (eval c_test) + rc = qmckl_context_destroy(context); + assert (rc == QMCKL_SUCCESS); + + return 0; +} + #+end_src + +*** Compute file names + #+begin_src emacs-lisp +; The following is required to compute the file names + +(setq pwd (file-name-directory buffer-file-name)) +(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) +(setq f (concat pwd name "_f.f90")) +(setq fh (concat pwd name "_fh.f90")) +(setq c (concat pwd name ".c")) +(setq h (concat name ".h")) +(setq h_private (concat name "_private.h")) +(setq c_test (concat pwd "test_" name ".c")) +(setq f_test (concat pwd "test_" name "_f.f90")) + +; Minted +(require 'ox-latex) +(setq org-latex-listings 'minted) +(add-to-list 'org-latex-packages-alist '("" "listings")) +(add-to-list 'org-latex-packages-alist '("" "color")) + + #+end_src + + +# -*- mode: org -*- +# vim: syntax=c + + diff --git a/org/table_of_contents b/org/table_of_contents index 37dcf21..4e0600c 100644 --- a/org/table_of_contents +++ b/org/table_of_contents @@ -6,6 +6,7 @@ qmckl_numprec.org qmckl_nucleus.org qmckl_electron.org qmckl_ao.org +qmckl_mo.org qmckl_jastrow.org qmckl_sherman_morrison_woodbury.org qmckl_distance.org From a52d6683cc39e9c1e9d2806761c4ea45f45bb0d6 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 17:18:43 +0200 Subject: [PATCH 06/23] Allow 0 beta electrons. #32 --- org/qmckl_electron.org | 4 +- org/qmckl_mo.org | 153 ++++++++++++++++++++++++++++++++++++++--- org/qmckl_tests.org | 2 +- 3 files changed, 145 insertions(+), 14 deletions(-) diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 945a214..9f2f89f 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -532,11 +532,11 @@ qmckl_set_electron_num(qmckl_context context, "up_num <= 0"); } - if (down_num <= 0) { + if (down_num < 0) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_set_electron_num", - "down_num <= 0"); + "down_num < 0"); } int32_t mask = 1 << 0; diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index c6a6b6b..52c829f 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -74,7 +74,6 @@ int main() { #include "qmckl_mo_private_func.h" #+end_src - * Context The following arrays are stored in the context: @@ -469,34 +468,34 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & do ielec = 1, elec_num ! Value info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & + ao_vgl(:, :, iwalk, 1), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,1),size(mo_vgl,1) * 1_8) + mo_vgl(:,:,iwalk,1),size(mo_vgl,1) * 1_8) ! Grad_x info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & + ao_vgl(:, :, iwalk, 2), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,2),size(mo_vgl,1) * 1_8) + mo_vgl(:,:,iwalk,2),size(mo_vgl,1) * 1_8) ! Grad_y info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & + ao_vgl(:, :, iwalk, 3), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,3),size(mo_vgl,1) * 1_8) + mo_vgl(:,:,iwalk,3),size(mo_vgl,1) * 1_8) ! Grad_z info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & + ao_vgl(:, :, iwalk, 4), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,4),size(mo_vgl,1) * 1_8) + mo_vgl(:,:,iwalk,4),size(mo_vgl,1) * 1_8) ! Lapl_z info_qmckl_dgemm_lap = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & + ao_vgl(:, :, iwalk, 5), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,5),size(mo_vgl,1) * 1_8) + mo_vgl(:,:,iwalk,5),size(mo_vgl,1) * 1_8) end do end do @@ -623,6 +622,138 @@ print ( "[4][1][15][14] : %25.15e"% lf(a,x,y)) *** Test #+begin_src c :tangle (eval c_test) :exports none +{ +#define walk_num chbrclf_walk_num +#define elec_num chbrclf_elec_num +#define shell_num chbrclf_shell_num +#define ao_num chbrclf_ao_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]); +const int64_t nucl_num = chbrclf_nucl_num; +const double* nucl_charge = chbrclf_charge; +const double* nucl_coord = &(chbrclf_nucl_coord[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); + +rc = qmckl_set_nucleus_num (context, nucl_num); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0])); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_nucleus_charge(context, nucl_charge); +assert(rc == QMCKL_SUCCESS); + +assert(qmckl_nucleus_provided(context)); + +//const int64_t shell_num = chbrclf_shell_num; +//const int64_t prim_num = chbrclf_prim_num; +//const int64_t ao_num = chbrclf_ao_num; +const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]); +const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]); +const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]); +const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]); +const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]); +const double * shell_factor = &(chbrclf_basis_shell_factor[0]); +const double * exponent = &(chbrclf_basis_exponent[0]); +const double * coefficient = &(chbrclf_basis_coefficient[0]); +const double * prim_factor = &(chbrclf_basis_prim_factor[0]); +const double * ao_factor = &(chbrclf_basis_ao_factor[0]); + +const char typ = 'G'; + +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_type (context, typ); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_shell_factor (context, shell_factor); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_exponent (context, exponent); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_coefficient (context, coefficient); +assert(rc == QMCKL_SUCCESS); +assert(!qmckl_ao_basis_provided(context)); + +rc = qmckl_set_ao_basis_prim_factor (context, prim_factor); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_ao_basis_ao_factor (context, ao_factor); +assert(rc == QMCKL_SUCCESS); + +assert(qmckl_ao_basis_provided(context)); + + +double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num]; + +rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_get_mo_basis_vgl(context, &(ao_vgl[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +//printf("\n"); +//printf(" ao_vgl ao_vgl[0][0][26][219] %25.15e\n", ao_vgl[0][0][26][219]); +//printf(" ao_vgl ao_vgl[1][0][26][219] %25.15e\n", ao_vgl[1][0][26][219]); +//printf(" ao_vgl ao_vgl[0][0][26][220] %25.15e\n", ao_vgl[0][0][26][220]); +//printf(" ao_vgl ao_vgl[1][0][26][220] %25.15e\n", ao_vgl[1][0][26][220]); +//printf(" ao_vgl ao_vgl[0][0][26][221] %25.15e\n", ao_vgl[0][0][26][221]); +//printf(" ao_vgl ao_vgl[1][0][26][221] %25.15e\n", ao_vgl[1][0][26][221]); +//printf(" ao_vgl ao_vgl[0][0][26][222] %25.15e\n", ao_vgl[0][0][26][222]); +//printf(" ao_vgl ao_vgl[1][0][26][222] %25.15e\n", ao_vgl[1][0][26][222]); +//printf(" ao_vgl ao_vgl[0][0][26][223] %25.15e\n", ao_vgl[0][0][26][223]); +//printf(" ao_vgl ao_vgl[1][0][26][223] %25.15e\n", ao_vgl[1][0][26][223]); +//printf(" ao_vgl ao_vgl[0][0][26][224] %25.15e\n", ao_vgl[0][0][26][224]); +//printf(" ao_vgl ao_vgl[1][0][26][224] %25.15e\n", ao_vgl[1][0][26][224]); +//printf("\n"); +} #+end_src diff --git a/org/qmckl_tests.org b/org/qmckl_tests.org index 8e08c6a..da89fbb 100644 --- a/org/qmckl_tests.org +++ b/org/qmckl_tests.org @@ -530,7 +530,7 @@ F 1 int64_t chbrclf_basis_nucleus_index[chbrclf_nucl_num] = {0, 14, 23, 37, 53}; -int64_t chbrclf_basis_nucleus_shell_num[chbrclf_nucl_num] = {14, 9, 14, 16, 19}; +int64_t chbrclf_basis_nucleus_shell_num[chbrclf_nucl_num] = {14, 9, 14, 16, 19}; // C, H, F, Cl, Br int32_t chbrclf_basis_shell_ang_mom[chbrclf_shell_num] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 0, 1, 1, 1, 2, 2, 0, From ff5e7882d077e4b6ff64d3241fba3fd368574648 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 18:16:51 +0200 Subject: [PATCH 07/23] Working on adding TransA and TransB to DGEMM. #32 --- org/qmckl_blas.org | 51 +++++++++++++++++++++++++++++++++++++--------- 1 file changed, 41 insertions(+), 10 deletions(-) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index b05cac5..4bb7384 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -95,11 +95,36 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, real*8 , intent(in) :: B(k,n) integer*8 , intent(in) :: ldc real*8 , intent(out) :: C(m,n) + real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:) - integer*8 :: i,j,l + integer*8 :: i,j,l, LDA_2, LDB_2 info = QMCKL_SUCCESS + if (TransA) then + allocate(AT(k,m)) + do i = 1, m + do j = 1, k + AT(j,i) = A(i,j) + end do + end do + LDA_2 = M + else + LDA_2 = LDA + endif + + if (TransB) then + allocate(BT(n,k)) + do i = 1, k + do j = 1, n + BT(j,i) = B(i,j) + end do + end do + LDB_2 = K + else + LDB_2 = LDB + endif + if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return @@ -120,12 +145,12 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, return endif - if (LDA .ne. m) then + if (LDA_2 .ne. m) then info = QMCKL_INVALID_ARG_9 return endif - if (LDB .ne. k) then + if (LDB_2 .ne. k) then info = QMCKL_INVALID_ARG_10 return endif @@ -135,7 +160,13 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, return endif - C = matmul(A,B) + if (TransA) then + C = matmul(AT,B) + else if (TransB) then + C = matmul(A,BT) + else + C = matmul(A,B) + endif end function qmckl_dgemm_f #+end_src @@ -153,8 +184,8 @@ end function qmckl_dgemm_f implicit none integer (c_int64_t) , intent(in) , value :: context - logical , intent(in) , value :: TransA - logical , intent(in) , value :: TransB + logical*8 , intent(in) , value :: TransA + logical*8 , intent(in) , value :: TransB integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: k @@ -188,8 +219,8 @@ end function qmckl_dgemm_f implicit none integer (c_int64_t) , intent(in) , value :: context - logical , intent(in) , value :: TransA - logical , intent(in) , value :: TransB + logical*8 , intent(in) , value :: TransA + logical*8 , intent(in) , value :: TransB integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: k @@ -216,7 +247,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:) integer*8 :: m, n, k, LDA, LDB, LDC integer*8 :: i,j,l - logical :: TransA, TransB + logical*8 :: TransA, TransB double precision :: x, alpha, beta TransA = .False. @@ -271,7 +302,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) deallocate(A,B,C,D) end function test_qmckl_dgemm #+end_src - + #+begin_src c :comments link :tangle (eval c_test) qmckl_exit_code test_qmckl_dgemm(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); From a0557322270d02b4de9f0db0760913a2468e4b28 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 18:17:21 +0200 Subject: [PATCH 08/23] Fixed bug and working on tests for MOs. #32 --- org/qmckl_mo.org | 78 +++++++++++++++++++++++++++++++++++------------- 1 file changed, 57 insertions(+), 21 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 52c829f..8171d5b 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -306,6 +306,9 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v rc = qmckl_provide_ao_vgl(context); if (rc != QMCKL_SUCCESS) return rc; + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); @@ -333,11 +336,11 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v *** Provide #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_provide_mo_basis_vgl(qmckl_context context); +qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_mo_basis_vgl(qmckl_context context) +qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { @@ -439,7 +442,7 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & double precision , intent(in) :: ao_vgl(ao_num,elec_num,walk_num,5) double precision , intent(in) :: coef_normalized(mo_num,ao_num) double precision , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5) - logical :: TransA, TransB + logical*8 :: TransA, TransB double precision :: alpha, beta integer :: info_qmckl_dgemm_value integer :: info_qmckl_dgemm_Gx @@ -456,46 +459,61 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & beta = 0.0d0 info = QMCKL_SUCCESS + info_qmckl_dgemm_value = QMCKL_SUCCESS + info_qmckl_dgemm_Gx = QMCKL_SUCCESS + info_qmckl_dgemm_Gy = QMCKL_SUCCESS + info_qmckl_dgemm_Gz = QMCKL_SUCCESS + info_qmckl_dgemm_lap = QMCKL_SUCCESS ! Don't compute exponentials when the result will be almost zero. ! TODO : Use numerical precision here cutoff = -dlog(1.d-15) - M = mo_num * 1_8 - N = ao_num * 1_8 - K = 1_8 + M = 1_8 + N = 1_8 + K = mo_num * 1_8 do iwalk = 1, walk_num do ielec = 1, elec_num ! Value + TransA = .True. + TransB = .True. info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 1), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,1),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,1),1_8) ! Grad_x + TransA = .True. + TransB = .True. info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 2), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,2),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,2),1_8) ! Grad_y + TransA = .True. + TransB = .True. info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 3), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,3),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,3),1_8) ! Grad_z + TransA = .True. + TransB = .True. info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 4), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,4),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,4),1_8) ! Lapl_z - info_qmckl_dgemm_lap = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 5), size(ao_vgl,1) * 1_8, & + TransA = .True. + TransB = .True. + info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,5),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,5),1_8) end do end do @@ -657,9 +675,6 @@ assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); -//const int64_t shell_num = chbrclf_shell_num; -//const int64_t prim_num = chbrclf_prim_num; -//const int64_t ao_num = chbrclf_ao_num; const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]); const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]); const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]); @@ -736,9 +751,30 @@ double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num]; rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0])); assert (rc == QMCKL_SUCCESS); -rc = qmckl_get_mo_basis_vgl(context, &(ao_vgl[0][0][0][0])); +/* Set up MO data */ +rc = qmckl_set_mo_basis_type(context, typ); assert (rc == QMCKL_SUCCESS); +const int64_t mo_num = chbrclf_ao_num; +rc = qmckl_set_mo_basis_mo_num(context, mo_num); +assert (rc == QMCKL_SUCCESS); + +double mo_coefficient[chbrclf_ao_num][mo_num]; +const int64_t min_num = (((chbrclf_ao_num) < (mo_num)) ? (chbrclf_ao_num) : (mo_num)); +for (int i=0; i Date: Mon, 27 Sep 2021 18:20:41 +0200 Subject: [PATCH 09/23] Deactivate mo test atm. --- org/qmckl_mo.org | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 8171d5b..31da1cd 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -770,9 +770,9 @@ assert (rc == QMCKL_SUCCESS); assert(qmckl_ao_basis_provided(context)); -double mo_vgl[5][walk_num][elec_num][chbrclf_ao_num]; -rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +//double mo_vgl[5][walk_num][elec_num][chbrclf_ao_num]; +//rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); +//assert (rc == QMCKL_SUCCESS); //printf("\n"); From c1c5a4a8aaa5021b6eb46044d57ad07445c0dc72 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 18:30:59 +0200 Subject: [PATCH 10/23] Provided function to check for MO basis. #32 --- org/qmckl_mo.org | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 31da1cd..e13cbc5 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -180,6 +180,19 @@ int64_t qmckl_get_mo_basis_mo_num (const qmckl_context context) { return ctx->mo_basis.mo_num; } +bool qmckl_mo_basis_provided(const 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); + + return ctx->mo_basis.provided; +} + + #+end_src ** Initialization functions @@ -768,7 +781,7 @@ for (int i=0; i Date: Mon, 27 Sep 2021 18:39:45 +0200 Subject: [PATCH 11/23] Check whether MO basis has been initialized before computations. #32 --- org/qmckl_mo.org | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index e13cbc5..2914c62 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -377,6 +377,13 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) NULL); } + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + /* Compute if necessary */ if (ctx->electron.coord_new_date > ctx->mo_basis.mo_vgl_date) { From 922c52c36a206dcc3bc6478f89c6270ca47e2ad7 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 18:41:01 +0200 Subject: [PATCH 12/23] Compiles locally. #32 --- org/qmckl_mo.org | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 2914c62..4c8be42 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -790,8 +790,8 @@ assert (rc == QMCKL_SUCCESS); assert(qmckl_mo_basis_provided(context)); -//double mo_vgl[5][walk_num][elec_num][chbrclf_ao_num]; -//rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); +double mo_vgl[5][walk_num][elec_num][chbrclf_ao_num]; +rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); //assert (rc == QMCKL_SUCCESS); From d3311ec2d95506436f285eaab92c32777cfb11cd Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 22:01:13 +0200 Subject: [PATCH 13/23] Testing provide mo. --- org/qmckl_mo.org | 66 ++++++++++++++++++++++++------------------------ 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 4c8be42..80139b1 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -405,24 +405,24 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) } qmckl_exit_code rc; - if (ctx->mo_basis.type == 'G') { - rc = qmckl_compute_mo_basis_gaussian_vgl(context, - ctx->ao_basis.ao_num, - ctx->mo_basis.mo_num, - ctx->electron.num, - ctx->electron.walk_num, - ctx->mo_basis.coefficient, - ctx->ao_basis.ao_vgl, - ctx->mo_basis.mo_vgl); - } else { - return qmckl_failwith( context, - QMCKL_FAILURE, - "compute_mo_basis_vgl", - "Not yet implemented"); - } - if (rc != QMCKL_SUCCESS) { - return rc; - } + //if (ctx->mo_basis.type == 'G') { + // rc = qmckl_compute_mo_basis_gaussian_vgl(context, + // ctx->ao_basis.ao_num, + // ctx->mo_basis.mo_num, + // ctx->electron.num, + // ctx->electron.walk_num, + // ctx->mo_basis.coefficient, + // ctx->ao_basis.ao_vgl, + // ctx->mo_basis.mo_vgl); + //} else { + // return qmckl_failwith( context, + // QMCKL_FAILURE, + // "compute_mo_basis_vgl", + // "Not yet implemented"); + //} + //if (rc != QMCKL_SUCCESS) { + // return rc; + //} ctx->mo_basis.mo_vgl_date = ctx->date; } @@ -656,7 +656,7 @@ print ( "[4][1][15][14] : %25.15e"% lf(a,x,y)) #+end_src - + *** Test #+begin_src c :tangle (eval c_test) :exports none @@ -795,20 +795,20 @@ rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); //assert (rc == QMCKL_SUCCESS); -//printf("\n"); -//printf(" ao_vgl ao_vgl[0][0][26][219] %25.15e\n", ao_vgl[0][0][26][219]); -//printf(" ao_vgl ao_vgl[1][0][26][219] %25.15e\n", ao_vgl[1][0][26][219]); -//printf(" ao_vgl ao_vgl[0][0][26][220] %25.15e\n", ao_vgl[0][0][26][220]); -//printf(" ao_vgl ao_vgl[1][0][26][220] %25.15e\n", ao_vgl[1][0][26][220]); -//printf(" ao_vgl ao_vgl[0][0][26][221] %25.15e\n", ao_vgl[0][0][26][221]); -//printf(" ao_vgl ao_vgl[1][0][26][221] %25.15e\n", ao_vgl[1][0][26][221]); -//printf(" ao_vgl ao_vgl[0][0][26][222] %25.15e\n", ao_vgl[0][0][26][222]); -//printf(" ao_vgl ao_vgl[1][0][26][222] %25.15e\n", ao_vgl[1][0][26][222]); -//printf(" ao_vgl ao_vgl[0][0][26][223] %25.15e\n", ao_vgl[0][0][26][223]); -//printf(" ao_vgl ao_vgl[1][0][26][223] %25.15e\n", ao_vgl[1][0][26][223]); -//printf(" ao_vgl ao_vgl[0][0][26][224] %25.15e\n", ao_vgl[0][0][26][224]); -//printf(" ao_vgl ao_vgl[1][0][26][224] %25.15e\n", ao_vgl[1][0][26][224]); -//printf("\n"); +printf("\n"); +printf(" mo_vgl mo_vgl[0][0][26][219] %25.15e\n", mo_vgl[0][0][2][3]); +printf(" mo_vgl mo_vgl[1][0][26][219] %25.15e\n", mo_vgl[1][0][2][3]); +printf(" mo_vgl mo_vgl[0][0][26][220] %25.15e\n", mo_vgl[0][0][2][3]); +printf(" mo_vgl mo_vgl[1][0][26][220] %25.15e\n", mo_vgl[1][0][2][3]); +printf(" mo_vgl mo_vgl[0][0][26][221] %25.15e\n", mo_vgl[0][0][2][3]); +printf(" mo_vgl mo_vgl[1][0][26][221] %25.15e\n", mo_vgl[1][0][2][3]); +printf(" mo_vgl mo_vgl[0][0][26][222] %25.15e\n", mo_vgl[0][0][2][3]); +printf(" mo_vgl mo_vgl[1][0][26][222] %25.15e\n", mo_vgl[1][0][2][3]); +printf(" mo_vgl mo_vgl[0][0][26][223] %25.15e\n", mo_vgl[0][0][2][3]); +printf(" mo_vgl mo_vgl[1][0][26][223] %25.15e\n", mo_vgl[1][0][2][3]); +printf(" mo_vgl mo_vgl[0][0][26][224] %25.15e\n", mo_vgl[0][0][2][3]); +printf(" mo_vgl mo_vgl[1][0][26][224] %25.15e\n", mo_vgl[1][0][2][3]); +printf("\n"); } #+end_src From 5220df5185a699d409d34847a22b9bae951e9473 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 22:06:49 +0200 Subject: [PATCH 14/23] Checking ao -> mo. --- org/qmckl_mo.org | 144 +++++++++++++++++++++++------------------------ 1 file changed, 72 insertions(+), 72 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 80139b1..cfc4c14 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -405,24 +405,24 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) } qmckl_exit_code rc; - //if (ctx->mo_basis.type == 'G') { - // rc = qmckl_compute_mo_basis_gaussian_vgl(context, - // ctx->ao_basis.ao_num, - // ctx->mo_basis.mo_num, - // ctx->electron.num, - // ctx->electron.walk_num, - // ctx->mo_basis.coefficient, - // ctx->ao_basis.ao_vgl, - // ctx->mo_basis.mo_vgl); - //} else { - // return qmckl_failwith( context, - // QMCKL_FAILURE, - // "compute_mo_basis_vgl", - // "Not yet implemented"); - //} - //if (rc != QMCKL_SUCCESS) { - // return rc; - //} + if (ctx->mo_basis.type == 'G') { + rc = qmckl_compute_mo_basis_gaussian_vgl(context, + ctx->ao_basis.ao_num, + ctx->mo_basis.mo_num, + ctx->electron.num, + ctx->electron.walk_num, + ctx->mo_basis.coefficient, + ctx->ao_basis.ao_vgl, + ctx->mo_basis.mo_vgl); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_mo_basis_vgl", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } ctx->mo_basis.mo_vgl_date = ctx->date; } @@ -445,7 +445,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) | ~int64_t~ | ~elec_num~ | in | Number of electrons | | ~int64_t~ | ~walk_num~ | in | Number of walkers | | ~double~ | ~coef_normalized[mo_num][ao_num]~ | in | AO to MO transformation matrix | - | ~double~ | ~ao_vgl[5][walk_num][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the AOs | + | ~double~ | ~ao_vgl[5][walk_num][elec_num][ao_num]~ | in | Value, gradients and Laplacian of the AOs | | ~double~ | ~mo_vgl[5][walk_num][elec_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs | #+begin_src f90 :comments org :tangle (eval f) :noweb yes @@ -460,7 +460,7 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: walk_num double precision , intent(in) :: ao_vgl(ao_num,elec_num,walk_num,5) - double precision , intent(in) :: coef_normalized(mo_num,ao_num) + double precision , intent(in) :: coef_normalized(ao_num,mo_num) double precision , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5) logical*8 :: TransA, TransB double precision :: alpha, beta @@ -492,50 +492,50 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & N = 1_8 K = mo_num * 1_8 - do iwalk = 1, walk_num - do ielec = 1, elec_num - ! Value - TransA = .True. - TransB = .True. - info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & - coef_normalized,size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,1),1_8) - ! Grad_x - TransA = .True. - TransB = .True. - info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & - coef_normalized,size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,2),1_8) - ! Grad_y - TransA = .True. - TransB = .True. - info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & - coef_normalized,size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,3),1_8) - ! Grad_z - TransA = .True. - TransB = .True. - info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & - coef_normalized,size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,4),1_8) - ! Lapl_z - TransA = .True. - TransB = .True. - info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & - coef_normalized,size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,5),1_8) - end do - end do + !do iwalk = 1, walk_num + ! do ielec = 1, elec_num + ! ! Value + ! TransA = .True. + ! TransB = .True. + ! info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ! ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & + ! coef_normalized,size(coef_normalized,1) * 1_8, & + ! beta, & + ! mo_vgl(:,ielec,iwalk,1),1_8) + ! ! Grad_x + ! TransA = .True. + ! TransB = .True. + ! info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ! ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & + ! coef_normalized,size(coef_normalized,1) * 1_8, & + ! beta, & + ! mo_vgl(:,ielec,iwalk,2),1_8) + ! ! Grad_y + ! TransA = .True. + ! TransB = .True. + ! info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ! ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & + ! coef_normalized,size(coef_normalized,1) * 1_8, & + ! beta, & + ! mo_vgl(:,ielec,iwalk,3),1_8) + ! ! Grad_z + ! TransA = .True. + ! TransB = .True. + ! info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ! ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & + ! coef_normalized,size(coef_normalized,1) * 1_8, & + ! beta, & + ! mo_vgl(:,ielec,iwalk,4),1_8) + ! ! Lapl_z + ! TransA = .True. + ! TransB = .True. + ! info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & + ! ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & + ! coef_normalized,size(coef_normalized,1) * 1_8, & + ! beta, & + ! mo_vgl(:,ielec,iwalk,5),1_8) + ! end do + !end do if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. & info_qmckl_dgemm_Gx .eq. QMCKL_SUCCESS .and. & @@ -553,7 +553,7 @@ end function qmckl_compute_mo_basis_gaussian_vgl_f #+CALL: generate_c_header(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl")) #+RESULTS: - #+BEGIN_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_compute_mo_basis_gaussian_vgl ( const qmckl_context context, const int64_t ao_num, @@ -563,13 +563,13 @@ end function qmckl_compute_mo_basis_gaussian_vgl_f const double* coef_normalized, const double* ao_vgl, double* const mo_vgl ); - #+END_src + #+end_src #+CALL: generate_c_interface(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl")) #+RESULTS: - #+BEGIN_src f90 :tangle (eval f) :comments org :exports none + #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_mo_basis_gaussian_vgl & (context, ao_num, mo_num, elec_num, walk_num, coef_normalized, ao_vgl, mo_vgl) & bind(C) result(info) @@ -583,7 +583,7 @@ end function qmckl_compute_mo_basis_gaussian_vgl_f integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: walk_num real (c_double ) , intent(in) :: coef_normalized(ao_num,mo_num) - real (c_double ) , intent(in) :: ao_vgl(mo_num,elec_num,walk_num,5) + real (c_double ) , intent(in) :: ao_vgl(ao_num,elec_num,walk_num,5) real (c_double ) , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5) integer(c_int32_t), external :: qmckl_compute_mo_basis_gaussian_vgl_f @@ -591,9 +591,12 @@ end function qmckl_compute_mo_basis_gaussian_vgl_f (context, ao_num, mo_num, elec_num, walk_num, coef_normalized, ao_vgl, mo_vgl) end function qmckl_compute_mo_basis_gaussian_vgl - #+END_src + #+end_src + +*** Test + #+begin_src python :results output :exports none import numpy as np @@ -656,9 +659,6 @@ print ( "[4][1][15][14] : %25.15e"% lf(a,x,y)) #+end_src - -*** Test - #+begin_src c :tangle (eval c_test) :exports none { #define walk_num chbrclf_walk_num From 97749325b1f19d23d6fc905b559157da7b1f17ff Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 22:37:51 +0200 Subject: [PATCH 15/23] Removed TransA and TransB. --- org/qmckl_mo.org | 78 +++++++++++++++++++++--------------------------- 1 file changed, 34 insertions(+), 44 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index cfc4c14..ab311e3 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -492,50 +492,40 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & N = 1_8 K = mo_num * 1_8 - !do iwalk = 1, walk_num - ! do ielec = 1, elec_num - ! ! Value - ! TransA = .True. - ! TransB = .True. - ! info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ! ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & - ! coef_normalized,size(coef_normalized,1) * 1_8, & - ! beta, & - ! mo_vgl(:,ielec,iwalk,1),1_8) - ! ! Grad_x - ! TransA = .True. - ! TransB = .True. - ! info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ! ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & - ! coef_normalized,size(coef_normalized,1) * 1_8, & - ! beta, & - ! mo_vgl(:,ielec,iwalk,2),1_8) - ! ! Grad_y - ! TransA = .True. - ! TransB = .True. - ! info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ! ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & - ! coef_normalized,size(coef_normalized,1) * 1_8, & - ! beta, & - ! mo_vgl(:,ielec,iwalk,3),1_8) - ! ! Grad_z - ! TransA = .True. - ! TransB = .True. - ! info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ! ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & - ! coef_normalized,size(coef_normalized,1) * 1_8, & - ! beta, & - ! mo_vgl(:,ielec,iwalk,4),1_8) - ! ! Lapl_z - ! TransA = .True. - ! TransB = .True. - ! info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & - ! ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & - ! coef_normalized,size(coef_normalized,1) * 1_8, & - ! beta, & - ! mo_vgl(:,ielec,iwalk,5),1_8) - ! end do - !end do + do iwalk = 1, walk_num + do ielec = 1, elec_num + ! Value + info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & + coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,1),1_8) + ! Grad_x + info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & + coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,2),1_8) + ! Grad_y + info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & + coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,3),1_8) + ! Grad_z + info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & + coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,4),1_8) + ! Lapl_z + info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & + coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,5),1_8) + end do + end do if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. & info_qmckl_dgemm_Gx .eq. QMCKL_SUCCESS .and. & From 25bb2fd97eb0673c696f095616f35d286731b2db Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 23:36:11 +0200 Subject: [PATCH 16/23] Fixed dimensions in call to DGEMM. #32 --- org/qmckl_blas.org | 2 +- org/qmckl_mo.org | 30 ++++++++++++++++-------------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 4bb7384..9cc1df1 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -252,7 +252,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) TransA = .False. TransB = .False. - m = 5_8 + m = 1_8 k = 4_8 n = 6_8 LDA = m diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index ab311e3..66c7b33 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -489,41 +489,44 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & ! TODO : Use numerical precision here cutoff = -dlog(1.d-15) M = 1_8 - N = 1_8 - K = mo_num * 1_8 + N = mo_num * 1_8 + K = ao_num * 1_8 + LDA = M + LDB = K + LDC = M do iwalk = 1, walk_num do ielec = 1, elec_num ! Value info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 1), LDA, & coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,1),1_8) + mo_vgl(:,ielec,iwalk,1),LDC) ! Grad_x info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 2), LDA, & coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,2),1_8) + mo_vgl(:,ielec,iwalk,2),LDC) ! Grad_y info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 3), LDA, & coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,3),1_8) + mo_vgl(:,ielec,iwalk,3),LDC) ! Grad_z info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 4), LDA, & coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,4),1_8) + mo_vgl(:,ielec,iwalk,4),LDC) ! Lapl_z info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 5), LDA, & coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,ielec,iwalk,5),1_8) + mo_vgl(:,ielec,iwalk,5),LDC) end do end do @@ -782,8 +785,7 @@ assert(qmckl_mo_basis_provided(context)); double mo_vgl[5][walk_num][elec_num][chbrclf_ao_num]; rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); -//assert (rc == QMCKL_SUCCESS); - +assert (rc == QMCKL_SUCCESS); printf("\n"); printf(" mo_vgl mo_vgl[0][0][26][219] %25.15e\n", mo_vgl[0][0][2][3]); From 296a2886ff77b535e025b68a98a6b28b4b62a16a Mon Sep 17 00:00:00 2001 From: v1j4y Date: Wed, 29 Sep 2021 18:12:41 +0200 Subject: [PATCH 17/23] Added reading actual mo coef for testing. #32 --- org/qmckl_mo.org | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 66c7b33..b923f24 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -768,17 +768,17 @@ assert (rc == QMCKL_SUCCESS); rc = qmckl_set_mo_basis_type(context, typ); assert (rc == QMCKL_SUCCESS); -const int64_t mo_num = chbrclf_ao_num; +const int64_t mo_num = chbrclf_mo_num; rc = qmckl_set_mo_basis_mo_num(context, mo_num); assert (rc == QMCKL_SUCCESS); -double mo_coefficient[chbrclf_ao_num][mo_num]; -const int64_t min_num = (((chbrclf_ao_num) < (mo_num)) ? (chbrclf_ao_num) : (mo_num)); -for (int i=0; i Date: Wed, 29 Sep 2021 18:55:02 +0200 Subject: [PATCH 18/23] Fixed bug in mo_dims. --- org/qmckl_mo.org | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index b923f24..571c5ff 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -773,26 +773,32 @@ rc = qmckl_set_mo_basis_mo_num(context, mo_num); assert (rc == QMCKL_SUCCESS); const double * mo_coefficient = &(chbrclf_mo_coef[0]); -//const int64_t min_num = (((chbrclf_ao_num) < (mo_num)) ? (chbrclf_ao_num) : (mo_num)); -//for (int i=0; i Date: Fri, 1 Oct 2021 00:23:07 +0200 Subject: [PATCH 19/23] Fixed bug in mo_coef dimensions. --- org/qmckl_mo.org | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 571c5ff..f6de762 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -469,7 +469,7 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & integer :: info_qmckl_dgemm_Gy integer :: info_qmckl_dgemm_Gz integer :: info_qmckl_dgemm_lap - integer*8 :: M, N, K, LDA, LDB, LDC + integer*8 :: M, N, K, LDA, LDB, LDC, i,j integer*8 :: inucl, iprim, iwalk, ielec, ishell double precision :: x, y, z, two_a, ar2, r2, v, cutoff @@ -500,31 +500,31 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & ! Value info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & ao_vgl(:, ielec, iwalk, 1), LDA, & - coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & beta, & mo_vgl(:,ielec,iwalk,1),LDC) ! Grad_x info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & ao_vgl(:, ielec, iwalk, 2), LDA, & - coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & beta, & mo_vgl(:,ielec,iwalk,2),LDC) ! Grad_y info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & ao_vgl(:, ielec, iwalk, 3), LDA, & - coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & beta, & mo_vgl(:,ielec,iwalk,3),LDC) ! Grad_z info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & ao_vgl(:, ielec, iwalk, 4), LDA, & - coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & beta, & mo_vgl(:,ielec,iwalk,4),LDC) ! Lapl_z info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & ao_vgl(:, ielec, iwalk, 5), LDA, & - coef_normalized(1:mo_num,1:ao_num),size(coef_normalized,1) * 1_8, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & beta, & mo_vgl(:,ielec,iwalk,5),LDC) end do @@ -783,7 +783,7 @@ double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num]; rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); assert (rc == QMCKL_SUCCESS); -// Overlap of MOs +//// Overlap of MOs //double ovl[mo_num][mo_num]; //for (int i=0; i Date: Mon, 4 Oct 2021 09:08:38 +0200 Subject: [PATCH 20/23] Added test for MO but deactivated because slow. #32 --- org/qmckl_mo.org | 52 +++++++++++++++++++++++++++++++++--------------- 1 file changed, 36 insertions(+), 16 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index f6de762..79d031f 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -120,7 +120,7 @@ typedef struct qmckl_mo_basis_struct { this mechanism. ** Access functions - + #+begin_src c :comments org :tangle (eval h_private_func) :exports none char qmckl_get_mo_basis_type (const qmckl_context context); int64_t qmckl_get_mo_basis_mo_num (const qmckl_context context); @@ -277,7 +277,6 @@ qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const dou qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double); double* new_array = (double*) qmckl_malloc(context, mem_info); - if (new_array == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, @@ -783,23 +782,44 @@ double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num]; rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); assert (rc == QMCKL_SUCCESS); -//// Overlap of MOs -//double ovl[mo_num][mo_num]; -//for (int i=0; i Date: Mon, 4 Oct 2021 09:50:16 +0200 Subject: [PATCH 21/23] Fixed test. #32 --- org/qmckl_mo.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 79d031f..a884ed6 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -815,7 +815,7 @@ assert (rc == QMCKL_SUCCESS); // double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num]; // rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); // assert (rc == QMCKL_SUCCESS); -// ovlmo1 += mo_vgl[0][0][0][0]; +// ovlmo1 += mo_vgl[0][0][0][0]*mo_vgl[0][0][0][0]*dr3; // } // } //} From f3a3c8f0b159000719bb18c3d819febc5563eb49 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 4 Oct 2021 15:43:41 +0200 Subject: [PATCH 22/23] Improved doc for context. #32 --- org/qmckl_mo.org | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index a884ed6..331be69 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -80,13 +80,8 @@ int main() { |---------------------+--------------------+--------------------------------------------------------------| | ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) | - | ~nucleus_index~ | ~[nucl_num]~ | Index of the first shell of each nucleus | - | ~nucleus_shell_num~ | ~[nucl_num]~ | Number of shells per nucleus | - | ~ao_num~ | | Number of AOs | - | ~ao_cartesian~ | | If true, use polynomials. Otherwise, use spherical harmonics | - | ~ao_factor~ | ~[ao_num]~ | Normalization factor of the AO | - | ~ao_shell~ | ~[ao_num]~ | For each AO, specify to which shell it belongs | - | ~coefficient~ | ~[mo_num, ao_num]~ | Orbital coefficients | + | ~mo_num~ | | Number of MOs | + | ~coefficient~ | ~[ao_num, mo_num]~ | Orbital coefficients | Computed data: From 5746487c2e021df80008168c4ee75d2137cbd4e5 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 4 Oct 2021 15:46:47 +0200 Subject: [PATCH 23/23] Fixed ordering in Doc. #32 --- org/qmckl_mo.org | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 331be69..f521b4d 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -81,7 +81,7 @@ int main() { |---------------------+--------------------+--------------------------------------------------------------| | ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) | | ~mo_num~ | | Number of MOs | - | ~coefficient~ | ~[ao_num, mo_num]~ | Orbital coefficients | + | ~coefficient~ | ~[mo_num, ao_num]~ | Orbital coefficients | Computed data: