Verified Local energy with QMC=Chem for Be2. #41

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.
@@ -90,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
@@ -114,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
@@ -161,26 +162,75 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
   endif
 
    if (TransA) then
-
+     info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, B, LDB_2, beta, c, LDC)
-     if (alpha == 1.0d0 .and. beta == 0.0d0) then
-       C = matmul(AT,B)
-     else
-       C = beta*C + alpha*matmul(AT,B)
-     endif
    else if (TransB) then
-     if (alpha == 1.0d0 .and. beta == 0.0d0) 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
-     else
+     info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, BT, LDB_2, beta, c, LDC)
-       C = beta*C + alpha*matmul(A,BT)
-     endif
    else
-     if (alpha == 1.0d0 .and. beta == 0.0d0) then
+     info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC)
-       C = matmul(A,B)
-     else
-       C = beta*C + alpha*matmul(A,B)
-     endif
    endif
 end function qmckl_dgemm_f
+
+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:
@@ -393,6 +443,7 @@ integer function qmckl_invert_f(context, ma, na, LDA, A, det_l) &
   case (4)
 !DIR$ forceinline
     call invert4(a,LDA,na,det_l)
+    
   case (3)
 !DIR$ forceinline
     call invert3(a,LDA,na,det_l)
@@ -414,11 +465,184 @@ subroutine invert1(a,LDA,na,det_l)
  integer*8, intent(in)             :: na
  double precision, intent(inout) :: det_l
 
+ call cofactor1(a,LDA,na,det_l)
+end
+
+subroutine invert2(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 invert3(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 invert4(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 invert5(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 invert2(a,LDA,na,det_l)
+subroutine cofactor2(a,LDA,na,det_l)
  implicit none
  double precision :: a (LDA,na)
  integer*8        :: LDA
@@ -437,7 +661,7 @@ subroutine invert2(a,LDA,na,det_l)
  a(2,2) = b(1,1)
 end
 
-subroutine invert3(a,LDA,na,det_l)
+subroutine cofactor3(a,LDA,na,det_l)
  implicit none
  double precision, intent(inout) :: a (LDA,na)
  integer*8, intent(in)             :: LDA
@@ -468,7 +692,7 @@ subroutine invert3(a,LDA,na,det_l)
 
 end
 
-subroutine invert4(a,LDA,na,det_l)
+subroutine cofactor4(a,LDA,na,det_l)
  implicit none
  double precision, intent(inout) :: a (LDA,na)
  integer*8, intent(in)             :: LDA
@@ -518,7 +742,7 @@ subroutine invert4(a,LDA,na,det_l)
 
 end
 
-subroutine invert5(a,LDA,na,det_l)
+subroutine cofactor5(a,LDA,na,det_l)
  implicit none
  double precision, intent(inout) :: a (LDA,na)
  integer*8, intent(in)             :: LDA

org/qmckl_determinant.org

@@ -1040,19 +1040,19 @@ integer function qmckl_compute_det_vgl_beta_f(context, &
       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, ielec, 1)
+        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, ielec, 2)
+        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, ielec, 3)
+        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, ielec, 4)
+        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, ielec, 5)
+        det_vgl_beta(imo, ielec, 5, iwalk, idet) = mo_vgl(mo_id, alpha_num + ielec, 5)
       end do
     end do
   end do
@@ -1726,7 +1726,7 @@ integer function qmckl_compute_det_inv_matrix_alpha_f(context, &
   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
+  integer*8 :: idet, iwalk, ielec, mo_id, imo, LDA, res, i, j
 
   allocate(matA(alpha_num, alpha_num))
 
@@ -1756,7 +1756,7 @@ integer function qmckl_compute_det_inv_matrix_alpha_f(context, &
   do idet = 1, det_num_alpha
   do iwalk = 1, walk_num
     ! Value
-    matA = det_vgl_alpha(1:alpha_num, 1:alpha_num, 1, iwalk, idet)
+    matA(1:alpha_num,1:alpha_num) = det_vgl_alpha(1:alpha_num, 1:alpha_num, 1, iwalk, idet)
     res = qmckl_invert(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
@@ -1887,7 +1887,7 @@ integer function qmckl_compute_det_inv_matrix_beta_f(context, &
   do idet = 1, det_num_beta
   do iwalk = 1, walk_num
     ! Value
-    matA = det_vgl_beta(1:beta_num, 1:beta_num, 1, iwalk, idet)
+    matA(1:beta_num,1:beta_num) = det_vgl_beta(1:beta_num, 1:beta_num, 1, iwalk, idet)
     res = qmckl_invert(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

org/qmckl_electron.org

@@ -1727,7 +1727,7 @@ integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, &
 
   ee_pot = 0.0d0
   do nw=1,walk_num
-    do j=1,elec_num
+    do j=2,elec_num
       do i=1,j-1
         ee_pot(nw) = ee_pot(nw) + 1.0d0/(ee_distance(i,j,nw))
       end do

org/qmckl_local_energy.org

@@ -245,6 +245,7 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context);
     #+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;
   }
@@ -286,6 +287,21 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
                            "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) {
@@ -306,7 +322,6 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
       ctx->local_energy.e_kin = e_kin;
     }
 
-    qmckl_exit_code rc; 
     if (ctx->det.type == 'G') {
       rc = qmckl_compute_kinetic_energy(context,
                                  ctx->det.walk_num,
@@ -319,8 +334,10 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
                                  ctx->det.mo_index_beta,
                                  ctx->mo_basis.mo_num,
                                  ctx->mo_basis.mo_vgl,
-                                 ctx->det.det_adj_matrix_alpha,
+                                 ctx->det.det_value_alpha,
-                                 ctx->det.det_adj_matrix_beta,
+                                 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,
@@ -358,14 +375,16 @@ qmckl_exit_code qmckl_provide_kinetic_energy(qmckl_context context) {
    | ~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_adj_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in  | Value, gradients and Laplacian of the Det |
+   | ~double~        | ~det_value_alpha[det_num_alpha][walk_num]~                            | in  | Det of wavefunction                       |
-   | ~double~        | ~det_adj_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~     | in  | Value, gradients and Laplacian of the Det |
+   | ~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_adj_matrix_alpha, det_adj_matrix_beta, e_kin) &
+     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
@@ -380,9 +399,12 @@ integer function qmckl_compute_kinetic_energy_f(context, walk_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_adj_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha)
+  double precision, intent(in)      :: det_value_alpha(walk_num, det_num_alpha)
-  double precision, intent(in)      :: det_adj_matrix_beta(beta_num, beta_num, walk_num, det_num_beta)
+  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
@@ -416,20 +438,34 @@ integer function qmckl_compute_kinetic_energy_f(context, walk_num, &
   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(ielec, iwalk, idet)
+      mo_id = mo_index_alpha(imo, iwalk, idet)
-      e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_adj_matrix_alpha(imo, ielec, 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(ielec, iwalk, idet)
+      mo_id = mo_index_beta(imo, iwalk, idet)
-      e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_adj_matrix_beta(imo, ielec, iwalk, idet) * &
+      e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * &
-                                     mo_vgl(mo_id, ielec, 5)
+                                     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
@@ -453,8 +489,10 @@ end function qmckl_compute_kinetic_energy_f
 	 const int64_t* mo_index_beta,
 	 const int64_t mo_num,
 	 const double* mo_vgl,
-	 const double* det_adj_matrix_alpha,
+	 const double* det_value_alpha,
-	 const double* det_adj_matrix_beta,
+	 const double* det_value_beta,
+	 const double* det_inv_matrix_alpha,
+	 const double* det_inv_matrix_beta,
 	 double* const e_kin ); 
    #+end_src
 
@@ -474,8 +512,10 @@ end function qmckl_compute_kinetic_energy_f
 		 mo_index_beta, &
 		 mo_num, &
 		 mo_vgl, &
-		 det_adj_matrix_alpha, &
+		 det_value_alpha, &
-		 det_adj_matrix_beta, &
+		 det_value_beta, &
+		 det_inv_matrix_alpha, &
+		 det_inv_matrix_beta, &
 		 e_kin) &
 	bind(C) result(info)
 
@@ -493,8 +533,10 @@ end function qmckl_compute_kinetic_energy_f
       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_adj_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha)
+      real    (c_double ) , intent(in)          :: det_value_alpha(walk_num,det_num_alpha)
-      real    (c_double ) , intent(in)          :: det_adj_matrix_beta(beta_num,beta_num,walk_num,det_num_beta)
+      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
@@ -510,8 +552,10 @@ end function qmckl_compute_kinetic_energy_f
 		 mo_index_beta, &
 		 mo_num, &
 		 mo_vgl, &
-		 det_adj_matrix_alpha, &
+		 det_value_alpha, &
-		 det_adj_matrix_beta, &
+		 det_value_beta, &
+		 det_inv_matrix_alpha, &
+		 det_inv_matrix_beta, &
 		 e_kin)
 
     end function qmckl_compute_kinetic_energy
@@ -830,14 +874,6 @@ qmckl_exit_code qmckl_provide_potential_energy(qmckl_context context) {
                            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_nucleus_repulsion(context);
   if (rc != QMCKL_SUCCESS) {
     return qmckl_failwith( context,
@@ -958,6 +994,7 @@ integer function qmckl_compute_potential_energy_f(context, walk_num, &
   endif
 
   e_pot = 0.0d0 + repulsion
+  print *,repulsion
   do iwalk = 1, walk_num
     e_pot(iwalk) = e_pot(iwalk) + ee_pot(iwalk) + en_pot(iwalk)
   end do

org/qmckl_mo.org

@@ -511,7 +511,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
    | ~int64_t~       | ~ao_num~                          | in  | Number of AOs                             |
    | ~int64_t~       | ~mo_num~                          | in  | Number of MOs                             |
    | ~int64_t~       | ~elec_num~                        | in  | Number of electrons                       |
-   | ~double~        | ~coef_normalized[ao_num][mo_num]~ | in  | AO to MO transformation matrix            |
+   | ~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 |
 
@@ -527,22 +527,26 @@ integer function qmckl_compute_mo_basis_vgl_f(context, &
   integer*8             , intent(in)  :: ao_num, mo_num
   integer*8             , intent(in)  :: elec_num
   double precision      , intent(in)  :: ao_vgl(ao_num,elec_num,5)
-  double precision      , intent(in)  :: coef_normalized(mo_num,ao_num)
+  double precision      , intent(in)  :: coef_normalized(ao_num,mo_num)
   double precision      , intent(out) :: mo_vgl(mo_num,elec_num,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*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
 
   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 = .False.
+  TransA = .True.
   TransB = .False. 
   alpha = 1.0d0
   beta  = 0.0d0
@@ -556,15 +560,29 @@ integer function qmckl_compute_mo_basis_vgl_f(context, &
   M = mo_num
   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
 
   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,    &
+                                 coef_normalized,size(coef_normalized,1)*1_8,    &
-                                 ao_vgl_big, LDB, &
+                                 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/))
@@ -605,8 +623,7 @@ end function qmckl_compute_mo_basis_vgl_f
       integer (c_int64_t) , intent(in)  , value :: ao_num
       integer (c_int64_t) , intent(in)  , value :: mo_num
       integer (c_int64_t) , intent(in)  , value :: elec_num
-      real    (c_double ) , intent(in)          :: coef_normalized(mo_num,ao_num)
+      real    (c_double ) , intent(in)          :: coef_normalized(ao_num,mo_num)
-
       real    (c_double ) , intent(in)          :: ao_vgl(ao_num,elec_num,5)
       real    (c_double ) , intent(out)         :: mo_vgl(mo_num,elec_num,5)