diff --git a/docs/src/manual/operations.md b/docs/src/manual/operations.md
index 2d0a26c29..87f76aa87 100644
--- a/docs/src/manual/operations.md
+++ b/docs/src/manual/operations.md
@@ -117,7 +117,7 @@ solver (including SCS and Mosek).
 | `matrixfrac(x, P)`                | $x^TP^{-1}x$                       | convex       | not monotonic | IC: P is positive semidefinite |
 | `sumlargesteigs(x, k)`            | sum of top $k$ eigenvalues of $x$  | convex       | not monotonic | IC: P symmetric        |
 | `T in GeomMeanHypoCone(A, B, t)`  | $T \preceq A \#_t B = A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}$ | concave | increasing    | IC: $A \succeq 0$, $B \succeq 0$, $t \in [0,1]$ |
-| `T in GeomMeanEpiCone(A, B, t)`   | $T \succeq A \#_t B = A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}$ | convex  | not monotonic | IC: $A \succeq 0$, $B \succeq 0$, $t \in [-1, 0] \cup [1, 2]$ |
+| `Convex.GenericConstraint((T, A, B), GeometricMeanEpiConeSquare(t, size(T, 1)))`   | $T \succeq A \#_t B = A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}$ | convex  | not monotonic | IC: $A \succeq 0$, $B \succeq 0$, $t \in [-1, 0] \cup [1, 2]$ |
 | `quantum_entropy(X)`              | $-\textrm{Tr}(X \log X)$           | concave      | not monotonic | IC: $X \succeq 0$; uses natural log |
 | `quantum_relative_entropy(A, B)`  | $\textrm{Tr}(A \log A - A \log B)$ | convex       | not monotonic | IC: $A \succeq 0$, $B \succeq 0$; uses natural log |
 | `trace_logm(X, C)`                | $\textrm{Tr}(C \log X)$            | concave in X | not monotonic | IC: $X \succeq 0$, $C \succeq 0$, $C$ constant; uses natural log |
diff --git a/src/Convex.jl b/src/Convex.jl
index 3797f13b7..a4aff3901 100644
--- a/src/Convex.jl
+++ b/src/Convex.jl
@@ -45,7 +45,7 @@ export conv,
     sumsmallest,
     sumsquares,
     GeomMeanHypoCone,
-    GeomMeanEpiCone,
+    GeometricMeanEpiConeSquare,
     RelativeEntropyEpiCone,
     quantum_relative_entropy,
     quantum_entropy,
diff --git a/src/atoms/TraceMpowerAtom.jl b/src/atoms/TraceMpowerAtom.jl
index 5ae57c4fe..f0440ab78 100644
--- a/src/atoms/TraceMpowerAtom.jl
+++ b/src/atoms/TraceMpowerAtom.jl
@@ -95,7 +95,13 @@ function new_conic_form!(context::Context{T}, atom::TraceMpowerAtom) where {T}
     if 0 <= atom.t <= 1
         add_constraint!(context, tmp in GeomMeanHypoCone(I, A, atom.t, false))
     else
-        add_constraint!(context, tmp in GeomMeanEpiCone(I, A, atom.t, false))
+        add_constraint!(
+            context,
+            Convex.GenericConstraint(
+                (tmp, I, A),
+                GeometricMeanEpiConeSquare(atom.t, size(A, 1)),
+            ),
+        )
     end
     # It's already a real mathematically, but Convex doesn't know it
     return conic_form!(context, real(LinearAlgebra.tr(atom.C * tmp)))
diff --git a/src/constraints/GeomMeanEpiCone.jl b/src/constraints/GeomMeanEpiCone.jl
deleted file mode 100644
index d2ecd14d5..000000000
--- a/src/constraints/GeomMeanEpiCone.jl
+++ /dev/null
@@ -1,166 +0,0 @@
-# Copyright (c) 2014: Madeleine Udell and contributors
-# Copyright (c) 2021: Hamza Fawzi
-#
-# Use of this source code is governed by an MIT-style license that can be found
-# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
-
-"""
-Constrains T to
-  A #_t B ⪯ T
-where:
-  * A #_t B is the t-weighted geometric mean of A and B:
-    A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}
-    Parameter t should be in [-1, 0] or [1, 2].
-  * Constraints A ⪰ 0, B ⪰ 0 are added.
-All expressions and atoms are subtypes of AbstractExpr.
-Please read expressions.jl first.
-
-Note on fullhyp: GeomMeanEpiCone will always return a full epigraph cone
-(unlike GeomMeanHypoCone) and so this parameter is not really used.  It is
-here just for consistency with the GeomMeanHypoCone function.
-
-REFERENCE
-  Ported from CVXQUAD which is based on the paper: "Lieb's concavity
-  theorem, matrix geometric means and semidefinite optimization" by Hamza
-  Fawzi and James Saunderson (arXiv:1512.03401)
-"""
-mutable struct GeomMeanEpiCone
-    A::AbstractExpr
-    B::AbstractExpr
-    t::Rational
-    size::Tuple{Int,Int}
-
-    function GeomMeanEpiCone(
-        A::AbstractExpr,
-        B::AbstractExpr,
-        t::Rational,
-        fullhyp::Bool = true,
-    )
-        if size(A) != size(B)
-            throw(DimensionMismatch("A and B must be the same size"))
-        end
-        n = size(A)[1]
-        if size(A) != (n, n)
-            throw(DimensionMismatch("A and B must be square"))
-        end
-        if t < -1 || (t > 0 && t < 1) || t > 2
-            throw(DomainError(t, "t must be in the range [-1, 0] or [1, 2]"))
-        end
-        return new(A, B, t, (n, n))
-    end
-
-    function GeomMeanEpiCone(
-        A::Value,
-        B::AbstractExpr,
-        t::Rational,
-        fullhyp::Bool = true,
-    )
-        return GeomMeanEpiCone(constant(A), B, t, fullhyp)
-    end
-
-    function GeomMeanEpiCone(
-        A::AbstractExpr,
-        B::Value,
-        t::Rational,
-        fullhyp::Bool = true,
-    )
-        return GeomMeanEpiCone(A, constant(B), t, fullhyp)
-    end
-
-    function GeomMeanEpiCone(
-        A::Value,
-        B::Value,
-        t::Rational,
-        fullhyp::Bool = true,
-    )
-        return GeomMeanEpiCone(constant(A), constant(B), t, fullhyp)
-    end
-
-    function GeomMeanEpiCone(
-        A::Union{AbstractExpr,Value},
-        B::Union{AbstractExpr,Value},
-        t::Integer,
-        fullhyp::Bool = true,
-    )
-        return GeomMeanEpiCone(A, B, t // 1, fullhyp)
-    end
-end
-
-mutable struct GeomMeanEpiConeConstraint <: Constraint
-    T::AbstractExpr
-    cone::GeomMeanEpiCone
-
-    function GeomMeanEpiConeConstraint(T::AbstractExpr, cone::GeomMeanEpiCone)
-        if size(T) != cone.size
-            throw(DimensionMismatch("T must be size $(cone.size)"))
-        end
-        return new(T, cone)
-    end
-
-    function GeomMeanEpiConeConstraint(T::Value, cone::GeomMeanEpiCone)
-        return GeomMeanEpiConeConstraint(constant(T), cone)
-    end
-end
-
-head(io::IO, ::GeomMeanEpiConeConstraint) = print(io, "∈(GeomMeanEpiCone)")
-
-Base.in(T, cone::GeomMeanEpiCone) = GeomMeanEpiConeConstraint(T, cone)
-
-function AbstractTrees.children(constraint::GeomMeanEpiConeConstraint)
-    return (
-        constraint.T,
-        constraint.cone.A,
-        constraint.cone.B,
-        "t=$(constraint.cone.t)",
-    )
-end
-
-# For t ∈ [-1, 0] ∪ [1, 2] the t-weighted matrix geometric mean is matrix convex
-# (arxiv:1512.03401).
-# So if A and B are convex sets, then T ⪰ A #_t B will be a convex set.
-function vexity(constraint::GeomMeanEpiConeConstraint)
-    A = vexity(constraint.cone.A)
-    B = vexity(constraint.cone.B)
-    T = vexity(constraint.T)
-
-    # NOTE: can't say A == NotDcp() because the NotDcp constructor prints a
-    # warning message.
-    if typeof(A) == ConcaveVexity || typeof(A) == NotDcp
-        return NotDcp()
-    elseif typeof(B) == ConcaveVexity || typeof(B) == NotDcp
-        return NotDcp()
-    end
-    vex = ConvexVexity() + (-T)
-    if vex == ConcaveVexity()
-        return NotDcp()
-    end
-    return vex
-end
-
-function _add_constraint!(
-    context::Context,
-    constraint::GeomMeanEpiConeConstraint,
-)
-    A = constraint.cone.A
-    B = constraint.cone.B
-    t = constraint.cone.t
-    T = constraint.T
-    is_complex =
-        sign(A) == ComplexSign() ||
-        sign(B) == ComplexSign() ||
-        sign(T) == ComplexSign()
-    Z = if is_complex
-        HermitianSemidefinite(size(A)[1])
-    else
-        Semidefinite(size(A)[1])
-    end
-    if t <= 0
-        add_constraint!(context, [T A; A Z] ⪰ 0)
-        add_constraint!(context, Z in GeomMeanHypoCone(A, B, -t, false))
-    else
-        @assert t >= 1 # range of t checked in GeomMeanEpiCone constructor
-        add_constraint!(context, [T B; B Z] ⪰ 0)
-        add_constraint!(context, Z in GeomMeanHypoCone(A, B, 2 - t, false))
-    end
-    return
-end
diff --git a/src/constraints/GeometricMeanEpiCone.jl b/src/constraints/GeometricMeanEpiCone.jl
new file mode 100644
index 000000000..3c7a06ad0
--- /dev/null
+++ b/src/constraints/GeometricMeanEpiCone.jl
@@ -0,0 +1,102 @@
+# Copyright (c) 2014: Madeleine Udell and contributors
+# Copyright (c) 2021: Hamza Fawzi
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+"""
+    GeometricMeanEpiConeSquare(t::Rational, side_dimension::Int)
+
+The constraint `(T, A, B) in GeometricMeanEpiConeSquare(t, side_dimension)`
+constrains T to
+
+  A #_t B ⪯ T
+
+where:
+
+  * A #_t B is the t-weighted geometric mean of A and B:
+    A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}
+  * Parameter t must be in [-1, 0] or [1, 2].
+  * Constraints A ⪰ 0, B ⪰ 0 are added.
+
+## Reference
+
+Ported from CVXQUAD which is based on the paper: "Lieb's concavity theorem,
+matrix geometric means and semidefinite optimization" by Hamza Fawzi and James
+Saunderson (arXiv:1512.03401)
+"""
+struct GeometricMeanEpiConeSquare <: MOI.AbstractVectorSet
+    t::Rational
+    side_dimension::Int
+
+    function GeometricMeanEpiConeSquare(t::Rational, side_dimension::Int)
+        if !(-1 <= t <= 0 || 1 <= t <= 2)
+            throw(DomainError(t, "t must be in the range [-1, 0] or [1, 2]"))
+        end
+        return new(t, side_dimension)
+    end
+end
+
+MOI.dimension(set::GeometricMeanEpiConeSquare) = 3 * set.side_dimension^2
+
+function head(io::IO, ::GeometricMeanEpiConeSquare)
+    return print(io, "GeometricMeanEpiConeSquare")
+end
+
+function GenericConstraint(func::Tuple, set::GeometricMeanEpiConeSquare)
+    for f in func
+        n = LinearAlgebra.checksquare(f)
+        if n != set.side_dimension
+            throw(
+                DimensionMismatch(
+                    "Matrix of side dimension `$n` does not match set of side dimension `$(set.side_dimension)`",
+                ),
+            )
+        end
+    end
+    return GenericConstraint(vcat(vec.(func)...), set)
+end
+
+function _get_matrices(c::GenericConstraint{GeometricMeanEpiConeSquare})
+    n = c.set.side_dimension
+    d = n^2
+    T = reshape(c.child[1:d], n, n)
+    A = reshape(c.child[d.+(1:d)], n, n)
+    B = reshape(c.child[2d.+(1:d)], n, n)
+    return T, A, B
+end
+
+# For t ∈ [-1, 0] ∪ [1, 2] the t-weighted matrix geometric mean is matrix convex
+# (arxiv:1512.03401).
+# So if A and B are convex sets, then T ⪰ A #_t B will be a convex set.
+function vexity(constraint::GenericConstraint{GeometricMeanEpiConeSquare})
+    T, A, B = _get_matrices(constraint)
+    if vexity(A) in (ConcaveVexity(), NotDcp()) ||
+       vexity(B) in (ConcaveVexity(), NotDcp())
+        return NotDcp()
+    end
+    return ConvexVexity() + -vexity(T)
+end
+
+function _add_constraint!(
+    context::Context,
+    constraint::GenericConstraint{GeometricMeanEpiConeSquare},
+)
+    T, A, B = _get_matrices(constraint)
+    t = constraint.set.t
+    Z = if sign(constraint.child) == ComplexSign()
+        HermitianSemidefinite(size(A, 1))
+    else
+        Semidefinite(size(A, 1))
+    end
+    if t <= 0
+        add_constraint!(context, [T A; A Z] ⪰ 0)
+        add_constraint!(context, Z in GeomMeanHypoCone(A, B, -t, false))
+    else
+        # range of t checked in GeometricMeanEpiConeSquare constructor.
+        @assert t >= 1
+        add_constraint!(context, [T B; B Z] ⪰ 0)
+        add_constraint!(context, Z in GeomMeanHypoCone(A, B, 2 - t, false))
+    end
+    return
+end
diff --git a/src/expressions.jl b/src/expressions.jl
index 958302f69..8ca8e5534 100644
--- a/src/expressions.jl
+++ b/src/expressions.jl
@@ -75,6 +75,8 @@ function vexity(x::AbstractExpr)
     return vex
 end
 
+vexity(::Value) = ConstVexity()
+
 evaluate(x) = output(x) # fallback
 
 Base.size(x::AbstractExpr) = x.size
diff --git a/src/problem_depot/problems/sdp.jl b/src/problem_depot/problems/sdp.jl
index 469dc4038..6f6b00943 100644
--- a/src/problem_depot/problems/sdp.jl
+++ b/src/problem_depot/problems/sdp.jl
@@ -739,79 +739,6 @@ end
     end
 end
 
-@add_problem sdp function sdp_geom_mean_epicone_argcheck(
-    handle_problem!,
-    ::Val{test},
-    atol,
-    rtol,
-    ::Type{T},
-) where {T,test}
-    if test
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            zeros(2, 3),
-            zeros(2, 3),
-            -1 // 2,
-        )
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            zeros(2, 3),
-            Variable(2, 3),
-            -1 // 2,
-        )
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            Variable(2, 3),
-            zeros(2, 3),
-            -1 // 2,
-        )
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            Variable(2, 3),
-            Variable(2, 3),
-            -1 // 2,
-        )
-
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            zeros(2, 2),
-            zeros(3, 3),
-            -1 // 2,
-        )
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            zeros(2, 2),
-            Variable(3, 3),
-            -1 // 2,
-        )
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            Variable(2, 2),
-            zeros(3, 3),
-            -1 // 2,
-        )
-        @test_throws DimensionMismatch GeomMeanEpiCone(
-            Variable(2, 2),
-            Variable(2, 3),
-            -1 // 2,
-        )
-
-        @test_throws DimensionMismatch Variable(2, 2) in GeomMeanEpiCone(
-            Variable(3, 3),
-            Variable(3, 3),
-            -1 // 2,
-        )
-        @test_throws DomainError GeomMeanEpiCone(
-            Variable(3, 3),
-            Variable(3, 3),
-            -3 // 2,
-        )
-        @test_throws DomainError GeomMeanEpiCone(
-            Variable(3, 3),
-            Variable(3, 3),
-            1 // 2,
-        )
-        @test_throws DomainError GeomMeanEpiCone(
-            Variable(3, 3),
-            Variable(3, 3),
-            5 // 2,
-        )
-    end
-end
-
 @add_problem sdp function sdp_relative_entropy_argcheck(
     handle_problem!,
     ::Val{test},
@@ -1245,7 +1172,10 @@ end
     B += 0.2 * LinearAlgebra.I # prevent numerical instability
     A = Variable(n, n)
 
-    c1 = eye(n) in GeomMeanEpiCone(A, B, -1)
+    c1 = Convex.GenericConstraint(
+        (eye(n), A, B),
+        GeometricMeanEpiConeSquare(-1 // 1, n),
+    )
     objective = tr(A)
     p = maximize(objective, c1; numeric_type = T)
 
@@ -1271,7 +1201,10 @@ end
     A /= tr(A) # solver has problems if B is large
     B = Variable(n, n)
 
-    c1 = eye(n) in GeomMeanEpiCone(A, B, -1)
+    c1 = Convex.GenericConstraint(
+        (eye(n), A, B),
+        GeometricMeanEpiConeSquare(-1 // 1, n),
+    )
     objective = tr(B)
     p = minimize(objective, c1; numeric_type = T)
 
@@ -1296,7 +1229,10 @@ end
     B += 0.2 * LinearAlgebra.I # prevent numerical instability
     A = Variable(n, n)
 
-    c1 = eye(n) in GeomMeanEpiCone(A, B, -3 // 5)
+    c1 = Convex.GenericConstraint(
+        (eye(n), A, B),
+        GeometricMeanEpiConeSquare(-3 // 5, n),
+    )
     objective = tr(A)
     p = maximize(objective, c1; numeric_type = T)
 
@@ -1322,7 +1258,10 @@ end
     A /= tr(A) # solver has problems if B is large
     B = Variable(n, n)
 
-    c1 = eye(n) in GeomMeanEpiCone(A, B, -3 // 5)
+    c1 = Convex.GenericConstraint(
+        (eye(n), A, B),
+        GeometricMeanEpiConeSquare(-3 // 5, n),
+    )
     objective = tr(B)
     p = minimize(objective, c1; numeric_type = T)
 
@@ -1348,7 +1287,10 @@ end
     B /= tr(B) # solver has problems if B is large
     A = Variable(n, n)
 
-    c1 = eye(n) in GeomMeanEpiCone(A, B, 8 // 5)
+    c1 = Convex.GenericConstraint(
+        (eye(n), A, B),
+        GeometricMeanEpiConeSquare(8 // 5, n),
+    )
     objective = tr(A)
     p = minimize(objective, c1; numeric_type = T)
 
@@ -1373,7 +1315,10 @@ end
     A += 0.2 * LinearAlgebra.I # prevent numerical instability
     B = Variable(n, n)
 
-    c1 = eye(n) in GeomMeanEpiCone(A, B, 8 // 5)
+    c1 = Convex.GenericConstraint(
+        (eye(n), A, B),
+        GeometricMeanEpiConeSquare(8 // 5, n),
+    )
     objective = tr(B)
     p = maximize(objective, c1; numeric_type = T)
 
diff --git a/src/reformulations/lieb_ando.jl b/src/reformulations/lieb_ando.jl
index cfd3b59ea..b61d5265b 100644
--- a/src/reformulations/lieb_ando.jl
+++ b/src/reformulations/lieb_ando.jl
@@ -97,7 +97,10 @@ function lieb_ando(
         # Convex function
         add_constraint!(
             T,
-            T in GeomMeanEpiCone(kron(A, Im), kron(In, conj(B)), t, false),
+            Convex.GenericConstraint(
+                (T, kron(A, Im), kron(In, conj(B))),
+                GeometricMeanEpiConeSquare(t, size(T, 1)),
+            ),
         )
         return real(LinearAlgebra.tr(KvKv * T))
     else
diff --git a/test/test_constraints.jl b/test/test_constraints.jl
index b6424425b..84eb59c54 100644
--- a/test/test_constraints.jl
+++ b/test/test_constraints.jl
@@ -297,6 +297,45 @@ function test_GenericConstraint_RelativeEntropyConeConstraint()
     return
 end
 
+### constraints/GenericConstraint_GeometricMeanEpiConeSquare
+
+function test_GeometricMeanEpiConeSquare()
+    T = Variable(2, 2)
+    A = Variable(2, 2)
+    B = Variable(2, 2)
+    C = [1 0; 0 1]
+    @test_throws DomainError GeometricMeanEpiConeSquare(-3 // 2, 3)
+    @test_throws DomainError GeometricMeanEpiConeSquare(1 // 2, 3)
+    @test_throws DomainError GeometricMeanEpiConeSquare(5 // 2, 3)
+    set = Convex.GeometricMeanEpiConeSquare(3 // 2, 2)
+    @test sprint(Convex.head, set) == "GeometricMeanEpiConeSquare"
+    @test MOI.dimension(set) == 12
+    for (f, dcp) in (
+        (T, A, B) => Convex.ConvexVexity(),
+        (T, A, C) => Convex.ConvexVexity(),
+        (T, C, B) => Convex.ConvexVexity(),
+        (C, A, B) => Convex.ConvexVexity(),
+        (T * T, A, B) => Convex.NotDcp(),
+        (T .^ 2, A, B) => Convex.NotDcp(),
+        (T, A * A, B) => Convex.NotDcp(),
+        (T, sqrt(A), B) => Convex.NotDcp(),
+        (T, A, B * B) => Convex.NotDcp(),
+        (T, A, sqrt(B)) => Convex.NotDcp(),
+        (sqrt(T), A, B) => Convex.ConvexVexity(),
+        (T + C, A, B) => Convex.ConvexVexity(),
+        (T * C, A, B) => Convex.ConvexVexity(),
+    )
+        c = Convex.GenericConstraint(f, set)
+        @test vexity(c) == dcp
+        @test sign(c.child) == Convex.NoSign()
+    end
+    Z = Variable(3, 3)
+    @test_throws DimensionMismatch Convex.GenericConstraint((Z, A, B), set)
+    @test_throws DimensionMismatch Convex.GenericConstraint((T, Z, B), set)
+    @test_throws DimensionMismatch Convex.GenericConstraint((T, A, Z), set)
+    return
+end
+
 end  # TestConstraints
 
 TestConstraints.runtests()