From 1869756ea48fc5dc8f82311b95e311104be29b82 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Tue, 26 Oct 2021 19:13:20 +0200 Subject: [PATCH] Verified Local energy with QMC=Chem for Be2. #41 --- org/qmckl_blas.org | 284 +++++++++++++++++++++++++++++++++---- org/qmckl_determinant.org | 16 +-- org/qmckl_electron.org | 2 +- org/qmckl_local_energy.org | 95 +++++++++---- org/qmckl_mo.org | 39 +++-- 5 files changed, 357 insertions(+), 79 deletions(-) diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index cd04b2e..918cd9b 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -22,7 +22,7 @@ int main() { * Matrix operations ** ~qmckl_dgemm~ - + Matrix multiply: $C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function. TODO: Add description about the external library dependence. @@ -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)) - do i = 1, m - do j = 1, k + allocate(AT(m,k)) + do i = 1, k + do j = 1, m AT(j,i) = A(i,j) end do end do @@ -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)) - do i = 1, k - do j = 1, n + allocate(BT(k,n)) + do i = 1, n + do j = 1, k BT(j,i) = B(i,j) end do end do @@ -161,26 +162,75 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, endif if (TransA) then - - if (alpha == 1.0d0 .and. beta == 0.0d0) then - C = matmul(AT,B) - else - C = beta*C + alpha*matmul(AT,B) - endif + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, B, LDB_2, beta, c, LDC) else if (TransB) then - if (alpha == 1.0d0 .and. beta == 0.0d0) then - C = matmul(A,BT) - else - C = beta*C + alpha*matmul(A,BT) - endif + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, BT, LDB_2, beta, c, LDC) + else if (TransA .and. TransB) then + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, BT, LDB_2, beta, c, LDC) else - if (alpha == 1.0d0 .and. beta == 0.0d0) then - C = matmul(A,B) - else - C = beta*C + alpha*matmul(A,B) - endif + info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC) endif end function qmckl_dgemm_f + +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 diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org index 6beba36..805cf34 100644 --- a/org/qmckl_determinant.org +++ b/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 diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index a26276c..931f259 100644 --- a/org/qmckl_electron.org +++ b/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 diff --git a/org/qmckl_local_energy.org b/org/qmckl_local_energy.org index fa350d9..b380719 100644 --- a/org/qmckl_local_energy.org +++ b/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_adj_matrix_beta, + ctx->det.det_value_alpha, + ctx->det.det_value_beta, + ctx->det.det_inv_matrix_alpha, + ctx->det.det_inv_matrix_beta, ctx->local_energy.e_kin); } else { return qmckl_failwith( context, @@ -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_adj_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det | + | ~double~ | ~det_value_alpha[det_num_alpha][walk_num]~ | in | Det of wavefunction | + | ~double~ | ~det_value_beta[det_num_beta][walk_num]~ | in | Det of wavefunction | + | ~double~ | ~det_inv_matrix_alpha[det_num_alpha][walk_num][alpha_num][alpha_num]~ | in | Value, gradients and Laplacian of the Det | + | ~double~ | ~det_inv_matrix_beta[det_num_beta][walk_num][beta_num][beta_num]~ | in | Value, gradients and Laplacian of the Det | | ~double~ | ~e_kin[walk_num]~ | out | Kinetic energy | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_kinetic_energy_f(context, walk_num, & det_num_alpha, det_num_beta, alpha_num, beta_num, elec_num, mo_index_alpha, mo_index_beta, & - mo_num, mo_vgl, det_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_adj_matrix_beta(beta_num, beta_num, walk_num, det_num_beta) + double precision, intent(in) :: det_value_alpha(walk_num, det_num_alpha) + double precision, intent(in) :: det_value_beta(walk_num, det_num_beta) + double precision, intent(in) :: det_inv_matrix_alpha(alpha_num, alpha_num, walk_num, det_num_alpha) + double precision, intent(in) :: det_inv_matrix_beta(beta_num, beta_num, walk_num, det_num_beta) double precision, intent(inout) :: e_kin(walk_num) + double precision :: tmp_e integer*8 :: idet, iwalk, ielec, mo_id, imo info = QMCKL_SUCCESS @@ -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) - e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_adj_matrix_alpha(imo, ielec, iwalk, idet) * & + mo_id = mo_index_alpha(imo, iwalk, idet) + e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * & mo_vgl(mo_id, ielec, 5) + !print *,"det alpha = ",det_inv_matrix_alpha(imo,ielec,iwalk,idet) + !print *,mo_vgl(mo_id,ielec,5) + !!print *," det val = ",det_value_alpha(iwalk,idet) + !tmp_e = tmp_e - 0.5d0 * det_inv_matrix_alpha(imo, ielec, iwalk, idet) * & + ! mo_vgl(mo_id, ielec, 5) end do + !print *,"e_kin = ",tmp_e end do ! Beta part + tmp_e = 0.0d0 do imo = 1, beta_num do ielec = 1, beta_num - mo_id = mo_index_beta(ielec, iwalk, idet) - e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_adj_matrix_beta(imo, ielec, iwalk, idet) * & - mo_vgl(mo_id, ielec, 5) + mo_id = mo_index_beta(imo, iwalk, idet) + e_kin(iwalk) = e_kin(iwalk) - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * & + mo_vgl(mo_id, alpha_num + ielec, 5) + !print *,"det beta = ",det_inv_matrix_beta(imo,ielec,iwalk,idet) + !print *,mo_vgl(mo_id,alpha_num+ielec,5) + !!print *," det val = ",det_value_alpha(iwalk,idet) + !tmp_e = tmp_e - 0.5d0 * det_inv_matrix_beta(imo, ielec, iwalk, idet) * & + ! mo_vgl(mo_id, alpha_num + ielec, 5) end do + !print *,"e_kin = ",tmp_e end do end do end do @@ -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_adj_matrix_beta, + const double* det_value_alpha, + const double* det_value_beta, + const double* det_inv_matrix_alpha, + const double* det_inv_matrix_beta, double* const e_kin ); #+end_src @@ -474,8 +512,10 @@ end function qmckl_compute_kinetic_energy_f mo_index_beta, & mo_num, & mo_vgl, & - det_adj_matrix_alpha, & - det_adj_matrix_beta, & + det_value_alpha, & + 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_adj_matrix_beta(beta_num,beta_num,walk_num,det_num_beta) + real (c_double ) , intent(in) :: det_value_alpha(walk_num,det_num_alpha) + real (c_double ) , intent(in) :: det_value_beta(walk_num,det_num_beta) + real (c_double ) , intent(in) :: det_inv_matrix_alpha(alpha_num,alpha_num,walk_num,det_num_alpha) + real (c_double ) , intent(in) :: det_inv_matrix_beta(beta_num,beta_num,walk_num,det_num_beta) real (c_double ) , intent(out) :: e_kin(walk_num) integer(c_int32_t), external :: qmckl_compute_kinetic_energy_f @@ -510,8 +552,10 @@ end function qmckl_compute_kinetic_energy_f mo_index_beta, & mo_num, & mo_vgl, & - det_adj_matrix_alpha, & - det_adj_matrix_beta, & + det_value_alpha, & + 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 diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index b9d220c..60f39f4 100644 --- a/org/qmckl_mo.org +++ b/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 - LDB = K - LDC = M + LDA = size(coef_normalized,1) + idx = 0 + !do j = 1,ao_num + !do i = 1,mo_num + ! idx = idx + 1 + ! coef_all(idx) = coef_normalized(i,j) + !end do + !end do + !idx = 0 + !do j = 1,mo_num + !do i = 1,ao_num + ! idx = idx + 1 + ! coef_trans(j,i) = coef_all(idx) + !end do + !end do ao_vgl_big = reshape(ao_vgl(:, :, :),(/ao_num, elec_num*5_8/)) + LDB = size(ao_vgl_big,1) + LDC = size(mo_vgl_big,1) info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - coef_normalized,size(coef_normalized,1) * 1_8, & - ao_vgl_big, LDB, & + coef_normalized,size(coef_normalized,1)*1_8, & + ao_vgl_big, size(ao_vgl_big,1)*1_8, & beta, & mo_vgl_big,LDC) mo_vgl = reshape(mo_vgl_big,(/mo_num,elec_num,5_8/)) @@ -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)