diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 3d59d7b..279fd80 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -22,7 +22,7 @@ int main() { * Matrix operations ** ~qmckl_dgemm~ - + Matrix multiply: $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. @@ -55,7 +55,7 @@ int main() { - ~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") @@ -79,7 +79,6 @@ int main() { 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) & @@ -91,21 +90,22 @@ 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,k) + real*8 , intent(in) :: A(lda,*) integer*8 , intent(in) :: ldb - real*8 , intent(in) :: B(k,n) + real*8 , intent(in) :: B(ldb,*) integer*8 , intent(in) :: ldc - real*8 , intent(out) :: C(m,n) + real*8 , intent(out) :: C(ldc,*) real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:) + integer*4 :: qmckl_dgemm_N_N_f 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 + allocate(AT(m,k)) + do i = 1, k + do j = 1, m AT(j,i) = A(i,j) end do end do @@ -115,9 +115,9 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, endif if (TransB) then - allocate(BT(n,k)) - do i = 1, k - do j = 1, n + allocate(BT(k,n)) + do i = 1, n + do j = 1, k BT(j,i) = B(i,j) end do end do @@ -162,27 +162,77 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, endif if (TransA) then - if (alpha == 1.d0 .and. beta == 0.d0) then - C = matmul(AT,B) - else - C = beta*C + alpha*matmul(AT,B) - endif + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, B, LDB_2, beta, c, LDC) else if (TransB) then - if (alpha == 1.d0 .and. beta == 0.d0) then - C = matmul(A,BT) - else - C = beta*C + alpha*matmul(A,BT) - endif + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, BT, LDB_2, beta, c, LDC) + else if (TransA .and. TransB) then + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, BT, LDB_2, beta, c, LDC) else - if (alpha == 1.d0 .and. beta == 0.d0) then - C = matmul(A,B) - else - C = beta*C + alpha*matmul(A,B) - endif + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC) endif end function qmckl_dgemm_f - #+end_src +integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8 , intent(in) :: m, n, k + real*8 , intent(in) :: alpha, beta + integer*8 , intent(in) :: lda + real*8 , intent(in) :: A(lda,k) + integer*8 , intent(in) :: ldb + real*8 , intent(in) :: B(ldb,n) + integer*8 , intent(in) :: ldc + real*8 , intent(out) :: C(ldc,n) + + integer*8 :: i,j,l, LDA_2, LDB_2 + + 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 /= k) then + info = QMCKL_INVALID_ARG_10 + return + endif + + if (LDC /= m) then + info = QMCKL_INVALID_ARG_13 + return + endif + + if (alpha == 1.0d0 .and. beta == 0.0d0) then + C = matmul(A,B) + else + C = beta*C + alpha*matmul(A,B) + endif +end function qmckl_dgemm_N_N_f + #+end_src + *** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") @@ -321,6 +371,687 @@ qmckl_exit_code test_qmckl_dgemm(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); #+end_src +** ~qmckl_adjoint~ + + Matrix adjoint. Given a matrix M, returns a matrix M⁻¹ such that: + + + \[ +M · M^{-1} = I +\] + + This is a native Fortran implementation hand written (by: A. Scemama) + only for small matrices (<=5x5). + + TODO: Add description about the external library dependence. + + #+NAME: qmckl_adjoint_args + | qmckl_context | context | in | Global state | + | int64_t | m | in | Number of rows of the input matrix | + | int64_t | n | in | Number of columns of the input matrix | + | int64_t | lda | in | Leading dimension of array ~A~ | + | double | A[][lda] | inout | Array containing the $m \times n$ matrix $A$ | + | double | det_l | inout | determinant of A | + +*** Requirements + + - ~context~ is not ~QMCKL_NULL_CONTEXT~ + - ~m > 0~ + - ~n > 0~ + - ~lda >= m~ + - ~A~ is allocated with at least $m \times n \times 8$ bytes + +*** C header + + #+CALL: generate_c_header(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint") + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_adjoint ( + const qmckl_context context, + const int64_t m, + const int64_t n, + const int64_t lda, + double* A, + double det_l ); + #+end_src + +*** Source + #+begin_src f90 :tangle (eval f) +integer function qmckl_adjoint_f(context, ma, na, LDA, A, det_l) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + double precision, intent(inout) :: A (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: ma + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + integer :: i,j + + info = QMCKL_SUCCESS + + select case (na) + case default +!DIR$ forceinline + print *," TODO: Implement general adjoint" + stop 0 + case (5) +!DIR$ forceinline + call adjoint5(a,LDA,na,det_l) + case (4) +!DIR$ forceinline + call adjoint4(a,LDA,na,det_l) + + case (3) +!DIR$ forceinline + call adjoint3(a,LDA,na,det_l) + case (2) +!DIR$ forceinline + call adjoint2(a,LDA,na,det_l) + case (1) +!DIR$ forceinline + call adjoint1(a,LDA,na,det_l) + case (0) + det_l=1.d0 + end select +end function qmckl_adjoint_f + +subroutine adjoint1(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + + call cofactor1(a,LDA,na,det_l) +end + +subroutine adjoint2(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + double precision :: b(2,2) + + call cofactor2(a,LDA,na,det_l) + + ! Calculate the transpose + b(1,1) = a(1,1) + b(1,2) = a(2,1) + b(2,1) = a(1,2) + b(2,2) = a(2,2) + a(1,1) = b(1,1) + a(1,2) = b(1,2) + a(2,1) = b(2,1) + a(2,2) = b(2,2) +end + +subroutine adjoint3(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + double precision :: b(3,3) + + call cofactor3(a,LDA,na,det_l) + + ! Calculate the transpose + b(1,1) = a(1,1) + b(1,2) = a(2,1) + b(1,3) = a(3,1) + b(2,1) = a(1,2) + b(2,2) = a(2,2) + b(2,3) = a(3,2) + b(3,1) = a(1,3) + b(3,2) = a(2,3) + b(3,3) = a(3,3) + ! copy + a(1,1) = b(1,1) + a(2,1) = b(2,1) + a(3,1) = b(3,1) + a(1,2) = b(1,2) + a(2,2) = b(2,2) + a(3,2) = b(3,2) + a(1,3) = b(1,3) + a(2,3) = b(2,3) + a(3,3) = b(3,3) +end + +subroutine adjoint4(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + double precision :: b(4,4) + + call cofactor4(a,LDA,na,det_l) + + ! Calculate the transpose + b(1,1) = a(1,1) + b(1,2) = a(2,1) + b(1,3) = a(3,1) + b(1,4) = a(4,1) + b(2,1) = a(1,2) + b(2,2) = a(2,2) + b(2,3) = a(3,2) + b(2,4) = a(4,2) + b(3,1) = a(1,3) + b(3,2) = a(2,3) + b(3,3) = a(3,3) + b(3,4) = a(4,3) + b(4,1) = a(1,4) + b(4,2) = a(2,4) + b(4,3) = a(3,4) + b(4,4) = a(4,4) + ! copy + a(1,1) = b(1,1) + a(2,1) = b(2,1) + a(3,1) = b(3,1) + a(4,1) = b(4,1) + a(1,2) = b(1,2) + a(2,2) = b(2,2) + a(3,2) = b(3,2) + a(4,2) = b(4,2) + a(1,3) = b(1,3) + a(2,3) = b(2,3) + a(3,3) = b(3,3) + a(4,3) = b(4,3) + a(1,4) = b(1,4) + a(2,4) = b(2,4) + a(3,4) = b(3,4) + a(4,4) = b(4,4) +end + +subroutine adjoint5(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + double precision :: b(5,5) + + call cofactor5(a,LDA,na,det_l) + + ! Calculate the transpose + b(1,1) = a(1,1) + b(1,2) = a(2,1) + b(1,3) = a(3,1) + b(1,4) = a(4,1) + b(1,5) = a(5,1) + b(2,1) = a(1,2) + b(2,2) = a(2,2) + b(2,3) = a(3,2) + b(2,4) = a(4,2) + b(2,5) = a(5,2) + b(3,1) = a(1,3) + b(3,2) = a(2,3) + b(3,3) = a(3,3) + b(3,4) = a(4,3) + b(3,5) = a(5,3) + b(4,1) = a(1,4) + b(4,2) = a(2,4) + b(4,3) = a(3,4) + b(4,4) = a(4,4) + b(4,5) = a(5,4) + b(5,1) = a(1,5) + b(5,2) = a(2,5) + b(5,3) = a(3,5) + b(5,4) = a(4,5) + b(5,5) = a(5,5) + ! copy + a(1,1) = b(1,1) + a(2,1) = b(2,1) + a(3,1) = b(3,1) + a(4,1) = b(4,1) + a(5,1) = b(5,1) + a(1,2) = b(1,2) + a(2,2) = b(2,2) + a(3,2) = b(3,2) + a(4,2) = b(4,2) + a(5,2) = b(5,2) + a(1,3) = b(1,3) + a(2,3) = b(2,3) + a(3,3) = b(3,3) + a(4,3) = b(4,3) + a(5,3) = b(5,3) + a(1,4) = b(1,4) + a(2,4) = b(2,4) + a(3,4) = b(3,4) + a(4,4) = b(4,4) + a(5,4) = b(5,4) + a(1,5) = b(1,5) + a(2,5) = b(2,5) + a(3,5) = b(3,5) + a(4,5) = b(4,5) + a(5,5) = b(5,5) +end + +subroutine cofactor1(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + + det_l = a(1,1) + a(1,1) = 1.d0 +end + +subroutine cofactor2(a,LDA,na,det_l) + implicit none + double precision :: a (LDA,na) + integer*8 :: LDA + integer*8 :: na + double precision :: det_l + double precision :: b(2,2) + + b(1,1) = a(1,1) + b(2,1) = a(2,1) + b(1,2) = a(1,2) + b(2,2) = a(2,2) + det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1) + a(1,1) = b(2,2) + a(2,1) = -b(2,1) + a(1,2) = -b(1,2) + a(2,2) = b(1,1) +end + +subroutine cofactor3(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + double precision :: b(4,3) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b + integer :: i + det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) & + -a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) & + +a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1)) + do i=1,4 + b(i,1) = a(i,1) + b(i,2) = a(i,2) + b(i,3) = a(i,3) + enddo + a(1,1) = b(2,2)*b(3,3) - b(2,3)*b(3,2) + a(2,1) = b(2,3)*b(3,1) - b(2,1)*b(3,3) + a(3,1) = b(2,1)*b(3,2) - b(2,2)*b(3,1) + + a(1,2) = b(1,3)*b(3,2) - b(1,2)*b(3,3) + a(2,2) = b(1,1)*b(3,3) - b(1,3)*b(3,1) + a(3,2) = b(1,2)*b(3,1) - b(1,1)*b(3,2) + + a(1,3) = b(1,2)*b(2,3) - b(1,3)*b(2,2) + a(2,3) = b(1,3)*b(2,1) - b(1,1)*b(2,3) + a(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1) + +end + +subroutine cofactor4(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + double precision :: b(4,4) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b + integer :: i,j + det_l = a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) & + -a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) & + +a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) & + -a(1,2)*(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) & + -a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) & + +a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))) & + +a(1,3)*(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) & + -a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) & + +a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) & + -a(1,4)*(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)) & + -a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)) & + +a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) + do i=1,4 + b(1,i) = a(1,i) + b(2,i) = a(2,i) + b(3,i) = a(3,i) + b(4,i) = a(4,i) + enddo + + a(1,1) = b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)) + a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)) + a(3,1) = b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) + a(4,1) = -b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))+b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))-b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) + + a(1,2) = -b(1,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(1,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(1,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)) + a(2,2) = b(1,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(1,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(1,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)) + a(3,2) = -b(1,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(1,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(1,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) + a(4,2) = b(1,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(1,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(1,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) + + a(1,3) = b(1,2)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))-b(1,3)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))+b(1,4)*(b(2,2)*b(4,3)-b(2,3)*b(4,2)) + a(2,3) = -b(1,1)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))+b(1,3)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))-b(1,4)*(b(2,1)*b(4,3)-b(2,3)*b(4,1)) + a(3,3) = b(1,1)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))-b(1,2)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))+b(1,4)*(b(2,1)*b(4,2)-b(2,2)*b(4,1)) + a(4,3) = -b(1,1)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))+b(1,2)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))-b(1,3)*(b(2,1)*b(4,2)-b(2,2)*b(4,1)) + + a(1,4) = -b(1,2)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))+b(1,3)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))-b(1,4)*(b(2,2)*b(3,3)-b(2,3)*b(3,2)) + a(2,4) = b(1,1)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))-b(1,3)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))+b(1,4)*(b(2,1)*b(3,3)-b(2,3)*b(3,1)) + a(3,4) = -b(1,1)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))+b(1,2)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))-b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1)) + a(4,4) = b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))-b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))+b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1)) + +end + +subroutine cofactor5(a,LDA,na,det_l) + implicit none + double precision, intent(inout) :: a (LDA,na) + integer*8, intent(in) :: LDA + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + double precision :: b(5,5) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b + integer :: i,j + det_l = a(1,1)*(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( & + a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))- & + a(2,3)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)- & + a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(2,4)*(a(3,2)*( & + a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+ & + a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,5)*(a(3,2)*(a(4,3)*a(5,4)- & + a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)* & + a(5,3)-a(4,3)*a(5,2))))-a(1,2)*(a(2,1)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)* & + a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)- & + a(4,4)*a(5,3)))-a(2,3)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( & + a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+ & + a(2,4)*(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)- & + a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(2,5)*(a(3,1)*( & + a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+ & + a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))+a(1,3)*(a(2,1)*(a(3,2)*(a(4,4)* & + a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*( & + a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)* & + a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)- & + a(4,4)*a(5,1)))+a(2,4)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*( & + a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))- & + a(2,5)*(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)- & + a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))-a(1,4)*(a(2,1)*( & + a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)* & + a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,3)* & + a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*( & + a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)* & + a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)- & + a(4,2)*a(5,1)))-a(2,5)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*( & + a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))+ & + a(1,5)*(a(2,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)* & + a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*( & + a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)* & + a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)* & + a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*( & + a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,4)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)* & + a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)- & + a(4,2)*a(5,1)))) + + do i=1,5 + b(1,i) = a(1,i) + b(2,i) = a(2,i) + b(3,i) = a(3,i) + b(4,i) = a(4,i) + b(5,i) = a(5,i) + enddo + + a(1,1) = & + (b(2,2)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))-b(2,3)* & + (b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))+b(2,4)* & + (b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(2,5)* & + (b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))) + a(2,1) = & + (-b(2,1)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))+b(2,3)* & + (b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))-b(2,4)* & + (b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(2,5)* & + (b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))) + a(3,1) = & + (b(2,1)*(b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))-b(2,2)* & + (b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))+b(2,4)* & + (b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(2,5)* & + (b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + a(4,1) = & + (-b(2,1)*(b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(2,2)* & + (b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(2,3)* & + (b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(2,5)* & + (b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + a(5,1) = & + (b(2,1)*(b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(2,2)* & + (b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(2,3)* & + (b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(2,4)* & + (b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + + a(1,2) = & + (-b(1,2)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))+b(1,3)* & + (b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))-b(1,4)* & + (b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(1,5)* & + (b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))) + a(2,2) = & + (b(1,1)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))-b(1,3)* & + (b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))+b(1,4)* & + (b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(1,5)* & + (b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))) + a(3,2) = & + (-b(1,1)*(b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))+b(1,2)* & + (b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))-b(1,4)* & + (b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(1,5)* & + (b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + a(4,2) = & + (b(1,1)*(b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(1,2)* & + (b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(1,3)* & + (b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(1,5)* & + (b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + a(5,2) = & + (-b(1,1)*(b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(1,2)* & + (b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(1,3)* & + (b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(1,4)* & + (b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + + a(1,3) = & + (b(1,2)*(b(2,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(2,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))-b(1,3)* & + (b(2,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))+b(1,4)* & + (b(2,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(1,5)* & + (b(2,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(2,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))) + a(2,3) = & + (-b(1,1)*(b(2,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(2,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))+b(1,3)* & + (b(2,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))-b(1,4)* & + (b(2,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(1,5)* & + (b(2,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))) + a(3,3) = & + (b(1,1)*(b(2,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))-b(1,2)* & + (b(2,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))+b(1,4)* & + (b(2,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(2,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(1,5)* & + (b(2,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(2,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + a(4,3) = & + (-b(1,1)*(b(2,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(1,2)* & + (b(2,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(1,3)* & + (b(2,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(2,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(1,5)* & + (b(2,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(2,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(2,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + a(5,3) = & + (b(1,1)*(b(2,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(2,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(1,2)* & + (b(2,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(1,3)* & + (b(2,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(2,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(1,4)* & + (b(2,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(2,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(2,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))) + + a(1,4) = & + (-b(1,2)*(b(2,3)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))+b(2,5)*(b(3,3)*b(5,4)-b(3,4)*b(5,3)))+b(1,3)* & + (b(2,2)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,4)-b(3,4)*b(5,2)))-b(1,4)* & + (b(2,2)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,3)-b(3,3)*b(5,2)))+b(1,5)* & + (b(2,2)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))+b(2,4)*(b(3,2)*b(5,3)-b(3,3)*b(5,2)))) + a(2,4) = & + (b(1,1)*(b(2,3)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))+b(2,5)*(b(3,3)*b(5,4)-b(3,4)*b(5,3)))-b(1,3)* & + (b(2,1)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,4)-b(3,4)*b(5,1)))+b(1,4)* & + (b(2,1)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,3)-b(3,3)*b(5,1)))-b(1,5)* & + (b(2,1)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,3)-b(3,3)*b(5,1)))) + a(3,4) = & + (-b(1,1)*(b(2,2)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,4)-b(3,4)*b(5,2)))+b(1,2)* & + (b(2,1)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,4)-b(3,4)*b(5,1)))-b(1,4)* & + (b(2,1)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))-b(2,2)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))+b(1,5)* & + (b(2,1)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))-b(2,2)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))) + a(4,4) = & + (b(1,1)*(b(2,2)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,3)-b(3,3)*b(5,2)))-b(1,2)* & + (b(2,1)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,3)-b(3,3)*b(5,1)))+b(1,3)* & + (b(2,1)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))-b(2,2)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))-b(1,5)* & + (b(2,1)*(b(3,2)*b(5,3)-b(3,3)*b(5,2))-b(2,2)*(b(3,1)*b(5,3)-b(3,3)*b(5,1))+b(2,3)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))) + a(5,4) = & + (-b(1,1)*(b(2,2)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))+b(2,4)*(b(3,2)*b(5,3)-b(3,3)*b(5,2)))+b(1,2)* & + (b(2,1)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,3)-b(3,3)*b(5,1)))-b(1,3)* & + (b(2,1)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))-b(2,2)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))+b(1,4)* & + (b(2,1)*(b(3,2)*b(5,3)-b(3,3)*b(5,2))-b(2,2)*(b(3,1)*b(5,3)-b(3,3)*b(5,1))+b(2,3)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))) + + a(1,5) = & + (b(1,2)*(b(2,3)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))+b(2,5)*(b(3,3)*b(4,4)-b(3,4)*b(4,3)))-b(1,3)* & + (b(2,2)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,4)-b(3,4)*b(4,2)))+b(1,4)* & + (b(2,2)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)))-b(1,5)* & + (b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)))) + a(2,5) = & + (-b(1,1)*(b(2,3)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))+b(2,5)*(b(3,3)*b(4,4)-b(3,4)*b(4,3)))+b(1,3)* & + (b(2,1)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,4)-b(3,4)*b(4,1)))-b(1,4)* & + (b(2,1)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)))+b(1,5)* & + (b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)))) + a(3,5) = & + (b(1,1)*(b(2,2)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,4)-b(3,4)*b(4,2)))-b(1,2)* & + (b(2,1)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,4)-b(3,4)*b(4,1)))+b(1,4)* & + (b(2,1)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))-b(2,2)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))-b(1,5)* & + (b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))) + a(4,5) = & + (-b(1,1)*(b(2,2)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)))+b(1,2)* & + (b(2,1)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)))-b(1,3)* & + (b(2,1)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))-b(2,2)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))+b(1,5)* & + (b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))) + a(5,5) = & + (b(1,1)*(b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)))-b(1,2)* & + (b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)))+b(1,3)* & + (b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))-b(1,4)* & + (b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))) + +end + #+end_src + +*** C interface :noexport: + + #+CALL: generate_c_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint") + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_adjoint & + (context, m, n, lda, A, det_l) & + 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 :: m + integer (c_int64_t) , intent(in) , value :: n + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(inout) :: A(lda,*) + real (c_double ) , intent(inout) :: det_l + + integer(c_int32_t), external :: qmckl_adjoint_f + info = qmckl_adjoint_f & + (context, m, n, lda, A, det_l) + + end function qmckl_adjoint + #+end_src + + #+CALL: generate_f_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint") + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_adjoint & + (context, m, n, lda, A, det_l) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(inout) :: A(lda,*) + real (c_double ) , intent(inout) :: det_l + + end function qmckl_adjoint + end interface + #+end_src + +*** Test :noexport: + #+begin_src f90 :tangle (eval f_test) +integer(qmckl_exit_code) function test_qmckl_adjoint(context) bind(C) + use qmckl + implicit none + integer(qmckl_context), intent(in), value :: context + + double precision, allocatable :: A(:,:), C(:,:) + integer*8 :: m, n, k, LDA, LDB, LDC + integer*8 :: i,j,l + double precision :: x, det_l, det_l_ref + + m = 4_8 + k = 4_8 + LDA = m + LDB = m + LDC = m + + allocate( A(LDA,k), C(LDC,k)) + + A = 0.10d0 + C = 0.d0 + A(1,1) = 1.0d0; + A(2,2) = 2.0d0; + A(3,3) = 3.0d0; + A(4,4) = 4.0d0; + + ! Exact inverse (Mathematica) + C(1,1) = 1.0102367161391992d0 + C(2,2) = 0.5036819224578257d0 + C(3,3) = 0.33511197860555897d0 + C(4,4) = 0.2510382472105688d0 + C(1,2) = -0.047782608144589914d0 + C(1,3) = -0.031305846715420985d0 + C(1,4) = -0.023278706531979707d0 + C(2,3) = -0.014829085286252043d0 + C(2,4) = -0.011026755725674596d0 + C(3,4) = -0.007224426165097149d0 + C(2,1) = -0.047782608144589914d0 + C(3,1) = -0.031305846715420985d0 + C(4,1) = -0.023278706531979707d0 + C(3,2) = -0.014829085286252043d0 + C(4,2) = -0.011026755725674596d0 + C(4,3) = -0.007224426165097149d0 + det_l_ref = 23.6697d0 + + test_qmckl_adjoint = qmckl_adjoint(context, m, k, LDA, A, det_l) + + if (test_qmckl_adjoint /= QMCKL_SUCCESS) return + + test_qmckl_adjoint = QMCKL_FAILURE + + x = 0.d0 + do j=1,m + do i=1,k + x = x + (A(i,j) - (C(i,j) * det_l_ref))**2 + end do + end do + + if (dabs(x) <= 1.d-15 .and. (dabs(det_l_ref - det_l)) <= 1.d-15) then + test_qmckl_adjoint = QMCKL_SUCCESS + endif + + deallocate(A,C) +end function test_qmckl_adjoint + #+end_src + + #+begin_src c :comments link :tangle (eval c_test) +qmckl_exit_code test_qmckl_adjoint(qmckl_context context); +assert(QMCKL_SUCCESS == test_qmckl_adjoint(context)); + #+end_src + * End of files :noexport: #+begin_src c :comments link :tangle (eval c_test) diff --git a/org/qmckl_context.org b/org/qmckl_context.org index bd1a0ce..7b26a7f 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -33,10 +33,14 @@ int main() { #include "qmckl_ao_private_type.h" #include "qmckl_mo_private_type.h" #include "qmckl_jastrow_private_type.h" +#include "qmckl_determinant_private_type.h" +#include "qmckl_local_energy_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" +#include "qmckl_determinant_private_func.h" +#include "qmckl_local_energy_private_func.h" #+end_src #+begin_src c :tangle (eval c) @@ -118,15 +122,15 @@ typedef struct qmckl_context_struct { uint64_t date; /* -- Molecular system -- */ - qmckl_nucleus_struct nucleus; - qmckl_electron_struct electron; - qmckl_ao_basis_struct ao_basis; - qmckl_mo_basis_struct mo_basis; - qmckl_jastrow_struct jastrow; + qmckl_nucleus_struct nucleus; + qmckl_electron_struct electron; + qmckl_ao_basis_struct ao_basis; + qmckl_mo_basis_struct mo_basis; + qmckl_jastrow_struct jastrow; + qmckl_determinant_struct det; + qmckl_local_energy_struct local_energy; /* To be implemented: - qmckl_mo_struct mo; - qmckl_determinant_struct det; ,*/ } qmckl_context_struct; @@ -240,6 +244,12 @@ qmckl_context qmckl_context_create() { rc = qmckl_init_ao_basis(context); assert (rc == QMCKL_SUCCESS); + + rc = qmckl_init_mo_basis(context); + assert (rc == QMCKL_SUCCESS); + + rc = qmckl_init_determinant(context); + assert (rc == QMCKL_SUCCESS); } /* Allocate qmckl_memory_struct */ diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org new file mode 100644 index 0000000..1e8e559 --- /dev/null +++ b/org/qmckl_determinant.org @@ -0,0 +1,2073 @@ +#+TITLE: Slater Determinant +#+SETUPFILE: ../tools/theme.setup +#+INCLUDE: ../tools/lib.org + +The slater deteminant is required for the calculation of the +wavefunction, gradient, and derivatives. These quantities will be used +to calculate the local Energy (\[E_L\]). + +ψ(x) = det|ϕ₁(x₁)...ϕᵢ(yᵢ)...ϕₙ(xₙ)| + +The above slater-matrix is also required and is denoted by Dᵢⱼ(x) such that: + +ψ(x) = det|Dᵢⱼ(x)| + +We also require the inverse of the slater-matrix which is denoted by D⁻¹ᵢⱼ(x). +Using this notation, the acceptance probability which is proportional to +ψ(y)/ψ(x) can be calculated as follows: + +ψ(yᵢ)/ψ(xᵢ) = ∑ⱼDᵢⱼ(y)D⁻¹ⱼᵢ(x) + +Concerning the gradient and laplacian, in fact what is actually +calculated is the ratio of the gradient/laplacian and the determinant +of the slater matrix: + +∇ψ(x)/ψ(x) + +and + +∇²ψ(x)/ψ(x) + +This avoids the unnecessary multiplication and division of by the +determinant ψ(x). + + +* 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_DETERMINANT_HPT +#define QMCKL_DETERMINANT_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" +#include "qmckl_determinant_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" +#include "qmckl_determinant_private_type.h" +#include "qmckl_determinant_private_func.h" + #+end_src + +* Context + + The following arrays are stored in the context: + + |------------------+------------------------------------------------+------------------------------------| + | ~type~ | ~char~ | α (~'A'~) or β (~'B'~) determinant | + | ~walk_num~ | ~int64_t~ | Number of walkers | + | ~det_num_alpha~ | ~int64_t~ | Number of determinants per walker | + | ~det_num_beta~ | ~int64_t~ | Number of determinants per walker | + | ~mo_index_alpha~ | ~mo_index[det_num_alpha][walk_num][alpha_num]~ | Index of MOs for each walker | + | ~mo_index_beta~ | ~mo_index[det_num_beta][walk_num][beta_num]~ | Index of MOs for each walker | + + Computed data: + + |-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------| + | ~up_num~ | ~int64_t~ | Number of number of α electrons | + | ~donwn_num~ | ~int64_t~ | Number of number of β electrons | + | ~det_value_alpha~ | ~[det_num_alpha][walk_num]~ | The α slater matrix for each determinant of each walker. | + | ~det_value_alpha_date~ | ~int64_t~ | Date of The α slater matrix for each determinant of each walker. | + | ~det_value_beta~ | ~[det_num_beta][walk_num]~ | The β slater matrix for each determinant of each walker. | + | ~det_value_beta_date~ | ~int64_t~ | Date of The β slater matrix for each determinant of each walker. | + | ~det_adj_matrix_alpha~ | ~[det_num_alpha][walk_num][alpha_num][alpha_num]~ | Adjoint of the α slater matrix for each determinant of each walker. | + | ~det_adj_matrix_alpha_date~ | ~int64_t~ | Date of the Adjoint of the α slater matrix for each determinant of each walker. | + | ~det_adj_matrix_beta~ | ~[det_num_beta][walk_num][beta_num][beta_num]~ | Adjoint of the β slater matrix for each determinant of each walker. | + | ~det_adj_matrix_beta_date~ | ~int64_t~ | Date of the Adjoint of the β slater matrix for each determinant of each walker. | + |-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------| + | ~det_vgl_alpha~ | ~[5][det_num_alpha][walk_num][alpha_num][alpha_num]~ | Value, gradients, Laplacian of Dᵅᵢⱼ(x) at electron positions | + | ~det_vgl_alpha_date~ | ~int64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions | + | ~det_vgl_beta~ | ~[5][det_num_beta][walk_num][beta_num][beta_num]~ | Value, gradients, Laplacian of Dᵝᵢⱼ(x) at electron positions | + | ~det_vgl_beta_date~ | ~int64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions | + | ~det_inv_matrix_alpha~ | ~[det_num_alpha][walk_num][alpha_num][alpha_num]~ | Inverse of the α electron slater matrix for each determinant of each walker. | + | ~det_inv_matrix_alpha_date~ | ~int64_t~ | Date for the Inverse of the α electron slater matrix for each determinant of each walker. | + | ~det_inv_matrix_beta~ | ~[det_num_beta][walk_num][beta_num][beta_num]~ | Inverse of the β electron slater matrix for each determinant of each walker. | + | ~det_inv_matrix_beta_date~ | ~int64_t~ | Date for the Inverse of the β electron slater matrix for each determinant of each walker. | + |-----------------------------+------------------------------------------------------+-------------------------------------------------------------------------------------------| + +** Data structure + + #+begin_src c :comments org :tangle (eval h_private_type) +typedef struct qmckl_determinant_struct { + char type; + int64_t walk_num; + int64_t det_num_alpha; + int64_t det_num_beta ; + int64_t up_num; + int64_t down_num; + int64_t* mo_index_alpha; + int64_t* mo_index_beta; + + double * det_value_alpha; + double * det_value_beta; + double * det_vgl_alpha; + double * det_adj_matrix_alpha; + double * det_inv_matrix_alpha; + double * det_vgl_beta; + double * det_adj_matrix_beta; + double * det_inv_matrix_beta; + int64_t det_value_alpha_date; + int64_t det_vgl_alpha_date; + int64_t det_adj_matrix_alpha_date; + int64_t det_inv_matrix_alpha_date; + int64_t det_value_beta_date; + int64_t det_vgl_beta_date; + int64_t det_adj_matrix_beta_date; + int64_t det_inv_matrix_beta_date; + + int32_t uninitialized; + bool provided; +} qmckl_determinant_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. + + #+begin_src c :comments org :tangle (eval h_private_func) +qmckl_exit_code qmckl_init_determinant(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) +qmckl_exit_code qmckl_init_determinant(qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return false; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + ctx->det.uninitialized = (1 << 6) - 1; + + return QMCKL_SUCCESS; +} + #+end_src + +** Access functions + + #+begin_src c :comments org :tangle (eval h_private_func) :exports none +char qmckl_get_determinant_type (const qmckl_context context); +int64_t qmckl_get_determinant_walk_num (const qmckl_context context); +int64_t qmckl_get_determinant_det_num_alpha (const qmckl_context context); +int64_t qmckl_get_determinant_det_num_beta (const qmckl_context context); +int64_t* qmckl_get_determinant_mo_index_alpha (const qmckl_context context); +int64_t* qmckl_get_determinant_mo_index_beta (const qmckl_context context); + #+end_src + + When all the data for the slater determinants have been provided, the following + function returns ~true~. + + #+begin_src c :comments org :tangle (eval h_func) +bool qmckl_determinant_provided (const qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +bool qmckl_determinant_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->det.provided; +} + #+end_src + + #+NAME:post + #+begin_src c :exports none +if ( (ctx->det.uninitialized & mask) != 0) { + return NULL; +} + #+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +char qmckl_get_determinant_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->det.uninitialized & mask) != 0) { + return (char) 0; + } + + assert (ctx->det.type != (char) 0); + return ctx->det.type; +} + +int64_t qmckl_get_determinant_walk_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->det.uninitialized & mask) != 0) { + return (int64_t) 0; + } + + assert (ctx->det.walk_num > (int64_t) 0); + return ctx->det.walk_num; +} + +int64_t qmckl_get_determinant_det_num_alpha (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 << 2; + + if ( (ctx->det.uninitialized & mask) != 0) { + return (int64_t) 0; + } + + assert (ctx->det.det_num_alpha > (int64_t) 0); + return ctx->det.det_num_alpha; +} + +int64_t qmckl_get_determinant_det_num_beta (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 << 3; + + if ( (ctx->det.uninitialized & mask) != 0) { + return (int64_t) 0; + } + + assert (ctx->det.det_num_beta > (int64_t) 0); + return ctx->det.det_num_beta; +} + +int64_t* qmckl_get_determinant_mo_index_alpha (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 << 4; + + if ( (ctx->det.uninitialized & mask) != 0) { + return (int64_t) 0; + } + + assert (ctx->det.mo_index_alpha != NULL); + return ctx->det.mo_index_alpha; +} + +int64_t* qmckl_get_determinant_mo_index_beta (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 << 5; + + if ( (ctx->det.uninitialized & mask) != 0) { + return (int64_t) 0; + } + + assert (ctx->det.mo_index_beta != NULL); + return ctx->det.mo_index_beta; +} + + #+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_determinant_type (const qmckl_context context, const char t); +qmckl_exit_code qmckl_set_determinant_walk_num (const qmckl_context context, const int64_t walk_num); +qmckl_exit_code qmckl_set_determinant_det_num_alpha (const qmckl_context context, const int64_t det_num_alpha); +qmckl_exit_code qmckl_set_determinant_det_num_beta (const qmckl_context context, const int64_t det_num_beta); +qmckl_exit_code qmckl_set_determinant_mo_index_alpha (const qmckl_context context, const int64_t* mo_index_alpha); +qmckl_exit_code qmckl_set_determinant_mo_index_beta (const qmckl_context context, const int64_t* mo_index_beta); + #+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->det.uninitialized &= ~mask; +ctx->det.provided = (ctx->det.uninitialized == 0); +if (ctx->det.provided) { + qmckl_exit_code rc_ = qmckl_finalize_determinant(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_determinant_type(qmckl_context context, const char t) { + <> + + if (t != 'G' && t != 'S') { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_determinant_type", + NULL); + } + + int32_t mask = 1; + ctx->det.type = t; + + <> +} + +qmckl_exit_code qmckl_set_determinant_walk_num(qmckl_context context, const int64_t walk_num) { + <> + + if (walk_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_determinant_walk_num", + "walk_num <= 0"); + } + + int32_t mask = 1 << 1; + ctx->det.walk_num = walk_num; + + <> +} + +qmckl_exit_code qmckl_set_determinant_det_num_alpha(qmckl_context context, const int64_t det_num_alpha) { + <> + + if (det_num_alpha <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_determinant_det_num_alpha", + "det_num_alpha <= 0"); + } + + int32_t mask = 1 << 2; + ctx->det.det_num_alpha = det_num_alpha; + + <> +} + +qmckl_exit_code qmckl_set_determinant_det_num_beta(qmckl_context context, const int64_t det_num_beta) { + <> + + if (det_num_beta <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_determinant_det_num_beta", + "det_num_beta <= 0"); + } + + int32_t mask = 1 << 3; + ctx->det.det_num_beta = det_num_beta; + + <> +} + +qmckl_exit_code qmckl_set_determinant_mo_index_alpha(qmckl_context context, const int64_t* mo_index_alpha) { + <> + + int32_t mask = 1 << 4; + + if (ctx->det.mo_index_alpha != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->det.mo_index_alpha); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_determinant_mo_index_alpha", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha * + ctx->electron.up_num * sizeof(int64_t); + int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info); + if (new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_determinant_mo_index_alpha", + NULL); + } + + memcpy(new_array, mo_index_alpha, mem_info.size); + + ctx->det.mo_index_alpha = new_array; + + <> +} + +qmckl_exit_code qmckl_set_determinant_mo_index_beta(qmckl_context context, const int64_t* mo_index_beta) { + <> + + int32_t mask = 1 << 5; + + if (ctx->det.mo_index_beta != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->det.mo_index_beta); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_determinant_mo_index_beta", + NULL); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta * + ctx->electron.down_num * sizeof(int64_t); + int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info); + if (new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_determinant_mo_index_beta", + NULL); + } + + memcpy(new_array, mo_index_beta, mem_info.size); + + ctx->det.mo_index_beta = new_array; + + <> +} + + #+end_src + + When the basis set is completely entered, other data structures are + computed to accelerate the calculations. + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_determinant(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_finalize_determinant", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + qmckl_exit_code rc; + rc = qmckl_provide_det_vgl_alpha(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_finalize_determinant", + NULL); + } + rc = qmckl_provide_det_vgl_beta(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_finalize_determinant", + NULL); + } + rc = qmckl_provide_det_inv_matrix_alpha(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_finalize_determinant", + NULL); + } + rc = qmckl_provide_det_inv_matrix_beta(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_finalize_determinant", + NULL); + } +} + #+end_src + +** Fortran Interfaces +** Test +* Computation +** Determinant matrix + :PROPERTIES: + :Name: qmckl_compute_det_vgl + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double* const det_vgl_alpha); +qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double* const det_vgl_beta); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_det_vgl_alpha(qmckl_context context, double * const det_vgl_alpha) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_alpha(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->det.det_num_alpha * ctx->det.walk_num * + ctx->electron.up_num * ctx->electron.up_num * sizeof(double); + memcpy(det_vgl_alpha, ctx->det.det_vgl_alpha, sze); + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_det_vgl_beta(qmckl_context context, double * const det_vgl_beta) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_beta(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->det.det_num_beta * ctx->det.walk_num * + ctx->electron.down_num * ctx->electron.down_num * sizeof(double); + memcpy(det_vgl_beta, ctx->det.det_vgl_beta, sze); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context); +qmckl_exit_code qmckl_provide_det_vgl_beta(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_det_vgl_alpha(qmckl_context context) { + + qmckl_exit_code rc; + 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->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_determinant", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->det.det_vgl_alpha_date) { + + /* Allocate array */ + if (ctx->det.det_vgl_alpha == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 5 * ctx->det.walk_num * ctx->det.det_num_alpha * + ctx->electron.up_num * ctx->electron.up_num * sizeof(double); + double* det_vgl_alpha = (double*) qmckl_malloc(context, mem_info); + + if (det_vgl_alpha == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_vgl_alpha", + NULL); + } + ctx->det.det_vgl_alpha = det_vgl_alpha; + } + + if (ctx->det.type == 'G') { + rc = qmckl_compute_det_vgl_alpha(context, + ctx->det.det_num_alpha, + ctx->det.walk_num, + ctx->electron.up_num, + ctx->electron.down_num, + ctx->electron.num, + ctx->det.mo_index_alpha, + ctx->mo_basis.mo_num, + ctx->mo_basis.mo_vgl, + ctx->det.det_vgl_alpha); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_det_vgl_alpha", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->det.det_vgl_alpha_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_provide_det_vgl_beta(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->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->det.det_vgl_beta_date) { + + /* Allocate array */ + if (ctx->det.det_vgl_beta == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 5 * ctx->det.walk_num * ctx->det.det_num_beta * + ctx->electron.down_num * ctx->electron.down_num * sizeof(double); + double* det_vgl_beta = (double*) qmckl_malloc(context, mem_info); + + if (det_vgl_beta == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_vgl_beta", + NULL); + } + ctx->det.det_vgl_beta = det_vgl_beta; + } + + qmckl_exit_code rc; + if (ctx->det.type == 'G') { + rc = qmckl_compute_det_vgl_beta(context, + ctx->det.det_num_beta, + ctx->det.walk_num, + ctx->electron.up_num, + ctx->electron.down_num, + ctx->electron.num, + ctx->det.mo_index_beta, + ctx->mo_basis.mo_num, + ctx->mo_basis.mo_vgl, + ctx->det.det_vgl_beta); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_det_vgl_beta", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->det.det_vgl_beta_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src +*** Compute alpha + :PROPERTIES: + :Name: qmckl_compute_det_vgl_alpha + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_compute_det_vgl_alpha_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~det_num_alpha~ | in | Number of determinants | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~int64_t~ | ~alpha_num~ | in | Number of electrons | + | ~int64_t~ | ~beta_num~ | in | Number of electrons | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~mo_index_alpha[det_num_alpha][walk_num][alpha_num]~ | in | MO indices for electrons | + | ~int64_t~ | ~mo_num~ | in | Number of MOs | + | ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs | + | ~double~ | ~det_vgl_alpha[det_num_alpha][walk_num][5][alpha_num][alpha_num]~ | out | Value, gradients and Laplacian of the Det | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_det_vgl_alpha_f(context, & + det_num_alpha, walk_num, alpha_num, beta_num, elec_num, & + mo_index_alpha, mo_num, mo_vgl, det_vgl_alpha) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: det_num_alpha + integer*8, intent(in) :: walk_num + integer*8, intent(in) :: alpha_num + integer*8, intent(in) :: beta_num + integer*8, intent(in) :: elec_num + integer*8, intent(in) :: mo_num + integer*8, intent(in) :: mo_index_alpha(alpha_num, walk_num, det_num_alpha) + double precision, intent(in) :: mo_vgl(mo_num, elec_num, 5) + double precision, intent(inout) :: det_vgl_alpha(alpha_num, alpha_num, 5, walk_num, det_num_alpha) + integer*8 :: idet, iwalk, ielec, mo_id, imo + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (alpha_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + do idet = 1, det_num_alpha + do iwalk = 1, walk_num + do ielec = 1, alpha_num + do imo = 1, alpha_num + mo_id = mo_index_alpha(imo,iwalk,idet) + ! Value + det_vgl_alpha(imo, ielec, 1, iwalk, idet) = mo_vgl(mo_id, ielec, 1) + + ! Grad_x + det_vgl_alpha(imo, ielec, 2, iwalk, idet) = mo_vgl(mo_id, ielec, 2) + + ! Grad_y + det_vgl_alpha(imo, ielec, 3, iwalk, idet) = mo_vgl(mo_id, ielec, 3) + + ! Grad_z + det_vgl_alpha(imo, ielec, 4, iwalk, idet) = mo_vgl(mo_id, ielec, 4) + + ! Lap + det_vgl_alpha(imo, ielec, 5, iwalk, idet) = mo_vgl(mo_id, ielec, 5) + end do + end do + end do + end do + +end function qmckl_compute_det_vgl_alpha_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_compute_det_vgl_alpha_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_vgl_alpha")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_det_vgl_alpha ( + const qmckl_context context, + const int64_t det_num_alpha, + const int64_t walk_num, + const int64_t alpha_num, + const int64_t beta_num, + const int64_t elec_num, + const int64_t* mo_index_alpha, + const int64_t mo_num, + const double* mo_vgl, + double* const det_vgl_alpha ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_compute_det_vgl_alpha_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_vgl_alpha")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_det_vgl_alpha & + (context, & + det_num_alpha, & + walk_num, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_alpha, & + mo_num, & + mo_vgl, & + det_vgl_alpha) & + 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 :: det_num_alpha + integer (c_int64_t) , intent(in) , value :: walk_num + integer (c_int64_t) , intent(in) , value :: alpha_num + integer (c_int64_t) , intent(in) , value :: beta_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) :: mo_index_alpha(alpha_num,walk_num,det_num_alpha) + integer (c_int64_t) , intent(in) , value :: mo_num + real (c_double ) , intent(in) :: mo_vgl(mo_num,elec_num,5) + real (c_double ) , intent(out) :: det_vgl_alpha(alpha_num,alpha_num,5,walk_num,det_num_alpha) + + integer(c_int32_t), external :: qmckl_compute_det_vgl_alpha_f + info = qmckl_compute_det_vgl_alpha_f & + (context, & + det_num_alpha, & + walk_num, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_alpha, & + mo_num, & + mo_vgl, & + det_vgl_alpha) + + end function qmckl_compute_det_vgl_alpha + #+end_src + +*** Compute beta + :PROPERTIES: + :Name: qmckl_compute_det_vgl_beta + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_compute_det_vgl_beta_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~det_num_beta~ | in | Number of determinants | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~int64_t~ | ~alpha_num~ | in | Number of electrons | + | ~int64_t~ | ~beta_num~ | in | Number of electrons | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~mo_index_beta[det_num_beta][walk_num][beta_num]~ | in | Number of electrons | + | ~int64_t~ | ~mo_num~ | in | Number of MOs | + | ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs | + | ~double~ | ~det_vgl_beta[det_num_beta][walk_num][5][beta_num][beta_num]~ | out | Value, gradients and Laplacian of the Det | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_det_vgl_beta_f(context, & + det_num_beta, walk_num, alpha_num, beta_num, elec_num, & + mo_index_beta, mo_num, mo_vgl, det_vgl_beta) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: det_num_beta + integer*8, intent(in) :: walk_num + integer*8, intent(in) :: alpha_num + integer*8, intent(in) :: beta_num + integer*8, intent(in) :: elec_num + integer*8, intent(in) :: mo_num + integer*8, intent(in) :: mo_index_beta(beta_num,walk_num,det_num_beta) + double precision, intent(in) :: mo_vgl(mo_num, elec_num, 5) + double precision, intent(inout) :: det_vgl_beta(beta_num, beta_num, 5, walk_num, det_num_beta) + integer*8 :: idet, iwalk, ielec, mo_id, imo + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (beta_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + do idet = 1, det_num_beta + do iwalk = 1, walk_num + do ielec = 1, beta_num + do imo = 1, beta_num + mo_id = mo_index_beta(imo, iwalk, idet) + ! Value + det_vgl_beta(imo, ielec, 1, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 1) + + ! Grad_x + det_vgl_beta(imo, ielec, 2, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 2) + + ! Grad_y + det_vgl_beta(imo, ielec, 3, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 3) + + ! Grad_z + det_vgl_beta(imo, ielec, 4, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 4) + + ! Lap + det_vgl_beta(imo, ielec, 5, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 5) + end do + end do + end do + end do + +end function qmckl_compute_det_vgl_beta_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_compute_det_vgl_beta_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_vgl_beta")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_det_vgl_beta ( + const qmckl_context context, + const int64_t det_num_beta, + const int64_t walk_num, + const int64_t alpha_num, + const int64_t beta_num, + const int64_t elec_num, + const int64_t* mo_index_beta, + const int64_t mo_num, + const double* mo_vgl, + double* const det_vgl_beta ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_compute_det_vgl_beta_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_vgl_beta")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_det_vgl_beta & + (context, & + det_num_beta, & + walk_num, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_beta, & + mo_num, & + mo_vgl, & + det_vgl_beta) & + 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 :: det_num_beta + integer (c_int64_t) , intent(in) , value :: walk_num + integer (c_int64_t) , intent(in) , value :: alpha_num + integer (c_int64_t) , intent(in) , value :: beta_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) :: mo_index_beta(beta_num,walk_num,det_num_beta) + integer (c_int64_t) , intent(in) , value :: mo_num + real (c_double ) , intent(in) :: mo_vgl(mo_num,elec_num,5) + real (c_double ) , intent(out) :: det_vgl_beta(beta_num,beta_num,5,walk_num,det_num_beta) + + integer(c_int32_t), external :: qmckl_compute_det_vgl_beta_f + info = qmckl_compute_det_vgl_beta_f & + (context, & + det_num_beta, & + walk_num, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_beta, & + mo_num, & + mo_vgl, & + det_vgl_beta) + + end function qmckl_compute_det_vgl_beta + #+end_src + +*** 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 * 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); + +/* Set up MO data */ +const int64_t mo_num = chbrclf_mo_num; +rc = qmckl_set_mo_basis_mo_num(context, mo_num); +assert (rc == QMCKL_SUCCESS); + +const double * mo_coefficient = &(chbrclf_mo_coef[0]); + +rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient); +assert (rc == QMCKL_SUCCESS); + +assert(qmckl_mo_basis_provided(context)); + +double mo_vgl[5][elec_num][chbrclf_mo_num]; +rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0])); +assert (rc == QMCKL_SUCCESS); + +/* Set up determinant data */ + +const int64_t det_num_alpha = 1; +const int64_t det_num_beta = 1; +int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num]; +int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num]; + +int i, j, k; +for(k = 0; k < det_num_alpha; ++k) + for(i = 0; i < walk_num; ++i) + for(j = 0; j < elec_up_num; ++j) + mo_index_alpha[k][i][j] = j + 1; +for(k = 0; k < det_num_beta; ++k) + for(i = 0; i < walk_num; ++i) + for(j = 0; j < elec_up_num; ++j) + mo_index_beta[k][i][j] = j + 1; + +rc = qmckl_set_determinant_type (context, typ); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_walk_num (context, walk_num); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_det_num_beta (context, det_num_beta); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0])); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0])); +assert (rc == QMCKL_SUCCESS); + +// Get slater-determinant + +double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num]; +double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num]; + +rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + + #+end_src + +** Inverse of Determinant matrix + :PROPERTIES: + :Name: qmckl_compute_det_inv_matrix + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double* const det_inv_matrix_alpha); +qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double* const det_inv_matrix_beta); +qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double* const det_adj_matrix_alpha); +qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double* const det_adj_matrix_beta); +qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double* const det_adj_matrix_alpha); +qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double* const det_adj_matrix_beta); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_det_inv_matrix_alpha(qmckl_context context, double * const det_inv_matrix_alpha) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_alpha(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_inv_matrix_alpha(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num; + memcpy(det_inv_matrix_alpha, ctx->det.det_inv_matrix_alpha, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_det_inv_matrix_beta(qmckl_context context, double * const det_inv_matrix_beta) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_beta(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_inv_matrix_beta(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num; + memcpy(det_inv_matrix_beta, ctx->det.det_inv_matrix_beta, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_det_adj_matrix_alpha(qmckl_context context, double * const det_adj_matrix_alpha) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_alpha(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_inv_matrix_alpha(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.up_num * ctx->electron.up_num; + memcpy(det_adj_matrix_alpha, ctx->det.det_adj_matrix_alpha, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_det_adj_matrix_beta(qmckl_context context, double * const det_adj_matrix_beta) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_beta(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_inv_matrix_beta(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num * ctx->electron.down_num * ctx->electron.down_num; + memcpy(det_adj_matrix_beta, ctx->det.det_adj_matrix_beta, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_det_alpha(qmckl_context context, double * const det_value_alpha) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_alpha(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_inv_matrix_alpha(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num; + memcpy(det_value_alpha, ctx->det.det_value_alpha, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_get_det_beta(qmckl_context context, double * const det_value_beta) { + + 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; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_vgl_beta(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_det_inv_matrix_beta(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->det.det_num_alpha * ctx->det.walk_num; + memcpy(det_value_beta, ctx->det.det_value_beta, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(qmckl_context context); +qmckl_exit_code qmckl_provide_det_inv_matrix_beta(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_det_inv_matrix_alpha(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->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->det.det_inv_matrix_alpha_date) { + + /* Allocate array */ + if (ctx->det.det_inv_matrix_alpha == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha * + ctx->electron.up_num * ctx->electron.up_num * sizeof(double); + double* det_inv_matrix_alpha = (double*) qmckl_malloc(context, mem_info); + + if (det_inv_matrix_alpha == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_inv_matrix_alpha", + NULL); + } + ctx->det.det_inv_matrix_alpha = det_inv_matrix_alpha; + } + + if (ctx->det.det_adj_matrix_alpha == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha * + ctx->electron.up_num * ctx->electron.up_num * sizeof(double); + double* det_adj_matrix_alpha = (double*) qmckl_malloc(context, mem_info); + + if (det_adj_matrix_alpha == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_adj_matrix_alpha", + NULL); + } + ctx->det.det_adj_matrix_alpha = det_adj_matrix_alpha; + } + + if (ctx->det.det_value_alpha == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_alpha * sizeof(double); + double* det_value_alpha = (double*) qmckl_malloc(context, mem_info); + + if (det_value_alpha == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_value_alpha", + NULL); + } + ctx->det.det_value_alpha = det_value_alpha; + } + + qmckl_exit_code rc; + if (ctx->det.type == 'G') { + rc = qmckl_compute_det_inv_matrix_alpha(context, + ctx->det.det_num_alpha, + ctx->det.walk_num, + ctx->electron.up_num, + ctx->det.det_vgl_alpha, + ctx->det.det_value_alpha, + ctx->det.det_adj_matrix_alpha, + ctx->det.det_inv_matrix_alpha); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_det_inv_matrix_alpha", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->det.det_value_alpha_date = ctx->date; + ctx->det.det_adj_matrix_alpha_date = ctx->date; + ctx->det.det_inv_matrix_alpha_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_provide_det_inv_matrix_beta(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->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->det.det_inv_matrix_beta_date) { + + /* Allocate array */ + if (ctx->det.det_inv_matrix_beta == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta * + ctx->electron.down_num * ctx->electron.down_num * sizeof(double); + double* det_inv_matrix_beta = (double*) qmckl_malloc(context, mem_info); + + if (det_inv_matrix_beta == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_inv_matrix_beta", + NULL); + } + ctx->det.det_inv_matrix_beta = det_inv_matrix_beta; + } + + if (ctx->det.det_adj_matrix_beta == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta * + ctx->electron.down_num * ctx->electron.down_num * sizeof(double); + double* det_adj_matrix_beta = (double*) qmckl_malloc(context, mem_info); + + if (det_adj_matrix_beta == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_adj_matrix_beta", + NULL); + } + ctx->det.det_adj_matrix_beta = det_adj_matrix_beta; + } + + if (ctx->det.det_value_beta == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->det.walk_num * ctx->det.det_num_beta * sizeof(double); + double* det_value_beta = (double*) qmckl_malloc(context, mem_info); + + if (det_value_beta == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_det_value_beta", + NULL); + } + ctx->det.det_value_beta = det_value_beta; + } + + qmckl_exit_code rc; + if (ctx->det.type == 'G') { + rc = qmckl_compute_det_inv_matrix_beta(context, + ctx->det.det_num_beta, + ctx->det.walk_num, + ctx->electron.down_num, + ctx->det.det_vgl_beta, + ctx->det.det_value_beta, + ctx->det.det_adj_matrix_beta, + ctx->det.det_inv_matrix_beta); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_det_inv_matrix_beta", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->det.det_value_beta_date = ctx->date; + ctx->det.det_adj_matrix_beta_date = ctx->date; + ctx->det.det_inv_matrix_beta_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute alpha + :PROPERTIES: + :Name: qmckl_compute_det_inv_matrix_alpha + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_det_inv_matrix_alpha_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~det_num_alpha~ | in | Number of determinants | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~int64_t~ | ~alpha_num~ | in | Number of electrons | + | ~double~ | ~det_vgl_alpha[det_num_alpha][walk_num][5][alpha_num][alpha_num]~ | in | determinant matrix Value, gradients and Laplacian of the MOs | + | ~double~ | ~det_value_alpha[det_num_alpha][walk_num]~ | out | value of determinant matrix | + | ~double~ | ~det_adj_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | out | adjoint of determinant matrix | + | ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | out | inverse of determinant matrix | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_det_inv_matrix_alpha_f(context, & + det_num_alpha, walk_num, alpha_num, det_vgl_alpha, det_value_alpha, det_adj_matrix_alpha, det_inv_matrix_alpha) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: det_num_alpha + integer*8, intent(in) :: walk_num + integer*8, intent(in) :: alpha_num + double precision, intent(in) :: det_vgl_alpha(alpha_num, alpha_num, 5, walk_num, det_num_alpha) + double precision, intent(inout) :: det_value_alpha(walk_num, det_num_alpha) + double precision, intent(inout) :: det_adj_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha) + double precision, intent(inout) :: det_inv_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha) + double precision,dimension(:,:),allocatable :: matA + double precision :: det_l + integer*8 :: idet, iwalk, ielec, mo_id, imo, LDA, res, i, j + + allocate(matA(alpha_num, alpha_num)) + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (det_num_alpha <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (alpha_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + LDA = alpha_num + do idet = 1, det_num_alpha + do iwalk = 1, walk_num + ! Value + matA(1:alpha_num,1:alpha_num) = det_vgl_alpha(1:alpha_num, 1:alpha_num, 1, iwalk, idet) + res = qmckl_adjoint(context, alpha_num, alpha_num, LDA, matA, det_l) + det_adj_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA + det_inv_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA/det_l + det_value_alpha(iwalk, idet) = det_l + end do + end do + + deallocate(matA) +end function qmckl_compute_det_inv_matrix_alpha_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_det_inv_matrix_alpha_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_inv_matrix_alpha")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_det_inv_matrix_alpha ( + const qmckl_context context, + const int64_t det_num_alpha, + const int64_t walk_num, + const int64_t alpha_num, + const double* det_vgl_alpha, + double* const det_value_alpha, + double* const det_adj_matrix_alpha, + double* const det_inv_matrix_alpha ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_det_inv_matrix_alpha_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_inv_matrix_alpha")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_det_inv_matrix_alpha & + (context, & + det_num_alpha, & + walk_num, & + alpha_num, & + det_vgl_alpha, & + det_value_alpha, & + det_adj_matrix_alpha, & + det_inv_matrix_alpha) & + 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 :: det_num_alpha + integer (c_int64_t) , intent(in) , value :: walk_num + integer (c_int64_t) , intent(in) , value :: alpha_num + real (c_double ) , intent(in) :: det_vgl_alpha(alpha_num,alpha_num,5,walk_num,det_num_alpha) + real (c_double ) , intent(out) :: det_value_alpha(walk_num,det_num_alpha) + real (c_double ) , intent(out) :: det_adj_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha) + real (c_double ) , intent(out) :: det_inv_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha) + + integer(c_int32_t), external :: qmckl_compute_det_inv_matrix_alpha_f + info = qmckl_compute_det_inv_matrix_alpha_f & + (context, & + det_num_alpha, & + walk_num, & + alpha_num, & + det_vgl_alpha, & + det_value_alpha, & + det_adj_matrix_alpha, & + det_inv_matrix_alpha) + + end function qmckl_compute_det_inv_matrix_alpha + #+end_src + +*** Compute beta + :PROPERTIES: + :Name: qmckl_compute_det_inv_matrix_beta + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_det_inv_matrix_beta_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~det_num_beta~ | in | Number of determinants | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~int64_t~ | ~beta_num~ | in | Number of electrons | + | ~double~ | ~det_vgl_beta[det_num_beta][walk_num][5][beta_num][beta_num]~ | in | determinant matrix Value, gradients and Laplacian of the MOs | + | ~double~ | ~det_value_beta[det_num_beta][walk_num]~ | out | value of determinant matrix | + | ~double~ | ~det_adj_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | out | adjoint of determinant matrix | + | ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | out | inverse of determinant matrix | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_det_inv_matrix_beta_f(context, & + det_num_beta, walk_num, beta_num, det_vgl_beta, det_value_beta, det_adj_matrix_beta, det_inv_matrix_beta) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: det_num_beta + integer*8, intent(in) :: walk_num + integer*8, intent(in) :: beta_num + double precision, intent(in) :: det_vgl_beta(beta_num, beta_num, 5, walk_num, det_num_beta) + double precision, intent(inout) :: det_value_beta(walk_num, det_num_beta) + double precision, intent(inout) :: det_adj_matrix_beta(beta_num, beta_num, walk_num, det_num_beta) + double precision, intent(inout) :: det_inv_matrix_beta(beta_num, beta_num, walk_num, det_num_beta) + double precision,dimension(:,:),allocatable :: matA + double precision :: det_l + integer*8 :: idet, iwalk, ielec, mo_id, imo, LDA, res + + allocate(matA(beta_num, beta_num)) + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (det_num_beta <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (beta_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + LDA = beta_num + do idet = 1, det_num_beta + do iwalk = 1, walk_num + ! Value + matA(1:beta_num,1:beta_num) = det_vgl_beta(1:beta_num, 1:beta_num, 1, iwalk, idet) + res = qmckl_adjoint(context, beta_num, beta_num, LDA, matA, det_l) + det_adj_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA + det_inv_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA/det_l + det_value_beta(iwalk, idet) = det_l + end do + end do + + deallocate(matA) +end function qmckl_compute_det_inv_matrix_beta_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_det_inv_matrix_beta_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_inv_matrix_beta")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_det_inv_matrix_beta ( + const qmckl_context context, + const int64_t det_num_beta, + const int64_t walk_num, + const int64_t beta_num, + const double* det_vgl_beta, + double* const det_value_beta, + double* const det_adj_matrix_beta, + double* const det_inv_matrix_beta ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_det_inv_matrix_beta_args,rettyp=get_value("CRetType"),fname="qmckl_compute_det_inv_matrix_beta")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_det_inv_matrix_beta & + (context, & + det_num_beta, & + walk_num, & + beta_num, & + det_vgl_beta, & + det_value_beta, & + det_adj_matrix_beta, & + det_inv_matrix_beta) & + 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 :: det_num_beta + integer (c_int64_t) , intent(in) , value :: walk_num + integer (c_int64_t) , intent(in) , value :: beta_num + real (c_double ) , intent(in) :: det_vgl_beta(beta_num,beta_num,5,walk_num,det_num_beta) + real (c_double ) , intent(out) :: det_value_beta(walk_num,det_num_beta) + real (c_double ) , intent(out) :: det_adj_matrix_beta(beta_num,beta_num,walk_num,det_num_beta) + real (c_double ) , intent(out) :: det_inv_matrix_beta(beta_num,beta_num,walk_num,det_num_beta) + + integer(c_int32_t), external :: qmckl_compute_det_inv_matrix_beta_f + info = qmckl_compute_det_inv_matrix_beta_f & + (context, & + det_num_beta, & + walk_num, & + beta_num, & + det_vgl_beta, & + det_value_beta, & + det_adj_matrix_beta, & + det_inv_matrix_beta) + + end function qmckl_compute_det_inv_matrix_beta + #+end_src + + +*** Test + #+begin_src c :tangle (eval c_test) :exports none +// Get adjoint of the slater-determinant + +double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num]; +double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num]; + +rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + + #+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/qmckl_electron.org b/org/qmckl_electron.org index 3104331..aff2002 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -68,29 +68,38 @@ int main() { The following data stored in the context: - | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | - | ~num~ | ~int64_t~ | Total number of electrons | - | ~up_num~ | ~int64_t~ | Number of up-spin electrons | - | ~down_num~ | ~int64_t~ | Number of down-spin electrons | - | ~walk_num~ | ~int64_t~ | Number of walkers | - | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | - | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | - | ~provided~ | ~bool~ | If true, ~electron~ is valid | - | ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates | - | ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates | - | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | - | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | - | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances | - | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | - | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | - | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | - | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + |---------------------------+----------------------------+-------------------------------------------| + | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | + | ~num~ | ~int64_t~ | Total number of electrons | + | ~up_num~ | ~int64_t~ | Number of up-spin electrons | + | ~down_num~ | ~int64_t~ | Number of down-spin electrons | + | ~walk_num~ | ~int64_t~ | Number of walkers | + | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | + | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | + | ~provided~ | ~bool~ | If true, ~electron~ is valid | + | ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates | + | ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates | + | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | + + Computed data: + + |-------------------------------------+--------------------------------------+----------------------------------------------------------------------| + | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | + | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances | + | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | + | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~ee_pot~ | ~double[walk_num]~ | Electron-electron rescaled distances derivatives | + | ~ee_pot_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~en_pot~ | double[walk_num] | Electron-nucleus potential energy | + | ~en_pot_date~ | int64_t | Date when the electron-nucleus potential energy was computed | + | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | + | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | ** Data structure @@ -105,6 +114,8 @@ typedef struct qmckl_electron_struct { int64_t coord_new_date; int64_t ee_distance_date; int64_t en_distance_date; + int64_t ee_pot_date; + int64_t en_pot_date; int64_t ee_distance_rescaled_date; int64_t ee_distance_rescaled_deriv_e_date; int64_t en_distance_rescaled_date; @@ -113,6 +124,8 @@ typedef struct qmckl_electron_struct { double* coord_old; double* ee_distance; double* en_distance; + double* ee_pot; + double* en_pot; double* ee_distance_rescaled; double* ee_distance_rescaled_deriv_e; double* en_distance_rescaled; @@ -188,7 +201,7 @@ if ( (ctx->electron.uninitialized & mask) != 0) { return NULL; } #+end_src - + *** Number of electrons #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -1383,7 +1396,7 @@ qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context co #+end_src *** Provide :noexport: - + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context); #+end_src @@ -1570,6 +1583,205 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal #+end_src +** Electron-electron potential + + ~ee_pot~ calculates the ~ee~ potential energy. + + \[ + \mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j>i}^{N_e}\frac{1}{r_{ij}} + \] + + where \(\mathcal{V}_{ee}\) is the ~ee~ potential and \[r_{ij}\] the ~ee~ + distance. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_ee_potential(qmckl_context context, double* const ee_pot) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ee_potential(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.walk_num * sizeof(double); + memcpy(ee_pot, ctx->electron.ee_pot, sze); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide :noexport: + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_potential(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->electron.provided) return QMCKL_NOT_PROVIDED; + + qmckl_exit_code rc = qmckl_provide_ee_distance(context); + if (rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->electron.ee_pot_date) { + + /* Allocate array */ + if (ctx->electron.ee_pot == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* ee_pot = (double*) qmckl_malloc(context, mem_info); + + if (ee_pot == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_ee_potential", + NULL); + } + ctx->electron.ee_pot = ee_pot; + } + + qmckl_exit_code rc = + qmckl_compute_ee_potential(context, + ctx->electron.num, + ctx->electron.walk_num, + ctx->electron.ee_distance, + ctx->electron.ee_pot); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->electron.ee_pot_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_ee_potential + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_ee_potential_args + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | int64_t | walk_num | in | Number of walkers | + | double | ee_distance[walk_num][elec_num][elec_num] | in | Electron-electron rescaled distances | + | double | ee_pot[walk_num] | out | Electron-electron potential | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, & + ee_distance, ee_pot) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + double precision , intent(out) :: ee_pot(walk_num) + + integer*8 :: nw, i, j + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + ee_pot = 0.0d0 + do nw=1,walk_num + do j=2,elec_num + do i=1,j-1 + if (dabs(ee_distance(i,j,nw)) > 1e-5) then + ee_pot(nw) = ee_pot(nw) + 1.0d0/(ee_distance(i,j,nw)) + endif + end do + end do + end do + +end function qmckl_compute_ee_potential_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_ee_potential ( + const qmckl_context context, + const int64_t elec_num, + const int64_t walk_num, + const double* ee_distance, + double* const ee_pot ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ee_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_ee_potential & + (context, elec_num, walk_num, ee_distance, ee_pot) & + 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 :: elec_num + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) + real (c_double ) , intent(out) :: ee_pot(walk_num) + + integer(c_int32_t), external :: qmckl_compute_ee_potential_f + info = qmckl_compute_ee_potential_f & + (context, elec_num, walk_num, ee_distance, ee_pot) + + end function qmckl_compute_ee_potential + #+end_src + +*** Test + #+begin_src c :tangle (eval c_test) +double ee_pot[walk_num]; + +rc = qmckl_get_electron_ee_potential(context, &(ee_pot[0])); +assert (rc == QMCKL_SUCCESS); + #+end_src ** Electron-nucleus distances *** Get @@ -2407,6 +2619,216 @@ assert (rc == QMCKL_SUCCESS); #+end_src +** Electron-nucleus potential + ~en_potential~ stores the ~en~ potential energy + + \[ + \mathcal{V}_{en} = -\sum_{i=1}^{N_e}\sum_{A=1}^{N_n}\frac{Z_A}{r_{iA}} + \] + + where \(\mathcal{V}_{en}\) is the ~en~ potential, \[r_{iA}\] the ~en~ + distance and \[Z_A\] is the nuclear charge. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_en_potential(qmckl_context context, double* const en_pot) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_en_potential(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.walk_num * sizeof(double); + memcpy(en_pot, ctx->electron.en_pot, sze); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide :noexport: + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_potential(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_potential(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->electron.provided) return QMCKL_NOT_PROVIDED; + if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; + + qmckl_exit_code rc = qmckl_provide_en_distance(context); + if (rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->electron.en_pot_date) { + + /* Allocate array */ + if (ctx->electron.en_pot == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* en_pot = (double*) qmckl_malloc(context, mem_info); + + if (en_pot == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_en_potential", + NULL); + } + ctx->electron.en_pot = en_pot; + } + + qmckl_exit_code rc = + qmckl_compute_en_potential(context, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.walk_num, + ctx->nucleus.charge, + ctx->electron.en_distance, + ctx->electron.en_pot); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->electron.en_pot_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_en_potential + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_en_potential_args + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | int64_t | nucl_num | in | Number of nucleii | + | int64_t | walk_num | in | Number of walkers | + | double | charge[nucl_num] | in | charge of nucleus | + | double | en_distance[walk_num][nucl_num][elec_num] | in | Electron-electron rescaled distances | + | double | en_pot[walk_num] | out | Electron-electron potential | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_en_potential_f(context, elec_num, nucl_num, walk_num, & + charge, en_distance, en_pot) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: nucl_num + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: charge(nucl_num) + double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + double precision , intent(out) :: en_pot(walk_num) + + integer*8 :: nw, i, j + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + en_pot = 0.0d0 + do nw=1,walk_num + do j=1,nucl_num + do i=1,elec_num + if (dabs(en_distance(i,j,nw)) > 1e-5) then + en_pot(nw) = en_pot(nw) - charge(j)/(en_distance(i,j,nw)) + endif + end do + end do + end do + +end function qmckl_compute_en_potential_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_en_potential ( + const qmckl_context context, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t walk_num, + const double* charge, + const double* en_distance, + double* const en_pot ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_en_potential_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_en_potential & + (context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot) & + 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 :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: charge(nucl_num) + real (c_double ) , intent(in) :: en_distance(elec_num,nucl_num,walk_num) + real (c_double ) , intent(out) :: en_pot(walk_num) + + integer(c_int32_t), external :: qmckl_compute_en_potential_f + info = qmckl_compute_en_potential_f & + (context, elec_num, nucl_num, walk_num, charge, en_distance, en_pot) + + end function qmckl_compute_en_potential + #+end_src + +*** Test + #+begin_src c :tangle (eval c_test) +double en_pot[walk_num]; + +rc = qmckl_get_electron_en_potential(context, &(en_pot[0])); +assert (rc == QMCKL_SUCCESS); + #+end_src + * End of files :noexport: #+begin_src c :tangle (eval h_private_type) diff --git a/org/qmckl_local_energy.org b/org/qmckl_local_energy.org new file mode 100644 index 0000000..a386d01 --- /dev/null +++ b/org/qmckl_local_energy.org @@ -0,0 +1,1705 @@ +#+TITLE: Local Energy +#+SETUPFILE: ../tools/theme.setup +#+INCLUDE: ../tools/lib.org + + +Here we calculate the final expectation value of the +local energy \[E_L\] as the sum of the kinetic energy and potential energy. + +\[ +E_L = KE + PE +\] + +Where the kinetic energy is given as: + +\[ +KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi} +\] + +The laplacian of the wavefunction in the single-determinant +case is given as follows: + +\[ +\frac{\bigtriangleup \Psi(r)}{\Psi(r)} = \sum_{j=1}^{N_e} \bigtriangleup \Phi_j(r_i) D_{ji}^{-1}(r) +\] + +The potential energy is the sum of all the following terms + +\[ +PE = \mathcal{V}_{ee} + \mathcal{V}_{en} + \mathcal{V}_{nn} +\] + +The potential for is calculated as the sum of single electron +contributions. + +\[ +\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{j + #+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" +#include "qmckl_determinant_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" +#include "qmckl_determinant_private_type.h" +#include "qmckl_determinant_private_func.h" + #+end_src + +* Context + + The following arrays are stored in the context: + + | | + + Computed data: + + |--------------+---------------------------+----------------------------| + | ~e_kin~ | ~[walk_num]~ | total kinetic energy | + | ~e_pot~ | ~[walk_num]~ | total potential energy | + | ~e_local~ | ~[walk_num]~ | local energy | + | ~r_drift~ | ~[3][walk_num][elec_num]~ | The drift vector | + | ~y_move~ | ~[3][walk_num]~ | The diffusion move | + | ~accep_prob~ | ~[walk_num]~ | The acceptance probability | + |--------------+---------------------------+----------------------------| + +** Data structure + + #+begin_src c :comments org :tangle (eval h_private_type) +typedef struct qmckl_local_energy_struct { + double * e_kin; + double * e_pot; + double * e_local; + double * accep_prob; + double * r_drift; + double * y_move; + int64_t e_kin_date; + int64_t e_pot_date; + int64_t e_local_date; + int64_t accep_prob_date; + int64_t r_drift_date; + int64_t y_move_date; + + int32_t uninitialized; + bool provided; +} qmckl_local_energy_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. + +* Computation +** Kinetic energy + :PROPERTIES: + :Name: qmckl_compute_kinetic_energy + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + +Where the kinetic energy is given as: + +\[ +KE = -\frac{1}{2}\frac{\bigtriangleup \Psi}{\Psi} +\] + +The laplacian of the wavefunction in the single-determinant +case is given as follows: + +\[ +\frac{\bigtriangleup \Psi(r)}{\Psi(r)} = \sum_{j=1}^{N_e} \bigtriangleup \Phi_j(r_i) D_{ji}^{-1}(r) +\] + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double* const kinetic_energy); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_kinetic_energy(qmckl_context context, double * const kinetic_energy) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED; + + if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED; + + rc = qmckl_provide_ao_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_kinetic_energy(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.walk_num * sizeof(double); + memcpy(kinetic_energy, ctx->local_energy.e_kin, sze); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) { + + qmckl_exit_code rc; + 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->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + rc = qmckl_provide_det_inv_matrix_alpha(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_det_inv_matrix_alpha", + NULL); + } + + rc = qmckl_provide_det_inv_matrix_beta(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_det_inv_matrix_beta", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->local_energy.e_kin_date) { + + /* Allocate array */ + if (ctx->local_energy.e_kin == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* e_kin = (double*) qmckl_malloc(context, mem_info); + + if (e_kin == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_e_kin", + NULL); + } + ctx->local_energy.e_kin = e_kin; + } + + if (ctx->det.type == 'G') { + rc = qmckl_compute_kinetic_energy(context, + ctx->det.walk_num, + ctx->det.det_num_alpha, + ctx->det.det_num_beta, + ctx->electron.up_num, + ctx->electron.down_num, + ctx->electron.num, + ctx->det.mo_index_alpha, + ctx->det.mo_index_beta, + ctx->mo_basis.mo_num, + ctx->mo_basis.mo_vgl, + ctx->det.det_value_alpha, + ctx->det.det_value_beta, + ctx->det.det_inv_matrix_alpha, + ctx->det.det_inv_matrix_beta, + ctx->local_energy.e_kin); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_kinetic_energy", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->local_energy.e_kin_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute kinetic enregy + :PROPERTIES: + :Name: qmckl_compute_kinetic_energy + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_compute_kinetic_energy_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~int64_t~ | ~det_num_alpha~ | in | Number of determinants | + | ~int64_t~ | ~det_num_beta~ | in | Number of determinants | + | ~int64_t~ | ~alpha_num~ | in | Number of electrons | + | ~int64_t~ | ~beta_num~ | in | Number of electrons | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~mo_index_alpha[det_num_alpha][walk_num][alpha_num]~ | in | MO indices for electrons | + | ~int64_t~ | ~mo_index_beta[det_num_beta][walk_num][beta_num]~ | in | MO indices for electrons | + | ~int64_t~ | ~mo_num~ | in | Number of MOs | + | ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs | + | ~double~ | ~det_value_alpha[det_num_alpha][walk_num]~ | in | Det of wavefunction | + | ~double~ | ~det_value_beta[det_num_beta][walk_num]~ | in | Det of wavefunction | + | ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in | Value, gradients and Laplacian of the Det | + | ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det | + | ~double~ | ~e_kin[walk_num]~ | out | Kinetic energy | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_kinetic_energy_f(context, walk_num, & + det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, & + mo_num, mo_vgl, det_value_alpha, det_value_beta, det_inv_matrix_alpha, det_inv_matrix_beta, e_kin) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: walk_num + integer*8, intent(in) :: det_num_alpha + integer*8, intent(in) :: det_num_beta + integer*8, intent(in) :: alpha_num + integer*8, intent(in) :: beta_num + integer*8, intent(in) :: elec_num + integer*8, intent(in) :: mo_num + integer*8, intent(in) :: mo_index_alpha(alpha_num, walk_num, det_num_alpha) + integer*8, intent(in) :: mo_index_beta(beta_num, walk_num, det_num_beta) + double precision, intent(in) :: mo_vgl(mo_num, elec_num, 5) + double precision, intent(in) :: det_value_alpha(walk_num, det_num_alpha) + double precision, intent(in) :: det_value_beta(walk_num, det_num_beta) + double precision, intent(in) :: det_inv_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha) + double precision, intent(in) :: det_inv_matrix_beta(beta_num, beta_num, walk_num, det_num_beta) + double precision, intent(inout) :: e_kin(walk_num) + double precision :: tmp_e + integer*8 :: idet, iwalk, ielec, mo_id, imo + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (alpha_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (beta_num < 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + e_kin = 0.0d0 + do idet = 1, det_num_alpha + do iwalk = 1, walk_num + ! Alpha part + tmp_e = 0.0d0 + do imo = 1, alpha_num + do ielec = 1, alpha_num + mo_id = mo_index_alpha(imo, iwalk, idet) + e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, ielec, 5) + !print *,"det alpha = ",det_inv_matrix_alpha(imo,ielec,iwalk,idet) + !print *,mo_vgl(mo_id,ielec,5) + !!print *," det val = ",det_value_alpha(iwalk,idet) + !tmp_e = tmp_e - 0.5d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * & + ! mo_vgl(mo_id, ielec, 5) + end do + !print *,"e_kin = ",tmp_e + end do + ! Beta part + tmp_e = 0.0d0 + do imo = 1, beta_num + do ielec = 1, beta_num + mo_id = mo_index_beta(imo, iwalk, idet) + e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, alpha_num + ielec, 5) + !print *,"det beta = ",det_inv_matrix_beta(imo,ielec,iwalk,idet) + !print *,mo_vgl(mo_id,alpha_num+ielec,5) + !!print *," det val = ",det_value_alpha(iwalk,idet) + !tmp_e = tmp_e - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * & + ! mo_vgl(mo_id, alpha_num + ielec, 5) + end do + !print *,"e_kin = ",tmp_e + end do + end do + end do + +end function qmckl_compute_kinetic_energy_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_compute_kinetic_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_kinetic_energy")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_kinetic_energy ( + const qmckl_context context, + const int64_t walk_num, + const int64_t det_num_alpha, + const int64_t det_num_beta, + const int64_t alpha_num, + const int64_t beta_num, + const int64_t elec_num, + const int64_t* mo_index_alpha, + const int64_t* mo_index_beta, + const int64_t mo_num, + const double* mo_vgl, + const double* det_value_alpha, + const double* det_value_beta, + const double* det_inv_matrix_alpha, + const double* det_inv_matrix_beta, + double* const e_kin ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_compute_kinetic_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_kinetic_energy")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_kinetic_energy & + (context, & + walk_num, & + det_num_alpha, & + det_num_beta, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_alpha, & + mo_index_beta, & + mo_num, & + mo_vgl, & + det_value_alpha, & + det_value_beta, & + det_inv_matrix_alpha, & + det_inv_matrix_beta, & + e_kin) & + 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 :: walk_num + integer (c_int64_t) , intent(in) , value :: det_num_alpha + integer (c_int64_t) , intent(in) , value :: det_num_beta + integer (c_int64_t) , intent(in) , value :: alpha_num + integer (c_int64_t) , intent(in) , value :: beta_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) :: mo_index_alpha(alpha_num,walk_num,det_num_alpha) + integer (c_int64_t) , intent(in) :: mo_index_beta(beta_num,walk_num,det_num_beta) + integer (c_int64_t) , intent(in) , value :: mo_num + real (c_double ) , intent(in) :: mo_vgl(mo_num,elec_num,5) + real (c_double ) , intent(in) :: det_value_alpha(walk_num,det_num_alpha) + real (c_double ) , intent(in) :: det_value_beta(walk_num,det_num_beta) + real (c_double ) , intent(in) :: det_inv_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha) + real (c_double ) , intent(in) :: det_inv_matrix_beta(beta_num,beta_num,walk_num,det_num_beta) + real (c_double ) , intent(out) :: e_kin(walk_num) + + integer(c_int32_t), external :: qmckl_compute_kinetic_energy_f + info = qmckl_compute_kinetic_energy_f & + (context, & + walk_num, & + det_num_alpha, & + det_num_beta, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_alpha, & + mo_index_beta, & + mo_num, & + mo_vgl, & + det_value_alpha, & + det_value_beta, & + det_inv_matrix_alpha, & + det_inv_matrix_beta, & + e_kin) + + end function qmckl_compute_kinetic_energy + #+end_src + +*** 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 * 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); + +/* Set up MO data */ +const int64_t mo_num = chbrclf_mo_num; +rc = qmckl_set_mo_basis_mo_num(context, mo_num); +assert (rc == QMCKL_SUCCESS); + +const double * mo_coefficient = &(chbrclf_mo_coef[0]); + +rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient); +assert (rc == QMCKL_SUCCESS); + +assert(qmckl_mo_basis_provided(context)); + +double mo_vgl[5][elec_num][chbrclf_mo_num]; +rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0])); +assert (rc == QMCKL_SUCCESS); + +/* Set up determinant data */ + +const int64_t det_num_alpha = 1; +const int64_t det_num_beta = 1; +int64_t mo_index_alpha[det_num_alpha][walk_num][elec_up_num]; +int64_t mo_index_beta[det_num_alpha][walk_num][elec_dn_num]; + +int i, j, k; +for(k = 0; k < det_num_alpha; ++k) + for(i = 0; i < walk_num; ++i) + for(j = 0; j < elec_up_num; ++j) + mo_index_alpha[k][i][j] = j + 1; +for(k = 0; k < det_num_beta; ++k) + for(i = 0; i < walk_num; ++i) + for(j = 0; j < elec_up_num; ++j) + mo_index_beta[k][i][j] = j + 1; + +rc = qmckl_set_determinant_type (context, typ); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_walk_num (context, walk_num); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_det_num_beta (context, det_num_beta); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0])); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0])); +assert (rc == QMCKL_SUCCESS); + +// Get alpha determinant + +double det_vgl_alpha[det_num_alpha][walk_num][5][elec_up_num][elec_up_num]; +double det_vgl_beta[det_num_beta][walk_num][5][elec_dn_num][elec_dn_num]; + +rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +// Get adjoint of the slater-determinant + +double det_inv_matrix_alpha[det_num_alpha][walk_num][elec_up_num][elec_up_num]; +double det_inv_matrix_beta[det_num_beta][walk_num][elec_dn_num][elec_dn_num]; + +rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0])); +assert (rc == QMCKL_SUCCESS); + +// Calculate the Kinetic energy + +double kinetic_energy[walk_num]; + +rc = qmckl_get_kinetic_energy(context, &(kinetic_energy[0])); +assert (rc == QMCKL_SUCCESS); + + #+end_src + +** Potential energy + :PROPERTIES: + :Name: qmckl_compute_potential_energy + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + +The potential energy is the sum of all the following terms + +\[ +PE = \mathcal{V}_{ee} + \mathcal{V}_{en} + \mathcal{V}_{nn} +\] + +The potential for is calculated as the sum of single electron +contributions. + +\[ +\mathcal{V}_{ee} = \sum_{i=1}^{N_e}\sum_{jelectron.walk_num * sizeof(double); + memcpy(potential_energy, ctx->local_energy.e_pot, sze); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_potential_energy(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); + + qmckl_exit_code rc; + if(!(ctx->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + rc = qmckl_provide_nucleus_repulsion(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_nucleus_repulsion", + NULL); + } + + rc = qmckl_provide_ee_potential(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ee_potential", + NULL); + } + + rc = qmckl_provide_en_potential(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_en_potential", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->local_energy.e_pot_date) { + + /* Allocate array */ + if (ctx->local_energy.e_pot == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* e_pot = (double*) qmckl_malloc(context, mem_info); + + if (e_pot == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_e_pot", + NULL); + } + ctx->local_energy.e_pot = e_pot; + } + + if (ctx->det.type == 'G') { + rc = qmckl_compute_potential_energy(context, + ctx->det.walk_num, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.ee_pot, + ctx->electron.en_pot, + ctx->nucleus.repulsion, + ctx->local_energy.e_pot); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_potential_energy", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->local_energy.e_pot_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute potential enregy + :PROPERTIES: + :Name: qmckl_compute_potential_energy + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_compute_potential_energy_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~nucl_num~ | in | Number of MOs | + | ~double~ | ~ee_pot[walk_num]~ | in | ee potential | + | ~double~ | ~en_pot[walk_num]~ | in | en potential | + | ~double~ | ~repulsion~ | in | en potential | + | ~double~ | ~e_pot[walk_num]~ | out | Potential energy | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_potential_energy_f(context, walk_num, & + elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: walk_num + integer*8, intent(in) :: elec_num + integer*8, intent(in) :: nucl_num + double precision, intent(in) :: ee_pot(walk_num) + double precision, intent(in) :: en_pot(walk_num) + double precision, intent(in) :: repulsion + double precision, intent(inout) :: e_pot(walk_num) + integer*8 :: idet, iwalk, ielec, mo_id, imo + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + e_pot = 0.0d0 + repulsion + do iwalk = 1, walk_num + e_pot(iwalk) = e_pot(iwalk) + ee_pot(iwalk) + en_pot(iwalk) + end do + +end function qmckl_compute_potential_energy_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_compute_potential_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_potential_energy")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_potential_energy ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const double* ee_pot, + const double* en_pot, + const double repulsion, + double* const e_pot ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_compute_potential_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_potential_energy")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_potential_energy & + (context, walk_num, elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) & + 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 :: walk_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + real (c_double ) , intent(in) :: ee_pot(walk_num) + real (c_double ) , intent(in) :: en_pot(walk_num) + real (c_double ) , intent(in) , value :: repulsion + real (c_double ) , intent(out) :: e_pot(walk_num) + + integer(c_int32_t), external :: qmckl_compute_potential_energy_f + info = qmckl_compute_potential_energy_f & + (context, walk_num, elec_num, nucl_num, ee_pot, en_pot, repulsion, e_pot) + + end function qmckl_compute_potential_energy + #+end_src + +*** Test + #+begin_src c :tangle (eval c_test) :exports none +// Calculate the Potential energy + +double potential_energy[walk_num]; + +rc = qmckl_get_potential_energy(context, &(potential_energy[0])); +assert (rc == QMCKL_SUCCESS); + + #+end_src + +** Local energy + :PROPERTIES: + :Name: qmckl_compute_local_energy + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + +The local energy is the sum of kinetic and potential energies. + +\[ +E_L = KE + PE +\] + + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double* const local_energy); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_local_energy(qmckl_context context, double * const local_energy) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED; + + if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED; + + rc = qmckl_provide_ao_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_local_energy(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.walk_num * sizeof(double); + memcpy(local_energy, ctx->local_energy.e_local, sze); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_local_energy(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_local_energy(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->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + qmckl_exit_code rc; + rc = qmckl_provide_kinetic_energy(context); + if(rc != QMCKL_SUCCESS){ + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_kinetic_energy", + NULL); + } + + rc = qmckl_provide_potential_energy(context); + if(rc != QMCKL_SUCCESS){ + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_potential_energy", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->local_energy.e_local_date) { + + /* Allocate array */ + if (ctx->local_energy.e_local == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* local_energy = (double*) qmckl_malloc(context, mem_info); + + if (local_energy == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_local_energy", + NULL); + } + ctx->local_energy.e_local = local_energy; + } + + if (ctx->det.type == 'G') { + rc = qmckl_compute_local_energy(context, + ctx->det.walk_num, + ctx->local_energy.e_kin, + ctx->local_energy.e_pot, + ctx->local_energy.e_local); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_local_energy", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->local_energy.e_local_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute local enregy + :PROPERTIES: + :Name: qmckl_compute_local_energy + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_compute_local_energy_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~double~ | ~e_kin[walk_num]~ | in | e kinetic | + | ~double~ | ~e_pot[walk_num]~ | in | e potential | + | ~double~ | ~e_local[walk_num]~ | out | local energy | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_local_energy_f(context, walk_num, & + e_kin, e_pot, e_local) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: walk_num + double precision, intent(in) :: e_kin(walk_num) + double precision, intent(in) :: e_pot(walk_num) + double precision, intent(inout) :: e_local(walk_num) + integer*8 :: idet, iwalk, ielec, mo_id, imo + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + e_local = 0.0d0 + do iwalk = 1, walk_num + e_local(iwalk) = e_local(iwalk) + e_kin(iwalk) + e_pot(iwalk) + end do + +end function qmckl_compute_local_energy_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_compute_local_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_local_energy")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_local_energy ( + const qmckl_context context, + const int64_t walk_num, + const double* e_kin, + const double* e_pot, + double* const e_local ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_compute_local_energy_args,rettyp=get_value("CRetType"),fname="qmckl_compute_local_energy")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_local_energy & + (context, walk_num, e_kin, e_pot, e_local) & + 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 :: walk_num + real (c_double ) , intent(in) :: e_kin(walk_num) + real (c_double ) , intent(in) :: e_pot(walk_num) + real (c_double ) , intent(out) :: e_local(walk_num) + + integer(c_int32_t), external :: qmckl_compute_local_energy_f + info = qmckl_compute_local_energy_f & + (context, walk_num, e_kin, e_pot, e_local) + + end function qmckl_compute_local_energy + #+end_src + +*** Test + #+begin_src c :tangle (eval c_test) :exports none +// Calculate the Local energy + +double local_energy[walk_num]; + +rc = qmckl_get_local_energy(context, &(local_energy[0])); +assert (rc == QMCKL_SUCCESS); + + #+end_src + +** Drift vector + :PROPERTIES: + :Name: qmckl_compute_drift_vector + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + +The drift vector is calculated as the ration of the gradient +with the determinant of the wavefunction. + +\[ +\mathbf{F} = 2 \frac{\nabla \Psi}{\Psi} +\] + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double* const drift_vector); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_drift_vector(qmckl_context context, double * const drift_vector) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + if(!qmckl_electron_provided(context)) return QMCKL_NOT_PROVIDED; + + if(!qmckl_nucleus_provided(context)) return QMCKL_NOT_PROVIDED; + + rc = qmckl_provide_ao_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_drift_vector(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.walk_num * ctx->electron.num * 3 * sizeof(double); + memcpy(drift_vector, ctx->local_energy.r_drift, sze); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_drift_vector(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_drift_vector(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->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if(!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + if (!ctx->det.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->local_energy.r_drift_date) { + + /* Allocate array */ + if (ctx->local_energy.r_drift == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * ctx->electron.num * 3 * sizeof(double); + double* r_drift = (double*) qmckl_malloc(context, mem_info); + + if (r_drift == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_r_drift", + NULL); + } + ctx->local_energy.r_drift = r_drift; + } + + qmckl_exit_code rc; + if (ctx->det.type == 'G') { + rc = qmckl_compute_drift_vector(context, + ctx->det.walk_num, + ctx->det.det_num_alpha, + ctx->det.det_num_beta, + ctx->electron.up_num, + ctx->electron.down_num, + ctx->electron.num, + ctx->det.mo_index_alpha, + ctx->det.mo_index_beta, + ctx->mo_basis.mo_num, + ctx->mo_basis.mo_vgl, + ctx->det.det_inv_matrix_alpha, + ctx->det.det_inv_matrix_beta, + ctx->local_energy.r_drift); + } else { + return qmckl_failwith( context, + QMCKL_FAILURE, + "compute_drift_vector", + "Not yet implemented"); + } + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->local_energy.r_drift_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute drift vector + :PROPERTIES: + :Name: qmckl_compute_drift_vector + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_compute_drift_vector_args + | ~qmckl_context~ | ~context~ | in | Global state | + | ~int64_t~ | ~walk_num~ | in | Number of walkers | + | ~int64_t~ | ~det_num_alpha~ | in | Number of determinants | + | ~int64_t~ | ~det_num_beta~ | in | Number of determinants | + | ~int64_t~ | ~alpha_num~ | in | Number of electrons | + | ~int64_t~ | ~beta_num~ | in | Number of electrons | + | ~int64_t~ | ~elec_num~ | in | Number of electrons | + | ~int64_t~ | ~mo_index_alpha[det_num_alpha][walk_num][alpha_num]~ | in | MO indices for electrons | + | ~int64_t~ | ~mo_index_beta[det_num_beta][walk_num][beta_num]~ | in | MO indices for electrons | + | ~int64_t~ | ~mo_num~ | in | Number of MOs | + | ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | in | Value, gradients and Laplacian of the MOs | + | ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in | Value, gradients and Laplacian of the Det | + | ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det | + | ~double~ | ~r_drift[walk_num][elec_num][3]~ | out | Kinetic energy | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_drift_vector_f(context, walk_num, & + det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, & + mo_num, mo_vgl, det_inv_matrix_alpha, det_inv_matrix_beta, r_drift) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8, intent(in) :: walk_num + integer*8, intent(in) :: det_num_alpha + integer*8, intent(in) :: det_num_beta + integer*8, intent(in) :: alpha_num + integer*8, intent(in) :: beta_num + integer*8, intent(in) :: elec_num + integer*8, intent(in) :: mo_num + integer*8, intent(in) :: mo_index_alpha(alpha_num, walk_num, det_num_alpha) + integer*8, intent(in) :: mo_index_beta(beta_num, walk_num, det_num_beta) + double precision, intent(in) :: mo_vgl(mo_num, elec_num, 5) + double precision, intent(in) :: det_inv_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha) + double precision, intent(in) :: det_inv_matrix_beta(beta_num, beta_num, walk_num, det_num_beta) + double precision, intent(inout) :: r_drift(3,elec_num,walk_num) + integer*8 :: idet, iwalk, ielec, mo_id, imo + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (alpha_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + if (beta_num < 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + r_drift = 0.0d0 + do idet = 1, det_num_alpha + do iwalk = 1, walk_num + ! Alpha part + do imo = 1, alpha_num + do ielec = 1, alpha_num + mo_id = mo_index_alpha(imo, iwalk, idet) + r_drift(1,ielec,iwalk) = r_drift(1,ielec,iwalk) + 2.0d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, ielec, 2) + r_drift(2,ielec,iwalk) = r_drift(2,ielec,iwalk) + 2.0d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, ielec, 3) + r_drift(3,ielec,iwalk) = r_drift(3,ielec,iwalk) + 2.0d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, ielec, 4) + end do + end do + ! Beta part + do imo = 1, beta_num + do ielec = 1, beta_num + mo_id = mo_index_beta(imo, iwalk, idet) + r_drift(1,alpha_num + ielec,iwalk) = r_drift(1,alpha_num + ielec,iwalk) + & + 2.0d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, alpha_num + ielec, 2) + r_drift(2,alpha_num + ielec,iwalk) = r_drift(2,alpha_num + ielec,iwalk) + & + 2.0d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, alpha_num + ielec, 3) + r_drift(3,alpha_num + ielec,iwalk) = r_drift(3,alpha_num + ielec,iwalk) + & + 2.0d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, alpha_num + ielec, 4) + end do + end do + end do + end do + +end function qmckl_compute_drift_vector_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_compute_drift_vector_args,rettyp=get_value("CRetType"),fname="qmckl_compute_drift_vector")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_drift_vector ( + const qmckl_context context, + const int64_t walk_num, + const int64_t det_num_alpha, + const int64_t det_num_beta, + const int64_t alpha_num, + const int64_t beta_num, + const int64_t elec_num, + const int64_t* mo_index_alpha, + const int64_t* mo_index_beta, + const int64_t mo_num, + const double* mo_vgl, + const double* det_inv_matrix_alpha, + const double* det_inv_matrix_beta, + double* const r_drift ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_compute_drift_vector_args,rettyp=get_value("CRetType"),fname="qmckl_compute_drift_vector")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_drift_vector & + (context, & + walk_num, & + det_num_alpha, & + det_num_beta, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_alpha, & + mo_index_beta, & + mo_num, & + mo_vgl, & + det_inv_matrix_alpha, & + det_inv_matrix_beta, & + r_drift) & + 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 :: walk_num + integer (c_int64_t) , intent(in) , value :: det_num_alpha + integer (c_int64_t) , intent(in) , value :: det_num_beta + integer (c_int64_t) , intent(in) , value :: alpha_num + integer (c_int64_t) , intent(in) , value :: beta_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) :: mo_index_alpha(alpha_num,walk_num,det_num_alpha) + integer (c_int64_t) , intent(in) :: mo_index_beta(beta_num,walk_num,det_num_beta) + integer (c_int64_t) , intent(in) , value :: mo_num + real (c_double ) , intent(in) :: mo_vgl(mo_num,elec_num,5) + real (c_double ) , intent(in) :: det_inv_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha) + real (c_double ) , intent(in) :: det_inv_matrix_beta(beta_num,beta_num,walk_num,det_num_beta) + real (c_double ) , intent(out) :: r_drift(3,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_drift_vector_f + info = qmckl_compute_drift_vector_f & + (context, & + walk_num, & + det_num_alpha, & + det_num_beta, & + alpha_num, & + beta_num, & + elec_num, & + mo_index_alpha, & + mo_index_beta, & + mo_num, & + mo_vgl, & + det_inv_matrix_alpha, & + det_inv_matrix_beta, & + r_drift) + + end function qmckl_compute_drift_vector + #+end_src + +*** Test + #+begin_src c :tangle (eval c_test) :exports none +// Calculate the Drift vector + +double drift_vector[walk_num][3]; + +rc = qmckl_get_drift_vector(context, &(drift_vector[0][0])); +assert (rc == QMCKL_SUCCESS); + #+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/qmckl_mo.org b/org/qmckl_mo.org index 61e619f..60f39f4 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -83,6 +83,7 @@ int main() { The following arrays are stored in the context: + |---------------+--------------------+----------------------| | ~mo_num~ | | Number of MOs | | ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients | @@ -117,6 +118,26 @@ typedef struct qmckl_mo_basis_struct { Some values are initialized by default, and are not concerned by this mechanism. + #+begin_src c :comments org :tangle (eval h_private_func) +qmckl_exit_code qmckl_init_mo_basis(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) +qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return false; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + ctx->mo_basis.uninitialized = (1 << 2) - 1; + + return QMCKL_SUCCESS; +} + #+end_src + ** Access functions #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -258,10 +279,9 @@ qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; 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); + qmckl_exit_code rc_ = qmckl_finalize_mo_basis(context); if (rc_ != QMCKL_SUCCESS) return rc_; - } - +} return QMCKL_SUCCESS; #+end_src @@ -318,6 +338,34 @@ qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const dou When the basis set is completely entered, other data structures are computed to accelerate the calculations. + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_finalize_mo_basis", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + qmckl_exit_code rc; + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_finalize_mo_basis", + NULL); + } + return rc; +} + #+end_src + * Computation ** Computation of MOs @@ -377,6 +425,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context); qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) { + qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } @@ -391,6 +440,14 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) NULL); } + rc = qmckl_provide_ao_vgl(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + if(!(ctx->electron.provided)) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, @@ -424,7 +481,6 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) ctx->mo_basis.mo_vgl = mo_vgl; } - qmckl_exit_code rc; rc = qmckl_compute_mo_basis_vgl(context, ctx->ao_basis.ao_num, ctx->mo_basis.mo_num, @@ -450,7 +506,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) :FRetType: qmckl_exit_code :END: - #+NAME: qmckl_mo_basis_gaussian_vgl_args + #+NAME: qmckl_mo_basis_vgl_args | ~qmckl_context~ | ~context~ | in | Global state | | ~int64_t~ | ~ao_num~ | in | Number of AOs | | ~int64_t~ | ~mo_num~ | in | Number of MOs | @@ -458,6 +514,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) | ~double~ | ~coef_normalized[mo_num][ao_num]~ | in | AO to MO transformation matrix | | ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | in | Value, gradients and Laplacian of the AOs | | ~double~ | ~mo_vgl[5][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_vgl_f(context, & @@ -473,85 +530,70 @@ integer function qmckl_compute_mo_basis_vgl_f(context, & double precision , intent(in) :: coef_normalized(ao_num,mo_num) double precision , intent(out) :: mo_vgl(mo_num,elec_num,5) logical*8 :: TransA, TransB + double precision,dimension(:,:),allocatable :: mo_vgl_big + double precision,dimension(:,:),allocatable :: ao_vgl_big + !double precision,dimension(:,:),allocatable :: coef_trans + !double precision,dimension(:),allocatable :: coef_all 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, i,j + integer*8 :: M, N, K, LDA, LDB, LDC, i,j, idx integer*8 :: inucl, iprim, iwalk, ielec, ishell double precision :: x, y, z, two_a, ar2, r2, v, cutoff - TransA = .False. + + allocate(mo_vgl_big(mo_num,elec_num*5)) + allocate(ao_vgl_big(ao_num,elec_num*5)) + !allocate(coef_all(mo_num*ao_num)) + !allocate(coef_trans(mo_num,ao_num)) + + TransA = .True. TransB = .False. alpha = 1.0d0 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 = 1_8 - N = mo_num * 1_8 + M = mo_num + N = elec_num*5 K = ao_num * 1_8 - LDA = M - LDB = K - LDC = M + LDA = size(coef_normalized,1) + idx = 0 + !do j = 1,ao_num + !do i = 1,mo_num + ! idx = idx + 1 + ! coef_all(idx) = coef_normalized(i,j) + !end do + !end do + !idx = 0 + !do j = 1,mo_num + !do i = 1,ao_num + ! idx = idx + 1 + ! coef_trans(j,i) = coef_all(idx) + !end do + !end do - do ielec = 1, elec_num - ! Value - info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, 1), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,1),LDC) - ! Grad_x - info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, 2), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,2),LDC) - ! Grad_y - info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, 3), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,3),LDC) - ! Grad_z - info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, 4), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,4),LDC) - ! Lapl_z - info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, 5), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,5),LDC) - end do + ao_vgl_big = reshape(ao_vgl(:, :, :),(/ao_num, elec_num*5_8/)) + LDB = size(ao_vgl_big,1) + LDC = size(mo_vgl_big,1) + + info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + coef_normalized,size(coef_normalized,1)*1_8, & + ao_vgl_big, size(ao_vgl_big,1)*1_8, & + beta, & + mo_vgl_big,LDC) + mo_vgl = reshape(mo_vgl_big,(/mo_num,elec_num,5_8/)) + + deallocate(mo_vgl_big) + deallocate(ao_vgl_big) - 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_vgl_f #+end_src - #+CALL: generate_c_header(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl")) + #+CALL: generate_c_header(table=qmckl_mo_basis_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org @@ -566,7 +608,7 @@ end function qmckl_compute_mo_basis_vgl_f #+end_src - #+CALL: generate_c_interface(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl")) + #+CALL: generate_c_interface(table=qmckl_mo_basis_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_vgl")) #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none @@ -785,10 +827,10 @@ rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0])); assert (rc == QMCKL_SUCCESS); // Test overlap of MO -//double point_x[100]; -//double point_y[100]; -//double point_z[100]; -//int32_t npoints=100; +//double point_x[10]; +//double point_y[10]; +//double point_z[10]; +//int32_t npoints=10; //// obtain points //double dr = 20./(npoints-1); //double dr3 = dr*dr*dr; @@ -802,10 +844,11 @@ assert (rc == QMCKL_SUCCESS); //double ovlmo1 = 0.0; //// Calculate overlap //for (int i=0;i 1e-5) then + energy = energy + charge(i) * charge(j) / nn_distance(i,j) + endif end do end do diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org index 9662938..4db1837 100644 --- a/org/qmckl_trexio.org +++ b/org/qmckl_trexio.org @@ -704,6 +704,42 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file) } #+end_src +*** AO Normalization + #+begin_src c :tangle (eval c) + { + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ao_num * sizeof(double); + + double* ao_normalization = (double*) qmckl_malloc(context, mem_info); + + if (ao_normalization == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_trexio_read_ao_normalization_X", + NULL); + } + + assert (ao_normalization != NULL); + + rcio = trexio_read_ao_normalization_64(file, ao_normalization); + if (rcio != TREXIO_SUCCESS) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "trexio_read_ao_normalization", + trexio_string_of_error(rcio)); + } + + rc = qmckl_set_ao_basis_ao_factor(context, ao_normalization); + + qmckl_free(context, ao_normalization); + ao_normalization = NULL; + + if (rc != QMCKL_SUCCESS) + return rc; + } + #+end_src + + #+begin_src c :tangle (eval c) @@ -711,7 +747,6 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file) } #endif #+end_src - ** Molecular orbitals In this section we read the MO coefficients. @@ -880,6 +915,7 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name) #+begin_src c :tangle (eval c_test) #ifdef HAVE_TREXIO +#define walk_num 2 qmckl_exit_code rc; char fname[256]; @@ -892,6 +928,8 @@ char message[256]; strncpy(fname, QMCKL_TEST_DIR,255); strncat(fname, "/chbrclf", 255); printf("Test file: %s\n", fname); + +rc = qmckl_set_electron_walk_num(context, walk_num); rc = qmckl_trexio_read(context, fname); if (rc != QMCKL_SUCCESS) { diff --git a/org/table_of_contents b/org/table_of_contents index 2a74909..4551964 100644 --- a/org/table_of_contents +++ b/org/table_of_contents @@ -1,17 +1,20 @@ qmckl.org -qmckl_error.org +qmckl_ao.org +qmckl_blas.org qmckl_context.org +qmckl_determinant.org +qmckl_distance.org +qmckl_electron.org +qmckl_error.org +qmckl_jastrow.org +qmckl_local_energy.org qmckl_memory.org +qmckl_mo.org qmckl_numprec.org qmckl_distance.org qmckl_nucleus.org -qmckl_electron.org -qmckl_ao.org -qmckl_mo.org -qmckl_jastrow.org qmckl_sherman_morrison_woodbury.org qmckl_utils.org -qmckl_blas.org qmckl_trexio.org qmckl_verificarlo.org qmckl_tests.org