Merge branch 'release/0.1.5'

MolarVerse · Feb 12, 2024 · 4bef0ac · 4bef0ac · galjos · Feb 18, 2024
2 parents 86f7d57 + c34f027
commit 4bef0ac
Show file tree

Hide file tree

Showing 27 changed files with 223 additions and 115 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,12 +1,12 @@
 name = "VibrationalAnalysis"
 uuid = "04082e5c-1d57-4577-9a0f-61e6cd56dade"
 authors = ["Josef M. Gallmetzer"]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-julia = "1.9"
+julia = "1.10"
diff --git a/README.md b/README.md
@@ -1,5 +1,6 @@
-[![codecov](https://codecov.io/gh/galjos/VibrationalAnalysis.jl/graph/badge.svg?token=PIG1D1QIEE)](https://codecov.io/gh/galjos/VibrationalAnalysis.jl)
-[![Build CI](https://github.com/galjos/VibrationalAnalysis.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/galjos/VibrationalAnalysis.jl/actions/workflows/CI.yml)
+[![codecov](https://codecov.io/gh/MolarVerse/VibrationalAnalysis.jl/graph/badge.svg?token=kESDHEzXcY)](https://codecov.io/gh/MolarVerse/VibrationalAnalysis.jl)
+[![CI](https://github.com/MolarVerse/VibrationalAnalysis.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/MolarVerse/VibrationalAnalysis.jl/actions/workflows/CI.yml)
+[![TagBot](https://github.com/MolarVerse/VibrationalAnalysis.jl/actions/workflows/TagBot.yml/badge.svg)](https://github.com/MolarVerse/VibrationalAnalysis.jl/actions/workflows/TagBot.yml)
 
 # VibrationalAnalysis.jl
 
@@ -13,7 +14,8 @@ julia> ] add VibrationalAnalysis
 ## Usage
 ```julia-repl
 julia> using VibrationalAnalysis
-julia> calculate("restart.rst", "hessian.dat", "moldescriptor.dat")
+julia> read_calculate("restart.rst", "hessian.dat", "moldescriptor.dat")
+julia> calculate(atom_masses, atom_coords, atom_charges, hessian)
 ```
 
 ## Documentation

diff --git a/script/dicts__massesdict.jl b/script/dicts__massesdict.jl
@@ -0,0 +1,112 @@
+masses = Dict(
+    "h" => 1.00794,
+    "d" => 2.014101778,
+    "t" => 3.0160492675,
+    "he" => 4.002602,
+    "li" => 6.941,
+    "be" => 9.012182,
+    "b" => 10.811,
+    "c" => 12.0107,
+    "n" => 14.0067,
+    "o" => 15.9994,
+    "f" => 18.9984032,
+    "ne" => 20.1797,
+    "na" => 22.989770,
+    "mg" => 24.3050,
+    "al" => 26.981538,
+    "si" => 28.0855,
+    "p" => 30.973761,
+    "s" => 32.065,
+    "cl" => 35.453,
+    "ar" => 39.948,
+    "k" => 39.0983,
+    "ca" => 40.078,
+    "sc" => 44.955910,
+    "ti" => 47.880,
+    "v" => 50.9415,
+    "cr" => 51.9961,
+    "mn" => 54.938049,
+    "fe" => 55.845,
+    "co" => 58.933200,
+    "ni" => 58.6934,
+    "cu" => 63.546,
+    "zn" => 65.399,
+    "ga" => 69.723,
+    "ge" => 72.64,
+    "as" => 74.92160,
+    "se" => 78.96,
+    "br" => 79.904,
+    "kr" => 83.798,
+    "rb" => 85.4678,
+    "sr" => 87.62,
+    "y" => 88.90585,
+    "zr" => 91.224,
+    "nb" => 92.90638,
+    "mo" => 95.94,
+    "tc" => 98.9063,
+    "ru" => 101.07,
+    "rh" => 102.9055,
+    "pd" => 106.42,
+    "ag" => 107.8682,
+    "cd" => 112.411,
+    "in" => 114.818,
+    "sn" => 118.71,
+    "sb" => 121.76,
+    "te" => 127.6,
+    "i" => 126.90447,
+    "xe" => 131.293,
+    "cs" => 132.90546,
+    "ba" => 137.327,
+    "la" => 138.9055,
+    "ce" => 140.116,
+    "pr" => 140.90765,
+    "nd" => 144.24,
+    "pm" => 146.9151,
+    "sm" => 150.36,
+    "lr" => 260.1053,
+    "eu" => 151.964,
+    "gd" => 157.25,
+    "tb" => 158.92534,
+    "dy" => 162.5,
+    "ho" => 164.93032,
+    "er" => 167.259,
+    "tm" => 168.93421,
+    "yb" => 173.04,
+    "lu" => 174.967,
+    "hf" => 178.49,
+    "ta" => 180.9479,
+    "w" => 183.84,
+    "re" => 186.207,
+    "os" => 190.23,
+    "ir" => 192.217,
+    "pt" => 195.078,
+    "au" => 196.96655,
+    "hg" => 200.59,
+    "tl" => 204.3833,
+    "pb" => 207.2,
+    "bi" => 208.98038,
+    "po" => 208.9824,
+    "at" => 209.9871,
+    "rn" => 222.0176,
+    "fr" => 223.0197,
+    "ra" => 226.0254,
+    "ac" => 227.0278,
+    "th" => 232.0381,
+    "pa" => 231.03588,
+    "u" => 238.0289,
+    "np" => 237.0482,
+    "pu" => 244.0642,
+    "am" => 243.0614,
+    "cm" => 247.0703,
+    "bk" => 247.0703,
+    "cf" => 251.0796,
+    "es" => 252.0829,
+    "fm" => 257.0951,
+    "md" => 258.0986,
+    "no" => 259.1009,
+    "q" => 999.00000,
+    "x" => 999.00000,
+    "cav" => 1000.00000,
+    "sup" => 1000000.0,
+    "dum" => 1.0
+)
diff --git a/script/test_script.jl b/script/test_script.jl
@@ -5,15 +5,15 @@ c = 2.99792458e10 # in cm/s https://discourse.julialang.org/t/the-speed-of-light
 
 include("dicts__massesdict.jl")
 
-rst = readdlm("ceo2_min.rst")
-atomNames = rst[:,1][1:end-1]
+rst = readdlm("../test/data/test_h2o.rst")
+atomNames = rst[:,1][3:end]
 nAtoms = size(atomNames)[1]
 
 masses = (x->masses[lowercase(x)]).(atomNames)
 masses_repeat = repeat(masses, inner=3)
 masses_mat = (masses_repeat * masses_repeat').^(1/2)
 
-hess = readdlm("hess.mat")
+hess = readdlm("../test/data/test_hessian_h2o.dat")
 
 H = hess' * hess
 
@@ -24,7 +24,7 @@ H_sym = V * Diagonal(sqrt.(abs.(v))) * V'
 H_sym_w = H_sym ./ masses_mat
 
 ### center of mass
-coord = rst[:, 4:6][1:end-1, :]
+coord = rst[:, 4:6][3:end, :]
 Rcm = sum(masses .* coord, dims = 1) ./ sum(masses)
 
 ### Inertia tensor
@@ -108,7 +108,7 @@ k = ω.^2 .* red_mass' / 6.022 / 1E28
 # writedlm(stdout, k)
 
 ### Charge HARD CODED ATTENTION!!!!!!! 
-charge = Dict("ce" => 1.796238775774, "o" => -0.898119387887)
+charge = Dict("o" => -0.65966, "h" => 0.32983)
 q = (x->charge[lowercase(x)]).(atomNames)
 
 intensities = []
@@ -121,8 +121,10 @@ for i in 1:size(Inormal)[1]
 end
 
 #### Print intensities
-println("Intensities")
-writedlm(stdout, hcat(ν_tilde[7:end], intensities[7:end]))
+println("Intensities: ", intensities)
+println("Wavenumber:", ν_tilde)
+println("Force Constants:", k)
+println("Reduced Masses:", red_mass)
 
 
 # Print all modes
@@ -157,49 +159,4 @@ writedlm(stdout, hcat(ν_tilde[7:end], intensities[7:end]))
 
 # println("Quadrupole moment: ", quad)
 
-using Plots
-using Trapz
-using LaTeXStrings
-
-# add a broadening function
-function add_broadening(list_ex_energy, list_osci_strength, line_profile::String="Lorentzian", line_param::Float64=10.0, step::Float64=0.1) where T<:Real
-    x_min = 0
-    x_max = 1250
-    x = x_min:step:x_max
-    y = zeros(length(x))
-
-    for xp in 1:length(x)
-        for (e, f) in zip(list_ex_energy, list_osci_strength)
-            if line_profile == "Gaussian"
-                y[xp] += f * exp(-((e - x[xp]) / line_param)^2)
-            elseif line_profile == "Lorentzian"
-                y[xp] += 0.5 * line_param * f / (pi * ((x[xp] - e)^2 + 0.25 * line_param^2))
-            end
-        end
-    end
-
-    # # normalize the lineshape so that its integral is 1
-    # integral = trapz(y, x)
-    # y /= integral
-
-    # # scale the lineshape to retain the same peak height
-    # max_y = maximum(y)
-    # y *= maximum(list_osci_strength) / max_y
-
-    return x, y
-end
-
-# plot the spectrum
-function plotSpectrum(ν, I)
-    # scatter(ν, I, label="Spectrum")
-    x,y = add_broadening(ν, I)
-    plot!(x, y, label="Spectrum", xlabel=L"wavenumber in cm${^-1}$", ylabel=L"Intesity in km mol$^{-1}$")
-end
-
-file = open("spectrum_ceo2.txt", "w")
-
-# write the spectrum to a file
-writedlm(file, hcat(add_broadening(real.(ν_tilde), intensities)...))
 
-# close the file
-close(file)
diff --git a/src/VibrationalAnalysis.jl b/src/VibrationalAnalysis.jl
@@ -2,11 +2,27 @@
 # VibrationalAnalysis.jl 
 
 This module contains functions to perform vibrational analysis on a QMCFC output.
+
+## Usage
+
+Can read directly from rst-file, moldescriptor-file and hessian-file and perform vibrational analysis on the system.
+
+```julia-repl
+julia> read_calculate("restart.rst", "hessian.dat", "moldescriptor.dat")
+```
+
+Or directly perform a vibrational analysis with atom masses, atom coordinates, atom charges and hessian of the system.
+
+```julia-repl
+julia> read_calculate(atom_masses, atom_coords, atom_charges, hessian)
+```
+
 """
 module VibrationalAnalysis
 
 using LinearAlgebra
 
+include("massesdict.jl")
 include("read_files.jl")
 include("coordinate_transform.jl")
 include("symmetrize.jl")

diff --git a/src/calculate.jl b/src/calculate.jl
@@ -1,7 +1,7 @@
-export calculate
+export calculate, read_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}
+    read_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.
 
@@ -12,26 +12,51 @@ Reads the restart file, the hessian and the atom charges and calculates the wave
 
 # Example
 ```julia-repl
-julia> calculate("restart.rst", "hessian.dat", "moldescriptor.dat")
+julia> read_calculate("restart.rst", "hessian.dat", "moldescriptor.dat")
 ```
 """
-function calculate(rst_file::String, hessian_file::String, moldescriptor_file::String)
+function read_calculate(rst_file::String, hessian_file::String, moldescriptor_file::String)
 
     # Read restart restart file 
     atom_names, atom_masses, atom_coords, atom_types = read_rst(rst_file)
 
     # Read the hessian
     hessian = read_hessian(hessian_file)
 
+    # Read the atom charges
+    atom_charges = read_moldescriptor(moldescriptor_file, atom_names, atom_types)
+
+    wavenumbers, intensities, force_constants, reduced_masses = calculate(atom_masses, atom_coords, atom_charges, hessian)
+
+    return wavenumbers, intensities, force_constants, reduced_masses
+end
+
+"""
+    calculate(atom_masses::Vector{Float64}, atom_coords::Matrix{Float64}, atom_charges::Vector{Float64}, hessian::Matrix{Float64}) -> wavenumbers::Vector{Float64}, intensities::Vector{Float64}, force_constants::Vector{Float64}, reduced_masses::Vector{Float64}
+
+Calculates the wavenumbers, intensities, force constants and reduced masses from the atom masses, atom coordinates, atom charges and the hessian.
+
+# Arguments
+- `atom_masses::Vector{Float64}`: Vector of atom_masses (n)
+- `atom_coords::Matrix{Float64}`: Matrix of atom coordinates (3xn)
+- `atom_charges::Vector{Float64}`: Vector of atom charges (n)
+- `hessian::Matrix{Float64}`: Matrix of the hessian (3nx3n)
+
+# Example
+```julia-repl
+julia> calculate(atom_masses, atom_coords, atom_charges, hessian)
+```
+
+"""
+
+function calculate(atom_masses::Vector{Float64}, atom_coords::Matrix{Float64}, atom_charges::Vector{Float64}, hessian::Matrix{Float64})
+
     # Symmetrize the hessian and mass weighting
     hessian_mw = mass_weighted_hessian(hessian, atom_masses)
 
-    # Transforme in internal coordinates
+    # Transform in internal coordinates
     eigenvalues, eigenvectors_internal_normalized, normalization = internal_coordinates(atom_coords, atom_masses, hessian_mw)
 
-    # Read the atom charges
-    atom_charges = read_moldescriptor(moldescriptor_file, atom_names, atom_types)
-
     # Calculate observables
     wavenumbers, omega = wavenumber_kcal(eigenvalues)
     reduced_masses = reduced_mass(normalization)

diff --git a/src/coordinate_transform.jl b/src/coordinate_transform.jl
@@ -1,5 +1,3 @@
-export center_to_com, inertia_tensor
-
 """
     center_to_com(coord::Matrix{Float64}, masses::Vector{Float64})
 

diff --git a/src/observables.jl b/src/observables.jl
@@ -1,6 +1,4 @@
 
-export wavenumber_kcal, wavenumber_dftb, reduced_mass, force_constant, infrared_intensity
-
 """
     wavenumber_dftb(eigenvalues::Vector{Float64}) -> wavenumbers::Vector{Float64}, omega::Vector{Float64}
 
@@ -76,7 +74,7 @@ Calculate the force constant from the wavenumbers and 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
+    force_const =  omega.^2 .* reduced_mass / 6.022 / 1E28
     return force_const[:] # convert to vector
 end
 

diff --git a/src/read_files.jl b/src/read_files.jl
@@ -1,9 +1,5 @@
-# Read Restart File
-
 export read_rst, read_hessian, read_moldescriptor
 
-include("massesdict.jl")
-
 """
     read_rst(rst_file::String) -> atom_names::Vector{String}, atom_masses::Vector{Float64}, atom_coords::Matrix{Float64}, atom_types::Vector{Int64}
 

diff --git a/src/symmetrize.jl b/src/symmetrize.jl
@@ -1,7 +1,3 @@
-using LinearAlgebra
-
-export mass_weighted_hessian, mass_weighted_hessian_add, symmetrize_addition, symmetrize_multiplication
-
 """
     symmetrize_addition(hessian::Matrix{Float64}) -> hessian_sym::Matrix{Float64}
 

diff --git a/src/transformation.jl b/src/transformation.jl
@@ -1,6 +1,3 @@
-
-export internal_coordinates, transformation_matrix, translational_modes, rotational_modes
-
 """
 	translational_modes(atom_masses::Vector{Float64}) -> translation::Matrix{Float64}