diff --git a/src/CommonOPF.jl b/src/CommonOPF.jl
index 65cd26d..194ab9b 100644
--- a/src/CommonOPF.jl
+++ b/src/CommonOPF.jl
@@ -53,13 +53,14 @@ export
     next_bus_above_with_outdegree_more_than_one,
     paths_between,
 
-    # impedance
+    # impedance and shunt admittance
     rij,
     rij_per_unit,
     xij,
     xij_per_unit,
     zij,
     zij_per_unit,
+    yj,
 
     # io
     dss_to_Network,
@@ -121,6 +122,7 @@ include("bounds.jl")
 
 include("network.jl")
 include("edges/impedances.jl")  # Network type in signatures, move the struct to types?
+include("busses/shunt_admittances.jl")
 include("network_reduction.jl")
 include("decomposition.jl")
 
diff --git a/src/busses/busses.jl b/src/busses/busses.jl
index a764c35..e1fa4b1 100644
--- a/src/busses/busses.jl
+++ b/src/busses/busses.jl
@@ -19,7 +19,7 @@ function build_busses(
         d[:bus] = string(d[:bus])
     end
     busses = ConcreteBusType[ConcreteBusType(;bdict...) for bdict in dicts]
-    check_busses!(busses)  # dispatch on Vector{}
+    check_busses!(busses)  # dispatch on AbstractVector{ConcreteBusType}
     return busses
 end
 
diff --git a/src/busses/shunt_admittances.jl b/src/busses/shunt_admittances.jl
new file mode 100644
index 0000000..4ced045
--- /dev/null
+++ b/src/busses/shunt_admittances.jl
@@ -0,0 +1,24 @@
+"""
+    yj(j::AbstractString, net::Network{SinglePhase})::ComplexF64
+
+Shunt admittance of bus `j`
+"""
+function yj(j::AbstractString, net::Network{SinglePhase})::ComplexF64
+    if :ShuntAdmittance in keys(net[j])
+        return net[j][:ShuntAdmittance].g + im * net[j][:ShuntAdmittance].b
+    end
+    return 0.0 + im * 0.0
+end
+
+
+"""
+    yj(j::AbstractString, net::Network{MultiPhase})::Matrix{ComplexF64}
+
+Shunt admittance of bus `j`
+"""
+function yj(j::AbstractString, net::Network{MultiPhase})::Matrix{ComplexF64}
+    if :ShuntAdmittance in keys(net[j])
+        return net[j][:ShuntAdmittance].gmatrix + im * net[j][:ShuntAdmittance].bmatrix
+    end
+    return zeros(3,3) + im * zeros(3,3)
+end
diff --git a/src/busses/shunts.jl b/src/busses/shunts.jl
index b423ec7..b5b9216 100644
--- a/src/busses/shunts.jl
+++ b/src/busses/shunts.jl
@@ -9,7 +9,8 @@ Required fields:
 @with_kw struct ShuntAdmittance <: AbstractBus
     # required values
     bus::String
-    g::Real
-    b::Real
-    # TODO matrices for MultiPhase models ?
+    g::Real = 0.0
+    b::Real = 0.0
+    gmatrix::AbstractArray = zeros(3,3)
+    bmatrix::AbstractArray = zeros(3,3)
 end
diff --git a/src/edges/conductors.jl b/src/edges/conductors.jl
index 43aefde..bbe5406 100644
--- a/src/edges/conductors.jl
+++ b/src/edges/conductors.jl
@@ -93,6 +93,9 @@ Conductor:
     The order of the `phases` is assumed to match the order of the `rmatrix` and `xmatrix`. For
     example using the example just above the 3x3 `rmatrix` looks like 
     ``[0.31, 0, 0.15; 0, 0, 0; 0.15, 0, 0.32]``
+
+Conductors also have a `cmatrix` attribute that is used when parsing OpenDSS models. The `cmatrix`
+is used to define `ShuntAdmittance` values for busses.
 """
 @with_kw mutable struct Conductor <: AbstractEdge
     # mutable because we set the rmatrix and xmatrix later in some cases
@@ -106,8 +109,10 @@ Conductor:
     x0::Union{Real, Missing} = missing
     r1::Union{Real, Missing} = missing
     x1::Union{Real, Missing} = missing
+    c1::Union{Real, Missing} = missing
     rmatrix::Union{AbstractArray, Missing} = missing
     xmatrix::Union{AbstractArray, Missing} = missing
+    cmatrix::Union{AbstractArray, Missing} = missing
     length::Union{Real, Missing} = missing
     amps_limit::Union{Real, Missing} = missing
 end
@@ -266,6 +271,7 @@ function validate_multiphase_edges!(conds::AbstractVector{Conductor})::Bool
             end
             c.xmatrix = template_cond.xmatrix
             c.rmatrix = template_cond.rmatrix
+            c.cmatrix = template_cond.cmatrix
         end
     end
 
diff --git a/src/io.jl b/src/io.jl
index 479b83c..71fcc82 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -63,16 +63,18 @@ end
 resize OpenDSS matrices to 3x3 and sort values by numerical phase order
 """
 function dss_impedance_matrices_to_three_phase(
-    dss_rmatrix::AbstractMatrix{Float64}, 
-    dss_xmatrix::AbstractMatrix{Float64}, 
+    dss_rmatrix::AbstractMatrix{Float64},
+    dss_xmatrix::AbstractMatrix{Float64},
+    dss_cmatrix::AbstractMatrix{Float64}, 
     phases::AbstractVector{Int}
-    )::Tuple{Matrix{Float64}, Matrix{Float64}}
-    r, x = zeros(3,3), zeros(3,3)
+    )::Tuple{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}
+    r, x, c = zeros(3,3), zeros(3,3), zeros(3,3)
     for (i, phs1) in enumerate(phases), (j, phs2) in enumerate(phases)
         r[phs1, phs2] = dss_rmatrix[i, j]
         x[phs1, phs2] = dss_xmatrix[i, j]
+        c[phs1, phs2] = dss_cmatrix[i, j]
     end
-    return r, x
+    return r, x, c
 end
 
 
@@ -220,6 +222,7 @@ function dss_to_Network(dssfilepath::AbstractString)::Network
     node_order = [lowercase(s) for s in OpenDSS.Circuit.YNodeOrder()]
     num_phases = opendss_circuit_num_phases()
     
+    conductors = opendss_lines(num_phases)
     net_dict = Dict{Symbol, Any}(
         :Network => Dict(
             :substation_bus => opendss_source_bus(),
@@ -227,10 +230,11 @@ function dss_to_Network(dssfilepath::AbstractString)::Network
             :Vbase => OpenDSS.Vsources.BasekV() * 1e3,
             # :Sbase => Tranformer value?,
         ),
-        :Conductor => opendss_lines(num_phases),
+        :Conductor => conductors,
         :Load => load_dicts,
         :Transformer => opendss_transformers(num_phases, Y, node_order),
         :VoltageRegulator => opendss_regulators(num_phases, Y, node_order),
+        :ShuntAdmittance => opendss_shunts(num_phases, conductors)
     )
 
     Network(net_dict)
@@ -254,9 +258,11 @@ function opendss_lines(num_phases::Int)::Vector{Dict{Symbol, Any}}
         
         # NodeOrder is list of node names in the order the nodes appear in the Y matrix.
         phases = OpenDSS.CktElement.NodeOrder()[1:OpenDSS.CktElement.NumPhases()]
-        rmatrix, xmatrix = dss_impedance_matrices_to_three_phase(
-            OpenDSS.Lines.RMatrix(), OpenDSS.Lines.XMatrix(), phases
+        rmatrix, xmatrix, cmatrix = dss_impedance_matrices_to_three_phase(
+            OpenDSS.Lines.RMatrix(), OpenDSS.Lines.XMatrix(), OpenDSS.Lines.CMatrix(), phases
         )
+        # TODO cmatrix into ShuntAdmittance? cmatrix is line capacitive susceptance -> divide it by
+        # 2 to get bus shunt susceptance.
 
         if num_phases == 1
             i = collect(phases)[1]
@@ -266,6 +272,7 @@ function opendss_lines(num_phases::Int)::Vector{Dict{Symbol, Any}}
                 :phases => phases,
                 :r1 => rmatrix[i, i],
                 :x1 => xmatrix[i, i],
+                :c1 => cmatrix[i, i],
                 :length => OpenDSS.Lines.Length()
             ))
         else
@@ -275,6 +282,7 @@ function opendss_lines(num_phases::Int)::Vector{Dict{Symbol, Any}}
                 :phases => phases,
                 :rmatrix => rmatrix,
                 :xmatrix => xmatrix,
+                :cmatrix => cmatrix,
                 :length => OpenDSS.Lines.Length()
             ))
         end
@@ -285,6 +293,35 @@ function opendss_lines(num_phases::Int)::Vector{Dict{Symbol, Any}}
 end
 
 
+"""
+
+Parse the line cmatrix values into shunt susceptance values (placing )
+"""
+function opendss_shunts(num_phases::Int, conductors::Vector{Dict{Symbol, Any}})::Vector{Dict{Symbol, Any}}
+    # for now using the line into the bus to get the shunt susceptance (not clear if we should
+    # average together the other line susceptances that are connected to the bus)
+    shunt_dicts = Dict{Symbol, Any}[]
+    f = 2π*60 * 1e-9  # nanofarads to siemens
+    for c in conductors
+        bus = c[:busses][2]
+
+        if num_phases == 1
+            push!(shunt_dicts, Dict(
+                :bus => bus,
+                :b => ismissing(c[:c1]) ? 0.0 : c[:c1] * c[:length] / 2 * f
+            ))
+        else
+            push!(shunt_dicts, Dict(
+                :bus => bus,
+                :bmatrix => ismissing(c[:cmatrix]) ? zeros(3,3) : c[:cmatrix] * c[:length] / 2 * f
+            ))
+        end
+    end
+    
+    return shunt_dicts
+end
+
+
 """
     opendss_source_bus()::String
 
@@ -345,7 +382,7 @@ function transformer_impedance(Y::Matrix{ComplexF64}, node_order::Vector{String}
     end
     # NOTE ignoring off diagonal terms since they cause numerical issues in IEEE13 source transformer
     Z = inv(Diagonal(Y_trfx)) * kV1 / kV2
-    r, x = dss_impedance_matrices_to_three_phase(abs.(real(Z)), -1*imag(Z), phases)
+    r, x, _ = dss_impedance_matrices_to_three_phase(abs.(real(Z)), -1*imag(Z), zeros(3,3), phases)
     return r, x, kV1, kV2
 end
 
diff --git a/src/network.jl b/src/network.jl
index 7b64d05..75bd221 100644
--- a/src/network.jl
+++ b/src/network.jl
@@ -123,11 +123,13 @@ Construct a `Network` from a yaml at the file path `fp`.
 """
 function Network(fp::String)
     # parse inputs
-    if endswith(lowercase(fp), "yaml") ||  endswith(lowercase(fp), "yml")
+    if endswith(lowercase(fp), ".yaml") ||  endswith(lowercase(fp), ".yml")
         d = load_yaml(fp)
+    elseif endswith(lowercase(fp), ".dss")
+        return CommonOPF.dss_to_Network(fp)
     else
         # TODO json
-        throw(error("Only parsing yaml (or yml) files so far."))
+        throw(error("Only parsing yaml (or yml) and dss files so far."))
     end
     Network(d)
 end
@@ -229,13 +231,15 @@ voltage_regulator_edges(net::AbstractNetwork) = collect(e for e in edges(net) if
 
 real_load_busses(net::Network{SinglePhase}) = collect(b for b in load_busses(net) if !ismissing(net[b][:Load].kws1))
 real_load_busses(net::Network{MultiPhase}) = collect(
-    b for b in load_busses(net) if !ismissing(net[b][:Load].kws1) || !ismissing(net[b][:Load].kws2) || !ismissing(net[b][:Load].kws3)
+    b for b in load_busses(net) 
+    if !ismissing(net[b][:Load].kws1) || !ismissing(net[b][:Load].kws2) || !ismissing(net[b][:Load].kws3)
 )
 
 
 reactive_load_busses(net::Network{SinglePhase}) = collect(b for b in load_busses(net) if !ismissing(net[b][:Load].kvars1))
 reactive_load_busses(net::Network{MultiPhase}) = collect(
-    b for b in load_busses(net) if !ismissing(net[b][:Load].kvars1) || !ismissing(net[b][:Load].kvars2) || !ismissing(net[b][:Load].kvars3)
+    b for b in load_busses(net) 
+        if !ismissing(net[b][:Load].kvars1) || !ismissing(net[b][:Load].kvars2) || !ismissing(net[b][:Load].kvars3)
 )
 
 total_load_kw(net::Network{SinglePhase}) = sum(net[load_bus][:Load].kws1 for load_bus in real_load_busses(net))
@@ -258,7 +262,9 @@ function leaf_busses(net::Network)
 end
 
 
-conductors(net::AbstractNetwork) = collect( net[ekey] for ekey in edges(net) if net[ekey] isa CommonOPF.Conductor )
+conductors(net::AbstractNetwork) = collect(
+    net[ekey] for ekey in edges(net) if net[ekey] isa CommonOPF.Conductor
+)
 
 
 function conductors_with_attribute_value(net::AbstractNetwork, attr::Symbol, val::Any)::AbstractVector{CommonOPF.Conductor}
@@ -274,8 +280,7 @@ end
 """
     fill_node_attributes!(g::MetaGraphsNext.AbstractGraph, vals::AbstractVector{<:AbstractBus})
 
-For each concrete bus in `vals` store a dict of key,val pairs for the bus attributes in a symbol key
-like :Load for CommonOPF.Load TODO just store the type
+For each concrete bus in `vals` store the concrete bus in the graph at `concrete_bus.node`.
 """
 function fill_node_attributes!(g::MetaGraphsNext.AbstractGraph, vals::AbstractVector{<:AbstractBus})
     for node in vals
@@ -285,7 +290,6 @@ function fill_node_attributes!(g::MetaGraphsNext.AbstractGraph, vals::AbstractVe
                 "You will have to manually add bus $(node.bus) if you want it in the graph."
             continue
         end
-        node_fieldnames = fieldnames(typeof(node))
         type = split(string(typeof(node)), ".")[end]  # e.g. "CommonOPF.Load" -> "Load"
         if !isempty( get(g[node.bus], Symbol(type), []) )
             @warn "Replacing existing attributes $(g[node.bus][Symbol(type)]) in node $(node.bus)"
diff --git a/test/data/case3_unbalanced.dss b/test/data/case3_unbalanced.dss
new file mode 100644
index 0000000..95f9f95
--- /dev/null
+++ b/test/data/case3_unbalanced.dss
@@ -0,0 +1,33 @@
+Clear  ! borrowed from OpenDSSDirect.jl
+New Circuit.3Bus_example
+!  define a really stiff source
+~ basekv=0.4   pu=0.9959  MVAsc1=1e6  MVAsc3=1e6 basemva=0.5
+
+
+New linecode.556MCM nphases=3 basefreq=50  ! ohms per 5 mile
+~ rmatrix = ( 0.1000 | 0.0400    0.1000 |  0.0400    0.0400    0.1000)
+~ xmatrix = ( 0.0583 |  0.0233    0.0583 | 0.0233    0.0233    0.0583)
+~ cmatrix = (50  | -0  50 | -0 -0 50  ) ! small capacitance
+
+
+New linecode.4/0QUAD nphases=3 basefreq=50  ! ohms per 100ft
+~ rmatrix = ( 0.1167 | 0.0467    0.1167 | 0.0467    0.0467    0.1167)
+~ xmatrix = (0.0667  |  0.0267    0.0667  |  0.0267    0.0267    0.0667 )
+~ cmatrix = (40  | -0  40 | -0 -0 40  )  ! small capacitance
+
+
+New Line.OHLine  bus1=sourcebus.1.2.3  Primary.1.2.3  linecode = 556MCM   length=1  ! 5 mile line
+New Line.Quad    Bus1=Primary.1.2.3  loadbus.1.2.3  linecode = 4/0QUAD  length=1   ! 100 ft
+
+
+New Load.L1 phases=1  loadbus.1.0   ( 0.4 3 sqrt / )   kW=9   kvar=3  model=1
+New Load.L2 phases=1  loadbus.2.0   ( 0.4 3 sqrt / )   kW=6   kvar=3  model=1
+New Load.L3 phases=1  loadbus.3.0   ( 0.4 3 sqrt / )   kW=6   kvar=3  model=1
+
+
+Set voltagebases=[0.4]
+Set tolerance=0.000001
+set defaultbasefreq=50
+Calcvoltagebases
+
+Solve
diff --git a/test/runtests.jl b/test/runtests.jl
index 38e9522..35c6642 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -47,6 +47,8 @@ with_logger(test_logger) do
 
     include("test_results.jl")
 
+    include("test_shunts.jl")
+
 end # outer-most testset
 end # with_logger
 
diff --git a/test/test_io.jl b/test/test_io.jl
index a268d97..6cff928 100644
--- a/test/test_io.jl
+++ b/test/test_io.jl
@@ -71,8 +71,14 @@
         @test CommonOPF.OpenDSS.CktElement.BusNames() == ["632.2.3", "645.2.3"]
         dss_rmatrix = CommonOPF.OpenDSS.Lines.RMatrix()
         dss_xmatrix = CommonOPF.OpenDSS.Lines.XMatrix()
+        dss_cmatrix = CommonOPF.OpenDSS.Lines.CMatrix()
         phases = CommonOPF.OpenDSS.CktElement.NodeOrder()[1:CommonOPF.OpenDSS.CktElement.NumPhases()]
-        r, x = CommonOPF.dss_impedance_matrices_to_three_phase(dss_rmatrix, dss_xmatrix, phases)
+        r, x = CommonOPF.dss_impedance_matrices_to_three_phase(
+            dss_rmatrix, 
+            dss_xmatrix, 
+            dss_cmatrix, 
+            phases
+        )
         @test r[2,2] == dss_rmatrix[1,1]
         @test x[3,2] == dss_xmatrix[2,1]
 
diff --git a/test/test_shunts.jl b/test/test_shunts.jl
new file mode 100644
index 0000000..dbab9c2
--- /dev/null
+++ b/test/test_shunts.jl
@@ -0,0 +1,18 @@
+@testset "shunt inputs and methods" begin
+
+    @testset "single phase" begin
+        net = Network_Papavasiliou_2018()
+        @test real(yj("1", net)) ≈ 0.0
+        @test imag(yj("1", net)) ≈ 0.0011
+
+    end
+
+    @testset "multi phase" begin
+        fp = joinpath("data", "case3_unbalanced.dss")
+        net = Network(fp)
+        f = 2π*60 * 1e-9  # nanofarads to siemens
+        @test all(imag(yj("primary", net)) .≈ f * [25.0 0 0; 0 25 0; 0 0 25])
+        @test all(imag(yj("loadbus", net)) .≈ f * [20.0 0 0; 0 20 0; 0 0 20])
+    end
+
+end