Merge branch 'release/0.1.4'

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"

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
+```

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

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

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

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,

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,43 +53,57 @@ 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
     return force_const[:] # convert to vector
 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
-
-

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})
 

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)