Prepare 0.4.0 structure input release

CHANGELOG.md

@@ -0,0 +1,7 @@
+# Changelog
+
+## 0.4.0 - 2026-06-09
+
+- add `read_structure(...)` as the primary public structure reader for restart and XYZ inputs
+- align the CLI and documentation on neutral `structure` input naming instead of restart-specific terminology
+- define XYZ plus moldescriptor compatibility explicitly: only single-molecule moldescriptor files are accepted because XYZ files do not encode molecule-type assignments

CITATION.cff

@@ -15,5 +15,5 @@ license:
 - mit
 title: VibrationalAnalysis.jl
 type: software
-version: 0.3.3
+version: 0.4.0
 url: https://github.com/MolarVerse/VibrationalAnalysis.jl

Project.toml

@@ -1,7 +1,7 @@
 name = "VibrationalAnalysis"
 uuid = "04082e5c-1d57-4577-9a0f-61e6cd56dade"
 authors = ["Josef M. Gallmetzer"]
-version = "0.3.3"
+version = "0.4.0"
 
 [deps]
 Comonicon = "863f3e99-da2a-4334-8734-de3dacbe5542"

README.md

@@ -30,7 +30,9 @@ Make sure you have `~/.julia/bin` in your PATH.
 shell> vibrationalanalysis -h
 ```
 
-The CLI accepts either restart files or standard single-structure XYZ files as the geometry input.
+The CLI accepts either restart files or standard single-structure XYZ files as the structure input.
+
+The recommended programmatic entrypoint is `read_structure(...)`, which dispatches to restart or XYZ parsing.
 
 # Acknowledgements
 This package was developed as part of the [MolarVerse](https://github.com/MolarVerse) organization. Significant contributions were made by:

docs/src/cli.md

@@ -29,10 +29,10 @@ Run the following command to see the available options:
 vibrationalanalysis -h
 ```
 
-Run a calculation with restart and Hessian files:
+Run a calculation with structure and Hessian files:
 
 ```bash
-vibrationalanalysis restart.rst hessian.dat --unit kcal --output wavenumbers.dat
+vibrationalanalysis structure.rst hessian.dat --unit kcal --output wavenumbers.dat
 ```
 
 Run the same calculation from a standard single-structure XYZ file:
@@ -44,13 +44,13 @@ vibrationalanalysis structure.xyz hessian.dat --unit kcal --output wavenumbers.d
 Include intensities by passing a moldescriptor file:
 
 ```bash
-vibrationalanalysis restart.rst hessian.dat --moldescriptor moldescriptor.dat --output wavenumbers.dat
+vibrationalanalysis structure.rst hessian.dat --moldescriptor moldescriptor.dat --output wavenumbers.dat
 ```
 
-When an XYZ file is used, atom types default to `1` for all atoms. Multi-molecule moldescriptor files therefore still require restart input with explicit atom-type information.
+When an XYZ file is used, atom types default to `1` for all atoms. If `--moldescriptor` is used with an XYZ file, the moldescriptor file must contain exactly one molecule definition because XYZ files do not encode molecule-type assignments.
 
 Write normal modes in matrix notation or xyz trajectories:
 
 ```bash
-vibrationalanalysis restart.rst hessian.dat --normal-modes normal_modes.dat --modes
+vibrationalanalysis structure.rst hessian.dat --normal-modes normal_modes.dat --modes
 ```

docs/src/examples.md

@@ -2,11 +2,11 @@
 
 Here are some examples of how to use the `VibrationalAnalysis` package.
 
-Calculate the vibrational analysis of a molecule from the output of a PQ calculation.
+Calculate the vibrational analysis of a molecule from a structure file.
 ```julia
 using VibrationalAnalysis
 
-atom_names, atom_masses, atom_coords, atom_types = read_rst("path/to/file.rst")
+atom_names, atom_masses, atom_coords, atom_types = read_structure("path/to/file.rst")
 hessian = read_hessian("path/to/file.hessian")
 atom_charges = read_moldescriptor("path/to/file.moldescriptor", atom_names, atom_types)
 
@@ -16,12 +16,12 @@ wavenumbers, intensities, force_constants, reduced_masses, normal_modes =
 write_calculate_output(wavenumbers, intensities, force_constants, reduced_masses)
 ```
 
-The same workflow works with a standard single-structure XYZ file when all atoms belong to the same molecule type.
+The same workflow works with a standard single-structure XYZ file. If intensities are needed, the moldescriptor file must contain exactly one molecule definition.
 
 ```julia
 using VibrationalAnalysis
 
-atom_names, atom_masses, atom_coords, atom_types = read_xyz("path/to/file.xyz")
+atom_names, atom_masses, atom_coords, atom_types = read_structure("path/to/file.xyz")
 hessian = read_hessian("path/to/file.hessian")
 
 wavenumbers, force_constants, reduced_masses, normal_modes =

src/VibrationalAnalysis.jl

@@ -12,8 +12,8 @@ This module contains functions to perform vibrational analysis on a QMCFC output
 Can read directly from restart files, XYZ files, moldescriptor files and hessian files and perform vibrational analysis on the system.
 
 ```julia-repl
-julia> atom_names, atom_masses, atom_coords, atom_types = read_rst("restart.rst")
-julia> atom_names, atom_masses, atom_coords, atom_types = read_xyz("structure.xyz")
+julia> atom_names, atom_masses, atom_coords, atom_types = read_structure("structure.rst")
+julia> atom_names, atom_masses, atom_coords, atom_types = read_structure("structure.xyz")
 julia> hessian = read_hessian("hessian.dat")
 julia> atom_charges = read_moldescriptor("moldescriptor.dat", atom_names, atom_types)
 julia> calculate(atom_masses, atom_coords, atom_charges, hessian)

src/main.jl

@@ -9,7 +9,7 @@ Calculate the wavenumbers, intensities, force constants, reduced masses and eige
 
 # Args
 
-- `restart`: The restart or XYZ structure file.
+- `structure`: The restart or XYZ structure file.
 - `hessian`: The hessian file.
 
 # Options
@@ -24,13 +24,15 @@ Calculate the wavenumbers, intensities, force constants, reduced masses and eige
 - `--modes`: Write the modes in xyz format. If not specified, the modes will not be written.
 
 """
-Comonicon.@main function vibrationalanalysis(restart::String, hessian::String; unit = "kcal", moldescriptor = nothing, output = nothing, normal_modes = nothing, modes::Bool = false)
+Comonicon.@main function vibrationalanalysis(structure::String, hessian::String; unit = "kcal", moldescriptor = nothing, output = nothing, normal_modes = nothing, modes::Bool = false)
 
 	# Check if the unit is valid
 	wavenumber = check_unit(unit)
 
+	input_format = structure_format(structure)
+
 	# Read structure file
-	atom_names, atom_masses, atom_coords, atom_types = read_structure(restart)
+	atom_names, atom_masses, atom_coords, atom_types = read_structure(structure, format = input_format)
 
 	# Read the hessian
 	hessian = read_hessian(hessian)
@@ -44,6 +46,9 @@ Comonicon.@main function vibrationalanalysis(restart::String, hessian::String; u
 		# write wavenumbers, force constants, reduced masses
 		write_calculate_output(wavenumbers, force_constants, reduced_masses, filename = output)
 	else
+		if input_format == :xyz && moldescriptor_molecule_count(moldescriptor) != 1
+			error("XYZ input requires a moldescriptor file with exactly one molecule type because XYZ files do not encode molecule-type assignments.")
+		end
 
 		# Read the atom charges
 		atom_charges = read_moldescriptor(moldescriptor, atom_names, atom_types)

src/read_files.jl

@@ -1,4 +1,4 @@
-export read_rst, read_xyz, read_hessian, read_moldescriptor
+export read_structure, read_rst, read_xyz, read_hessian, read_moldescriptor
 
 """
 	read_rst(rst_file::String) -> atom_names::Vector{String}, atom_masses::Vector{Float64}, atom_coords::Matrix{Float64}, atom_types::Vector{Int64}
@@ -132,22 +132,85 @@ function read_xyz(xyz_file::String)
 	return atom_names, atom_masses, Matrix(hcat(atom_coords...)'), ones(Int64, number_atoms)
 end
 
+function structure_format(structure_file::String; format::Symbol = :auto)
+	if format == :auto
+		return lowercase(splitext(structure_file)[2]) == ".xyz" ? :xyz : :rst
+	elseif format in (:rst, :xyz)
+		return format
+	end
+
+	error("Unsupported structure format. Use :auto, :rst or :xyz.")
+end
+
 """
-	read_structure(structure_file::String) -> atom_names::Vector{String}, atom_masses::Vector{Float64}, atom_coords::Matrix{Float64}, atom_types::Vector{Int64}
+	read_structure(structure_file::String; format::Symbol = :auto) -> atom_names::Vector{String}, atom_masses::Vector{Float64}, atom_coords::Matrix{Float64}, atom_types::Vector{Int64}
 
 Reads either a restart file or a single-structure XYZ file.
 
 # Arguments
 - `structure_file::String`: The restart or XYZ file.
+- `format::Symbol`: The structure format. Use `:auto`, `:rst` or `:xyz`.
 """
-function read_structure(structure_file::String)
-	if lowercase(splitext(structure_file)[2]) == ".xyz"
+function read_structure(structure_file::String; format::Symbol = :auto)
+	if structure_format(structure_file, format = format) == :xyz
 		return read_xyz(structure_file)
 	end
 
 	return read_rst(structure_file)
 end
 
+function moldescriptor_types(moldescriptor_file::String)
+
+	# Check if the file exists and
+	if !isfile(moldescriptor_file)
+		error("The moldescriptor file does not exist.")
+	end
+
+	# Check if the file is empty
+	if filesize(moldescriptor_file) == 0
+		error("The moldescriptor file is empty.")
+	end
+
+	# Delete empty lines, lines containing whitespaces and comments
+	moldescriptor_lines = filter(x -> x != "" && occursin(r"\w+", x) && !occursin(r"#", x), readlines(moldescriptor_file))
+
+	# Strip the moldescriptor lines that are not lines of 3 strings
+	moldescriptor_lines = filter(x -> length(split(x)) == 3, moldescriptor_lines)
+
+	# Check if the moldescriptor lines are empty
+	if length(moldescriptor_lines) == 0
+		error("The moldescriptor file is not a moldescriptor file.")
+	end
+
+	# Check if the moldescriptor file contains only lines of 3 strings
+	# First is a string, second is an integer and third is a float
+	for line in moldescriptor_lines
+		if !(occursin(r"\w+", split(line)[1]) && occursin(r"\d+", split(line)[2]) && occursin(r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?", split(line)[3]))
+			error("The moldescriptor file contains lines that are not of the form: string integer float.")
+		end
+	end
+
+	mol_types = Vector{Dict{String, Float64}}()
+	i = 1
+
+	while i <= length(moldescriptor_lines)
+		number_atoms = parse(Int64, split(moldescriptor_lines[i])[2])
+		molecule = moldescriptor_lines[i+1:i+number_atoms]
+		dict = Dict{String, Float64}()
+		for atom in molecule
+			dict[lowercase(split(atom)[1])] = parse(Float64, split(atom)[3])
+		end
+		push!(mol_types, dict)
+		i = i + number_atoms + 1
+	end
+
+	return mol_types
+end
+
+function moldescriptor_molecule_count(moldescriptor_file::String)
+	return length(moldescriptor_types(moldescriptor_file))
+end
+
 """
 	read_hessian(hessian_file::String) -> hessian::Matrix{Float64}
 
@@ -203,47 +266,10 @@ Reads a moldescriptor file and returns atom_charges.
 """
 function read_moldescriptor(moldescriptor_file::String, atom_names::Vector{String}, atom_types::Vector{Int64})
 
-	# Check if the file exists and
-	if !isfile(moldescriptor_file)
-		error("The moldescriptor file does not exist.")
-	end
-
-	# Check if the file is empty
-	if filesize(moldescriptor_file) == 0
-		error("The moldescriptor file is empty.")
-	end
-
-	# Delete empty lines, lines containing whitespaces and comments
-	moldescriptor_lines = filter(x -> x != "" && occursin(r"\w+", x) && !occursin(r"#", x), readlines(moldescriptor_file))
-
-	# Strip the moldescriptor lines that are not lines of 3 strings
-	moldescriptor_lines = filter(x -> length(split(x)) == 3, moldescriptor_lines)
-
-	# Check if the moldescriptor lines are empty
-	if length(moldescriptor_lines) == 0
-		error("The moldescriptor file is not a moldescriptor file.")
-	end
-
-	# Check if the moldescriptor file contains only lines of 3 strings
-	# First is a string, second is an integer and third is a float
-	for line in moldescriptor_lines
-		if !(occursin(r"\w+", split(line)[1]) && occursin(r"\d+", split(line)[2]) && occursin(r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?", split(line)[3]))
-			error("The moldescriptor file contains lines that are not of the form: string integer float.")
-		end
-	end
-
-	mol_types = Vector()
-	i = 1
+	mol_types = moldescriptor_types(moldescriptor_file)
 
-	while i <= length(moldescriptor_lines)
-		number_atoms = parse(Int64, split(moldescriptor_lines[i])[2])
-		molecule = moldescriptor_lines[i+1:i+number_atoms]
-		dict = Dict()
-		for atom in molecule
-			dict[lowercase(split(atom)[1])] = parse(Float64, split(atom)[3])
-		end
-		push!(mol_types, dict)
-		i = i + number_atoms + 1
+	if any(atom_type < 1 || atom_type > length(mol_types) for atom_type in atom_types)
+		error("The atom types are incompatible with the moldescriptor file.")
 	end
 
 	atom_charges = [mol_types[atom_types[i]][lowercase(atom_names[i])] for i in 1:length(atom_names)]

test/data/test_h2o_moldescriptor.dat

@@ -0,0 +1,5 @@
+# Molecule 1
+H2O 3 0.0
+O 0 -0.65966
+H 1 0.32983
+H 1 0.32983

test/main.jl

@@ -0,0 +1,27 @@
+using Test
+
+@testset "CLI structure input - xyz" begin
+    temp = mktempdir()
+    xyz_file = joinpath(@__DIR__, "data", "test_h2o.xyz")
+    hessian_file = joinpath(@__DIR__, "data", "test_hessian_h2o.dat")
+    moldescriptor_file = joinpath(@__DIR__, "data", "test_h2o_moldescriptor.dat")
+
+    cd(temp) do
+        VibrationalAnalysis.vibrationalanalysis(xyz_file, hessian_file, output = "xyz_wavenumbers.dat")
+        wavenumber_lines = readlines("xyz_wavenumbers.dat")
+        @test wavenumber_lines[1] == "# Wavenumbers (cm-1)  Force constants (mdyn Å-1)  Reduced masses (amu)"
+        @test length(wavenumber_lines) == 10
+
+        VibrationalAnalysis.vibrationalanalysis(xyz_file, hessian_file, moldescriptor = moldescriptor_file, output = "xyz_intensities.dat")
+        intensity_lines = readlines("xyz_intensities.dat")
+        @test intensity_lines[1] == "# Wavenumbers (cm-1)  Intensities (km mol-1)  Force constants (mdyn Å-1)  Reduced masses (amu)"
+        @test length(intensity_lines) == 10
+    end
+end
+
+@testset "CLI structure input - xyz moldescriptor restrictions" begin
+    xyz_file = joinpath(@__DIR__, "data", "test_h2o.xyz")
+    hessian_file = joinpath(@__DIR__, "data", "test_hessian_h2o.dat")
+    moldescriptor_file = joinpath(@__DIR__, "data", "test_moldescriptor.dat")
+    @test_throws ErrorException VibrationalAnalysis.vibrationalanalysis(xyz_file, hessian_file, moldescriptor = moldescriptor_file)
+end

test/read_files.jl

@@ -1,4 +1,4 @@
-using VibrationalAnalysis: read_rst, read_xyz, read_hessian, read_moldescriptor
+using VibrationalAnalysis: read_structure, read_rst, read_xyz, read_hessian, read_moldescriptor
 using Test
 
 # Test reading RST files src/read_files.jl
@@ -57,6 +57,16 @@ end
     @test atom_types == ones(Int64, 3)
 end
 
+@testset "Read Structure File" begin
+    rst = read_structure("data/test_h2o.rst")
+    xyz = read_structure("data/test_h2o.xyz")
+    xyz_explicit = read_structure("data/test_h2o.xyz", format = :xyz)
+    @test rst == read_rst("data/test_h2o.rst")
+    @test xyz == read_xyz("data/test_h2o.xyz")
+    @test xyz_explicit == read_xyz("data/test_h2o.xyz")
+    @test_throws ErrorException read_structure("data/test_h2o.rst", format = :invalid)
+end
+
 @testset "Read Hessian File Exceptions" begin
     # Test Hessian File Exceptions
     @test_throws ErrorException read_hessian("data/not_a_file.dat")

test/runtests.jl

@@ -13,7 +13,9 @@ include("transformation.jl")
 include("observables.jl")
 # Test Calculate Functions
 include("calculate.jl")
+# Test CLI entrypoint
+include("main.jl")
 # Test Write Functions
 include("write.jl")
 # Test Check Functions
-include("check.jl")
+include("check.jl")