Reading data and writing calculation inputs

Choose a reader based on what the file represents. State holds a molecular geometry and optional calculation metadata; Space holds multiple states; Data stores a numerical array and metadata. Readers do not run Gaussian, CP2K, GROMACS, ORCA, or any other simulation program.

File or task	Recommended entry point	Result
Single-frame or multi-frame XYZ	Space.from_xyz(path)	One State per frame
XYZ containing carbohydrates	Space.from_xyz(path, glycans=True)	One Glycan per frame
Last carbohydrate frame	Glycan.from_file(path, software="xyz")	A Glycan with sugar assignments
Gaussian output	State.from_file(path, software="gaussian")	Geometry and fields found in the log
CP2K distance/dihedral trajectory	Data.from_file(path, software="cp2k", colvar=..., timestep=...)	Time and measured coordinate columns
CP2K ring puckering	Data.from_file(..., observable="cremer_pople")	Named time/pucker columns
Numeric CSV/TSV	Data.from_file(path, software="csv", ...)	NumPy array in .data
GROMACS XVG	Data.from_file(path, software="gromacs")	Numeric data and an inferred time unit
Gaussian input	Space(..., output_path=..., theory=...).to_file()	One input.com per state
XYZ output	InputBuilder.build(space=space.states, ..., software="xyz")	One geometry.xyz per state

Read every XYZ frame

from matmalib import Space

space = Space.from_xyz("data/glucose_puckers.xyz")
print("Frames:", len(space.states))
first, last = space.states[0], space.states[-1]
print("First frame:", first.id, first.NAtoms, first.xyz.shape)
print("Last frame:", last.id)
print("First three elements:", [a.type for a in first.atoms[:3]])

Each state receives a string frame ID ("0", "1", …), a path, atoms, and NAtoms. Atom IDs match their zero-based positions. Coordinates retain the file's units, normally angstroms. XYZ comment lines are not parsed for energies or times. Plain XYZ states do not automatically receive a bond graph; call connectivity_matrix() when a workflow needs one, or use carbohydrate detection through glycans=True.

Try changing: load data/glucose.xyz for one frame or data/glycan.xyz with glycans=True. Access space.states[0].rings after choosing the glycan option. An empty trajectory raises ValueError; the XYZ parser rejects malformed frame counts, missing rows, and nonnumeric coordinate values. Plain XYZ loading does not check coordinate finiteness; downstream conformer and trajectory analysis perform additional validation.

State.from_file(..., software="xyz") is not a supported dispatch. Use Space.from_xyz(...).states[0] for a plain single geometry. For the final saved frame, use .states[-1].

Read a numeric table

The CSV adapter wraps numpy.loadtxt. Supply skiprows for a header, delimiter for the separator, and dtype for the element type. It does not infer column names or units.

from pathlib import Path
from matmalib import Data

path = Path("measurements.csv")
path.write_text("time_fs,distance_A\n0,1.42\n5,1.44\n10,1.41\n")
result = Data.from_file(str(path), software="csv", skiprows=1)
result.columns = ("time_fs", "distance_A")
result.time_unit = "fs"

print(result.data)
print("Mean distance:", round(result.data[:, 1].mean(), 3), "angstrom")
print("Columns:", result.columns)

The result is a (3, 2) array and the mean is approximately 1.423 Å. Data(**kwargs) simply stores attributes and repr(data) displays those attributes. A one-row or one-column input follows loadtxt's usual dimension squeezing; use np.atleast_2d if your application requires a two-dimensional table.

Try changing: use a tab-separated string and delimiter="\t", or add another sample. For quoting, missing values, mixed fields, or named headers, parse with an appropriate table library and wrap the numeric result in Data(data=array, ...).

Extract a CP2K-style collective variable

For an XYZ trajectory, colvar selects two atoms for a distance or four atoms for a dihedral. Pass the time between saved frames in femtoseconds, not necessarily the integrator timestep. The first frame is assigned time zero.

The bundled glucose_puckers.xyz file contains 38 reference conformers, not a sampled molecular dynamics trajectory. Its generated time values below demonstrate the reader's mechanics and have no physical time interpretation.

from matmalib import Data

trajectory = Data.from_file(
    "data/glucose_puckers.xyz",
    software="cp2k",
    colvar=[0, 1],
    timestep=0.5,
)
print("Shape:", trajectory.data.shape)
print("Time unit:", trajectory.time_unit)
print("First five [time, distance] rows:\n", trajectory.data[:5])

Column zero contains time in fs; column one contains distance in the input length units. Four-atom selections produce degrees. The legacy parser collects atoms in their file order, so use ascending, unique indices. It does not preserve an arbitrary requested torsion ordering. A three-atom selection currently returns an (angle, axis) tuple from the geometry helper and cannot be assembled into the intended numeric table; compute angles frame by frame with Space.from_xyz instead.

Try changing: double timestep to 1.0. The first column doubles while the distances stay fixed. Both missing colvar and missing timestep raise ValueError on this XYZ path.

For ring trajectories, use the structured puckering interface:

from matmalib import Data

result = Data.from_file(
    "data/glucose_puckers.xyz",
    software="cp2k",
    observable="cremer_pople",
    timestep=0.5,
    classify=True,
)
print("Columns:", result.columns)
print("First numeric row:", result.data[0])
print("First pucker labels:", result.pucker_labels[0])

This path supports rings, conn_mat, distXX, distXH, degrees, and classify; automatic detection uses the first frame. With degrees=True (the default), angular columns use degrees and amplitudes use input length units. With degrees=False, angular coordinate columns use radians; classification distances retain their documented degree-based convention. Read trajectory analysis for column names, frame validation, and fixed ring ordering.

For a .metadynLog path, the adapter simply calls numpy.loadtxt and labels the time unit fs. It does not interpret, reorder, or convert the log's columns. Verify their meaning against the producing calculation.

Read GROMACS XVG tables

from pathlib import Path
from matmalib import Data

path = Path("distance.xvg")
path.write_text(
    '# Example XVG header\n'
    '@ xaxis label "Time (ps)"\n'
    '@ yaxis label "Distance (nm)"\n'
    '0 0.142\n1 0.145\n2 0.143\n'
)
result = Data.from_file(str(path), software="gromacs")
print("Time unit:", result.time_unit)
print("Unconverted numeric table:\n", result.data)

The reader skips the number of lines containing # or @, and infers time_unit from an xaxis line's final token. This assumes a contiguous header followed by a numeric table. It does not convert GROMACS distances from nm into Å, label dependent columns, or support multi-dataset separators. For this example, result.data[:, 1] * 10 would convert nm to Å.

The same adapter contains an .xpm matrix reader. It assumes a square grid, uppercase single-character cell keys, and a particular numeric legend layout, and reverses the row direction. Treat it as a format-specific legacy parser: check its output against the source image and use a specialized parser for other XPM variants.

Read Gaussian calculations

from matmalib import State, Space

state = State.from_file("calculations/conf_001/input.log", software="gaussian")
print(state.id, state.NAtoms, state.E)
if hasattr(state, "F"):
    print("Thermal free energy (hartree):", state.F)

# Reads input.log in each immediate subdirectory.
ensemble = Space.from_dir("calculations", software="gaussian")
ensemble.sort("E")
print(ensemble.diff(quant="E"))  # kcal/mol above the lowest E

The Gaussian parser detects sp, opt, freq, optfreq, or nmr from route text. Depending on the job and available sections, it attaches:

Field	Meaning / units
E	SCF electronic energy, hartree
Ezpe	Electronic plus zero-point energy, hartree
H / F	Electronic plus thermal enthalpy / free energy, hartree
Freq	Frequencies as printed by Gaussian, normally cm⁻¹
Ints	IR intensities as printed, normally km/mol
NVibs, Vibs	Mode count and (n_modes, n_atoms, 3) printed displacement arrays
nmr	Lower-triangular nuclear spin coupling array; no symmetry completion

State.from_file calls StateBuilder.build, which infers connectivity and fragments, then assigns id from the containing folder (or filename when there is no folder). The parser expects particular Gaussian text layouts and six-decimal coordinate lines. It is not a general log validator, and missing sections can cause missing attributes or errors. Inspect successful termination, atom count, final geometry, mode count, and extracted energies for your own logs before analysis. Ordinary notebook/unit tests do not establish support for every Gaussian version or multi-job layout.

Space.from_dir reads only immediate subdirectories, assumes every one contains input.log, does not sort directory names, and currently supports only Gaussian. Missing or unrelated folders are not skipped.

Generate Gaussian input files

The following example writes a real input file into the browser's temporary filesystem and prints it. It does not submit a calculation. Use explicit theory settings; the library's historical defaults are not a recommendation for a scientific method.

from pathlib import Path
from matmalib import Space

space = Space.from_xyz("data/glucose.xyz")
space.states[0].id = "glucose"
space.theory = {
    "disp": "GD3BJ",
    "nprocs": 4,
    "mem": "4GB",
    "method": "PBE1PBE",
    "basis_set": "6-31G(d)",
    "jobtype": "opt freq",
    "charge": 0,
    "multiplicity": 1,
    "other_options": "SCF=Tight",
}
output = Path("gaussian_inputs")
output.mkdir(exist_ok=True)
space.output_path = str(output)
space.to_file(software="gaussian")
print((output / "glucose" / "input.com").read_text())

The writer creates <output_path>/<state.id>/input.com. The parent directory must exist; each ID should be a unique, safe directory name. Repeated runs overwrite the same file. Coordinates are rounded to three decimals. The charge and multiplicity are taken literally; the writer does not validate electronic structure or Gaussian syntax.

Try changing: change jobtype to "sp", or set disp=False, then compare the route line. Recognized dispersion values include True, "GD3BJ", "EmpiricalDispersion=GD3BJ", "GD3", and "EmpiricalDispersion=GD3". Other values produce a blank dispersion option. The writer mutates the theory["disp"] entry while normalizing it.

All shown dictionary keys except other_options are required by the Gaussian writer. Optional extra adds literal text after the coordinates, for example a basis or constraints section appropriate to your chosen calculation. Space.to_file returns None. InputBuilder.build(..., names=None) provides the lower-level entry point; assign state.id yourself because its non-None names branch currently fails.

Export XYZ geometries

The XYZ writer takes an iterable of states. Pass space.states directly to InputBuilder.build: Space.to_file(software="xyz") currently passes a non-iterable Space and fails.

from pathlib import Path
import numpy as np
from matmalib import Space
from matmalib.calc import InputBuilder

space = Space.from_xyz("data/glucose.xyz")
space.states[0].id = "glucose"
output = Path("xyz_exports")
output.mkdir(exist_ok=True)
InputBuilder().build(
    space=space.states, theory={}, software="xyz", path=str(output)
)
path = output / "glucose" / "geometry.xyz"
restored = Space.from_xyz(str(path)).states[0]
print("Written file:", path)
print("Restored atoms:", restored.NAtoms)
print("Maximum rounding difference:", np.abs(restored.xyz - space.states[0].xyz).max())

Each state becomes a separate geometry.xyz, with state.id as the comment. Coordinates are rounded to three decimals, so expect differences up to 0.0005 per coordinate in the input units. This writer does not produce a concatenated trajectory or preserve energy metadata.

Legacy readers and writers

The code exposes additional methods that should not be treated as working end-to-end workflows:

Method	Current constraint
StateBuilder.cp2k_geom()	Reads the final XYZ frame, but passes its atom index into the unused valence parameter, leaving Atom.id=None. The subsequent legacy connectivity step is unsafe. Use Space.from_xyz(path).states[-1] instead.
StateBuilder.orca()	The optimization branch constructs an invalid companion-file path and references uninitialized variables; the energy-only branch does not populate the geometry needed by build(). No complete ORCA import workflow is documented as supported.
InputBuilder.fhiaims()	Assumes a site-specific species-defaults directory and an iterable of states; also constructs an absolute /geometry.in path and formats Atom objects as strings. Do not use this writer as a working FHI-aims exporter.
InputBuilder.build(..., software=...)	Unknown software values silently do nothing, unlike reader builders, which raise ValueError.

The core reference notes document the lower-level builder methods and all public core APIs. Continue to ensemble energies, spectra, or trajectory plots.