diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml
index c43111d..0bc46ed 100644
--- a/SCF/HartreeFock.ml
+++ b/SCF/HartreeFock.ml
@@ -319,17 +319,27 @@ let make
       xt_o_x  (Mat.sub  fps  spf) m_X
     in
 
+    let diis, m_F_diis = 
+      let diis =
+        DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
+      in
+
+      try
+        let m_F_diis =
+          let x = 
+            Bigarray.genarray_of_array1 (DIIS.next diis)
+          in
+          Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
+        in
+        diis, m_F_diis
+
+      with Failure _ -> (* Failure in DIIS.next *)
+        DIIS.make (), m_F_ortho
+    in
     let diis = 
       DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
     in
 
-    let m_F_diis =
-      let x = 
-        Bigarray.genarray_of_array1 (DIIS.next diis)
-      in
-      Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
-    in
-
 
     (* MOs in orthogonal MO basis *)
     let m_C', _ =
@@ -515,15 +525,22 @@ let make
       xt_o_x  (Mat.sub  fps  spf) m_X
     in
 
-    let diis = 
-      DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
-    in
-
-    let m_F_diis =
-      let x = 
-        Bigarray.genarray_of_array1 (DIIS.next diis)
+    let diis, m_F_diis = 
+      let diis =
+        DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
       in
-      Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
+
+      try
+        let m_F_diis =
+          let x = 
+            Bigarray.genarray_of_array1 (DIIS.next diis)
+          in
+          Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
+        in
+        diis, m_F_diis
+
+      with Failure _ -> (* Failure in DIIS.next *)
+        DIIS.make (), m_F_ortho
     in
 
 
diff --git a/Utils/DIIS.mli b/Utils/DIIS.mli
index 9f1ce66..e658508 100644
--- a/Utils/DIIS.mli
+++ b/Utils/DIIS.mli
@@ -9,7 +9,7 @@ The DIIS approximate solution for iteration {% $m+1$ %} is given by
 
 {% \begin{align*}
 \mathbf{p}_{m+1} & = \sum_{i=1}^m c_i (\mathbf{p}^f + \mathbf{e}_i) \\
-                         & = \sum_{i=1}^m c_i \mathbf{p}^f + \sum_i c_i \mathbf{e}_i 
+                 & = \sum_{i=1}^m c_i \mathbf{p}^f + \sum_i c_i \mathbf{e}_i 
 \end{align*} %}
 
 where {% $\mathbf{p}^f$ %} is the exact solution. One wants to minimize the