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/Symmetry/AllowedAnisotropy.jl b/src/Symmetry/AllowedAnisotropy.jl
index b445cb721..d9f6f3eda 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.")
+        @warn "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.")
+        @warn "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..7894e3a34 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.""")
+        @warn """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
 
@@ -370,8 +369,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 +491,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.")
+        @warn "Atoms are being renumbered."
     end
 
     ret = Crystal(cryst.latvecs, cryst.prim_latvecs, new_positions, new_types, new_classes, new_sitesyms,
diff --git a/src/Symmetry/Printing.jl b/src/Symmetry/Printing.jl
index f04df3ff0..4c2b5607d 100644
--- a/src/Symmetry/Printing.jl
+++ b/src/Symmetry/Printing.jl
@@ -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.""")
+                    @warn """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/OnsiteCoupling.jl b/src/System/OnsiteCoupling.jl
index 087559e1f..ea70989b9 100644
--- a/src/System/OnsiteCoupling.jl
+++ b/src/System/OnsiteCoupling.jl
@@ -152,7 +152,7 @@ function set_onsite_coupling!(sys::System{N}, op::Matrix{ComplexF64}, i::Int) wh
     (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.")
+        @info "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..4695b5b78 100644
--- a/src/System/PairExchange.jl
+++ b/src/System/PairExchange.jl
@@ -99,7 +99,7 @@ function set_exchange!(sys::System{N}, J, bond::Bond; biquad=0.) where N
 
     # 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.")
+        @info "Overriding exchange for bond $bond."
     end
 
     for i in 1:natoms(sys.crystal)