diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 52c829f..8171d5b 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -306,6 +306,9 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v rc = qmckl_provide_ao_vgl(context); if (rc != QMCKL_SUCCESS) return rc; + rc = qmckl_provide_mo_vgl(context); + if (rc != QMCKL_SUCCESS) return rc; + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); @@ -333,11 +336,11 @@ qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_v *** Provide #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_provide_mo_basis_vgl(qmckl_context context); +qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_mo_basis_vgl(qmckl_context context) +qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { @@ -439,7 +442,7 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & double precision , intent(in) :: ao_vgl(ao_num,elec_num,walk_num,5) double precision , intent(in) :: coef_normalized(mo_num,ao_num) double precision , intent(out) :: mo_vgl(mo_num,elec_num,walk_num,5) - logical :: TransA, TransB + logical*8 :: TransA, TransB double precision :: alpha, beta integer :: info_qmckl_dgemm_value integer :: info_qmckl_dgemm_Gx @@ -456,46 +459,61 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & beta = 0.0d0 info = QMCKL_SUCCESS + info_qmckl_dgemm_value = QMCKL_SUCCESS + info_qmckl_dgemm_Gx = QMCKL_SUCCESS + info_qmckl_dgemm_Gy = QMCKL_SUCCESS + info_qmckl_dgemm_Gz = QMCKL_SUCCESS + info_qmckl_dgemm_lap = QMCKL_SUCCESS ! Don't compute exponentials when the result will be almost zero. ! TODO : Use numerical precision here cutoff = -dlog(1.d-15) - M = mo_num * 1_8 - N = ao_num * 1_8 - K = 1_8 + M = 1_8 + N = 1_8 + K = mo_num * 1_8 do iwalk = 1, walk_num do ielec = 1, elec_num ! Value + TransA = .True. + TransB = .True. info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 1), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 1), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,1),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,1),1_8) ! Grad_x + TransA = .True. + TransB = .True. info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 2), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 2), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,2),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,2),1_8) ! Grad_y + TransA = .True. + TransB = .True. info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 3), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 3), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,3),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,3),1_8) ! Grad_z + TransA = .True. + TransB = .True. info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 4), size(ao_vgl,1) * 1_8, & + ao_vgl(:, ielec, iwalk, 4), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,4),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,4),1_8) ! Lapl_z - info_qmckl_dgemm_lap = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, :, iwalk, 5), size(ao_vgl,1) * 1_8, & + TransA = .True. + TransB = .True. + info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 5), size(ao_vgl,1) * 1_8, & coef_normalized,size(coef_normalized,1) * 1_8, & beta, & - mo_vgl(:,:,iwalk,5),size(mo_vgl,1) * 1_8) + mo_vgl(:,ielec,iwalk,5),1_8) end do end do @@ -657,9 +675,6 @@ assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); -//const int64_t shell_num = chbrclf_shell_num; -//const int64_t prim_num = chbrclf_prim_num; -//const int64_t ao_num = chbrclf_ao_num; const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]); const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]); const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]); @@ -736,9 +751,30 @@ double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num]; rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0])); assert (rc == QMCKL_SUCCESS); -rc = qmckl_get_mo_basis_vgl(context, &(ao_vgl[0][0][0][0])); +/* Set up MO data */ +rc = qmckl_set_mo_basis_type(context, typ); assert (rc == QMCKL_SUCCESS); +const int64_t mo_num = chbrclf_ao_num; +rc = qmckl_set_mo_basis_mo_num(context, mo_num); +assert (rc == QMCKL_SUCCESS); + +double mo_coefficient[chbrclf_ao_num][mo_num]; +const int64_t min_num = (((chbrclf_ao_num) < (mo_num)) ? (chbrclf_ao_num) : (mo_num)); +for (int i=0; i