introduced doc version of qmckl_compute_forces_mo_g

Author: scemama

org/qmckl_ao.org

@@ -8229,7 +8229,7 @@ function qmckl_compute_ao_hessian_doc(context, &
      e_coord(1) = coord(ipoint,1)
      e_coord(2) = coord(ipoint,2)
      e_coord(3) = coord(ipoint,3)
-     !o_hessian(:,:,ipoint,:,:) = 0.d0
+     ao_hessian(:,:,ipoint,:,:) = 0.d0
      do inucl=1,nucl_num
         n_coord(1) = nucl_coord(inucl,1)
         n_coord(2) = nucl_coord(inucl,2)

org/qmckl_forces.org

@@ -6452,6 +6452,11 @@ qmckl_exit_code
 qmckl_get_forces_mo_g(qmckl_context context,
                           double* const forces_mo_g,
                           const int64_t size_max);
+
+qmckl_exit_code
+qmckl_get_forces_mo_g_inplace(qmckl_context context,
+                              double* const forces_mo_g,
+                              const int64_t size_max);
     #+end_src
 
     #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
@@ -6462,7 +6467,10 @@ qmckl_get_forces_mo_g(qmckl_context context,
 {
 
   if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
-    return QMCKL_NULL_CONTEXT;
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_get_forces_mo_g",
+                           NULL);
   }
 
   qmckl_exit_code rc;
@@ -6483,6 +6491,49 @@ qmckl_get_forces_mo_g(qmckl_context context,
 
   memcpy(forces_mo_g, ctx->forces.forces_mo_g, sze * sizeof(double));
 
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_forces_mo_g_inplace (qmckl_context context,
+                                     double* const forces_mo_g,
+                                     const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_get_forces_mo_g",
+                           NULL);
+  }
+
+  qmckl_exit_code rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  const int64_t sze = ctx->point.num * 3 * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_forces_mo_g",
+                           "input array too small");
+  }
+
+  rc = qmckl_context_touch(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  double* old_array =  ctx->forces.forces_mo_g;
+
+  ctx->forces.forces_mo_g = forces_mo_g;
+
+  rc = qmckl_provide_forces_mo_g(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  ctx->forces.forces_mo_g = old_array;
+
   return QMCKL_SUCCESS;
 }
     #+end_src
@@ -6500,6 +6551,19 @@ qmckl_get_forces_mo_g(qmckl_context context,
        integer (c_int64_t) , intent(in)  , value :: size_max
      end function qmckl_get_forces_mo_g
   end interface
+
+  interface
+     integer(qmckl_exit_code) function qmckl_get_forces_mo_g_inplace (context, &
+          forces_mo_g, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+
+       integer (c_int64_t) , intent(in)  , value :: context
+       real(c_double),       intent(out)         :: forces_mo_g(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_forces_mo_g_inplace
+  end interface
     #+end_src
 
 *** Provide :noexport:
@@ -6527,7 +6591,7 @@ qmckl_exit_code qmckl_provide_forces_mo_g(qmckl_context context)
   if (!ctx->mo_basis.provided) {
     return qmckl_failwith( context,
                            QMCKL_NOT_PROVIDED,
-                           "qmckl_provide_mo_basis_mo_vgl",
+                           "qmckl_provide_forces_mo_g",
                            NULL);
   }
 
@@ -6573,18 +6637,18 @@ qmckl_exit_code qmckl_provide_forces_mo_g(qmckl_context context)
                              NULL);
     }
 
-    rc = qmckl_compute_forces_mo_g_doc(context,
-                                         ctx->ao_basis.ao_num,
-                                         ctx->mo_basis.mo_num,
-                                         ctx->point.num,
-                                         ctx->nucleus.num,
-                                         ctx->ao_basis.shell_num,
-                                         ctx->ao_basis.nucleus_index,
-                                         ctx->ao_basis.nucleus_shell_num,
-                                         ctx->ao_basis.shell_ang_mom,
-                                         ctx->mo_basis.coefficient_t,
-                                         ctx->ao_basis.ao_hessian,
-                                         ctx->forces.forces_mo_g);
+    rc = qmckl_compute_forces_mo_g(context,
+                                   ctx->ao_basis.ao_num,
+                                   ctx->mo_basis.mo_num,
+                                   ctx->point.num,
+                                   ctx->nucleus.num,
+                                   ctx->ao_basis.shell_num,
+                                   ctx->ao_basis.nucleus_index,
+                                   ctx->ao_basis.nucleus_shell_num,
+                                   ctx->ao_basis.shell_ang_mom,
+                                   ctx->mo_basis.coefficient_t,
+                                   ctx->ao_basis.ao_hessian,
+                                   ctx->forces.forces_mo_g);
 
     if (rc != QMCKL_SUCCESS) {
       return rc;
@@ -6641,6 +6705,34 @@ integer function qmckl_compute_forces_mo_g_doc(context, &
   integer*8 :: i,j, m, n,a
   double precision :: c1
 
+
+  info = qmckl_dgemm(context,'N', 'N', mo_num, 3*point_num*3*nucl_num, ao_num, &
+                -1.d0, coefficient_t, mo_num, &
+                ao_hessian(1:ao_num, 1:3, 1:point_num, 1:3, 1:nucl_num), ao_num, &
+                0.d0, forces_mo_g(1:mo_num, 1:3, 1:point_num, 1:3, 1:nucl_num), mo_num)
+
+end function qmckl_compute_forces_mo_g_doc
+    #+end_src
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_forces_mo_g_hpc(context, &
+     ao_num, mo_num, point_num, nucl_num, &
+     shell_num, nucleus_index, nucleus_shell_num, shell_ang_mom, &
+     coefficient_t, ao_hessian, forces_mo_g) &
+     bind(C) result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in), value  :: context
+  integer (c_int64_t)   , intent(in), value  :: nucl_num, ao_num, mo_num, point_num, shell_num
+  integer (c_int64_t) ,   intent(in)         :: nucleus_index(nucl_num)
+  integer (c_int64_t) ,   intent(in)         :: nucleus_shell_num(nucl_num)
+  integer (c_int32_t) ,   intent(in)         :: shell_ang_mom(shell_num)
+  real    (c_double )   , intent(in)         :: coefficient_t(mo_num,ao_num)
+  real    (c_double )   , intent(in)         :: ao_hessian(ao_num,4,point_num,3,nucl_num)
+  real    (c_double )   , intent(out)        :: forces_mo_g(mo_num,3,point_num,3,nucl_num)
+  integer*8 :: i,j, m, n,a
+  double precision :: c1
+
   integer           :: l, il, k
   integer*8         :: ishell_start, ishell_end, ishell
   integer           :: lstart(0:20)
@@ -6673,18 +6765,14 @@ integer function qmckl_compute_forces_mo_g_doc(context, &
       do j=1,point_num
         do m = 1, 3
           do ishell = ishell_start, ishell_end
-            k = ao_index(ishell)
             l = shell_ang_mom(ishell)
-            do il = lstart(l), lstart(l+1)-1
+            il = lstart(l+1)-lstart(l) 
+            do k = ao_index(ishell), ao_index(ishell) + il - 1
               c1 = ao_hessian(k, m, j, n, a)
-              if (c1 == 0.d0) then
-                k = k + 1
-                cycle
-              end if
+              if (c1 == 0.d0) cycle
               do i=1,mo_num
                 forces_mo_g(i, m, j, n, a) = forces_mo_g(i, m, j, n, a) - coefficient_t(i,k) * c1
               end do
-              k = k+1
             end do
           end do
         end do
@@ -6695,26 +6783,39 @@ integer function qmckl_compute_forces_mo_g_doc(context, &
   deallocate(ao_index)
 
 
-
-
-  !do a=1,nucl_num
-  !  do n = 1,3
-  !    do m = 1, 3
-  !      info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, &
-  !              -1.d0, coefficient_t, mo_num, &
-  !              ao_hessian(:, m, :, n, a), ao_num, &
-  !              1.d0, forces_mo_g(:, m, :, n, a), mo_num)
-  !    end do
-  !  end do
-  !end do
-
-
-end function qmckl_compute_forces_mo_g_doc
+end function qmckl_compute_forces_mo_g_hpc
     #+end_src
 
-   #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
-   qmckl_exit_code qmckl_compute_forces_mo_g_doc (
+qmckl_exit_code qmckl_compute_forces_mo_g (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const int64_t nucl_num,
+         const int64_t shell_num,
+         const int64_t* nucleus_index,
+         const int64_t* nucleus_shell_num,
+         const int32_t* shell_ang_mom,
+         const double* coefficient_t,
+         const double* ao_hessian,
+         double* const forces_mo_g );
+
+qmckl_exit_code qmckl_compute_forces_mo_g_doc (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const int64_t nucl_num,
+         const int64_t shell_num,
+         const int64_t* nucleus_index,
+         const int64_t* nucleus_shell_num,
+         const int32_t* shell_ang_mom,
+         const double* coefficient_t,
+         const double* ao_hessian,
+         double* const forces_mo_g );
+
+qmckl_exit_code qmckl_compute_forces_mo_g_hpc (
          const qmckl_context context,
          const int64_t ao_num,
          const int64_t mo_num,
@@ -6728,6 +6829,31 @@ end function qmckl_compute_forces_mo_g_doc
          const double* ao_hessian,
          double* const forces_mo_g );
    #+end_src
+   
+    #+begin_src c :tangle (eval c) :comments org :exports none
+qmckl_exit_code qmckl_compute_forces_mo_g (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const int64_t nucl_num,
+         const int64_t shell_num,
+         const int64_t* nucleus_index,
+         const int64_t* nucleus_shell_num,
+         const int32_t* shell_ang_mom,
+         const double* coefficient_t,
+         const double* ao_hessian,
+         double* const forces_mo_g )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_forces_mo_g_hpc 
+#else
+  return qmckl_compute_forces_mo_g_doc 
+#endif
+   (context, ao_num, mo_num, point_num, nucl_num, shell_num, nucleus_index,
+   nucleus_shell_num, shell_ang_mom, coefficient_t, ao_hessian, forces_mo_g );
+}
+#+end_src
 
 *** Test :noexport:
 
@@ -6795,7 +6921,8 @@ for (int64_t a = 0; a < nucl_num; a++) {
             if (m == -4) {
               finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] = 0.0;
             }
-            finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] += coef[m + 4] * mo_output[i*mo_num*5 + (n+1)*mo_num + j]/delta_x;
+            finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] +=
+              coef[m + 4] * mo_output[i*mo_num*5 + (n+1)*mo_num + j]/delta_x;
           }
         }
       }
@@ -6817,15 +6944,19 @@ free(mo_output);
 
 
 
-for (int j = 0; j < mo_num; j++){
+for (int a = 0; a < nucl_num; a++) {
   for (int n = 0; n < 3; n++){
-    for (int i = 0; i < point_num; i++){
-      for (int a = 0; a < nucl_num; a++) {
-        for (int k = 0; k < 3; k++){
-          //printf("k=%i a=%i i=%i n=%i j=%i\n", k, a, i, n, j);
-          //printf("%.10f\t", finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j]);
-          //printf("%.10f\n", forces_mo_g[a*9*mo_num*point_num + k*3*mo_num*point_num + i*3*mo_num + n*mo_num + j]);
-          assert(fabs(finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] -  forces_mo_g[a*9*mo_num*point_num + k*3*mo_num*point_num + i*3*mo_num + n*mo_num + j]) < 1.e-9);
+    for (int j = 0; j < point_num; j++){
+      for (int m = 0; m < 3; m++){
+        for (int i = 0; i < mo_num; i++){
+          printf("a=%i n=%i j=%i m=%i i=%i\n", a, n, j, m, i);
+          printf("%.10f\t", finite_difference_force_mo_g[n*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + j*3*mo_num + m*mo_num + i]);
+          printf("%.10f\n", forces_mo_g[a*9*mo_num*point_num + n*3*mo_num*point_num + j*3*mo_num + m*mo_num + i]);
+          assert(fabs(finite_difference_force_mo_g[n*3*mo_num*point_num*nucl_num +
+                                                   a*3*mo_num*point_num + j*3*mo_num +
+                                                   m*mo_num + i] - 
+                      forces_mo_g[a*9*mo_num*point_num + n*3*mo_num*point_num +
+                                  j*3*mo_num + m*mo_num + i]) < 1.e-9);
         }
       }
     }