diff --git a/Project.toml b/Project.toml
index 2f6a447..495efb3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "VibrationalAnalysis"
 uuid = "04082e5c-1d57-4577-9a0f-61e6cd56dade"
 authors = ["Josef M. Gallmetzer"]
-version = "0.1.3"
+version = "0.1.4"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
diff --git a/README.md b/README.md
index b453c36..e728418 100644
--- a/README.md
+++ b/README.md
@@ -2,3 +2,22 @@
 [![Build CI](https://github.com/galjos/VibrationalAnalysis.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/galjos/VibrationalAnalysis.jl/actions/workflows/CI.yml)
 
 # VibrationalAnalysis.jl
+
+This package provides tools to perform [vibrational analysis](https://gaussian.com/vib/) of molecules. It is based on the output of the `QMCFC` molecular dynamics code developed in the [Hofer Lab](https://www.uibk.ac.at/en/aatc/ag-hofer/) of the [University of Innsbruck](https://www.uibk.ac.at/).
+
+## Installation
+```julia-repl
+julia> ] add VibrationalAnalysis
+```
+
+## Usage
+```julia-repl
+julia> using VibrationalAnalysis
+julia> calculate("restart.rst", "hessian.dat", "moldescriptor.dat")
+```
+
+## Documentation
+```julia-repl
+julia> using VibrationalAnalysis
+julia> ?VibrationalAnalysis
+```
diff --git a/data/test_hessian_h2o_kcal.dat b/data/test_hessian_h2o_kcal.dat
deleted file mode 100644
index 277361e..0000000
--- a/data/test_hessian_h2o_kcal.dat
+++ /dev/null
@@ -1,9 +0,0 @@
-	-1234.039840683900	272.742876207450	  0.000000000000	1101.440992683300	-87.766482380400	  0.000000000000	132.653144266450	-184.938341616250	  0.000000000000
-	272.742891449950	-1234.045484438450	  0.000000000000	-184.947127076600	132.655686883300	  0.000000000000	-87.757713722400	1101.444092755400	  0.000000000000
-	  0.000000000000	  0.000000000000	 -0.144077553600	  0.000000000000	  0.000000000000	  0.057823776500	  0.000000000000	  0.000000000000	  0.057753440200
-	1101.409093298700	-184.601642771250	  0.000000000000	-1087.384179219150	174.047538712650	  0.000000000000	-14.066329572050	 10.540419105500	  0.000000000000
-	-88.150018362800	132.633290377050	  0.000000000000	174.406707938500	-118.590048147650	  0.000000000000	-86.281056296900	-14.056122943500	  0.000000000000
-	  0.000000000000	  0.000000000000	  0.072073944700	  0.000000000000	  0.000000000000	 -0.056305990000	  0.000000000000	  0.000000000000	 -0.001517786900
-	132.630747385150	-88.141233436300	  0.000000000000	-14.056813464100	-86.281056332300	  0.000000000000	-118.586814694400	174.397922510750	  0.000000000000
-	-184.592873087100	1101.412194061400	  0.000000000000	 10.540419138100	-14.065638735650	  0.000000000000	174.038770019250	-1087.387969811900	  0.000000000000
-	  0.000000000000	  0.000000000000	  0.072003608900	  0.000000000000	  0.000000000000	 -0.001517786500	  0.000000000000	  0.000000000000	 -0.056235653400
diff --git a/src/VibrationalAnalysis.jl b/src/VibrationalAnalysis.jl
index 021980c..1e8a054 100644
--- a/src/VibrationalAnalysis.jl
+++ b/src/VibrationalAnalysis.jl
@@ -1,3 +1,8 @@
+""" 
+# VibrationalAnalysis.jl 
+
+This module contains functions to perform vibrational analysis on a QMCFC output.
+"""
 module VibrationalAnalysis
 
 using LinearAlgebra
diff --git a/src/calculate.jl b/src/calculate.jl
index 2d0c143..1843690 100644
--- a/src/calculate.jl
+++ b/src/calculate.jl
@@ -1,5 +1,20 @@
 export calculate
 
+"""
+    calculate(rst_file::String, hessian_file::String, moldescriptor_file::String) -> wavenumbers::Vector{Float64}, intensities::Vector{Float64}, force_constants::Vector{Float64}, reduced_masses::Vector{Float64}
+
+Reads the restart file, the hessian and the atom charges and calculates the wavenumbers, intensities, force constants and reduced masses.
+
+# Arguments
+- `rst_file::String`: The restart file.
+- `hessian_file::String`: The hessian file.
+- `moldescriptor_file::String`: The moldescriptor file.
+
+# Example
+```julia-repl
+julia> calculate("restart.rst", "hessian.dat", "moldescriptor.dat")
+```
+"""
 function calculate(rst_file::String, hessian_file::String, moldescriptor_file::String)
 
     # Read restart restart file 
@@ -14,7 +29,7 @@ function calculate(rst_file::String, hessian_file::String, moldescriptor_file::S
     # Transforme in internal coordinates
     eigenvalues, eigenvectors_internal_normalized, normalization = internal_coordinates(atom_coords, atom_masses, hessian_mw)
 
-    #Read the atom charges
+    # Read the atom charges
     atom_charges = read_moldescriptor(moldescriptor_file, atom_names, atom_types)
 
     # Calculate observables
diff --git a/src/coordinate_transform.jl b/src/coordinate_transform.jl
index 1c0a918..bf731d3 100644
--- a/src/coordinate_transform.jl
+++ b/src/coordinate_transform.jl
@@ -4,22 +4,41 @@ export center_to_com, inertia_tensor
     center_to_com(coord::Matrix{Float64}, masses::Vector{Float64})
 
 Translate the coordinates to the center of mass.
-"""
 
+# Arguments
+- `coord::Matrix{Float64}`: The coordinates.
+- `masses::Vector{Float64}`: The masses.
+
+# Example
+```julia-repl
+julia> center_to_com(coord, masses)
+```
+"""
 function center_to_com(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64})
+    
     # Calculate the center of mass
     com = sum(atom_coords .* atom_masses, dims=1) / sum(atom_masses)
+    
     # Translate the coordinates to the center of mass
-    atom_coords = atom_coords .- com
+    atom_coords = atom_coords .- com#
+
     return atom_coords
 end
 
 """
     inertia_tensor(coord::Matrix{Float64}, masses::Vector{Float64})
 
-Calculate the inertia tensor.    
-"""
+Calculate the inertia tensor of a molecule.
 
+# Arguments
+- `coord::Matrix{Float64}`: The coordinates.
+- `masses::Vector{Float64}`: The masses.
+
+# Example
+```julia-repl
+julia> inertia_tensor(coord, masses)
+```
+"""
 function inertia_tensor(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64})
     
     # center the coordinates to the center of mass
diff --git a/src/massesdict.jl b/src/massesdict.jl
index c24554e..f40a38b 100644
--- a/src/massesdict.jl
+++ b/src/massesdict.jl
@@ -1,3 +1,16 @@
+"""
+    masses::Dict{String, Float64}
+
+Atomic masses in amu. Letter all in lowercase.
+
+# Example
+```julia-repl
+julia> masses["h"]
+1.00794
+julia> masses["c"]
+12.0107
+```
+"""
 masses = Dict(
     "h"   => 1.00794,
     "d"   => 2.014101778,
diff --git a/src/observables.jl b/src/observables.jl
index ca28189..36ce943 100644
--- a/src/observables.jl
+++ b/src/observables.jl
@@ -2,12 +2,18 @@
 export wavenumber_kcal, wavenumber_dftb, reduced_mass, force_constant, infrared_intensity
 
 """
-    wavenumber_dftb(eigenvalues::Vector{Float64})
+    wavenumber_dftb(eigenvalues::Vector{Float64}) -> wavenumbers::Vector{Float64}, omega::Vector{Float64}
 
-Convert eigenvalues from Hartree Å^-2 g^-1 to wavenumbers in cm^-1.
+Convert `eigenvalues` from Hartree Å^-2 g^-1 to `wavenumbers` in cm^-1. Made for DFTB hessian files.
+Output include `omega` in s^-2.
 
-"""
+``ω = √v``
+
+``ν̃ = 1/(2π * c) * ω``
 
+# Arguments
+- `eigenvalues::Vector{Float64}`: The eigenvalues.
+"""
 function wavenumber_dftb(eigenvalues::Vector{Float64})
     
     # Conversion Hartree B^-2 g^-1 to s^-2: 
@@ -22,12 +28,18 @@ end
 
 
 """
-    wavenumber_kcal(eigenvalues::Vector{Float64})
+    wavenumber_kcal(eigenvalues::Vector{Float64}) -> wavenumbers::Vector{Float64}, omega::Vector{Float64}
 
-Convert eigenvalues from kcal Å^-2 g^-1 to wavenumbers in cm^-1.
+Convert `eigenvalues` from kcal Å^-2 g^-1 to `wavenumbers` in cm^-1. Made for QMCFC hessian files.
+Output include `omega` in s^-2.
 
-"""
+``ω = √v``
+
+``ν̃ = 1/(2π * c) * ω``
 
+# Arguments
+- `eigenvalues::Vector{Float64}`: The eigenvalues.
+"""
 function wavenumber_kcal(eigenvalues::Vector{Float64})
     
     # Convert eigenvalues in kcal Å^-2 g^-1 to wavenumbers in s^-2: 
@@ -41,22 +53,27 @@ function wavenumber_kcal(eigenvalues::Vector{Float64})
 end
 
 """
-    reduced_mass(normalization)
+    reduced_mass(normalization::Vector{Float64}) -> red_mass::Vector{Float64}
 
 Calculate the reduced mass from the normalization vector.
-"""
 
-function reduced_mass(normalization)
+# Arguments
+- `normalization::Vector{Float64}`: The normalization vector.
+"""
+function reduced_mass(normalization::Vector{Float64})
     red_mass = normalization .^ 2
     return red_mass[:] # convert to vector
 end
 
 """
-    force_constant(wavenumbers::Vector{Float64}, reduced_mass::Vector{Float64})
+    force_constant(wavenumbers::Vector{Float64}, reduced_mass::Vector{Float64}) -> force_const::Vector{Float64}
 
 Calculate the force constant from the wavenumbers and the reduced mass.
-"""
 
+# Arguments
+- `wavenumbers::Vector{Float64}`: The wavenumbers.
+- `reduced_mass::Vector{Float64}`: The reduced mass.
+"""
 function force_constant(omega::Vector{Float64}, reduced_mass::Vector{Float64})
     # Conversion g mol^-1 s^-2 to mdyn Å^-1: / 6.022 / 1E23 (mol) / 1E3 (kg/g) / 1E2 (mdyn/Å / N/m)
     force_const =  omega'.^2 .* reduced_mass / 6.022 / 1E28
@@ -64,20 +81,29 @@ function force_constant(omega::Vector{Float64}, reduced_mass::Vector{Float64})
 end
 
 """
-    infrared_intensity(eigenvectors_internal_normalized::Matrix{Float64}, coordinates::Matrix{Float64}, charges::Vector{Float64})
+    infrared_intensity(eigenvectors_internal_normalized::Matrix{Float64}, atom_charges::Vector{Float64}, reduced_mass::Vector{Float64}) -> intensities::Vector{Float64}
 
-Calculate the infrared intensity from the normalization eigen matrix, the coordinates and the charges.
-"""
+Calculate the infrared intensity in km mol^-1 from the normalization eigen matrix, the coordinates and the charges.
 
-function infrared_intensity(eigenvectors_internal_normalized::Matrix{Float64}, atom_charges::Vector{Float64}, reduced_masses)
+# Arguments
+- `eigenvectors_internal_normalized::Matrix{Float64}`: The normalized eigenvectors.
+- `atom_charges::Vector{Float64}`: The atom charges.
+- `reduced_mass::Vector{Float64}`: The reduced mass.
+"""
+function infrared_intensity(eigenvectors_internal_normalized::Matrix{Float64}, atom_charges::Vector{Float64}, reduced_masses::Vector{Float64})
+    
     intensities = []
+    
     for i in 1:size(eigenvectors_internal_normalized)[1]
+        
+        # Eigenvector of the i-th mode
         eigenvector = reshape(eigenvectors_internal_normalized[:, i], 3, :)'
-        # http://thiele.ruc.dk/~spanget/help/g09/k_constants.html
+
+        # Conversion factor taken from: http://thiele.ruc.dk/~spanget/help/g09/k_constants.html
         intensity = sum((sum(eigenvector .* atom_charges, dims=1) / 0.2081943  / norm(eigenvector)).^2 / reduced_masses[i] * 42.2561)
+        
         push!(intensities, intensity)
     end
+    
     return intensities
 end
-
-
diff --git a/src/read_files.jl b/src/read_files.jl
index e1a38bd..dacb82d 100644
--- a/src/read_files.jl
+++ b/src/read_files.jl
@@ -5,10 +5,12 @@ export read_rst, read_hessian, read_moldescriptor
 include("massesdict.jl")
 
 """
-    read_rst(rst_file::String)
+    read_rst(rst_file::String) -> atom_names::Vector{String}, atom_masses::Vector{Float64}, atom_coords::Matrix{Float64}, atom_types::Vector{Int64}
 
 Reads a restart file and returns a tuple of atom names, masses, and coordinates.
-    
+
+# Arguments
+- `rst_file::String`: The restart file.
 """
 function read_rst(rst_file::String)
 
@@ -48,10 +50,12 @@ function read_rst(rst_file::String)
 end
 
 """
-    read_hessian(hessian_file::String)
+    read_hessian(hessian_file::String) -> hessian::Matrix{Float64}
 
 Reads a hessian file and returns a nxn Matrix of the hessian.
 
+# Arguments
+- `hessian_file::String`: The hessian file.
 """
 function read_hessian(hessian_file::String)
     
@@ -91,10 +95,12 @@ function read_hessian(hessian_file::String)
 end
 
 """
-    read_moldescriptor(moldescriptor_file::String)
+    read_moldescriptor(moldescriptor_file::String) -> atom_charges::Vector{Float64}
 
 Reads a moldescriptor file and returns atom_charges.
-    
+
+# Arguments
+- `moldescriptor_file::String`: The moldescriptor file.
 """
 function read_moldescriptor(moldescriptor_file::String, atom_names::Vector{String}, atom_types::Vector{Int64})
 
diff --git a/src/symmetrize.jl b/src/symmetrize.jl
index bbc9d4c..86c79b6 100644
--- a/src/symmetrize.jl
+++ b/src/symmetrize.jl
@@ -1,23 +1,35 @@
 using LinearAlgebra
 
-export symmetrize_addition, symmetrize_multiplication, mass_weighted_hessian, mass_weighted_hessian_add
+export mass_weighted_hessian, mass_weighted_hessian_add, symmetrize_addition, symmetrize_multiplication
 
 """
-    symmetrize_addition(hessian::Matrix{Float64})
+    symmetrize_addition(hessian::Matrix{Float64}) -> hessian_sym::Matrix{Float64}
 
 Symmetrize a Hessian matrix by adding the transpose and dividing by 2.
-"""
 
+``H_sym = (H + H') / 2``
+
+# Arguments
+- `hessian::Matrix{Float64}`: The Hessian matrix.
+"""
 function symmetrize_addition(hessian::Matrix{Float64})
     return (hessian + hessian') / 2
 end
 
 """
-    symmetrize_multiplication(hessian::Matrix{Float64})
+    symmetrize_multiplication(hessian::Matrix{Float64}) -> hessian_sym::Matrix{Float64}
 
 Symmetrize a Hessian matrix by multiplying by the transpose and taking the square root of the absolute value of the eigenvalues.
-"""
 
+``H_2 = H * H'``
+
+``v, V = eigen(H_2)``
+
+``H_sym = V * √|v| * V'``
+
+# Arguments
+- `hessian::Matrix{Float64}`: The Hessian matrix.
+"""
 function symmetrize_multiplication(hessian::Matrix{Float64})
     _H = hessian' * hessian
     v, V = eigen(_H)
@@ -26,12 +38,14 @@ function symmetrize_multiplication(hessian::Matrix{Float64})
 end
 
 """
-    mass_weight_hessian(hessian::Matrix{Float64}, atom_masses::Vector{Float64})
+    mass_weight_hessian(hessian::Matrix{Float64}, atom_masses::Vector{Float64}) -> hessian::Matrix{Float64}
 
 Mass weight a Hessian matrix using matrix multiplication for symmetrization.
 
+# Arguments
+- `hessian::Matrix{Float64}`: The Hessian matrix.
+- `atom_masses::Vector{Float64}`: The masses of the atoms.
 """
-
 function mass_weighted_hessian(hessian::Matrix{Float64}, atom_masses::Vector{Float64})
     # Symmetrize the hessian
     hessian = symmetrize_multiplication(hessian)
@@ -52,8 +66,10 @@ end
 
 Mass weight a Hessian matrix using symmetrize by addition.
 
+# Arguments
+- `hessian::Matrix{Float64}`: The Hessian matrix.
+- `atom_masses::Vector{Float64}`: The masses of the atoms.
 """
-
 function mass_weighted_hessian_add(hessian::Matrix{Float64}, atom_masses::Vector{Float64})
     # Symmetrize the hessian
     hessian = symmetrize_addition(hessian)
diff --git a/src/transformation.jl b/src/transformation.jl
index 59ddf9f..3c807ba 100644
--- a/src/transformation.jl
+++ b/src/transformation.jl
@@ -2,11 +2,13 @@
 export internal_coordinates, transformation_matrix, translational_modes, rotational_modes
 
 """
-	translational_modes(atom_masses::Vector{Float64})
+	translational_modes(atom_masses::Vector{Float64}) -> translation::Matrix{Float64}
 	
-	Calculate the translational modes of a molecule.   
-"""
+Calculate the translational modes of a molecule. 
 
+# Arguments
+- `atom_masses::Vector{Float64}`: The masses of the atoms.
+"""
 function translational_modes(atom_masses::Vector{Float64})
 	# Create translation matrix 3N x 3
 	masses_diagonal = [diagm(0 => [i, i, i]) for i in atom_masses]
@@ -21,11 +23,14 @@ function translational_modes(atom_masses::Vector{Float64})
 end
 
 """
-	rotational_modes(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64})
+	rotational_modes(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64}) -> rotation::Matrix{Float64}
 	
-Calculate the rotational modes of a molecule.    
-"""
+Calculate the rotational modes of a molecule.
 
+# Arguments
+- `atom_coords::Matrix{Float64}`: The coordinates of the atoms.
+- `atom_masses::Vector{Float64}`: The masses of the atoms.
+"""
 function rotational_modes(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64})
 
 	# Translate the coordinates to the center of mass
@@ -54,11 +59,14 @@ function rotational_modes(atom_coords::Matrix{Float64}, atom_masses::Vector{Floa
 end
 
 """
-	transformation_matrix(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64})
+	transformation_matrix(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64}) -> transformation::Matrix{Float64}
 	
-	Calculate the transformation matrix.  
-"""
+Calculate the transformation matrix.
 
+# Arguments
+- `atom_coords::Matrix{Float64}`: The coordinates of the atoms.
+- `atom_masses::Vector{Float64}`: The masses of the atoms.
+"""
 function transformation_matrix(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64})
 
 	# Initialize transformation matrix 3N x 3N
@@ -77,12 +85,17 @@ function transformation_matrix(atom_coords::Matrix{Float64}, atom_masses::Vector
 end
 
 """
-	internal_coordinates(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64}, hessian_mw::Matrix{Float64})
+	internal_coordinates(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64}, hessian_mw::Matrix{Float64}) -> eigenvalues::Vector{Float64}, eigenvectors_internal_normalized::Matrix{Float64}, normalization::Vector{Float64}
 	
-	Calculate the internal coordinates of a molecule.    
-"""
+Calculate the internal coordinates of a molecule.
 
+# Arguments
+- `atom_coords::Matrix{Float64}`: The coordinates of the atoms.
+- `atom_masses::Vector{Float64}`: The masses of the atoms.
+- `hessian_mw::Matrix{Float64}`: The mass weighted hessian.
+"""
 function internal_coordinates(atom_coords::Matrix{Float64}, atom_masses::Vector{Float64}, hessian_mw::Matrix{Float64})
+	
 	# Calculate the transformation matrix
 	transformation = transformation_matrix(atom_coords, atom_masses)
 	# Gram-Schmidt orthogonalization
@@ -98,7 +111,7 @@ function internal_coordinates(atom_coords::Matrix{Float64}, atom_masses::Vector{
 	eigenvectors_internal = masses_matrix * transformation * eigenvectors
 
 	# Normalize Vectors
-	normalization = sqrt.(1 ./ sum(eigenvectors_internal .^ 2, dims = 1))
+	normalization = sqrt.(1 ./ sum(eigenvectors_internal .^ 2, dims = 1))[:] # convert to vector
 	# Normalize the eigenvectors
 	eigenvectors_internal_normalized = eigenvectors_internal .* normalization
 
diff --git a/src/write_out.jl b/src/write_out.jl
index 4e0eb0e..edee9e5 100644
--- a/src/write_out.jl
+++ b/src/write_out.jl
@@ -1,12 +1,21 @@
 
-export write_modes
+export write_modes, write_wavenumber_intensity
 
 """
-    write_modes(eigenvector_internal::Matrix{Float64}, coord::Matrix{Float64}, atom_names::Vector{String})
+    write_modes(eigenvector_internal::Matrix{Float64}, coord::Matrix{Float64}, atom_names::Vector{String}; filename="modes", amplitude=0.25, step=0.01) -> nothing
 
-Write the modes to a file.
-"""
+Write the modes to a file in xyz format.
+
+# Arguments
+- `eigenvector_internal::Matrix{Float64}`: The eigenvectors in internal coordinates.
+- `coord::Matrix{Float64}`: The coordinates of the atoms.
+- `atom_names::Vector{String}`: The names of the atoms.
 
+# Keyword Arguments
+- `filename::String`: The name of the file.
+- `amplitude::Float64`: The amplitude of the mode.
+- `step::Float64`: The step size of the mode.
+"""
 function write_modes(eigenvectors_internal_normalized::Matrix{Float64}, atom_coords::Matrix{Float64}, atom_names::Vector{String}; filename="modes", amplitude=0.25, step=0.01)
     for i in 1:size(eigenvectors_internal_normalized)[1]
 
@@ -33,11 +42,17 @@ function write_modes(eigenvectors_internal_normalized::Matrix{Float64}, atom_coo
 end
 
 """
-    write_wavenumber_frequency(wavenumbers::Vector{Float64}, intensities::Vector{Float64})
+    write_wavenumber_frequency(wavenumbers::Vector{Float64}, intensities::Vector{Float64}; filename=stdout) -> nothing
 
 Write the wavenumbers and intensities to a file.
-"""
 
+# Arguments
+- `wavenumbers::Vector{Float64}`: The wavenumbers.
+- `intensities::Vector{Float64}`: The intensities.
+
+# Keyword Arguments
+- `filename::String`: The name of the file.
+"""
 function write_wavenumber_intensity(wavenumbers, intensities; filename=stdout)
     
     # Open file
@@ -47,9 +62,9 @@ function write_wavenumber_intensity(wavenumbers, intensities; filename=stdout)
         file = filename
     end
 
-    println(file, "Wavenumbers (cm-1)    Intensities (km mol-1)")
+    println(file, "# Wavenumbers (cm-1)    Intensities (km mol-1)")
 
-    for i in 1:length(wavenumbers)
+    for i in eachindex(wavenumbers)
         println(file, wavenumbers[i], "    ", intensities[i])
     end