diff --git a/src/ASM/GenAlphaMats.C b/src/ASM/GenAlphaMats.C
deleted file mode 100644
index 46c567f7..00000000
--- a/src/ASM/GenAlphaMats.C
+++ /dev/null
@@ -1,102 +0,0 @@
-// $Id$
-//==============================================================================
-//!
-//! \file GenAlphaMats.C
-//!
-//! \date Jul 04 2013
-//!
-//! \author Knut Morten Okstad / SINTEF
-//!
-//! \brief Representation of the element matrices for a dynamic FEM problem.
-//!
-//==============================================================================
-
-#include "GenAlphaMats.h"
-
-
-GenAlphaMats::GenAlphaMats (double alpha, double a, double b) : NewmarkMats(a,b)
-{
-  beta = 0.25*(1.0-alpha)*(1.0-alpha);
-  gamma = 0.5 - alpha;
-}
-
-
-const Matrix& GenAlphaMats::getNewtonMatrix () const
-{
-  Matrix& N = const_cast<Matrix&>(A.front());
-  double alphaPlus1 = 1.5 - gamma; // = 1.0 + alpha
-
-  N = A[2];
-  N.multiply(alphaPlus1*(1.0 + alpha2*gamma/(beta*h)));
-  if (A.size() > 3)
-    N.add(A[3],alphaPlus1);
-  N.add(A[1],(alphaPlus1*gamma*alpha1 + 1.0/h)/(beta*h));
-
-#if SP_DEBUG > 2
-  std::cout <<"\nElement mass matrix"<< A[1];
-  if (A.size() > 3)
-    std::cout <<"Material stiffness matrix"<< A[2]
-              <<"Geometric stiffness matrix"<< A[3];
-  else
-    std::cout <<"Tangent stiffness matrix"<< A[2];
-  std::cout <<"Resulting Newton matrix"<< A[0];
-#endif
-
-  return A.front();
-}
-
-
-const Vector& GenAlphaMats::getRHSVector () const
-{
-  int ia = vec.size() - 1; // index to element acceleration vector (a)
-  int iv = vec.size() - 2; // index to element velocity vector (v)
-
-#if SP_DEBUG > 2
-  std::cout <<"GenAlphaMats::getRHSVector "
-            << (isPredictor ? "predictor" : "corrector") << std::endl;
-#endif
-
-  if (A.size() > 2 && b.size() > 1 && vec.size() > 2)
-  {
-    Vector& Fi = const_cast<Vector&>(b[1]);
-
-    // Find the actual inertia force from the dynamic equilibrium equation
-    // (Fi = Fext - Fs - Fd) and store it in the second RHS vector Fi=b[1]
-    Fi = b.front(); // Fi = Fext - Fs
-    if (alpha1 > 0.0 && iv > 0)
-      Fi.add(A[1]*vec[iv],-alpha1); // Fi -= alpha1*M*v
-    if (alpha2 > 0.0 && iv > 0)
-      Fi.add(A[2]*vec[iv],-alpha2); // Fi -= alpha2*K*v
-
-#if SP_DEBUG > 2
-    std::cout <<"Element inertia vector"<< Fi;
-#endif
-  }
-
-  if (A.size() > 2 && vec.size() > 2)
-  {
-    Vector& RHS = const_cast<Vector&>(b.front());
-#if SP_DEBUG > 2
-    std::cout <<"\nElement velocity vector"<< vec[iv];
-    std::cout <<"Element acceleration vector"<< vec[ia];
-    std::cout <<"S_ext - S_int"<< b.front();
-    std::cout <<"S_i"<< A[1]*vec[ia];
-#endif
-
-    // Find the right-hand-side vector
-    double alphaPlus1 = 1.5 - gamma;
-    RHS *= alphaPlus1;
-    if (alpha1 > 0.0)
-      RHS.add(A[1]*vec[iv], alpha1*(isPredictor ? alphaPlus1 : -alphaPlus1));
-    if (alpha2 > 0.0)
-      RHS.add(A[2]*vec[iv], alpha2*(isPredictor ? alphaPlus1 : -alphaPlus1));
-
-    RHS.add(A[1]*vec[ia], isPredictor ? 1.0 : -1.0);
-  }
-
-#if SP_DEBUG > 2
-  std::cout <<"\nElement right-hand-side vector"<< b.front();
-#endif
-
-  return b.front();
-}
diff --git a/src/ASM/GenAlphaMats.h b/src/ASM/GenAlphaMats.h
deleted file mode 100644
index b662342d..00000000
--- a/src/ASM/GenAlphaMats.h
+++ /dev/null
@@ -1,39 +0,0 @@
-// $Id$
-//==============================================================================
-//!
-//! \file GenAlphaMats.h
-//!
-//! \date Jul 04 2013
-//!
-//! \author Knut Morten Okstad / SINTEF
-//!
-//! \brief Representation of the element matrices for a dynamic FEM problem.
-//!
-//==============================================================================
-
-#ifndef _GEN_ALPHA_MATS_H
-#define _GEN_ALPHA_MATS_H
-
-#include "NewmarkMats.h"
-
-
-/*!
-  \brief Class representing the element matrices for a dynamic FEM problem
-  based on generalized alpha time integration.
-*/
-
-class GenAlphaMats : public NewmarkMats
-{
-public:
-  //! \brief The constructor initializes the time integration parameters.
-  GenAlphaMats(double alpha, double a, double b);
-  //! \brief Empty destructor.
-  virtual ~GenAlphaMats() {}
-
-  //! \brief Returns the element-level Newton matrix.
-  virtual const Matrix& getNewtonMatrix() const;
-  //! \brief Returns the element-level right-hand-side vector.
-  virtual const Vector& getRHSVector() const;
-};
-
-#endif
diff --git a/src/ASM/HHTMats.C b/src/ASM/HHTMats.C
index f1ae1f93..abc8d47b 100644
--- a/src/ASM/HHTMats.C
+++ b/src/ASM/HHTMats.C
@@ -14,10 +14,11 @@
 #include "HHTMats.h"
 
 
-HHTMats::HHTMats (double alpha, double a, double b) : NewmarkMats(a,b)
+HHTMats::HHTMats (double alpha, double a, double b, bool old) : NewmarkMats(a,b)
 {
   beta = 0.25*(1.0-alpha)*(1.0-alpha);
   gamma = 0.5 - alpha;
+  oldHHT = old;
 }
 
 
@@ -37,7 +38,7 @@ const Matrix& HHTMats::getNewtonMatrix () const
   N.add(A[1],(alphaPlus1*alpha1*gamma + 1.0/h)/(beta*h));
 
 #if SP_DEBUG > 2
-  std::cout <<"\nDynElmMats::getNewtonMatrix";
+  std::cout <<"\nHHTMats::getNewtonMatrix";
   std::cout <<"\nElement mass matrix"<< A[1];
   if (A.size() > 3)
     std::cout <<"Material stiffness matrix"<< A[2]
@@ -63,7 +64,7 @@ const Vector& HHTMats::getRHSVector () const
   int ia = vec.size() - 1; // index to element acceleration vector (a)
   int iv = vec.size() - 2; // index to element velocity vector (v)
 
-  if (isPredictor && b.size() > 1)
+  if ((isPredictor || oldHHT) && b.size() > 1)
   {
     Vector& Fia = const_cast<Vector&>(b[1]);
 
@@ -76,11 +77,12 @@ const Vector& HHTMats::getRHSVector () const
       Fia.add(A[1]*vec[iv],-alpha1); // Fia -= alpha1*M*v
     if (alpha2 > 0.0 && iv > 0)
       Fia.add(A[2]*vec[iv],-alpha2); // Fia -= alpha2*K*v
+#if SP_DEBUG > 2
+    std::cout <<"Element inertia vector"<< Fia;
+#endif
   }
 
 #if SP_DEBUG > 2
-  if (isPredictor && b.size() > 1)
-    std::cout <<"Element inertia vector"<< b[1];
   std::cout <<"Element velocity vector"<< vec[iv];
   std::cout <<"Element acceleration vector"<< vec[ia];
   if (b.size() > 2)
@@ -97,7 +99,9 @@ const Vector& HHTMats::getRHSVector () const
   // the predictor step force vector after the element assembly.
   Vector& RHS = const_cast<Vector&>(b.front());
   double alphaPlus1 = 1.5 - gamma;
-  if (isPredictor)
+  if (oldHHT)
+    RHS *= alphaPlus1; // RHS = -(1+alphaH)*S_int
+  else if (isPredictor)
   {
     int pa = iv - 1; // index to predicted element acceleration vector (A)
     A[1].multiply(vec[pa]-vec[ia],RHS); // RHS = M*(A-a)
@@ -121,9 +125,10 @@ const Vector& HHTMats::getRHSVector () const
     RHS.add(A[1]*vec[iv],alphaPlus1*alpha1); // RHS -= (1+alphaH)*alpha1*M*v
   if (alpha2 > 0.0)
     RHS.add(A[2]*vec[iv],alphaPlus1*alpha2); // RHS -= (1+alphaH)*alpha2*K*v
-
+  if (oldHHT)
+    RHS.add(A[1]*vec[ia], isPredictor ? 1.0 : -1.0); // RHS (+/-)= M*a
 #if SP_DEBUG > 2
-  std::cout <<"Element right-hand-side vector"<< b.front();
+  std::cout <<"\nElement right-hand-side vector"<< b.front();
 #endif
 
   return b.front();
diff --git a/src/ASM/HHTMats.h b/src/ASM/HHTMats.h
index 0958aa50..2e923f04 100644
--- a/src/ASM/HHTMats.h
+++ b/src/ASM/HHTMats.h
@@ -26,7 +26,7 @@ class HHTMats : public NewmarkMats
 {
 public:
   //! \brief The constructor initializes the time integration parameters.
-  HHTMats(double alpha, double a, double b);
+  HHTMats(double alpha, double a, double b, bool old = false);
   //! \brief Empty destructor.
   virtual ~HHTMats() {}
 
@@ -34,6 +34,9 @@ public:
   virtual const Matrix& getNewtonMatrix() const;
   //! \brief Returns the element-level right-hand-side vector.
   virtual const Vector& getRHSVector() const;
+
+private:
+  bool oldHHT; //!< If \e true, used toghether with NewmarkNLSIM
 };
 
 #endif