diff --git a/src/Operators/Symbolic.jl b/src/Operators/Symbolic.jl
index 2ec8eb841..8ec403757 100644
--- a/src/Operators/Symbolic.jl
+++ b/src/Operators/Symbolic.jl
@@ -67,7 +67,7 @@ function operator_to_matrix(p::DP.AbstractPolynomialLike; N)
         [stevens_operator_symbols[k] => stevens_matrices(k; N) for k=1:6]... 
     )
     if !(rep ≈ rep')
-        println("Warning: Symmetrizing non-Hermitian operator '$p'.")
+        @warn "Symmetrizing non-Hermitian operator '$p'."
     end
     # Symmetrize in any case for slightly more accuracy
     return (rep+rep')/2
diff --git a/src/Reshaping.jl b/src/Reshaping.jl
index fe4bbbb8b..1c3ed13c7 100644
--- a/src/Reshaping.jl
+++ b/src/Reshaping.jl
@@ -339,9 +339,9 @@ function suggest_magnetic_supercell(qs, latsize)
 
     qstrs = join(map(wavevec_str, qs), ", ")
 
-    println("Suggested magnetic supercell in multiples of lattice vectors:")
-    println()
-    println("    $best_A")
-    println()
-    println("for wavevectors [$qstrs].")
+    println("""Suggested magnetic supercell in multiples of lattice vectors:
+               
+                   $best_A
+               
+               for wavevectors [$qstrs].""")
 end
diff --git a/src/SunnyGfx/SunnyGfx.jl b/src/SunnyGfx/SunnyGfx.jl
index 58a5da558..db232ac73 100755
--- a/src/SunnyGfx/SunnyGfx.jl
+++ b/src/SunnyGfx/SunnyGfx.jl
@@ -23,14 +23,14 @@ directory `dir` is provided, an HTML file will be written at that location.
 """
 function browser(html_str::String; dir=Nothing)
     if Sys.islinux() && dir==Nothing
-        println("""If your browser cannot access `/tmp`, consider specifying a directory:
-                       browser(html_str; dir="~")""")
+        @info """If your browser cannot access `/tmp`, consider specifying a directory:
+                     browser(html_str; dir="~")"""
     end
 
     is_tempdir = (dir == Nothing)
     dir = is_tempdir ?  mktempdir() : expanduser(dir)
     path = joinpath(dir, "SunnyGfx.html")
-    !is_tempdir && println("Writing `$path`.")
+    !is_tempdir && @info "Writing `$path`."
     open(path, "w") do io
         write(io, html_str)
     end
diff --git a/src/Symmetry/AllowedAnisotropy.jl b/src/Symmetry/AllowedAnisotropy.jl
index b445cb721..28822df3c 100644
--- a/src/Symmetry/AllowedAnisotropy.jl
+++ b/src/Symmetry/AllowedAnisotropy.jl
@@ -135,7 +135,7 @@ function suggest_frame_for_atom(cryst::Crystal, i::Int)
     end
 
     if isempty(axes_counts)
-        println("Warning: Could not find a symmetry axis.")
+        @info "Could not find a symmetry axis."
         return Mat3(I)
     end
 
@@ -170,7 +170,7 @@ function suggest_frame_for_atom(cryst::Crystal, i::Int)
     orthogonal_axes_counts = filter(x -> abs(x[1]⋅z_dir) < 1e-12, axes_counts)
 
     if isempty(orthogonal_axes_counts)
-        println("Warning: Could not find a symmetry axis orthogonal to $z_dir.")
+        @info "Could not find a symmetry axis orthogonal to $z_dir."
         x_dir = (z_dir ≈ Vec3(1,0,0)) ? Vec3(0,0,1) : Vec3(1,0,0)
         x_dir = normalize(x_dir - (x_dir⋅z_dir)*z_dir)
     else
diff --git a/src/Symmetry/Crystal.jl b/src/Symmetry/Crystal.jl
index 4cec9700d..c4f9f1615 100644
--- a/src/Symmetry/Crystal.jl
+++ b/src/Symmetry/Crystal.jl
@@ -134,14 +134,13 @@ function Crystal(latvecs, positions, spacegroup_number::Int; types::Union{Nothin
 end
 
 function print_crystal_warnings(latvecs, positions)
-    det(latvecs) < 0 && println("Warning: Lattice vectors are not right-handed.")
+    det(latvecs) < 0 && @warn "Lattice vectors are not right-handed."
     if length(positions) >= 100
-        println("""
-            Warning: This a very large crystallographic cell, which Sunny does not handle well.
-            If the intention is to model chemical inhomogeneity, the recommended procedure is as
-            follows: First, create a small unit cell with an idealized structure. Next, create
-            a perfectly periodic `System` of the desired size. Finally, use `to_inhomogeneous`
-            and related functions to design a system with the desired inhomogeneities.""")
+        @info """This a very large crystallographic cell, which Sunny does not handle well.
+                 If the intention is to model chemical inhomogeneity, the recommended procedure is as
+                 follows: First, create a small unit cell with an idealized structure. Next, create
+                 a perfectly periodic `System` of the desired size. Finally, use `to_inhomogeneous`
+                 and related functions to design a system with the desired inhomogeneities."""
     end
 end
 
@@ -309,6 +308,7 @@ function crystal_from_symbol(latvecs::Mat3, positions::Vector{Vec3}, types::Vect
             end
         end
 
+        # TODO: Make this @warn
         println("The spacegroup '$symbol' allows for multiple settings!")
         println("Returning a list of the possible crystals:")
         for (i, (hall_number, c)) in enumerate(zip(hall_numbers, crysts))
@@ -370,8 +370,8 @@ function crystal_from_symops(latvecs::Mat3, positions::Vector{Vec3}, types::Vect
     end
 
     if !is_subgroup
-        println("""Warning: User provided symmetry operation could not be inferred by Spglib,
-                   which likely indicates a non-conventional unit cell.""")
+        @warn """User provided symmetry operation could not be inferred by Spglib,
+                 which likely indicates a non-conventional unit cell."""
     end
 
     # If the inferred symops match the provided ones, then we use the inferred
@@ -492,7 +492,7 @@ function subcrystal(cryst::Crystal, classes::Vararg{Int, N}) where N
     new_sitesyms = cryst.sitesyms[atoms]
 
     if atoms != 1:maximum(atoms)
-        println("Warning: atoms are being renumbered.")
+        @info "Atoms have been renumbered in subcrystal."
     end
 
     ret = Crystal(cryst.latvecs, cryst.prim_latvecs, new_positions, new_types, new_classes, new_sitesyms,
diff --git a/src/Symmetry/Parsing.jl b/src/Symmetry/Parsing.jl
index ddb175a90..adb63d830 100644
--- a/src/Symmetry/Parsing.jl
+++ b/src/Symmetry/Parsing.jl
@@ -108,16 +108,17 @@ function Crystal(filename::AbstractString; symprec=nothing)
         end
         (err, i) = findmax(errs)
         if err < 1e-12
-            println("Precision parameter is unspecified, but all coordinates seem to be simple fractions.")
-            println("Setting symprec=1e-12.")
+            @info """Precision parameter is unspecified, but all coordinates seem to be simple fractions.
+                     Setting symprec=1e-12."""
             symprec = 1e-12
         elseif 1e-12 < err < 1e-4
-            @printf "Precision parameter is unspecified, but coordinate string '%s' seems to have error %.1e.\n" strs[i] err
             symprec = 15err
-            @printf "Setting symprec=%.1e.\n" symprec
+            err_str = @sprintf "%.1e" err
+            symprec_str = @sprintf "%.1e" symprec
+            @info """Precision parameter is unspecified, but coordinate string '$s' seems to have error $err_str.
+                     Setting symprec=$symprec_str."""
         else
-            println("Error: Cannot infer precision. Please provide an explicit `symprec` parameter to load '$filename'")
-            return Nothing
+            error("Cannot infer precision. Please provide an explicit `symprec` parameter to load '$filename'")
         end
     end
 
diff --git a/src/Symmetry/Printing.jl b/src/Symmetry/Printing.jl
index f04df3ff0..6536e3a6c 100644
--- a/src/Symmetry/Printing.jl
+++ b/src/Symmetry/Printing.jl
@@ -216,9 +216,9 @@ function print_suggested_frame(cryst::Crystal, i::Int)
     R_strs = [number_to_math_string(x; digits=14, atol=1e-12) for x in R]
     R_strs = _add_padding_to_coefficients(R_strs)
 
-    println("R = [" * join(R_strs[1,:], " "))
-    println("     " * join(R_strs[2,:], " "))
-    println("     " * join(R_strs[3,:], " "), " ]")
+    println("""R = [$(join(R_strs[1,:], " "))
+                    $(join(R_strs[2,:], " "))
+                    $(join(R_strs[3,:], " "))]""")
 end
 
 
@@ -282,8 +282,8 @@ function print_allowed_anisotropy(cryst::Crystal, i::Int; R::Mat3, atol, digits,
             for b in reverse(collect(eachcol(B)))
 
                 if any(x -> 1e-12 < abs(x) < 1e-6, b)
-                    println("""Warning: Found a very small but nonzero expansion coefficient.
-                               This may indicate a slightly misaligned reference frame.""")
+                    @info """Found a very small but nonzero expansion coefficient.
+                             This may indicate a slightly misaligned reference frame."""
                 end
 
                 # rescale column by its minimum nonzero value
diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl
index 383e64c62..47ef6f097 100644
--- a/src/System/Interactions.jl
+++ b/src/System/Interactions.jl
@@ -4,6 +4,17 @@ function empty_interactions(na, N)
     end
 end
 
+# Warn up to `OverrideWarningMax` times about overriding a coupling
+OverrideWarningCnt::Int = 0
+OverrideWarningMax::Int = 5
+function warn_coupling_override(str)
+    global OverrideWarningCnt, OverrideWarningMax
+    OverrideWarningCnt < OverrideWarningMax && @info str
+    OverrideWarningCnt += 1
+    OverrideWarningCnt == OverrideWarningMax && @info "Suppressing future override notifications."
+end
+
+
 # Creates a clone of the lists of exchange interactions, which can be mutably
 # updated.
 function clone_interactions(ints::Interactions)
diff --git a/src/System/OnsiteCoupling.jl b/src/System/OnsiteCoupling.jl
index 087559e1f..54123ab77 100644
--- a/src/System/OnsiteCoupling.jl
+++ b/src/System/OnsiteCoupling.jl
@@ -144,15 +144,15 @@ function set_onsite_coupling!(sys::System{N}, op::Matrix{ComplexF64}, i::Int) wh
     iszero(op) && return 
 
     if !is_anisotropy_valid(sys.crystal, i, op)
-        println("Symmetry-violating anisotropy: $op.")
-        println("Use `print_site(crystal, $i)` for more information.")
+        @error """Symmetry-violating anisotropy: $op.
+                  Use `print_site(crystal, $i)` for more information."""
         error("Invalid anisotropy.")
     end
 
     (1 <= i <= natoms(sys.crystal)) || error("Atom index $i is out of range.")
 
     if !iszero(ints[i].onsite)
-        println("Warning: Overriding anisotropy for atom $i.")
+        warn_coupling_override("Overriding anisotropy for atom $i.")
     end
 
     onsite = OnsiteCoupling(sys, op, sys.Ns[1,1,1,i])
diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl
index 9260fc260..6a978da2f 100644
--- a/src/System/PairExchange.jl
+++ b/src/System/PairExchange.jl
@@ -92,14 +92,14 @@ function set_exchange!(sys::System{N}, J, bond::Bond; biquad=0.) where N
 
     # Verify that exchange is symmetry-consistent
     if !is_coupling_valid(sys.crystal, bond, J)
-        println("Symmetry-violating exchange: $J.")
-        println("Use `print_bond(crystal, $bond)` for more information.")
+        @error """Symmetry-violating exchange: $J.
+                  Use `print_bond(crystal, $bond)` for more information."""
         error("Interaction violates symmetry.")
     end
 
     # Print a warning if an interaction already exists for bond
     if any(x -> x.bond == bond, ints[bond.i].pair)
-        println("Warning: Overriding exchange for bond $bond.")
+        warn_coupling_override("Overriding coupling for $bond.")
     end
 
     for i in 1:natoms(sys.crystal)
@@ -124,8 +124,8 @@ function sites_to_internal_bond(sys::System{N}, site1::CartesianIndex{4}, site2:
         else
             cell1 = Tuple(to_cell(site1))
             cell2 = Tuple(to_cell(site2))
-            println("""Cells $cell1 and $cell2 are not compatible with the offset
-                       $n_ref for a system with lattice size $latsize.""")
+            @error """Cells $cell1 and $cell2 are not compatible with the offset
+                      $n_ref for a system with lattice size $latsize."""
             error("Incompatible displacement specified")
         end
     end
@@ -146,11 +146,10 @@ function sites_to_internal_bond(sys::System{N}, site1::CartesianIndex{4}, site2:
     if safety * distances[perm[1]] < distances[perm[2]] - 1e-12
         return bonds[perm[1]]
     else
-        println(distances[perm])
         n1 = bonds[perm[1]].n
         n2 = bonds[perm[2]].n
-        println("""Cannot find an obvious offset vector. Possibilities include $n1 and $n2.
-                   Try using a bigger system size, or pass an explicit offset vector.""")
+        @error """Cannot find an obvious offset vector. Possibilities include $n1 and $n2.
+                  Try using a bigger system size, or pass an explicit offset vector."""
         error("Ambiguous offset between sites.")
     end
 end