Merge branch 'master' into add_lds

org/qmckl_blas.org

@@ -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))
+  allocate(AT(m,k))
-    do i = 1, m
+    do i = 1, k
-      do j = 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))
+  allocate(BT(k,n))
-    do i = 1, k
+    do i = 1, n
-      do j = 1, n
+      do j = 1, k
         BT(j,i) = B(i,j)
       end do
     end do
@@ -162,25 +162,75 @@ 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
+     info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, B, LDB_2, beta, c, LDC)
-       C = matmul(AT,B)
-     else
-       C = beta*C + alpha*matmul(AT,B)
-     endif
    else if (TransB) then
-     if (alpha == 1.d0 .and. beta == 0.d0) then
+     info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, BT, LDB_2, beta, c, LDC)
-       C = matmul(A,BT)
+   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
-       C = beta*C + alpha*matmul(A,BT)
+     info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC)
    endif
-   else
+end function qmckl_dgemm_f
-     if (alpha == 1.d0 .and. beta == 0.d0) then
+
+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
-   endif
+end function qmckl_dgemm_N_N_f
-end function qmckl_dgemm_f
     #+end_src
     
 *** C interface                                                    :noexport:
@@ -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)

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)
@@ -123,10 +127,10 @@ typedef struct qmckl_context_struct {
   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 */

org/qmckl_electron.org

@@ -68,6 +68,7 @@ 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               |
@@ -79,6 +80,10 @@ int main() {
   | ~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                                           |
@@ -87,6 +92,10 @@ int main() {
   | ~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                     |
@@ -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;
@@ -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)

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                             |
@@ -459,6 +515,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
    | ~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, &
      ao_num, mo_num, elec_num, &
@@ -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*8                     :: M, N, K, LDA, LDB, LDC, i,j, idx
-  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 :: 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
+  M = mo_num
-  N = mo_num * 1_8
+  N = elec_num*5
   K = ao_num * 1_8
-  LDA = M
+  LDA = size(coef_normalized,1)
-  LDB = K
+  idx = 0
-  LDC = M
+  !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
+  ao_vgl_big = reshape(ao_vgl(:, :, :),(/ao_num, elec_num*5_8/))
-          ! Value
+  LDB = size(ao_vgl_big,1) 
-          info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     &
+  LDC = size(mo_vgl_big,1)
-                                         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
 
-  if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and.  &
+  info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     &
-     info_qmckl_dgemm_Gx    .eq. QMCKL_SUCCESS .and.  &
+                                 coef_normalized,size(coef_normalized,1)*1_8,    &
-     info_qmckl_dgemm_Gy    .eq. QMCKL_SUCCESS .and.  &
+                                 ao_vgl_big, size(ao_vgl_big,1)*1_8, &
-     info_qmckl_dgemm_Gz    .eq. QMCKL_SUCCESS .and.  &
+                                 beta,                                       &
-     info_qmckl_dgemm_lap   .eq. QMCKL_SUCCESS        ) then
+                                 mo_vgl_big,LDC)
-     info = QMCKL_SUCCESS
+  mo_vgl = reshape(mo_vgl_big,(/mo_num,elec_num,5_8/))
-   else
+
-     info = QMCKL_FAILURE
+  deallocate(mo_vgl_big)
-   end if
+  deallocate(ao_vgl_big)
 
 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_x[10];
-//double point_y[100];
+//double point_y[10];
-//double point_z[100];
+//double point_z[10];
-//int32_t npoints=100;
+//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<npoints;++i) {
-//  printf(".");
 //  fflush(stdout);
 //  for (int j=0;j<npoints;++j) {
+//  printf(" .. ");
 //    for (int k=0;k<npoints;++k) {
+//  printf(" . ");
 //      // Set point
 //      elec_coord[0] = point_x[i];
 //      elec_coord[1] = point_y[j];

org/qmckl_nucleus.org

@@ -67,18 +67,24 @@ int main() {
 
   The following data stored in the context:
 
+  |------------------------+----------------+-------------------------------------------|
   | ~uninitialized~        | int32_t        | Keeps bit set for uninitialized data      |
   | ~num~                  | int64_t        | Total number of nuclei                    |
   | ~provided~             | bool           | If true, ~nucleus~ is valid               |
   | ~charge~               | double[num]    | Nuclear charges                           |
   | ~coord~                | double[3][num] | Nuclear coordinates, in transposed format |
+  | ~coord_date~           | int64_t        | Nuclear coordinates, date if modified     |
+  | ~rescale_factor_kappa~ | double         | The distance scaling factor               |
+
+  Computed data:
+  
+  |-----------------------------+------------------+------------------------------------------------------------|
   | ~nn_distance~               | double[num][num] | Nucleus-nucleus distances                                  |
   | ~nn_distance_date~          | int64_t          | Date when Nucleus-nucleus distances were computed          |
   | ~nn_distance_rescaled~      | double[num][num] | Nucleus-nucleus rescaled distances                         |
   | ~nn_distance_rescaled_date~ | int64_t          | Date when Nucleus-nucleus rescaled distances were computed |
   | ~repulsion~                 | double           | Nuclear repulsion energy                                   |
   | ~repulsion_date~            | int64_t          | Date when the nuclear repulsion energy was computed        |
-   | ~rescale_factor_kappa~      | double           | The distance scaling factor                                |
 
 ** Data structure
 
@@ -88,6 +94,7 @@ typedef struct qmckl_nucleus_struct {
   int64_t   repulsion_date;
   int64_t   nn_distance_date;
   int64_t   nn_distance_rescaled_date;
+  int64_t   coord_date;
   double*   coord;
   double*   charge;
   double*   nn_distance;
@@ -130,8 +137,6 @@ qmckl_exit_code qmckl_init_nucleus(qmckl_context context) {
 }
    #+end_src
 
-
-
 ** Access functions
 
    #+begin_src c :comments org :tangle (eval h_func) :exports none
@@ -790,7 +795,6 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
 
      #+end_src
 
-
 ** Nucleus-nucleus rescaled distances
 
 *** Get
@@ -1083,7 +1087,9 @@ integer function qmckl_compute_nucleus_repulsion_f(context, nucl_num, charge, nn
   energy = 0.d0
   do j=2, nucl_num
      do i=1, j-1
+        if (dabs(nn_distance(i,j)) > 1e-5) then
           energy = energy + charge(i) * charge(j) / nn_distance(i,j)
+        endif
      end do
   end do
 

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) {

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