diff --git a/_setup/parameters.py b/_setup/parameters.py
index 5d607c73f52d5c1c1f9a0462147a18325a0eab86..f25dd477e204938188d0ecd7462f9ced5cc36942 100644
--- a/_setup/parameters.py
+++ b/_setup/parameters.py
@@ -19,7 +19,7 @@ class Parameters:
   '''
   
 
-  def __init__(self, A=None, Ms=None, K=0.0, D=0.0, dmi_axis='z', gamma0=2.212761569e+05, mu0=1.2566370614359173e-06, alpha=1.0, T_start=None, T_end=None, timeStepSize=None, H_ext=None, easy_axis=(1,0,0), g=2.0023193043617, muB=9.27400968e-24, e=1.6021766208e-19, P=None, xi=None, dmCoupling="none", maxwellCoupling="strayField", spintronicsCoupling="none", slonczewskiOutOfPlane_axis='z', hbar=1.054571800e-34, Je=None, flag_Je_stationary=False, d=None, p=(1,0,0), TOL=1.0e-8, preconditionerType="jacobi", alphaPreconditioner=None, QmatType="householder", QmatAxis="e3", guessType="old", theta=0.5, m0=None):
+  def __init__(self, A=None, Ms=None, K=0.0, D=0.0, dmi_axis='z', gamma0=2.212761569e+05, mu0=1.2566370614359173e-06, alpha=1.0, T_start=None, T_end=None, timeStepSize=None, H_ext=None, easy_axis=(1,0,0), g=2.0023193043617, muB=9.27400968e-24, e=1.6021766208e-19, P=None, xi=None, dmCoupling="none", dmImplicitTreatment=False, maxwellCoupling="strayField", spintronicsCoupling="none", slonczewskiOutOfPlane_axis='z', hbar=1.054571800e-34, Je=None, flag_Je_stationary=False, d=None, p=(1,0,0), TOL=1.0e-8, preconditionerType="jacobi", alphaPreconditioner=None, QmatType="householder", QmatAxis="e3", guessType="old", theta=0.5, m0=None):
     self.A = A
     self.Ms = Ms
     self.K = K
@@ -39,6 +39,7 @@ class Parameters:
     self.easy_axis = easy_axis
     
     self.dmCoupling = dmCoupling
+    self.dmImplicitTreatment = dmImplicitTreatment
     self.maxwellCoupling = maxwellCoupling
 #    self.sigma_in = 1.0 #TOECDO
 #    self.sigma_out = None #TOECDO
diff --git a/integrators/_forms/__init__.py b/integrators/_forms/__init__.py
index 85b57960b25926d2abd4335178cbeebfd8c1fb28..81234a555ffb38419b65d31696455141fe5120b1 100644
--- a/integrators/_forms/__init__.py
+++ b/integrators/_forms/__init__.py
@@ -6,6 +6,7 @@ from commix.integrators._forms.lf_common import LinearForms_Common
 from commix.integrators._forms.blf_magnetization import BilinearForms_Magnetization
 from commix.integrators._forms.lf_magnetization import LinearForms_Magnetization
 
+from commix.integrators._forms.blf_piImplicit import BilinearForms_piImplicit
 from commix.integrators._forms.lf_derivative import LinearForms_Derivative
 from commix.integrators._forms.lf_old import LinearForms_Old
 
@@ -13,7 +14,7 @@ from commix.integrators._forms.lf_eta import LinearForms_Eta
 
 __all__ = \
   [ 'BilinearForms_Common', 'LinearForms_Common' \
-  , 'BilinearForms_Magnetization', 'LinearForms_Magnetization' \
+  , 'BilinearForms_piImplicit', 'BilinearForms_Magnetization', 'LinearForms_Magnetization' \
   , 'LinearForms_Derivative', 'LinearForms_Old' \
   , 'LinearForms_Eta' \
   ]
diff --git a/integrators/_forms/_bilinearForms/__init__.py b/integrators/_forms/_bilinearForms/__init__.py
index e36e891185b2c2a67d088a4583881fe31398544f..6553ca28f6fa841048adcb06750d9770e08247e4 100644
--- a/integrators/_forms/_bilinearForms/__init__.py
+++ b/integrators/_forms/_bilinearForms/__init__.py
@@ -3,5 +3,7 @@ print("Importing module", __name__)
 from commix.integrators._forms._bilinearForms.mass import Mass
 from commix.integrators._forms._bilinearForms.stiff import Stiff
 from commix.integrators._forms._bilinearForms.skew import Skew
+from commix.integrators._forms._bilinearForms.bulkDMI import BulkDMI
 
-__all__ = ['Mass', 'Stiff', 'Skew']
+__all__ = ['Mass', 'Stiff', 'Skew' \
+  , 'BulkDMI']
diff --git a/integrators/_forms/_bilinearForms/bulkDMI.py b/integrators/_forms/_bilinearForms/bulkDMI.py
new file mode 100644
index 0000000000000000000000000000000000000000..03dc7cd5ef4c97511b17aa5fdf1b0d2e0195af0e
--- /dev/null
+++ b/integrators/_forms/_bilinearForms/bulkDMI.py
@@ -0,0 +1,12 @@
+
+
+################################################################################
+
+
+def BulkDMI(psi, phi, Dpsi, Dphi):
+  from commix._tools.algebra import InnerProd, CrossProd, Curl
+  
+  return InnerProd(Curl(Dpsi), phi) + InnerProd(phi, Curl(Dphi))
+
+
+################################################################################
diff --git a/integrators/_forms/blf_piImplicit.py b/integrators/_forms/blf_piImplicit.py
new file mode 100644
index 0000000000000000000000000000000000000000..cab30c92152d69c7114f9ca2beceba5280fd3b0e
--- /dev/null
+++ b/integrators/_forms/blf_piImplicit.py
@@ -0,0 +1,39 @@
+### collection of bilinear forms, stemming from implicit treatment of pi ###
+
+
+################################################################################
+
+
+class BilinearForms_piImplicit:
+  
+
+#------------------------------------------------------------------------------#
+
+
+  def _DerivLinearTerms_pi_Implicit_LHS(self):
+    return self._theta * self._k \
+      * self._DerivDzyaloshinskiiMoriya_Implicit_LHS_Term()
+
+
+#------------------------------------------------------------------------------#
+
+
+  def _DerivDzyaloshinskiiMoriya_Implicit_LHS_Term(self):
+    if not(self.parameters.dmImplicitTreatment):
+      return self._ZeroCF()
+    else:
+      return -0.5 * self._C_dm * self._DerivDzyaloshinskiiMoriya_Implicit_LHS()
+        
+        
+#------------------------------------------------------------------------------#
+
+
+  def _DerivDzyaloshinskiiMoriya_Implicit_LHS(self):
+    from . import _bilinearForms
+    return _bilinearForms.BulkDMI(self._psi, self._phi, self._Dpsi, self._Dphi)
+
+
+#------------------------------------------------------------------------------#
+
+
+################################################################################
diff --git a/integrators/_forms/lf_derivative.py b/integrators/_forms/lf_derivative.py
index 528308d6cf5f68009fb6f26c67d89063c1719274..11aab4ae6197c5e23b415c3a12aab91f866c50b2 100644
--- a/integrators/_forms/lf_derivative.py
+++ b/integrators/_forms/lf_derivative.py
@@ -17,6 +17,16 @@ class LinearForms_Derivative:
       + self._DerivDzyaloshinskiiMoriyaTerm() )
 
 
+#------------------------------------------------------------------------------#
+
+  
+  def _DerivLinearTerms_pi_Implicit_RHS(self):
+    return self._theta * self._k * \
+      ( self._DerivUniaxialAnisotropyTerm() \
+      + self._DerivStrayFieldTerm() \
+      + self._DerivDzyaloshinskiiMoriyaTerm_Implicit_RHS() )
+
+
 #------------------------------------------------------------------------------#
 
 
@@ -62,6 +72,16 @@ class LinearForms_Derivative:
     return -0.5 * self._C_dm * self._DerivDzyaloshinskiiMoriya()
 
 
+#------------------------------------------------------------------------------#
+
+
+  def _DerivDzyaloshinskiiMoriyaTerm_Implicit_RHS(self):
+    if self.parameters.dmImplicitTreatment:
+      return self._ZeroCF()
+    else:
+      return -0.5 * self._C_dm * self._DerivDzyaloshinskiiMoriya()
+
+
 #------------------------------------------------------------------------------#
 
 
diff --git a/integrators/_integrator.py b/integrators/_integrator.py
index 4d44b56341f04ac01e09560b141461b69e63bd6d..0d6676c1d090bad32ee18813e9d1223b063c2b02 100644
--- a/integrators/_integrator.py
+++ b/integrators/_integrator.py
@@ -102,9 +102,12 @@ class _Integrator( \
     import sys
     ## TODO search for all the prints and sys.exits --> replace with raise(...)
 
-    if(False): #TODO
-      print("Oh no, something went wrong")
-      sys.exit()
+    from . import tps2
+    if self.parameters.dmImplicitTreatment:
+      if not(isinstance(self, tps2.TPS2)):
+        raise ValueError("dmImplicitTreatment=True only available for TPS2")
+      if not(self.parameters.dmCoupling == "bulk"):
+        raise ValueError("dmImplicitTreatment=True only available for bulkDMI")
 
 
 #------------------------------------------------------------------------------#
diff --git a/integrators/tps2.py b/integrators/tps2.py
index b36d1cf5438e17e3e7c7119a0ae01354a673cf0d..a75010efa2fbc203ded8ea65dd918c017fc10743 100644
--- a/integrators/tps2.py
+++ b/integrators/tps2.py
@@ -6,7 +6,9 @@ from . import _forms
 ################################################################################
 
 
-class TPS2(_tpsx.TPSX, _forms.lf_derivative.LinearForms_Derivative):
+class TPS2(_tpsx.TPSX, \
+  _forms.lf_derivative.LinearForms_Derivative, \
+  _forms.blf_piImplicit.BilinearForms_piImplicit):
   # description of TPS2
 
 
@@ -141,7 +143,7 @@ class TPS2(_tpsx.TPSX, _forms.lf_derivative.LinearForms_Derivative):
     from ngsolve import TaskManager, TET, IntegrationRule, SymbolicLFI, LinearForm, SymbolicBFI, BilinearForm
 
     #blf
-    lhs_stat_cf = self._StiffTerm()
+    lhs_stat_cf = self._StiffTerm() + self._DerivLinearTerms_pi_Implicit_LHS()
     lhs_instat_cf = self._MassTerm() + self._SkewTerm()
     
     self._a_stat = BilinearForm(self._X)
@@ -169,7 +171,7 @@ class TPS2(_tpsx.TPSX, _forms.lf_derivative.LinearForms_Derivative):
     self._f_fix = LinearForm(self._X)
     self._f_fix += SymbolicLFI(rhs_fix_cf.Compile(realcompile=True, wait=True))
 
-    rhs_iter_cf = self._DerivLinearTerms_pi()
+    rhs_iter_cf = self._DerivLinearTerms_pi_Implicit_RHS()
     rhs_iter_cf += self._DerivNonLinearTerms_Pi()
     
     self._f_iter = LinearForm(self._X)