diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org
index 88fe69c..5d49030 100644
--- a/org/qmckl_mo.org
+++ b/org/qmckl_mo.org
@@ -666,25 +666,41 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
   integer*8 :: i,j,k
   double precision :: c1, c2, c3, c4, c5
 
-  do j=1,point_num
-     mo_vgl(:,:,j) = 0.d0
-     do k=1,ao_num
-        if (ao_vgl(k,1,j) /= 0.d0) then
-           c1 = ao_vgl(k,1,j)
-           c2 = ao_vgl(k,2,j)
-           c3 = ao_vgl(k,3,j)
-           c4 = ao_vgl(k,4,j)
-           c5 = ao_vgl(k,5,j)
-           do i=1,mo_num
-              mo_vgl(i,1,j) = mo_vgl(i,1,j) + coef_normalized_t(i,k) * c1
-              mo_vgl(i,2,j) = mo_vgl(i,2,j) + coef_normalized_t(i,k) * c2
-              mo_vgl(i,3,j) = mo_vgl(i,3,j) + coef_normalized_t(i,k) * c3
-              mo_vgl(i,4,j) = mo_vgl(i,4,j) + coef_normalized_t(i,k) * c4
-              mo_vgl(i,5,j) = mo_vgl(i,5,j) + coef_normalized_t(i,k) * c5
-           end do
-        end if
+  integer*8 :: LDA, LDB, LDC
+
+  info = QMCKL_SUCCESS
+  if (.False.)  then    ! fast algorithm
+     do j=1,point_num
+        mo_vgl(:,:,j) = 0.d0
+        do k=1,ao_num
+           if (ao_vgl(k,1,j) /= 0.d0) then
+              c1 = ao_vgl(k,1,j)
+              c2 = ao_vgl(k,2,j)
+              c3 = ao_vgl(k,3,j)
+              c4 = ao_vgl(k,4,j)
+              c5 = ao_vgl(k,5,j)
+              do i=1,mo_num
+                 mo_vgl(i,1,j) = mo_vgl(i,1,j) + coef_normalized_t(i,k) * c1
+                 mo_vgl(i,2,j) = mo_vgl(i,2,j) + coef_normalized_t(i,k) * c2
+                 mo_vgl(i,3,j) = mo_vgl(i,3,j) + coef_normalized_t(i,k) * c3
+                 mo_vgl(i,4,j) = mo_vgl(i,4,j) + coef_normalized_t(i,k) * c4
+                 mo_vgl(i,5,j) = mo_vgl(i,5,j) + coef_normalized_t(i,k) * c5
+              end do
+           end if
+        end do
      end do
-  end do
+     
+  else ! dgemm
+
+    LDA = size(coef_normalized_t,1)
+    LDB = size(ao_vgl,1) 
+    LDC = size(mo_vgl,1)
+
+    info = qmckl_dgemm(context,'N', 'N', mo_num, point_num*5_8, ao_num*1_8, 1.d0,     &
+                                    coef_normalized_t, LDA, ao_vgl, LDB, &
+                                    0.d0, mo_vgl, LDC)
+
+  end if
 
 end function qmckl_compute_mo_basis_mo_vgl_doc_f
     #+end_src