diff --git a/docs/src/matter.md b/docs/src/matter.md
index cf84fd018..088f715cf 100644
--- a/docs/src/matter.md
+++ b/docs/src/matter.md
@@ -66,9 +66,9 @@ length(coords)                 # number of particles, (5)
 
 ```jldoctest matter
 try
-  coords.x = rand(3, 5)
+    coords.x = rand(3, 5)
 catch ex
-  @error "Immutable!"
+    @error "Immutable!"
 end
 
 # output
diff --git a/src/PorousMaterials.jl b/src/PorousMaterials.jl
index e14b3f61d..294623e0c 100644
--- a/src/PorousMaterials.jl
+++ b/src/PorousMaterials.jl
@@ -1,7 +1,20 @@
 module PorousMaterials
 
-using CSV, DataFrames, Distributed, Graphs, JLD2, LinearAlgebra, OffsetArrays, Optim, 
-    Polynomials, Printf, ProgressMeter, Roots, SpecialFunctions, Statistics, StatsBase
+using CSV,
+    DataFrames,
+    Distributed,
+    Graphs,
+    JLD2,
+    LinearAlgebra,
+    OffsetArrays,
+    Optim,
+    Polynomials,
+    Printf,
+    ProgressMeter,
+    Roots,
+    SpecialFunctions,
+    Statistics,
+    StatsBase
 
 # extend Xtals
 using Reexport
@@ -34,8 +47,8 @@ function __init__()
     rc[:paths][:molecules] = ""
     rc[:paths][:grids] = ""
     rc[:paths][:simulations] = ""
-    set_paths(joinpath(pwd(), "data"), no_warn=true)
-    append_atomic_masses()
+    set_paths(joinpath(pwd(), "data"); no_warn=true)
+    return append_atomic_masses()
 end
 
 # overload unicode names w/ pure-ASCII forms
@@ -47,42 +60,82 @@ export
     fit_adsorption_isotherm,
 
     # molecule.jl
-    Molecule, translate_by!, translate_to!, random_rotation!, random_rotation_matrix, ion, distortion,
+    Molecule,
+    translate_by!,
+    translate_to!,
+    random_rotation!,
+    random_rotation_matrix,
+    ion,
+    distortion,
 
     # forcefield.jl
-    LJForceField, replication_factors, forcefield_coverage,
+    LJForceField,
+    replication_factors,
+    forcefield_coverage,
 
     # energy_utilities.jl
-    PotentialEnergy, SystemPotentialEnergy,
+    PotentialEnergy,
+    SystemPotentialEnergy,
 
     # vdw_energetics.jl
-    lennard_jones, vdw_energy, vdw_energy_no_PBC,
+    lennard_jones,
+    vdw_energy,
+    vdw_energy_no_PBC,
 
     # electrostatics.jl
-    electrostatic_potential_energy, precompute_kvec_wts,
-    setup_Ewald_sum, total, Eikr, total_electrostatic_potential_energy,
+    electrostatic_potential_energy,
+    precompute_kvec_wts,
+    setup_Ewald_sum,
+    total,
+    Eikr,
+    total_electrostatic_potential_energy,
 
     # mc_helpers.jl
-    random_insertion!, remove_molecule!, random_translation!, random_reinsertion!, needs_rotations,
+    random_insertion!,
+    remove_molecule!,
+    random_translation!,
+    random_reinsertion!,
+    needs_rotations,
     AdaptiveTranslationStepSize,
 
     # Grid.jl
     apply_periodic_boundary_condition!,
-    Grid, write_cube, read_cube, energy_grid, compute_accessibility_grid, accessible,
-    required_n_pts, xf_to_id, id_to_xf, update_density!, find_energy_minimum, origin,
- 
+    Grid,
+    write_cube,
+    read_cube,
+    energy_grid,
+    compute_accessibility_grid,
+    accessible,
+    required_n_pts,
+    xf_to_id,
+    id_to_xf,
+    update_density!,
+    find_energy_minimum,
+    origin,
+
     # EOS.jl
-    calculate_properties, PengRobinsonFluid, VdWFluid,
+    calculate_properties,
+    PengRobinsonFluid,
+    VdWFluid,
 
     # gcmc.jl
-    μVT_sim, muVT_sim, adsorption_isotherm, stepwise_adsorption_isotherm,
-    μVT_output_filename, muVT_output_filename, GCMCstats, MarkovCounts, isotherm_sim_results_to_dataframe,
- 
+    μVT_sim,
+    muVT_sim,
+    adsorption_isotherm,
+    stepwise_adsorption_isotherm,
+    μVT_output_filename,
+    muVT_output_filename,
+    GCMCstats,
+    MarkovCounts,
+    isotherm_sim_results_to_dataframe,
+
     # henry.jl
-    henry_coefficient, henry_result_savename,
+    henry_coefficient,
+    henry_result_savename,
 
     # energy_min.jl
-    find_energy_minimum, find_energy_minimum_gridsearch,
+    find_energy_minimum,
+    find_energy_minimum_gridsearch,
 
     # nvt.jl
     NVT_sim
diff --git a/src/muvt.jl b/src/muvt.jl
index f9dea5fd3..9bda91a7b 100644
--- a/src/muvt.jl
+++ b/src/muvt.jl
@@ -1,6 +1,15 @@
 module MuVT
 
-import ..Molecule, ..Frac, ..Crystal, ..LJForceField, ..@printf, ..@sprintf, ..@save, ..Cart, ..EwaldParams, ..Eikr
+import ..Molecule,
+    ..Frac,
+    ..Crystal,
+    ..LJForceField,
+    ..@printf,
+    ..@sprintf,
+    ..@save,
+    ..Cart,
+    ..EwaldParams,
+    ..Eikr
 
 ###
 #   Markov chain proposals
@@ -933,4 +942,4 @@ end
 
 export μVT_output_filename, μVT_sim
 
-end
\ No newline at end of file
+end
diff --git a/src/nvt.jl b/src/nvt.jl
index 3913abfba..d76a95933 100644
--- a/src/nvt.jl
+++ b/src/nvt.jl
@@ -1,31 +1,39 @@
 module NVT
 
-import ..Molecule, ..Frac, ..Crystal, ..LJForceField, ..@printf, ..@sprintf, ..@save, ..Cart, ..EwaldParams, ..Eikr
+import ..Molecule,
+    ..Frac,
+    ..Crystal,
+    ..LJForceField,
+    ..@printf,
+    ..@sprintf,
+    ..@save,
+    ..Cart,
+    ..EwaldParams,
+    ..Eikr
 
 const KB = 1.38064852e7 # Boltmann constant (Pa-m3/K --> Pa-A3/K)
 
 ###
 #   Markov chain proposals
 ###
-const PROPOSAL_ENCODINGS = Dict(1 => "translation",
-                                2 => "rotation",
-                                3 => "reinsertion"
-                                ) # helps with printing later
+const PROPOSAL_ENCODINGS = Dict(1 => "translation", 2 => "rotation", 3 => "reinsertion") # helps with printing later
 const N_PROPOSAL_TYPES = length(keys(PROPOSAL_ENCODINGS))
 # each proposal type gets an Int for clearer code
 const TRANSLATION = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["translation"]
-const ROTATION    = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["rotation"]
+const ROTATION = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["rotation"]
 const REINSERTION = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["reinsertion"]
 
 # TODO move this to MC helpers? but not sure if it will inline. so wait after test with @time
 # potential energy change after inserting/deleting/perturbing coordinates of molecules[molecule_id]
-@inline function potential_energy(molecule_id::Int,
-                                  molecules::Array{Molecule{Frac}, 1},
-                                  xtal::Crystal,
-                                  ljff::LJForceField,
-                                  eparams::EwaldParams,
-                                  eikr_gh::Eikr,
-                                  eikr_gg::Eikr)
+@inline function potential_energy(
+    molecule_id::Int,
+    molecules::Array{Molecule{Frac}, 1},
+    xtal::Crystal,
+    ljff::LJForceField,
+    eparams::EwaldParams,
+    eikr_gh::Eikr,
+    eikr_gg::Eikr
+)
     energy = SystemPotentialEnergy()
     # van der Waals interactions
     energy.gg.vdw = vdw_energy(molecule_id, molecules, ljff, xtal.box) # guest-guest
@@ -33,25 +41,41 @@ const REINSERTION = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["reinsertion
     # electrostatic interactions
     if has_charges(molecules[molecule_id])
         # guest-guest
-        energy.gg.es = total(electrostatic_potential_energy(molecules, molecule_id, eparams, xtal.box, eikr_gg))
+        energy.gg.es = total(
+            electrostatic_potential_energy(
+                molecules,
+                molecule_id,
+                eparams,
+                xtal.box,
+                eikr_gg
+            )
+        )
         # guest-host
         if has_charges(xtal)
-            energy.gh.es = total(electrostatic_potential_energy(xtal, molecules[molecule_id], eparams, eikr_gh))
+            energy.gh.es = total(
+                electrostatic_potential_energy(
+                    xtal,
+                    molecules[molecule_id],
+                    eparams,
+                    eikr_gh
+                )
+            )
         end
     end
     return energy
 end
 
-function NVT_sim(xtal::Crystal,
-                 molecule_template::Molecule{Cart},
-                 molecules::Array{Molecule{Cart}, 1},
-                 temperature::Float64,
-                 ljff::LJForceField;
-                 n_burn_cycles::Int=5000,
-                 n_sample_cycles::Int=5000,
-                 verbose=true,
-                 ewald_precision::Float64=1e-6
-                )
+function NVT_sim(
+    xtal::Crystal,
+    molecule_template::Molecule{Cart},
+    molecules::Array{Molecule{Cart}, 1},
+    temperature::Float64,
+    ljff::LJForceField;
+    n_burn_cycles::Int=5000,
+    n_sample_cycles::Int=5000,
+    verbose=true,
+    ewald_precision::Float64=1e-6
+)
     assert_P1_symmetry(xtal)
 
     start_time = time()
@@ -63,21 +87,21 @@ function NVT_sim(xtal::Crystal,
     #  NON-RIGOROUS PLACEHOLDER -- DELETE
     ###
     molecule = deepcopy(molecule_template)
-    
+
     ###
     #   replicate xtal so that nearest image convention can be applied for short-range interactions
     ###
     repfactors = replication_factors(xtal.box, ljff)
     xtal = replicate(xtal, repfactors) # frac coords still in [0, 1]
-    
+
     # convert molecules array to fractional using this box.
     molecules = Frac.(molecules, xtal.box)
 
-    if ! neutral(molecule.charges)
+    if !neutral(molecule.charges)
         error(@sprintf("Molecule %s is not charge neutral!\n", molecule.species))
     end
 
-    if ! (forcefield_coverage(xtal.atoms, ljff) & forcefield_coverage(molecule.atoms, ljff))
+    if !(forcefield_coverage(xtal.atoms, ljff) & forcefield_coverage(molecule.atoms, ljff))
         error("Missing atoms from forcefield.")
     end
 
@@ -85,9 +109,12 @@ function NVT_sim(xtal::Crystal,
     # pre-compute weights on k-vector contributions to long-rage interactions in
     #   Ewald summation for electrostatics
     #   allocate memory for exp^{i * n * k ⋅ r}
-    eparams = setup_Ewald_sum(xtal.box, sqrt(ljff.r²_cutoff),
-                        verbose=verbose & (has_charges(molecule) || has_charges(xtal)),
-                        ϵ=ewald_precision)
+    eparams = setup_Ewald_sum(
+        xtal.box,
+        sqrt(ljff.r²_cutoff);
+        verbose=verbose & (has_charges(molecule) || has_charges(xtal)),
+        ϵ=ewald_precision
+    )
     eikr_gh = Eikr(xtal, eparams)
     eikr_gg = Eikr(molecule, eparams)
 
@@ -103,19 +130,21 @@ function NVT_sim(xtal::Crystal,
             # assert that the molecules are inside the simulation box
             @assert inside(m) "initializing with molecules outside simulation box!"
             # ensure pair-wise bond distances match template
-            @assert ! distortion(m, Frac(molecule, xtal.box), xtal.box)
+            @assert !distortion(m, Frac(molecule, xtal.box), xtal.box)
         end
 
         system_energy.gh.vdw = total_vdw_energy(xtal, molecules, ljff)
         system_energy.gg.vdw = total_vdw_energy(molecules, ljff, xtal.box)
-        system_energy.gh.es  = total(total_electrostatic_potential_energy(xtal, molecules, eparams, eikr_gh))
-        system_energy.gg.es  = total(electrostatic_potential_energy(molecules, eparams, xtal.box, eikr_gg))
+        system_energy.gh.es =
+            total(total_electrostatic_potential_energy(xtal, molecules, eparams, eikr_gh))
+        system_energy.gg.es =
+            total(electrostatic_potential_energy(molecules, eparams, xtal.box, eikr_gg))
     end
 
     ####
     #   proposal probabilities
     ###
-    mc_proposal_probabilities = [0.0 for p = 1:N_PROPOSAL_TYPES]
+    mc_proposal_probabilities = [0.0 for p in 1:N_PROPOSAL_TYPES]
     # set defaults
     mc_proposal_probabilities[REINSERTION] = 0.05
     if needs_rotations(molecule)
@@ -130,19 +159,26 @@ function NVT_sim(xtal::Crystal,
     mc_proposal_probabilities = ProbabilityWeights(mc_proposal_probabilities)
     if verbose
         println("\tMarkov chain proposals:")
-        for p = 1:N_PROPOSAL_TYPES
-            @printf("\t\tprobability of %s: %f\n", PROPOSAL_ENCODINGS[p], mc_proposal_probabilities[p])
+        for p in 1:N_PROPOSAL_TYPES
+            @printf(
+                "\t\tprobability of %s: %f\n",
+                PROPOSAL_ENCODINGS[p],
+                mc_proposal_probabilities[p]
+            )
         end
     end
 
-    markov_counts = MarkovCounts(zeros(Int, length(PROPOSAL_ENCODINGS)), zeros(Int, length(PROPOSAL_ENCODINGS)))
+    markov_counts = MarkovCounts(
+        zeros(Int, length(PROPOSAL_ENCODINGS)),
+        zeros(Int, length(PROPOSAL_ENCODINGS))
+    )
 
     # (n_burn_cycles + n_sample_cycles) is number of outer cycles.
     #   for each outer cycle, peform max(20, # molecules in the system) MC proposals.
     markov_chain_time = 0
     outer_cycle_start = 1
-    for outer_cycle = outer_cycle_start:(n_burn_cycles + n_sample_cycles)
-        for inner_cycle = 1:max(20, length(molecules))
+    for outer_cycle in outer_cycle_start:(n_burn_cycles + n_sample_cycles)
+        for inner_cycle in 1:max(20, length(molecules))
             markov_chain_time += 1
 
             # choose proposed move randomly; keep track of proposals
@@ -154,14 +190,28 @@ function NVT_sim(xtal::Crystal,
                 molecule_id = rand(1:length(molecules))
 
                 # energy of the molecule before it was translated
-                energy_old = potential_energy(molecule_id, molecules, xtal, ljff,
-                                              eparams, eikr_gh, eikr_gg)
+                energy_old = potential_energy(
+                    molecule_id,
+                    molecules,
+                    xtal,
+                    ljff,
+                    eparams,
+                    eikr_gh,
+                    eikr_gg
+                )
 
                 old_molecule = random_translation!(molecules[molecule_id], xtal.box)
 
                 # energy of the molecule after it is translated
-                energy_new = potential_energy(molecule_id, molecules, xtal, ljff,
-                                              eparams, eikr_gh, eikr_gg)
+                energy_new = potential_energy(
+                    molecule_id,
+                    molecules,
+                    xtal,
+                    ljff,
+                    eparams,
+                    eikr_gh,
+                    eikr_gg
+                )
 
                 # Metropolis Hastings Acceptance for translation
                 if rand() < exp(-(sum(energy_new) - sum(energy_old)) / temperature)
@@ -178,8 +228,15 @@ function NVT_sim(xtal::Crystal,
                 molecule_id = rand(1:length(molecules))
 
                 # energy of the molecule before we rotate it
-                energy_old = potential_energy(molecule_id, molecules, xtal, ljff,
-                                              eparams, eikr_gh, eikr_gg)
+                energy_old = potential_energy(
+                    molecule_id,
+                    molecules,
+                    xtal,
+                    ljff,
+                    eparams,
+                    eikr_gh,
+                    eikr_gg
+                )
 
                 # store old molecule to restore old position in case move is rejected
                 old_molecule = deepcopy(molecules[molecule_id])
@@ -188,8 +245,15 @@ function NVT_sim(xtal::Crystal,
                 random_rotation!(molecules[molecule_id], xtal.box)
 
                 # energy of the molecule after it is translated
-                energy_new = potential_energy(molecule_id, molecules, xtal, ljff,
-                                              eparams, eikr_gh, eikr_gg)
+                energy_new = potential_energy(
+                    molecule_id,
+                    molecules,
+                    xtal,
+                    ljff,
+                    eparams,
+                    eikr_gh,
+                    eikr_gg
+                )
 
                 # Metropolis Hastings Acceptance for rotation
                 if rand() < exp(-(sum(energy_new) - sum(energy_old)) / temperature)
@@ -206,15 +270,29 @@ function NVT_sim(xtal::Crystal,
                 molecule_id = rand(1:length(molecules))
 
                 # compute the potential energy of the molecule we propose to re-insert
-                energy_old = potential_energy(molecule_id, molecules, xtal, ljff,
-                                             eparams, eikr_gh, eikr_gg)
+                energy_old = potential_energy(
+                    molecule_id,
+                    molecules,
+                    xtal,
+                    ljff,
+                    eparams,
+                    eikr_gh,
+                    eikr_gg
+                )
 
                 # reinsert molecule; store old configuration of the molecule in case proposal is rejected
                 old_molecule = random_reinsertion!(molecules[molecule_id], xtal.box)
 
                 # compute the potential energy of the molecule in its new configuraiton
-                energy_new = potential_energy(molecule_id, molecules, xtal, ljff,
-                                              eparams, eikr_gh, eikr_gg)
+                energy_new = potential_energy(
+                    molecule_id,
+                    molecules,
+                    xtal,
+                    ljff,
+                    eparams,
+                    eikr_gh,
+                    eikr_gg
+                )
 
                 # Metropolis Hastings Acceptance for reinsertion
                 if rand() < exp(-(sum(energy_new) - sum(energy_old)) / temperature)
@@ -233,20 +311,20 @@ function NVT_sim(xtal::Crystal,
 
     for m in molecules
         @assert inside(m) "molecule outside box!"
-        @assert ! distortion(m, Frac(molecule, xtal.box), xtal.box, atol=1e-10) "molecules distorted"
+        @assert !distortion(m, Frac(molecule, xtal.box), xtal.box; atol=1e-10) "molecules distorted"
     end
 
     # compute total energy, compare to `current_energy*` variables where were incremented
     system_energy_end = SystemPotentialEnergy()
     system_energy_end.gh.vdw = total_vdw_energy(xtal, molecules, ljff)
     system_energy_end.gg.vdw = total_vdw_energy(molecules, ljff, xtal.box)
-    system_energy_end.gh.es = total(total_electrostatic_potential_energy(xtal, molecules,
-                                                 eparams, eikr_gh))
-    system_energy_end.gg.es = total(total_electrostatic_potential_energy(molecules,
-                                        eparams, xtal.box, eikr_gg))
+    system_energy_end.gh.es =
+        total(total_electrostatic_potential_energy(xtal, molecules, eparams, eikr_gh))
+    system_energy_end.gg.es =
+        total(total_electrostatic_potential_energy(molecules, eparams, xtal.box, eikr_gg))
 
     # see Energetics_Util.jl for this function, overloaded isapprox to print mismatch
-    if ! isapprox(system_energy, system_energy_end, verbose=true, atol=0.01)
+    if !isapprox(system_energy, system_energy_end; verbose=true, atol=0.01)
         error("energy incremented improperly during simulation...")
     end
 
@@ -260,4 +338,4 @@ end # NVT_sim
 
 export NVT_sim
 
-end # module NVT
\ No newline at end of file
+end # module NVT