diff --git a/src/atoms/second_order_cone/EucNormAtom.jl b/src/atoms/second_order_cone/EucNormAtom.jl
deleted file mode 100644
index 7e635765b..000000000
--- a/src/atoms/second_order_cone/EucNormAtom.jl
+++ /dev/null
@@ -1,32 +0,0 @@
-mutable struct EucNormAtom <: AbstractExpr
-    children::Tuple{AbstractExpr}
-    size::Tuple{Int,Int}
-
-    EucNormAtom(x::AbstractExpr) = new((x,), (1, 1))
-end
-
-head(io::IO, ::EucNormAtom) = print(io, "norm2")
-
-Base.sign(::EucNormAtom) = Positive()
-
-monotonicity(x::EucNormAtom) = (sign(x.children[1]) * Nondecreasing(),)
-
-curvature(::EucNormAtom) = ConvexVexity()
-
-evaluate(x::EucNormAtom) = norm(vec(evaluate(x.children[1])))
-
-function new_conic_form!(context::Context{T}, A::EucNormAtom) where {T}
-    x = only(AbstractTrees.children(A))
-    t = conic_form!(context, Variable())
-    obj = conic_form!(context, x)
-    f = operate(vcat, T, sign(A), t, obj)
-    MOI_add_constraint(context.model, f, MOI.SecondOrderCone(length(x) + 1))
-    return t
-end
-
-function LinearAlgebra.norm2(x::AbstractExpr)
-    if sign(x) == ComplexSign()
-        return EucNormAtom([real(x); imag(x)])
-    end
-    return EucNormAtom(x)
-end
diff --git a/src/atoms/second_order_cone/EuclideanNormAtom.jl b/src/atoms/second_order_cone/EuclideanNormAtom.jl
new file mode 100644
index 000000000..48dd797a4
--- /dev/null
+++ b/src/atoms/second_order_cone/EuclideanNormAtom.jl
@@ -0,0 +1,32 @@
+mutable struct EuclideanNormAtom <: AbstractExpr
+    children::Tuple{AbstractExpr}
+    size::Tuple{Int,Int}
+
+    EuclideanNormAtom(x::AbstractExpr) = new((x,), (1, 1))
+end
+
+head(io::IO, ::EuclideanNormAtom) = print(io, "norm2")
+
+Base.sign(::EuclideanNormAtom) = Positive()
+
+monotonicity(x::EuclideanNormAtom) = (sign(x.children[1]) * Nondecreasing(),)
+
+curvature(::EuclideanNormAtom) = ConvexVexity()
+
+evaluate(x::EuclideanNormAtom) = norm(vec(evaluate(x.children[1])))
+
+function new_conic_form!(context::Context{T}, A::EuclideanNormAtom) where {T}
+    x = only(AbstractTrees.children(A))
+    t = conic_form!(context, Variable())
+    obj = conic_form!(context, x)
+    f = operate(vcat, T, sign(A), t, obj)
+    MOI_add_constraint(context.model, f, MOI.SecondOrderCone(length(x) + 1))
+    return t
+end
+
+function LinearAlgebra.norm2(x::AbstractExpr)
+    if sign(x) == ComplexSign()
+        return EuclideanNormAtom([real(x); imag(x)])
+    end
+    return EuclideanNormAtom(x)
+end
diff --git a/src/atoms/second_order_cone/GeoMeanAtom.jl b/src/atoms/second_order_cone/GeoMeanAtom.jl
index 5bcc1e707..f253cb648 100644
--- a/src/atoms/second_order_cone/GeoMeanAtom.jl
+++ b/src/atoms/second_order_cone/GeoMeanAtom.jl
@@ -1,16 +1,19 @@
+_to_expr(x::AbstractExpr) = x
+_to_expr(x::Value) = constant(x)
+
 mutable struct GeoMeanAtom <: AbstractExpr
-    children::NTuple{N,AbstractExpr} where {N}
+    children::Vector{AbstractExpr}
     size::Tuple{Int,Int}
 
     function GeoMeanAtom(args::AbstractExprOrValue...)
-        args = Tuple(arg isa Value ? constant(arg) : arg for arg in args)
-        sz = size(first(args))
-        if any(!=(sz), size.(args))
-            error("geo mean must take arguments of the same size")
-        elseif any(x -> sign(x) == ComplexSign(), args)
-            error("The arguments should be real, not complex")
+        new_args = AbstractExpr[_to_expr(arg) for arg in args]
+        sz = size(first(new_args))
+        if any(!=(sz), size.(new_args))
+            error("[GeoMeanAtom] geomean must take arguments of the same size")
+        elseif any(x -> sign(x) == ComplexSign(), new_args)
+            error("[GeoMeanAtom] the arguments must be real, not complex")
         end
-        return new(args, sz)
+        return new(new_args, sz)
     end
 end
 
@@ -19,19 +22,14 @@ head(io::IO, ::GeoMeanAtom) = print(io, "geomean")
 Base.sign(::GeoMeanAtom) = Positive()
 
 function monotonicity(q::GeoMeanAtom)
-    return fill(Nondecreasing(), length(q.children))
+    return ntuple(_ -> Nondecreasing(), length(q.children))
 end
 
 curvature(::GeoMeanAtom) = ConcaveVexity()
 
 _geomean(scalar_args...) = prod(scalar_args)^(1 / length(scalar_args))
 
-function evaluate(q::GeoMeanAtom)
-    n = length(q.children)
-    children = evaluate.(q.children)
-    c1 = first(children)
-    return [_geomean((children[i][I] for i in 1:n)...) for I in eachindex(c1)]
-end
+evaluate(q::GeoMeanAtom) = _geomean.(evaluate.(q.children)...)
 
 function new_conic_form!(context::Context{T}, q::GeoMeanAtom) where {T}
     n = length(q.children)
diff --git a/src/atoms/second_order_cone/HuberAtom.jl b/src/atoms/second_order_cone/HuberAtom.jl
index d433b0fbc..a7b8b9269 100644
--- a/src/atoms/second_order_cone/HuberAtom.jl
+++ b/src/atoms/second_order_cone/HuberAtom.jl
@@ -5,9 +5,9 @@ mutable struct HuberAtom <: AbstractExpr
 
     function HuberAtom(x::AbstractExpr, M::Real)
         if sign(x) == ComplexSign()
-            error("Argument must be real")
+            error("[HuberAtom] argument must be real")
         elseif M <= 0
-            error("Huber parameter must by a positive scalar")
+            error("[HuberAtom] parameter must by a positive scalar. Got `M=$M`")
         end
         return new((x,), x.size, M)
     end
@@ -38,7 +38,6 @@ function new_conic_form!(context::Context, x::HuberAtom)
     s = Variable(c.size)
     n = Variable(c.size)
     add_constraint!(context, c == s + n)
-    # objective given by s.^2 + 2 * M * |n|
     return conic_form!(context, square(s) + 2 * x.M * abs(n))
 end
 
diff --git a/src/atoms/second_order_cone/QolElemAtom.jl b/src/atoms/second_order_cone/QolElemAtom.jl
index d43d1f5c1..d20decb7e 100644
--- a/src/atoms/second_order_cone/QolElemAtom.jl
+++ b/src/atoms/second_order_cone/QolElemAtom.jl
@@ -5,7 +5,7 @@ mutable struct QolElemAtom <: AbstractExpr
     function QolElemAtom(x::AbstractExpr, y::AbstractExpr)
         if x.size != y.size
             error(
-                "elementwise quad over lin must take two arguments of the same size",
+                "[QolElemAtom] elementwise quad over lin must take two arguments of the same size",
             )
         end
         return new((x, y), x.size)
@@ -27,39 +27,32 @@ function evaluate(q::QolElemAtom)
 end
 
 function new_conic_form!(context::Context{T}, q::QolElemAtom) where {T}
-    sz = q.children[1].size
-    t = Variable(sz[1], sz[2])
-    t_obj = conic_form!(context, t)
     x, y = q.children
-    for i in 1:length(q.children[1])
+    t = Variable(x.size)
+    for i in 1:length(x)
         add_constraint!(
             context,
             SecondOrderConeConstraint(y[i] + t[i], y[i] - t[i], 2 * x[i]),
         )
     end
     add_constraint!(context, y >= 0)
-    return t_obj
+    return conic_form!(context, t)
 end
 
 qol_elementwise(x::AbstractExpr, y::AbstractExpr) = QolElemAtom(x, y)
 
-function Base.Broadcast.broadcasted(::typeof(^), x::AbstractExpr, k::Int)
-    if k != 2
-        error("raising variables to powers other than 2 is not implemented")
+function square(x::AbstractExpr)
+    if sign(x) == ComplexSign()
+        error(
+            "square of a complex number is not DCP. Did you mean square_modulus?",
+        )
     end
-    return QolElemAtom(x, constant(ones(x.size[1], x.size[2])))
+    return QolElemAtom(x, constant(ones(x.size)))
 end
 
-function Base.Broadcast.broadcasted(
-    ::typeof(Base.literal_pow),
-    ::typeof(^),
-    x::AbstractExpr,
-    ::Val{k},
-) where {k}
-    return Base.Broadcast.broadcasted(^, x, k)
-end
+sumsquares(x::AbstractExpr) = square(norm2(x))
 
-invpos(x::AbstractExpr) = QolElemAtom(constant(ones(x.size[1], x.size[2])), x)
+invpos(x::AbstractExpr) = QolElemAtom(constant(ones(x.size)), x)
 
 function Base.Broadcast.broadcasted(::typeof(/), x::Value, y::AbstractExpr)
     return dotmultiply(constant(x), invpos(y))
@@ -72,13 +65,18 @@ function Base.:/(x::Value, y::AbstractExpr)
     return MultiplyAtom(constant(x), invpos(y))
 end
 
-sumsquares(x::AbstractExpr) = square(norm2(x))
-
-function square(x::AbstractExpr)
-    if sign(x) == ComplexSign()
-        error(
-            "Square of complex number is not DCP. Did you mean square_modulus?",
-        )
+function Base.Broadcast.broadcasted(::typeof(^), x::AbstractExpr, k::Int)
+    if k != 2
+        error("raising variables to powers other than 2 is not implemented")
     end
-    return QolElemAtom(x, constant(ones(x.size[1], x.size[2])))
+    return square(x)
+end
+
+function Base.Broadcast.broadcasted(
+    ::typeof(Base.literal_pow),
+    ::typeof(^),
+    x::AbstractExpr,
+    ::Val{k},
+) where {k}
+    return Base.Broadcast.broadcasted(^, x, k)
 end
diff --git a/src/atoms/second_order_cone/QuadOverLinAtom.jl b/src/atoms/second_order_cone/QuadOverLinAtom.jl
index af314db8b..822b5cc88 100644
--- a/src/atoms/second_order_cone/QuadOverLinAtom.jl
+++ b/src/atoms/second_order_cone/QuadOverLinAtom.jl
@@ -3,8 +3,10 @@ mutable struct QuadOverLinAtom <: AbstractExpr
     size::Tuple{Int,Int}
 
     function QuadOverLinAtom(x::AbstractExpr, y::AbstractExpr)
-        if x.size[2] != 1 && y.size != (1, 1)
-            error("quad over lin arguments must be a vector and a scalar")
+        if x.size[2] != 1 || y.size != (1, 1)
+            error(
+                "[QuadOverLinAtom] quadoverlin arguments must be a vector and a scalar",
+            )
         end
         return new((x, y), (1, 1))
     end
diff --git a/src/variable.jl b/src/variable.jl
index 1405db5a6..73cebecbb 100644
--- a/src/variable.jl
+++ b/src/variable.jl
@@ -102,9 +102,10 @@ sign!(x::AbstractVariable, s::Sign) = x.sign = s
 Returns the current value of `x` if assigned; errors otherwise.
 """
 function evaluate(x::AbstractVariable)
-    return _value(x) === nothing ?
-           error("Value of the variable is yet to be calculated") :
-           output(_value(x))
+    if _value(x) === nothing
+        error("Value of the variable is yet to be calculated")
+    end
+    return output(copy(_value(x)))
 end
 
 """
diff --git a/test/test_atoms.jl b/test/test_atoms.jl
index 6ee2ef31d..d299fcb1b 100644
--- a/test/test_atoms.jl
+++ b/test/test_atoms.jl
@@ -984,6 +984,213 @@ function test_SumLargestAtom()
     return
 end
 
+### second_order_cone/EuclideanNormAtom
+
+function test_EuclideanNormAtom()
+    target = """
+    variables: t, x, y
+    minobjective: 1.0 * t
+    [1.0 * t, 1.0 * x, 1.0 * y] in SecondOrderCone(3)
+    """
+    _test_atom(target) do context
+        return norm2(Variable(2))
+    end
+    target = """
+    variables: t, x, y
+    minobjective: 1.0 * t
+    [1.0 * t, 1.0 * x, 2.0 * x] in SecondOrderCone(3)
+    """
+    _test_atom(target) do context
+        x = Variable()
+        return norm2((1 + 2im) * x)
+    end
+    return
+end
+
+### second_order_cone/GeoMeanAtom
+
+function test_GeoMeanAtom()
+    target = """
+    variables: t, x, y
+    minobjective: 1.0 * t
+    [1.0 * t, 1.0 * x, 1.0 * y] in GeometricMeanCone(3)
+    """
+    _test_atom(target) do context
+        return geomean(Variable(), Variable())
+    end
+    target = """
+    variables: t1, t2, x1, x2, y1, y2
+    minobjective: [1.0 * t1, 1.0 * t2]
+    [1.0 * t1, 1.0 * x1, 1.0 * y1] in GeometricMeanCone(3)
+    [1.0 * t2, 1.0 * x2, 1.0 * y2] in GeometricMeanCone(3)
+    """
+    _test_atom(target) do context
+        return geomean(Variable(2), Variable(2))
+    end
+    @test_throws(
+        ErrorException(
+            "[GeoMeanAtom] geomean must take arguments of the same size",
+        ),
+        geomean(Variable(2), Variable(3)),
+    )
+    @test_throws(
+        ErrorException("[GeoMeanAtom] the arguments must be real, not complex"),
+        geomean(Variable(), 2im),
+    )
+    x = Variable(2)
+    x.value = [2.0, 3.0]
+    y = Variable(2)
+    y.value = [4.0, 5.0]
+    atom = geomean(x, y)
+    @test evaluate(atom) ≈ [sqrt(8), sqrt(15)]
+    return
+end
+
+### second_order_cone/HuberAtom
+
+function test_HuberAtom()
+    target = """
+    variables: c, s, n, t, n_abs
+    minobjective: 1.0 * t + 4.0 * n_abs
+    [1.0 * c + -1.0 * s + -1.0 * n] in Zeros(1)
+    [-1.0 * n + 1.0 * n_abs] in Nonnegatives(1)
+    [1.0 * n + 1.0 * n_abs] in Nonnegatives(1)
+    [1.0 + 1.0 * t, 1.0 + -1.0 * t, 2.0 * s] in SecondOrderCone(3)
+    """
+    _test_atom(target) do context
+        return huber(Variable(), 2.0)
+    end
+    @test_throws(
+        ErrorException(
+            "[HuberAtom] parameter must by a positive scalar. Got `M=-2.0`",
+        ),
+        huber(Variable(2), -2.0),
+    )
+    @test_throws(
+        ErrorException("[HuberAtom] argument must be real"),
+        huber(im * Variable(2)),
+    )
+    x = Variable(2)
+    x.value = [2.0, 3.0]
+    atom = huber(x, 2.0)
+    @test evaluate(atom) ≈ [4.0, 8.0]
+    atom = huber(x)
+    @test evaluate(atom) ≈ [3.0, 5.0]
+    return
+end
+
+### second_order_cone/QolElemAtom
+
+function test_QolElemAtom()
+    target = """
+    variables: y1, y2, t1, t2, x1, x2
+    minobjective: [1.0 * t1, 1.0 * t2]
+    [1.0 * y1, 1.0 * y2] in Nonnegatives(2)
+    [1.0 * y1 + 1.0 * t1, 1.0 * y1 + -1.0 * t1, 2.0 * x1] in SecondOrderCone(3)
+    [1.0 * y2 + 1.0 * t2, 1.0 * y2 + -1.0 * t2, 2.0 * x2] in SecondOrderCone(3)
+    """
+    _test_atom(target) do context
+        return qol_elementwise(Variable(2), Variable(2))
+    end
+    target = """
+    variables: t1, t2, x1, x2
+    minobjective: [1.0 * t1, 1.0 * t2]
+    [1.0 + 1.0 * t1, 1.0 + -1.0 * t1, 2.0 * x1] in SecondOrderCone(3)
+    [1.0 + 1.0 * t2, 1.0 + -1.0 * t2, 2.0 * x2] in SecondOrderCone(3)
+    """
+    _test_atom(target) do context
+        return qol_elementwise(Variable(2), constant([1, 1]))
+    end
+    _test_atom(target) do context
+        return square(Variable(2))
+    end
+    _test_atom(target) do context
+        return Variable(2) .^ 2
+    end
+    _test_atom(target) do context
+        a = 2
+        return Variable(2) .^ a
+    end
+    target = """
+    variables: y1, y2, t1, t2
+    minobjective: [1.0 * t1, 1.0 * t2]
+    [1.0 * y1, 1.0 * y2] in Nonnegatives(2)
+    [1.0 * y1 + 1.0 * t1, 1.0 * y1 + -1.0 * t1, 2.0] in SecondOrderCone(3)
+    [1.0 * y2 + 1.0 * t2, 1.0 * y2 + -1.0 * t2, 2.0] in SecondOrderCone(3)
+    """
+    _test_atom(target) do context
+        return invpos(Variable(2))
+    end
+    _test_atom(target) do context
+        return 1 ./ Variable(2)
+    end
+    target = """
+    variables: y, t
+    minobjective: 3.0 * t
+    [1.0 * y] in Nonnegatives(1)
+    [1.0 * y + 1.0 * t, 1.0 * y + -1.0 * t, 2.0] in SecondOrderCone(3)
+    """
+    _test_atom(target) do context
+        return 3 / Variable()
+    end
+    y = Variable(2)
+    @test_throws(
+        ErrorException("cannot divide by a variable of size $(size(y))"),
+        1 / y,
+    )
+
+    @test_throws(
+        ErrorException(
+            "square of a complex number is not DCP. Did you mean square_modulus?",
+        ),
+        square(im * Variable()),
+    )
+    @test_throws(
+        ErrorException(
+            "raising variables to powers other than 2 is not implemented",
+        ),
+        Variable() .^ 3,
+    )
+    @test_throws(
+        ErrorException(
+            "[QolElemAtom] elementwise quad over lin must take two arguments of the same size",
+        ),
+        qol_elementwise(Variable(2), Variable())
+    )
+    x = Variable(2)
+    x.value = [2.0, 3.0]
+    y = Variable(2)
+    y.value = [4.0, 5.0]
+    atom = qol_elementwise(x, y)
+    @test evaluate(atom) ≈ [1.0, 9.0 / 5.0]
+    return
+end
+
+### second_order_cone/QuadOverLinAtom
+
+function test_QuadOverLinAtom()
+    target = """
+    variables: y, t, x1, x2
+    minobjective: 1.0 * t
+    [1.0 * y] in Nonnegatives(1)
+    [1.0 * y + 1.0 * t, 1.0 * y + -1.0 * t, 2.0 * x1, 2.0 * x2] in SecondOrderCone(4)
+    """
+    _test_atom(target) do context
+        return quadoverlin(Variable(2), Variable())
+    end
+    @test_throws(
+        ErrorException(
+            "[QuadOverLinAtom] quadoverlin arguments must be a vector and a scalar",
+        ),
+        quadoverlin(Variable(2), Variable(2))
+    )
+    x = Variable(2)
+    x.value = [2.0, 3.0]
+    atom = quadoverlin(x, constant(2.0))
+    @test evaluate(atom) ≈ 13 / 2
+    return
+end
+
 ### second_order_cone/RationalNormAtom
 
 function test_RationalNormAtom()