diff --git a/project_euromir/cvxpy_solver.py b/project_euromir/cvxpy_solver.py
index d51ad83..7f9a470 100644
--- a/project_euromir/cvxpy_solver.py
+++ b/project_euromir/cvxpy_solver.py
@@ -35,7 +35,7 @@
     from cvxpy.reductions.solution import Solution
     from cvxpy.reductions.solvers import utilities
     from cvxpy.reductions.solvers.conic_solvers.conic_solver import (
-        ConicSolver, NonNeg, Zero)
+        SOC, ConicSolver, NonNeg, Zero)
 except ImportError as exc: # pragma: no cover
     raise ImportError(
         "Can't use CVXPY interface if CVXPY is not installed!") from exc
@@ -73,7 +73,7 @@ class Solver(ConicSolver):
     """
 
     MIP_CAPABLE = False
-    SUPPORTED_CONSTRAINTS = [Zero, NonNeg]
+    SUPPORTED_CONSTRAINTS = [Zero, NonNeg, SOC]
     REQUIRES_CONSTR = False
 
     def name(self):
@@ -86,7 +86,7 @@ def solve_via_data(
 
         x_orig, y_orig, s_orig = solve(
             matrix=data['A'], b=data['b'], c=data['c'], zero=data['dims'].zero,
-            nonneg=data['dims'].nonneg)
+            nonneg=data['dims'].nonneg, soc=data['dims'].soc)
 
         # breakpoint()
 
diff --git a/project_euromir/loss_no_hsde.py b/project_euromir/loss_no_hsde.py
index 4fc53b9..aa39e8b 100644
--- a/project_euromir/loss_no_hsde.py
+++ b/project_euromir/loss_no_hsde.py
@@ -23,7 +23,7 @@
 # The loss function used in the main loop
 ###
 
-def create_workspace(m, n, zero):
+def create_workspace(m, n, zero, soc=()):
     workspace = {}
 
     # preallocate some variables
@@ -36,7 +36,7 @@ def create_workspace(m, n, zero):
     return workspace
 
 # variable is xy
-def loss_gradient(xy, m, n, zero, matrix, b, c, workspace):
+def loss_gradient(xy, m, n, zero, matrix, b, c, workspace, soc=()):
     """Function for LBFGS loop, used in line search as well."""
 
     x = xy[:n]
@@ -81,7 +81,7 @@ def loss_gradient(xy, m, n, zero, matrix, b, c, workspace):
 
     return loss, workspace['gradient']
 
-def hessian(xy, m, n, zero, matrix, b, c, workspace, regularizer = 0.):
+def hessian(xy, m, n, zero, matrix, b, c, workspace, soc=(), regularizer = 0.):
     """Hessian to use inside LBFGS loop."""
 
     x = xy[:n]
@@ -129,7 +129,7 @@ def _matvec(dxdy):
         matvec=_matvec
     )
 
-def residual(xy, m, n, zero, matrix, b, c):
+def residual(xy, m, n, zero, matrix, b, c, soc=()):
     """Residual function to use L-M approach instead."""
 
     x = xy[:n]
@@ -160,7 +160,7 @@ def residual(xy, m, n, zero, matrix, b, c):
 
     return res
 
-def Dresidual(xy, m, n, zero, matrix, b, c):
+def Dresidual(xy, m, n, zero, matrix, b, c, soc=()):
     """Linear operator to matrix multiply the residual derivative."""
 
     x = xy[:n]
diff --git a/project_euromir/refinement.py b/project_euromir/refinement.py
index 326a96f..1159474 100644
--- a/project_euromir/refinement.py
+++ b/project_euromir/refinement.py
@@ -22,7 +22,7 @@
 import scipy as sp
 
 
-def refine(z, matrix, b, c, zero, nonneg, max_iters=None):
+def refine(z, matrix, b, c, zero, nonneg, soc=(), max_iters=None):
     """All Python for now, will call all the rest.
 
     Equilibration is not done here, data must be already transformed.
diff --git a/project_euromir/solver_new.py b/project_euromir/solver_new.py
index 1eda864..e003f10 100644
--- a/project_euromir/solver_new.py
+++ b/project_euromir/solver_new.py
@@ -43,7 +43,7 @@
 
 QR_PRESOLVE = True
 
-def solve(matrix, b, c, zero, nonneg,
+def solve(matrix, b, c, zero, nonneg, soc=(),
         # xy = None, # need to import logic for equilibration
         ):
     "Main function."
@@ -65,7 +65,7 @@ def solve(matrix, b, c, zero, nonneg,
     d, e, sigma, rho, matrix_transf, b_transf, c_transf = \
     equilibrate.hsde_ruiz_equilibration(
             matrix, b, c, dimensions={
-                'zero': zero, 'nonneg': nonneg, 'second_order': ()},
+                'zero': zero, 'nonneg': nonneg, 'second_order': soc},
             max_iters=25)
 
     if QR_PRESOLVE:
@@ -82,18 +82,18 @@ def solve(matrix, b, c, zero, nonneg,
     def _local_loss(xy):
         return loss_gradient(
             xy, m=m, n=n, zero=zero, matrix=matrix_transf, b=b_transf,
-            c=c_transf, workspace=workspace)[0]
+            c=c_transf, workspace=workspace, soc=soc)[0]
 
     def _local_grad(xy):
         # very important, need to copy the output, to redesign
         return np.copy(loss_gradient(
             xy, m=m, n=n, zero=zero, matrix=matrix_transf, b=b_transf,
-            c=c_transf, workspace=workspace)[1])
+            c=c_transf, workspace=workspace, soc=soc)[1])
 
     def _local_hessian(xy):
         return hessian(
             xy, m=m, n=n, zero=zero, matrix=matrix_transf, b=b_transf,
-            c=c_transf, workspace=workspace)
+            c=c_transf, workspace=workspace, soc=soc)
 
     def _local_hessian_x_nogap(x):
         return hessian_x_nogap(
@@ -106,12 +106,12 @@ def _local_hessian_y_nogap(x):
     def _local_residual(xy):
         return residual(
             xy, m=m, n=n, zero=zero, matrix=matrix_transf, b=b_transf,
-            c=c_transf)
+            c=c_transf, soc=soc)
 
     def _local_derivative_residual(xy):
         return Dresidual(
             xy, m=m, n=n, zero=zero, matrix=matrix_transf, b=b_transf,
-            c=c_transf)
+            c=c_transf, soc=soc)
 
     xy = np.zeros(n+m)
     loss_xy = _local_loss(xy)
@@ -254,7 +254,7 @@ def _local_derivative_residual(xy):
     for _ in range(3):
         u, v = refine(
             z=u-v, matrix=matrix_transf, b=b_transf, c=c_transf, zero=zero,
-            nonneg=nonneg)
+            nonneg=nonneg, soc=soc)
 
     if u[-1] < 1e-8:
         raise FloatingPointError(