diff --git a/Project.toml b/Project.toml
index 3a0790f..e03176e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,6 +13,7 @@ MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 TravelingSalesmanHeuristics = "8c8f4381-2cdd-507c-846c-be2bcff6f45f"
 
 [compat]
@@ -22,6 +23,7 @@ MIToS = "2"
 MultivariateStats = "0.9, 0.10"
 Rotations = "1"
 StaticArrays = "1"
+StructArrays = "0.6"
 TravelingSalesmanHeuristics = "0.3"
 julia = "1.6"
 
diff --git a/src/GPCRAnalysis.jl b/src/GPCRAnalysis.jl
index 56ae8e1..1c635a9 100644
--- a/src/GPCRAnalysis.jl
+++ b/src/GPCRAnalysis.jl
@@ -13,6 +13,7 @@ using MultivariateStats
 using Distances
 using Hungarian
 using StaticArrays
+using StructArrays
 using TravelingSalesmanHeuristics
 
 export @res_str
@@ -21,7 +22,8 @@ export SequenceMapping
 export species, uniprotX
 export try_download_alphafold, download_alphafolds, getchain, findall_subseq
 export filter_species!, filter_long!, sortperm_msa, chimerax_script
-export project_sequences, columnwise_entropy, align, residue_centroid, residue_centroid_matrix, mapclosest, chargelocations, positive_locations, negative_locations
+export project_sequences, columnwise_entropy, align, residue_centroid, residue_centroid_matrix, mapclosest
+export chargelocations, positive_locations, negative_locations, charge_magnitudes, fields
 export StructAlign, residueindex, ismapped
 export BWScheme, lookupbw
 export aa_properties, aa_properties_zscored
diff --git a/src/analyze.jl b/src/analyze.jl
index d832859..24fc32d 100644
--- a/src/analyze.jl
+++ b/src/analyze.jl
@@ -146,3 +146,93 @@ end
 
 positive_locations(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] for loc in chargelocs if loc[3] ∈ ("ARG", "LYS", "HIS")]
 negative_locations(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] for loc in chargelocs if loc[3] ∈ ("ASP", "GLU")]
+
+_charge_magnitude(aa) = aa == "LYS" ? 1.0 :
+                        aa ∈ ("ARG", "HIS") ? 0.5 :
+                        aa ∈ ("ASP", "GLU") ? -0.5 : 0.0
+charge_magnitudes(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] => _charge_magnitude(loc[3]) for loc in chargelocs]
+
+struct Fields{T} <: FieldVector{2, T}
+    electrostatic::T
+    steric::T
+end
+
+gvwr(resname, a::PDBAtom) = get(vanderwaalsradius,
+                                (resname, a.atom),
+                                get(vanderwaalsradius, (resname, a.element), nothing))
+
+# Repulsive "meta-potentials" to estimate steric interactions
+gaussian(x::Coordinates, y::Coordinates, σ) = exp(-distance(x, y)^2 / (2 * σ^2))
+
+function lennardjones(x::Coordinates, y::Coordinates, σ)
+    # Lennard-Jones potential
+    r = distance(x, y)
+    k = σ / r
+    k6 = k^6
+    return iszero(r) ? Inf : k6^2 - k6
+end
+
+function steric(f::F, x::Coordinates, pdb::AbstractVector{PDBResidue}, σfac=1) where F
+    s = 0.0
+    for aa in pdb
+        resname = aa.id.name
+        for a in aa.atoms
+            vdwr = gvwr(resname, a)
+            vdwr === nothing && continue
+            s += f(x, a.coordinates, vdwr * σfac)
+        end
+    end
+    return s
+end
+function steric!(f::F, ss::AbstractArray, xs::AbstractArray{Coordinates}, pdb::AbstractVector{PDBResidue}, σfac=1) where F
+    # Performance optimization of `steric`
+    axes(ss) == axes(xs) || throw(DimensionMismatch("ss and xs must have commensurate axes, got $(axes(ss)) and $(axes(xs))"))
+    fill!(ss, 0.0)
+    for aa in pdb
+        resname = aa.id.name
+        for a in aa.atoms
+            vdwr = gvwr(resname, a)
+            vdwr === nothing && continue
+            for (i, x) in pairs(xs)
+                ss[i] += f(x, a.coordinates, vdwr * σfac)
+            end
+        end
+    end
+    return ss
+end
+
+const membrane_width = 34  # Å, hydrophobic thickness (see https://pubs.rsc.org/en/content/articlehtml/2016/sm/c6sm01777k)
+const λB1 = 568   # Bjerrum length for ϵr = 1, in Å
+
+function dielectric(density::Real, location::Coordinates)
+    extramembrane = 1 / (1 + exp(abs(2*location[3]/membrane_width)^8))
+    # @assert 0.0 <= extramembrane <= 1.0
+    ϵext = 30 * extramembrane + 2.2 * (1-extramembrane)
+    return 6.5 * density + ϵext * (1 - density)
+end
+
+function fields(pdb::AbstractVector{PDBResidue}, locations, r0=0.0)
+    sterics = similar(locations, Float64)
+    steric!(lennardjones, sterics, locations, pdb)
+    # Get the local dielectric constant: see
+    #   On the Dielectric "Constant" of Proteins: Smooth Dielectric Function for Macromolecular Modeling and Its Implementation in DelPhi,
+    #   Lin Li, Chuan Li, Zhe Zhang, Emil Alexov
+    #   J Chem Theory Comput. 2013 Apr 9;9(4):2126-2136. doi: 10.1021/ct400065j
+    # https://pubmed.ncbi.nlm.nih.gov/23585741/
+    # We use the "density" as an indicator of whether we are interior or exterior
+    density = similar(locations, Float64)
+    steric!(gaussian, density, locations, pdb)
+    dmin, dmax = extrema(density)
+    dielec = dielectric.((density .- dmin) ./ (dmax - dmin), locations)
+    cmags = charge_magnitudes(chargelocations(pdb))
+    electrostatic = similar(locations, Float64)
+    fill!(electrostatic, 0.0)
+    for (i, x) in pairs(locations)
+        λ = λB1 / dielec[i]
+        for (y, v) in cmags
+            r = distance(x, y) + r0
+            electrostatic[i] += v/r * exp(-r/λ)
+        end
+    end
+    return StructArray{Fields{Float64}}((; electrostatic, steric=sterics))
+end