[WIP]Trajectory -> Gaussview readable ".log" file

src/pymatgen/core/trajectory.py

@@ -475,6 +475,47 @@ def write_Xdatcar(
         with zopen(filename, mode="wt", encoding="utf-8") as file:
             file.write(xdatcar_str)  # type:ignore[arg-type]
 
+    # TODO: This needs testing. The actual reading can't be tested automatically, but writing different
+    # trajectories without errors can be tested.
+    def write_Gaussian_log(
+        self,
+        filename: PathLike = "calc.log",
+        site_property_map: dict[str, str] | None = None,
+        write_lattice: bool = True,
+    ) -> None:
+        """Write Trajectory object to a Gaussian-formatted log file.
+
+        The produced file will be readable by Gaussview.
+
+        Args:
+            filename: File to write. File does not need to end with '.log' to be readable by Gaussview.
+                The suffix ".logx" is recommended to distinguish from real Gaussian log files.
+            site_property_map: A mapping between implemented Gaussian properties and pymatgen properties.
+                i.e. providing {"mulliken charges": "charges"} will write the data stored in
+                `Trajectory[i].site_properties["charges"]` as Gaussian logs the Mulliken charges for each site.
+                The currently implemented property keys are "mulliken charges", "mulliken spin", "nbo charges",
+                "nbo spin", and "esp charges". By default all charge keys are mapped to "charges", and all spin
+                keys are mapped to "magmom". All keys must map to site properties of `list[float]`
+            write_lattice: Whether to write the lattice information in the log file. Setting this to false will
+                not write the lattice information in the log file, but will enable the "forces" site properties
+                to be readable by Gaussview.
+        """
+        # Import inside function as this is an expensive import
+        # (causes import time for Trajectory to surpass threshold for 3 test setups)
+        from pymatgen.io.gaussian import _traj_to_gaussian_log
+
+        # Each frame must have an energy to be read properly by Gaussview.
+        if hasattr(self, "frame_properties") and self.frame_properties is not None:
+            energies = [fp.get("energy", 0) for fp in self.frame_properties]
+        else:
+            energies = [0] * len(self)
+        self.to_positions()
+        _sp_map = site_property_map or {}
+        lines = _traj_to_gaussian_log(self, write_lattice, energies, _sp_map)
+        with open(filename, "w") as file:
+            file.write(lines)
+        file.close()
+
     def as_dict(self) -> dict:
         """Return the trajectory as a MSONable dict."""
         lat = self.lattice.tolist() if self.lattice is not None else None

src/pymatgen/io/gaussian.py

@@ -21,8 +21,11 @@
 if TYPE_CHECKING:
     from pathlib import Path
 
+    from numpy.typing import NDArray
     from typing_extensions import Self
 
+    from pymatgen.core.structure import Structure
+    from pymatgen.core.trajectory import Trajectory
     from pymatgen.util.typing import PathLike
 
 __author__ = "Shyue Ping Ong, Germain Salvato-Vallverdu, Xin Chen"
@@ -1331,3 +1334,110 @@ def to_input(
             link0_parameters=link0_parameters,
             dieze_tag=dieze_tag,
         )
+
+
+def _traj_to_gaussian_log(traj: Trajectory, do_cell: bool, energies: list[float], sp_map: dict | None) -> str:
+    dump_lines = ["", " Entering Link 1 ", " ", *_log_input_orientation(traj[0], do_cell=do_cell)]
+    for i in range(len(traj)):
+        dump_lines += [
+            *[*_log_input_orientation(traj[i], do_cell=do_cell), f"SCF Done:  E =  {energies[i]:.6f}  "],
+            *_log_mulliken(traj[i], "charges", "spin"),
+            *_log_esp(traj[i], "charges"),
+            *_log_nbo(traj[i], "charges", "spin"),
+            *[*_log_forces(traj[i]), "", " " + "Grad" * 18, "", f" Step number {i + 1}", "", " " + "Grad" * 18],
+        ]
+    dump_lines[-2:-2] = [" Optimization completed.", "    -- Stationary point found."]
+    dump_lines += [*_log_input_orientation(traj[-1], do_cell=do_cell), " Normal termination of Gaussian 16"]
+    return "\n".join(dump_lines)
+
+
+def _log_input_orientation(frame: Structure | Molecule, do_cell=True) -> list[str]:
+    at_ns = frame.atomic_numbers
+    dump_lines = [
+        *[" " * 24 + "Standard orientation:" + " " * 26, " " + "-" * 69],
+        " Center     Atomic      Atomic             Coordinates (Angstroms)",
+        *[" Number     Number       Type             X           Y           Z", " " + "-" * 69],
+        *[f"{i + 1} {at_ns[i]} 0 {p[0]:.6f} {p[1]:.6f} {p[2]:.6f} " for i, p in enumerate(frame.cart_coords)],
+    ]
+    if do_cell:
+        cell = frame.lattice.matrix
+        dump_lines += [f"{i} -2 0 {v[0]:.6f} {v[1]:.6f} {v[2]:.6f} " for i, v in enumerate(cell, start=len(at_ns) + 1)]
+    dump_lines.append(" ---------------------------------------------------------------------")
+    return dump_lines
+
+
+def _log_mulliken(frame: Structure | Molecule, charge_key: str, spin_key: str) -> list[str]:
+    dump_lines = []
+    if charge_key in frame.site_properties:
+        charges = np.array(frame.site_properties[charge_key])
+        dump_lines += ["Mulliken charges:", "1"]
+        spins: None | NDArray = None
+        if spin_key in frame.site_properties:
+            spins = np.array(frame.site_properties[spin_key])
+            dump_lines[-2] = dump_lines[-2].replace(":", " and spin densities:")
+            dump_lines[-1] += " 2"
+        for i, charge in enumerate(charges):
+            dump_lines.append(f"{i + 1} {frame.symbol_set[i]} {charge} " + "" or f"{spins[i]} ")
+        dump_lines.append(f"Sum of Mulliken charges = {np.sum(charges)}" + "" if spins is None else f" {np.sum(spins)}")
+    return dump_lines
+
+
+def _log_esp(frame: Structure | Molecule, charge_key: str) -> list[str]:
+    dump_lines = []
+    if charge_key in frame.site_properties:
+        charges = np.array(frame.site_properties[charge_key])
+        dump_lines += ["", " Charges from ESP fit", "  ESP charges:", "               1"]
+        dump_lines += [f"     {i + 1} {s} {charges[i]:.6f}" for i, s in enumerate(frame.symbol_set)]
+        dipole = np.sum(frame.cart_coords * charges[:, np.newaxis], axis=0)
+        tot = np.linalg.norm(dipole)
+        dump_lines.append(f" Charge= {np.sum(charges)} Dipole= {' '.join(str(v) for v in dipole)} Tot=   {tot}")
+    return dump_lines
+
+
+def _log_nbo_fields(frame: Structure | Molecule, charges: NDArray, totals: NDArray | None = None) -> list[str]:
+    _totals = totals if isinstance(totals, np.ndarray) else np.zeros_like(charges)
+    dump_lines = [
+        *[" Summary of Natural Population Analysis:  ", "", " " * 39 + "Natural Population "],
+        " " * 17 + "Natural " + "-" * 47 + " ",
+        *["    Atom  No    Charge         Core      Valence    Rydberg      Total", " " + "-" * 71],
+    ]
+    for i, s in enumerate(frame.symbol_set):
+        # Values for core, valence, and rydberg fields don't affect anything in gview.
+        dump_lines.append(f"{s} {i + 1} {charges[i]:.6f} {0.0:.6f} {0.0:.6f} {0.0:.6f} {_totals[i]:.6f}")
+    dump_lines += [
+        *[" " + "-" * 71, f"   * Total *   {np.sum(charges):.6f} " + f"{0.0:.6f} " * 3 + f"{np.sum(_totals):.6f}"],
+    ]
+    return dump_lines
+
+
+def _log_nbo(frame: Structure | Molecule, charge_key: str, spin_key: str) -> list[str]:
+    dump_lines = []
+    if charge_key in frame.site_properties:
+        nbo_charges = np.array(frame.site_properties[charge_key])
+        dump_lines += _log_nbo_fields(frame, nbo_charges)
+        if spin_key in frame.site_properties:  # NBO spins read as difference between alpha/beta totals
+            spins = np.array(frame.site_properties[spin_key])
+            alpha_spins = np.maximum(np.zeros(len(spins)), spins)
+            beta_spins = np.abs(np.minimum(np.zeros(len(spins)), spins))
+            for spin_type, spin_arr in zip(("Alpha", "Beta"), (alpha_spins, beta_spins), strict=False):
+                dump_lines += [f" {spin_type} spin orbitals "]
+                dump_lines += _log_nbo_fields(frame, np.zeros_like(spin_arr), totals=spin_arr)
+    return dump_lines
+
+
+def _log_forces(frame: Structure | Molecule, do_cell: bool = True) -> list[str]:
+    # TODO: Find the trick for getting Gaussview to read forces on PBC calculations
+    dump_lines = []
+    if "forces" in frame.site_properties:
+        forces = np.array(frame.site_properties["forces"])
+        dump_lines = [
+            *["-" * 67, " Center     Atomic" + " " * 19 + "Forces (Hartrees/Bohr)"],
+            *[" Number     " + (" " * 14).join(["Number", "X", "Y", "Z"]), "-" * 67],
+        ]
+        for i, n in enumerate(frame.atomic_numbers):
+            dump_lines.append(f" {i + 1} {n}\t" + "\t".join(f"{forces[i][j]:.9f}" for j in range(3)))
+        if do_cell:
+            dump_lines += [" " * 25 + "-2" + " " * 10 + "   ".join(f"{0.0:.9f}" for j in range(3)) for i in range(3)]
+        nforces = np.linalg.norm(forces, axis=1)
+        dump_lines += [" " + "-" * 67, f" Cartesian Forces:  Max {max(nforces):.9f} RMS {np.std(nforces):.9f}"]
+    return dump_lines