TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-18 17:03:43 +02:00
Code Issues Projects Releases Wiki Activity

Merge pull request #42 from v1j4y/wf_det_grad_cof

Browse Source 
Slater Determinant and Local Energy
This commit is contained in:
Anthony Scemama 2021-11-17 14:07:34 +01:00 committed by GitHub
commit ec5a6baf4a
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 5183 additions and 152 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

781
org/qmckl_blas.org
View File

@ -79,7 +79,6 @@ int main() {
const int64_t ldc ); const int64_t ldc );
#+END_src #+END_src
*** Source *** Source
#+begin_src f90 :tangle (eval f) #+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) & 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 integer*8 , intent(in) :: m, n, k
real*8 , intent(in) :: alpha, beta real*8 , intent(in) :: alpha, beta
integer*8 , intent(in) :: lda integer*8 , intent(in) :: lda
real*8 , intent(in) :: A(m,k) real*8 , intent(in) :: A(lda,*)
integer*8 , intent(in) :: ldb integer*8 , intent(in) :: ldb
real*8 , intent(in) :: B(k,n) real*8 , intent(in) :: B(ldb,*)
integer*8 , intent(in) :: ldc integer*8 , intent(in) :: ldc
real*8 , intent(out) :: C(m,n) real*8 , intent(out) :: C(ldc,*)
real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:) real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:)
integer*4 :: qmckl_dgemm_N_N_f
integer*8 :: i,j,l, LDA_2, LDB_2 integer*8 :: i,j,l, LDA_2, LDB_2
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
if (TransA) then 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) AT(j,i) = A(i,j)
end do end do
end do end do
@ -115,9 +115,9 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
endif endif
if (TransB) then 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) BT(j,i) = B(i,j)
end do end do
end do end do
@ -162,25 +162,75 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
endif endif
if (TransA) then 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 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
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 else
if (alpha == 1.d0 .and. beta == 0.d0) 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 endif
end function qmckl_dgemm_f 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 #+end_src
*** C interface :noexport: *** C interface :noexport:
@ -321,6 +371,687 @@ qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
#+end_src #+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: * End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test) #+begin_src c :comments link :tangle (eval c_test)

24
org/qmckl_context.org
View File

@ -33,10 +33,14 @@ int main() {
#include "qmckl_ao_private_type.h" #include "qmckl_ao_private_type.h"
#include "qmckl_mo_private_type.h" #include "qmckl_mo_private_type.h"
#include "qmckl_jastrow_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_nucleus_private_func.h"
#include "qmckl_electron_private_func.h" #include "qmckl_electron_private_func.h"
#include "qmckl_ao_private_func.h" #include "qmckl_ao_private_func.h"
#include "qmckl_mo_private_func.h" #include "qmckl_mo_private_func.h"
#include "qmckl_determinant_private_func.h"
#include "qmckl_local_energy_private_func.h"
#+end_src #+end_src
#+begin_src c :tangle (eval c) #+begin_src c :tangle (eval c)
@ -118,15 +122,15 @@ typedef struct qmckl_context_struct {
uint64_t date; uint64_t date;
/* -- Molecular system -- */ /* -- Molecular system -- */
qmckl_nucleus_struct nucleus; qmckl_nucleus_struct nucleus;
qmckl_electron_struct electron; qmckl_electron_struct electron;
qmckl_ao_basis_struct ao_basis; qmckl_ao_basis_struct ao_basis;
qmckl_mo_basis_struct mo_basis; qmckl_mo_basis_struct mo_basis;
qmckl_jastrow_struct jastrow; qmckl_jastrow_struct jastrow;
qmckl_determinant_struct det;
qmckl_local_energy_struct local_energy;
/* To be implemented: /* To be implemented:
qmckl_mo_struct mo;
qmckl_determinant_struct det;
,*/ ,*/
} qmckl_context_struct; } qmckl_context_struct;
@ -240,6 +244,12 @@ qmckl_context qmckl_context_create() {
rc = qmckl_init_ao_basis(context); rc = qmckl_init_ao_basis(context);
assert (rc == QMCKL_SUCCESS); 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 */ /* Allocate qmckl_memory_struct */

2073
org/qmckl_determinant.org Normal file
View File

File diff suppressed because it is too large Load Diff

468
org/qmckl_electron.org
View File

@ -68,29 +68,38 @@ int main() {
The following data stored in the context: The following data stored in the context:
| ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | |---------------------------+----------------------------+-------------------------------------------|
| ~num~ | ~int64_t~ | Total number of electrons | | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data |
| ~up_num~ | ~int64_t~ | Number of up-spin electrons | | ~num~ | ~int64_t~ | Total number of electrons |
| ~down_num~ | ~int64_t~ | Number of down-spin electrons | | ~up_num~ | ~int64_t~ | Number of up-spin electrons |
| ~walk_num~ | ~int64_t~ | Number of walkers | | ~down_num~ | ~int64_t~ | Number of down-spin electrons |
| ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | | ~walk_num~ | ~int64_t~ | Number of walkers |
| ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor |
| ~provided~ | ~bool~ | If true, ~electron~ is valid | | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor |
| ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates | | ~provided~ | ~bool~ | If true, ~electron~ is valid |
| ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates | | ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates |
| ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | | ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates |
| ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates |
| ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | Computed data:
| ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances | |-------------------------------------+--------------------------------------+----------------------------------------------------------------------|
| ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances |
| ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | | ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | | ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances |
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | | ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives |
| ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
| ~ee_pot~ | ~double[walk_num]~ | Electron-electron rescaled distances derivatives |
| ~ee_pot_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
| ~en_pot~ | double[walk_num] | Electron-nucleus potential energy |
| ~en_pot_date~ | int64_t | Date when the electron-nucleus potential energy was computed |
| ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances |
| ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances |
| ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives |
| ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives |
** Data structure ** Data structure
@ -105,6 +114,8 @@ typedef struct qmckl_electron_struct {
int64_t coord_new_date; int64_t coord_new_date;
int64_t ee_distance_date; int64_t ee_distance_date;
int64_t en_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_date;
int64_t ee_distance_rescaled_deriv_e_date; int64_t ee_distance_rescaled_deriv_e_date;
int64_t en_distance_rescaled_date; int64_t en_distance_rescaled_date;
@ -113,6 +124,8 @@ typedef struct qmckl_electron_struct {
double* coord_old; double* coord_old;
double* ee_distance; double* ee_distance;
double* en_distance; double* en_distance;
double* ee_pot;
double* en_pot;
double* ee_distance_rescaled; double* ee_distance_rescaled;
double* ee_distance_rescaled_deriv_e; double* ee_distance_rescaled_deriv_e;
double* en_distance_rescaled; double* en_distance_rescaled;
@ -1570,6 +1583,205 @@ rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescal
#+end_src #+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 ** Electron-nucleus distances
*** Get *** Get
@ -2407,6 +2619,216 @@ assert (rc == QMCKL_SUCCESS);
#+end_src #+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: * End of files :noexport:
#+begin_src c :tangle (eval h_private_type) #+begin_src c :tangle (eval h_private_type)

1705
org/qmckl_local_energy.org Normal file
View File

File diff suppressed because it is too large Load Diff

179
org/qmckl_mo.org
View File

@ -83,6 +83,7 @@ int main() {
The following arrays are stored in the context: The following arrays are stored in the context:
|---------------+--------------------+----------------------| |---------------+--------------------+----------------------|
| ~mo_num~ | | Number of MOs | | ~mo_num~ | | Number of MOs |
| ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients | | ~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 Some values are initialized by default, and are not concerned by
this mechanism. 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 ** Access functions
#+begin_src c :comments org :tangle (eval h_func) :exports none #+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.uninitialized &= ~mask;
ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0); ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0);
if (ctx->mo_basis.provided) { 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_; if (rc_ != QMCKL_SUCCESS) return rc_;
} }
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
#+end_src #+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 When the basis set is completely entered, other data structures are
computed to accelerate the calculations. 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
** Computation of MOs ** 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 qmckl_provide_mo_vgl(qmckl_context context)
{ {
qmckl_exit_code rc;
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT; return QMCKL_NULL_CONTEXT;
} }
@ -391,6 +440,14 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
NULL); 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)) { if(!(ctx->electron.provided)) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_NOT_PROVIDED, QMCKL_NOT_PROVIDED,
@ -424,7 +481,6 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
ctx->mo_basis.mo_vgl = mo_vgl; ctx->mo_basis.mo_vgl = mo_vgl;
} }
qmckl_exit_code rc;
rc = qmckl_compute_mo_basis_vgl(context, rc = qmckl_compute_mo_basis_vgl(context,
ctx->ao_basis.ao_num, ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num, ctx->mo_basis.mo_num,
@ -450,7 +506,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
:FRetType: qmckl_exit_code :FRetType: qmckl_exit_code
:END: :END:
#+NAME: qmckl_mo_basis_gaussian_vgl_args #+NAME: qmckl_mo_basis_vgl_args
| ~qmckl_context~ | ~context~ | in | Global state | | ~qmckl_context~ | ~context~ | in | Global state |
| ~int64_t~ | ~ao_num~ | in | Number of AOs | | ~int64_t~ | ~ao_num~ | in | Number of AOs |
| ~int64_t~ | ~mo_num~ | in | Number of MOs | | ~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~ | ~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 | | ~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 #+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_vgl_f(context, & integer function qmckl_compute_mo_basis_vgl_f(context, &
ao_num, mo_num, elec_num, & 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(in) :: coef_normalized(ao_num,mo_num)
double precision , intent(out) :: mo_vgl(mo_num,elec_num,5) double precision , intent(out) :: mo_vgl(mo_num,elec_num,5)
logical*8 :: TransA, TransB 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 double precision :: alpha, beta
integer :: info_qmckl_dgemm_value 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 integer*8 :: inucl, iprim, iwalk, ielec, ishell
double precision :: x, y, z, two_a, ar2, r2, v, cutoff 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. TransB = .False.
alpha = 1.0d0 alpha = 1.0d0
beta = 0.0d0 beta = 0.0d0
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
info_qmckl_dgemm_value = 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. ! Don't compute exponentials when the result will be almost zero.
! TODO : Use numerical precision here ! TODO : Use numerical precision here
cutoff = -dlog(1.d-15) 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 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 function qmckl_compute_mo_basis_vgl_f
#+end_src #+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: #+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org #+begin_src c :tangle (eval h_func) :comments org
@ -566,7 +608,7 @@ end function qmckl_compute_mo_basis_vgl_f
#+end_src #+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: #+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none #+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); assert (rc == QMCKL_SUCCESS);
// Test overlap of MO // 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 //// obtain points
//double dr = 20./(npoints-1); //double dr = 20./(npoints-1);
//double dr3 = dr*dr*dr; //double dr3 = dr*dr*dr;
@ -802,10 +844,11 @@ assert (rc == QMCKL_SUCCESS);
//double ovlmo1 = 0.0; //double ovlmo1 = 0.0;
//// Calculate overlap //// Calculate overlap
//for (int i=0;i<npoints;++i) { //for (int i=0;i<npoints;++i) {
// printf(".");
// fflush(stdout); // fflush(stdout);
// for (int j=0;j<npoints;++j) { // for (int j=0;j<npoints;++j) {
// printf(" .. ");
// for (int k=0;k<npoints;++k) { // for (int k=0;k<npoints;++k) {
// printf(" . ");
// // Set point // // Set point
// elec_coord[0] = point_x[i]; // elec_coord[0] = point_x[i];
// elec_coord[1] = point_y[j]; // elec_coord[1] = point_y[j];

38
org/qmckl_nucleus.org
View File

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

40
org/qmckl_trexio.org
View File

@ -704,6 +704,42 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
} }
#+end_src #+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) #+begin_src c :tangle (eval c)
@ -711,7 +747,6 @@ qmckl_trexio_read_ao_X(qmckl_context context, trexio_t* const file)
} }
#endif #endif
#+end_src #+end_src
** Molecular orbitals ** Molecular orbitals
In this section we read the MO coefficients. 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) #+begin_src c :tangle (eval c_test)
#ifdef HAVE_TREXIO #ifdef HAVE_TREXIO
#define walk_num 2
qmckl_exit_code rc; qmckl_exit_code rc;
char fname[256]; char fname[256];
@ -892,6 +928,8 @@ char message[256];
strncpy(fname, QMCKL_TEST_DIR,255); strncpy(fname, QMCKL_TEST_DIR,255);
strncat(fname, "/chbrclf", 255); strncat(fname, "/chbrclf", 255);
printf("Test file: %s\n", fname); printf("Test file: %s\n", fname);
rc = qmckl_set_electron_walk_num(context, walk_num);
rc = qmckl_trexio_read(context, fname); rc = qmckl_trexio_read(context, fname);
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {

15
org/table_of_contents
View File

@ -1,17 +1,20 @@
qmckl.org qmckl.org
qmckl_error.org qmckl_ao.org
qmckl_blas.org
qmckl_context.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_memory.org
qmckl_mo.org
qmckl_numprec.org qmckl_numprec.org
qmckl_distance.org qmckl_distance.org
qmckl_nucleus.org qmckl_nucleus.org
qmckl_electron.org
qmckl_ao.org
qmckl_mo.org
qmckl_jastrow.org
qmckl_sherman_morrison_woodbury.org qmckl_sherman_morrison_woodbury.org
qmckl_utils.org qmckl_utils.org
qmckl_blas.org
qmckl_trexio.org qmckl_trexio.org
qmckl_verificarlo.org qmckl_verificarlo.org
qmckl_tests.org qmckl_tests.org