From f81561bdf5d3f3932632c0f6011be0ea65629712 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Thu, 29 May 2025 15:08:36 -0600 Subject: [PATCH 01/18] Function stub --- src/pymatgen/core/trajectory.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index f37106406ef..2d276233b36 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -477,6 +477,21 @@ def write_Xdatcar( with zopen(filename, mode="wt", encoding="utf-8") as file: file.write(xdatcar_str) # type:ignore[arg-type] + def write_Gaussian_log( + self, + filename: PathLike = "calc.log", + ) -> 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. + """ + # TODO: Write me! :) + print(f"Writing {filename} ...") + def as_dict(self) -> dict: """Return the trajectory as a MSONable dict.""" lat = self.lattice.tolist() if self.lattice is not None else None From 48f5baa0f58ef9411f391658890108e7b58df78b Mon Sep 17 00:00:00 2001 From: benrich37 Date: Thu, 29 May 2025 16:06:51 -0600 Subject: [PATCH 02/18] stashing changes --- src/pymatgen/core/trajectory.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 2d276233b36..377eddeac35 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -491,6 +491,8 @@ def write_Gaussian_log( """ # TODO: Write me! :) print(f"Writing {filename} ...") + # Ensure trajectory is in position form + self.to_positions() def as_dict(self) -> dict: """Return the trajectory as a MSONable dict.""" From a5f4608b9b4aede9cdfe7f720ba8d26263351787 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Fri, 30 May 2025 16:08:08 -0600 Subject: [PATCH 03/18] Transfer from existing work --- src/pymatgen/core/trajectory.py | 203 +++++++++++++++++++++++++++++++- 1 file changed, 201 insertions(+), 2 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 377eddeac35..422ad338024 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -477,9 +477,13 @@ 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. @@ -488,11 +492,43 @@ def write_Gaussian_log( 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 {`gview_mulliken_charges_key`: `charges`} will write the data stored in + `Trajectory[i].site_properties["charges"]` as Gaussian logs the Mulliken charges for each site. + The keys of implemented Gaussian properties are stored in `pymatgen.core.trajectory` as constants + in the form of `gview__key`, e.g. `gview_mulliken_charges_key`. The current implemented + Gaussian properties and the data format they expect are: + - `gview_mulliken_charges_key`: list[float] + - `gview_esp_charges_key`: list[float] + - `gview_nbo_charges_key`: list[float] + If None (default), the mapping will be set to fill all currently implemented Gaussian properties + with the "charges" property from the site properties. + write_lattice: Whether to write the lattice information in the log file. This will enable the logged + forces to be readable by Gaussview. """ - # TODO: Write me! :) - print(f"Writing {filename} ...") + if site_property_map is None: + # Default mapping fills all currently implemented Gaussian properties + site_property_map = { + gview_mulliken_charges_key: "charges", + gview_esp_charges_key: "charges", + gview_nbo_charges_key: "charges", + } + # 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 = [] + for fp in self.frame_properties: + if "energy" in fp: + energies.append(fp["energy"]) + else: + energies.append(0) + else: + energies = [0] * len(self) # Ensure trajectory is in position form self.to_positions() + lines = _traj_to_logx_lines(self, write_lattice, energies, site_property_map, is_done=True) + with open(filename, "w") as file: + file.write("\n".join(lines)) + file.close() def as_dict(self) -> dict: """Return the trajectory as a MSONable dict.""" @@ -894,3 +930,166 @@ def to_ase( temp_file.close() return ase_traj + + +_logx_init_lines = ["", " Entering Link 1 ", " "] +_logx_finish_lines = [" Normal termination of Gaussian 16"] +gview_mulliken_charges_key = "mulliken charges" +gview_esp_charges_key = "esp charges" +gview_nbo_charges_key = "nbo charges" + + +def _traj_to_logx_lines( + traj: Trajectory, do_cell, energies: list[float], site_property_map: dict, is_done: bool = False +): + dump_lines = _logx_init_lines + _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) + dump_lines += [f" SCF Done: E = {energies[i]}", "", ""] + if gview_mulliken_charges_key in site_property_map: + dump_lines += _log_mulliken_charges(traj[i], site_property_map[gview_mulliken_charges_key]) + if gview_esp_charges_key in site_property_map: + dump_lines += _log_esp_charges(traj[i], site_property_map[gview_esp_charges_key]) + if gview_nbo_charges_key in site_property_map: + dump_lines += _log_nbo_charges(traj[i], site_property_map[gview_nbo_charges_key]) + dump_lines += _log_forces(traj[i]) + dump_lines += _opt_spacer(i, len(traj)) + if is_done: + dump_lines += _log_input_orientation(traj[-1]) + dump_lines += _logx_finish_lines + return dump_lines + + +def _opt_spacer(i, nSteps): + dump_lines = [ + "", + " GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad", + "", + f" Step number {i + 1}", + ] + if i + 1 == nSteps: + dump_lines += [" Optimization completed.", " -- Stationary point found."] + dump_lines += ["", " GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad"] + return dump_lines + + +def _log_mulliken_charges(structure: Structure, site_property: str = "charges"): + dump_lines = [] + if site_property in structure.site_properties: + charges = np.array(structure.site_properties[site_property]) + nAtoms = len(charges) + symbols = [site.specie.symbol for site in structure.sites] + dump_lines = [ + " **********************************************************************", + "", + " Population analysis using the SCF Density.", + "", + " **********************************************************************", + "", + " Mulliken charges:", + " 1", + ] + for i in range(nAtoms): + dump_lines.append(f"{int(i + 1)} {symbols[i]} {charges[i]} ") + dump_lines.append(f" Sum of Mulliken charges = {np.sum(charges)}") + return dump_lines + + +def _log_esp_charges(structure: Structure, site_property: str = "charges"): + dump_lines = [] + if site_property in structure.site_properties: + posns = structure.cart_coords + symbols = [site.specie.symbol for site in structure.sites] + charges = np.array(structure.site_properties[site_property]) + dump_lines = [ + " **********************************************************************", + " Electrostatic Properties Using The SCF Density.", + " **********************************************************************", + ] + for i, posn in enumerate(posns): + dump_lines.append(f" Atomic Center {i + 1} is at {' '.join(f'{x:.6f}' for x in posn)}") + dump_lines += ["", " Charges from ESP fit", " ESP charges:", " 1"] + for i, symbol in enumerate(symbols): + dump_lines.append(f" {i + 1} {symbol} {charges[i]:.6f}") + charge = np.sum(charges) + dipole = np.zeros(3) + for i, posn in enumerate(posns): + dipole += posn * charges[i] + tot = np.linalg.norm(dipole) + dump_lines.append(f" Charge= {charge:.6f} Dipole= {' '.join(str(v) for v in dipole)} Tot= {tot}") + return dump_lines + + +def _log_nbo_charges(structure: Structure, site_property: str = "charges"): + dump_lines = [] + if site_property in structure.site_properties: + dump_lines = [ + " Summary of Natural Population Analysis: ", + "", + " Natural Population ", + " Natural ----------------------------------------------- ", + " Atom No Charge Core Valence Rydberg Total", + " -----------------------------------------------------------------------", + ] + nbo_charges = np.array(structure.site_properties[site_property]) + nAtoms = len(nbo_charges) + symbols = [site.specie.symbol for site in structure.sites] + for i in range(nAtoms): + dump_lines.append(f"{symbols[i]} {i + 1} {nbo_charges[i]} {0.0} {0.0} {0.0} {0.0}") + dump_lines.append(" =======================================================================") + dump_lines.append(f" * Total * {np.sum(nbo_charges)} 0.00000 0.00000 0.00000 0.00000") + return dump_lines + + +def _log_input_orientation(structure: Structure, do_cell=False): + dump_lines = [ + " Standard orientation: ", + " ---------------------------------------------------------------------", + " Center Atomic Atomic Coordinates (Angstroms)", + " Number Number Type X Y Z", + " ---------------------------------------------------------------------", + ] + at_ns = [site.specie.number for site in structure.sites] + at_posns = structure.cart_coords + nAtoms = len(at_ns) + for i in range(nAtoms): + dump_lines.append(f" {i + 1} {at_ns[i]} 0 ") + for j in range(3): + dump_lines[-1] += f"{at_posns[i][j]} " + # TODO: Gaussian also logs the a/b/c angles and lattice vector lengths, but I don't know + # if those are read by Gaussview. + if do_cell: + cell = structure.lattice.matrix + for i in range(3): + dump_lines.append(f"{i + nAtoms + 1} -2 0 ") + for j in range(3): + dump_lines[-1] += f"{cell[i][j]} " + dump_lines.append(" ---------------------------------------------------------------------") + return dump_lines + + +def _log_forces(structure: Structure): + dump_lines = [ + "-------------------------------------------------------------------", + " Center Atomic Forces (Hartrees/Bohr)", + " Number Number X Y Z", + " -------------------------------------------------------------------", + ] + forces = [] + if "forces" in structure.site_properties: + forces = structure.site_properties["forces"] + atomic_numbers = [site.specie.number for site in structure.sites] + for i, number in enumerate(atomic_numbers): + add_str = f" {i + 1} {number}" + force = forces[i] + for j in range(3): + add_str += f"\t{force[j]:.9f}" + dump_lines.append(add_str) + lattice_forces = np.zeros((3, 3)) + for ilat in range(3): + dump_lines.append(" -2 " + " ".join(f"{lattice_forces[ilat][j]:.9f}" for j in range(3))) + dump_lines.append(" -------------------------------------------------------------------") + aforces = np.array(forces) + nforces = np.linalg.norm(aforces, axis=1) + dump_lines.append(f" Cartesian Forces: Max {max(nforces):.9f} RMS {np.std(nforces):.9f}") + return dump_lines From 01e5710a22bbab4a7086f385a38c40f266d09a9c Mon Sep 17 00:00:00 2001 From: benrich37 Date: Fri, 30 May 2025 16:12:24 -0600 Subject: [PATCH 04/18] Testing expansion to `Molecule` for a frame --- src/pymatgen/core/trajectory.py | 54 ++++++++++++++++----------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 422ad338024..f4b3a88fc93 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -960,6 +960,33 @@ def _traj_to_logx_lines( return dump_lines +def _log_input_orientation(frame: Structure | Molecule, do_cell=False): + dump_lines = [ + " Standard orientation: ", + " ---------------------------------------------------------------------", + " Center Atomic Atomic Coordinates (Angstroms)", + " Number Number Type X Y Z", + " ---------------------------------------------------------------------", + ] + at_ns = [site.specie.number for site in frame.sites] + at_posns = frame.cart_coords + nAtoms = len(at_ns) + for i in range(nAtoms): + dump_lines.append(f" {i + 1} {at_ns[i]} 0 ") + for j in range(3): + dump_lines[-1] += f"{at_posns[i][j]} " + # TODO: Gaussian also logs the a/b/c angles and lattice vector lengths, but I don't know + # if those are read by Gaussview. + if do_cell: + cell = frame.lattice.matrix + for i in range(3): + dump_lines.append(f"{i + nAtoms + 1} -2 0 ") + for j in range(3): + dump_lines[-1] += f"{cell[i][j]} " + dump_lines.append(" ---------------------------------------------------------------------") + return dump_lines + + def _opt_spacer(i, nSteps): dump_lines = [ "", @@ -1041,33 +1068,6 @@ def _log_nbo_charges(structure: Structure, site_property: str = "charges"): return dump_lines -def _log_input_orientation(structure: Structure, do_cell=False): - dump_lines = [ - " Standard orientation: ", - " ---------------------------------------------------------------------", - " Center Atomic Atomic Coordinates (Angstroms)", - " Number Number Type X Y Z", - " ---------------------------------------------------------------------", - ] - at_ns = [site.specie.number for site in structure.sites] - at_posns = structure.cart_coords - nAtoms = len(at_ns) - for i in range(nAtoms): - dump_lines.append(f" {i + 1} {at_ns[i]} 0 ") - for j in range(3): - dump_lines[-1] += f"{at_posns[i][j]} " - # TODO: Gaussian also logs the a/b/c angles and lattice vector lengths, but I don't know - # if those are read by Gaussview. - if do_cell: - cell = structure.lattice.matrix - for i in range(3): - dump_lines.append(f"{i + nAtoms + 1} -2 0 ") - for j in range(3): - dump_lines[-1] += f"{cell[i][j]} " - dump_lines.append(" ---------------------------------------------------------------------") - return dump_lines - - def _log_forces(structure: Structure): dump_lines = [ "-------------------------------------------------------------------", From 257d4a3f01eef854cb0f43f266e06dcae6bb636d Mon Sep 17 00:00:00 2001 From: benrich37 Date: Fri, 30 May 2025 16:19:57 -0600 Subject: [PATCH 05/18] Expansion to `Molecule` for a frame --- src/pymatgen/core/trajectory.py | 90 ++++++++++++++++----------------- 1 file changed, 45 insertions(+), 45 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index f4b3a88fc93..5de70ecc2b2 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -932,17 +932,15 @@ def to_ase( return ase_traj -_logx_init_lines = ["", " Entering Link 1 ", " "] -_logx_finish_lines = [" Normal termination of Gaussian 16"] gview_mulliken_charges_key = "mulliken charges" gview_esp_charges_key = "esp charges" gview_nbo_charges_key = "nbo charges" def _traj_to_logx_lines( - traj: Trajectory, do_cell, energies: list[float], site_property_map: dict, is_done: bool = False -): - dump_lines = _logx_init_lines + _log_input_orientation(traj[0], do_cell=do_cell) + traj: Trajectory, do_cell: bool, energies: list[float], site_property_map: dict, is_done: bool = False +) -> list[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) dump_lines += [f" SCF Done: E = {energies[i]}", "", ""] @@ -956,11 +954,11 @@ def _traj_to_logx_lines( dump_lines += _opt_spacer(i, len(traj)) if is_done: dump_lines += _log_input_orientation(traj[-1]) - dump_lines += _logx_finish_lines + dump_lines += [" Normal termination of Gaussian 16"] return dump_lines -def _log_input_orientation(frame: Structure | Molecule, do_cell=False): +def _log_input_orientation(frame: Structure | Molecule, do_cell=False) -> list[str]: dump_lines = [ " Standard orientation: ", " ---------------------------------------------------------------------", @@ -987,7 +985,7 @@ def _log_input_orientation(frame: Structure | Molecule, do_cell=False): return dump_lines -def _opt_spacer(i, nSteps): +def _opt_spacer(i, nSteps) -> list[str]: dump_lines = [ "", " GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad", @@ -1000,12 +998,12 @@ def _opt_spacer(i, nSteps): return dump_lines -def _log_mulliken_charges(structure: Structure, site_property: str = "charges"): +def _log_mulliken_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: dump_lines = [] - if site_property in structure.site_properties: - charges = np.array(structure.site_properties[site_property]) + if site_property in frame.site_properties: + charges = np.array(frame.site_properties[site_property]) nAtoms = len(charges) - symbols = [site.specie.symbol for site in structure.sites] + symbols = [site.specie.symbol for site in frame.sites] dump_lines = [ " **********************************************************************", "", @@ -1022,12 +1020,12 @@ def _log_mulliken_charges(structure: Structure, site_property: str = "charges"): return dump_lines -def _log_esp_charges(structure: Structure, site_property: str = "charges"): +def _log_esp_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: dump_lines = [] - if site_property in structure.site_properties: - posns = structure.cart_coords - symbols = [site.specie.symbol for site in structure.sites] - charges = np.array(structure.site_properties[site_property]) + if site_property in frame.site_properties: + posns = frame.cart_coords + symbols = [site.specie.symbol for site in frame.sites] + charges = np.array(frame.site_properties[site_property]) dump_lines = [ " **********************************************************************", " Electrostatic Properties Using The SCF Density.", @@ -1047,9 +1045,9 @@ def _log_esp_charges(structure: Structure, site_property: str = "charges"): return dump_lines -def _log_nbo_charges(structure: Structure, site_property: str = "charges"): +def _log_nbo_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: dump_lines = [] - if site_property in structure.site_properties: + if site_property in frame.site_properties: dump_lines = [ " Summary of Natural Population Analysis: ", "", @@ -1058,9 +1056,9 @@ def _log_nbo_charges(structure: Structure, site_property: str = "charges"): " Atom No Charge Core Valence Rydberg Total", " -----------------------------------------------------------------------", ] - nbo_charges = np.array(structure.site_properties[site_property]) + nbo_charges = np.array(frame.site_properties[site_property]) nAtoms = len(nbo_charges) - symbols = [site.specie.symbol for site in structure.sites] + symbols = [site.specie.symbol for site in frame.sites] for i in range(nAtoms): dump_lines.append(f"{symbols[i]} {i + 1} {nbo_charges[i]} {0.0} {0.0} {0.0} {0.0}") dump_lines.append(" =======================================================================") @@ -1068,28 +1066,30 @@ def _log_nbo_charges(structure: Structure, site_property: str = "charges"): return dump_lines -def _log_forces(structure: Structure): - dump_lines = [ - "-------------------------------------------------------------------", - " Center Atomic Forces (Hartrees/Bohr)", - " Number Number X Y Z", - " -------------------------------------------------------------------", - ] - forces = [] - if "forces" in structure.site_properties: - forces = structure.site_properties["forces"] - atomic_numbers = [site.specie.number for site in structure.sites] - for i, number in enumerate(atomic_numbers): - add_str = f" {i + 1} {number}" - force = forces[i] - for j in range(3): - add_str += f"\t{force[j]:.9f}" - dump_lines.append(add_str) - lattice_forces = np.zeros((3, 3)) - for ilat in range(3): - dump_lines.append(" -2 " + " ".join(f"{lattice_forces[ilat][j]:.9f}" for j in range(3))) - dump_lines.append(" -------------------------------------------------------------------") - aforces = np.array(forces) - nforces = np.linalg.norm(aforces, axis=1) - dump_lines.append(f" Cartesian Forces: Max {max(nforces):.9f} RMS {np.std(nforces):.9f}") +def _log_forces(frame: Structure | Molecule) -> list[str]: + dump_lines = [] + if "forces" in frame.site_properties: + dump_lines = [ + "-------------------------------------------------------------------", + " Center Atomic Forces (Hartrees/Bohr)", + " Number Number X Y Z", + " -------------------------------------------------------------------", + ] + forces = frame.site_properties["forces"] + atomic_numbers = [site.specie.number for site in frame.sites] + for i, number in enumerate(atomic_numbers): + add_str = f" {i + 1} {number}" + force = forces[i] + for j in range(3): + add_str += f"\t{force[j]:.9f}" + dump_lines.append(add_str) + lattice_forces = np.zeros((3, 3)) + for ilat in range(3): + dump_lines.append( + " -2 " + " ".join(f"{lattice_forces[ilat][j]:.9f}" for j in range(3)) + ) + dump_lines.append(" -------------------------------------------------------------------") + aforces = np.array(forces) + nforces = np.linalg.norm(aforces, axis=1) + dump_lines.append(f" Cartesian Forces: Max {max(nforces):.9f} RMS {np.std(nforces):.9f}") return dump_lines From 431f653f2ef7ed74017c1da546e99dba96aa6849 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Sat, 31 May 2025 16:40:31 -0600 Subject: [PATCH 06/18] preliminary spin dumping --- src/pymatgen/core/trajectory.py | 95 +++++++++++++++++++++++++++------ 1 file changed, 79 insertions(+), 16 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 5de70ecc2b2..0347bd683be 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -28,6 +28,7 @@ from collections.abc import Iterator, Sequence from typing import Any + from numpy.typing import NDArray from typing_extensions import Self from pymatgen.util.typing import PathLike, SitePropsType @@ -998,7 +999,9 @@ def _opt_spacer(i, nSteps) -> list[str]: return dump_lines -def _log_mulliken_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: +def _log_mulliken_charges( + frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" +) -> list[str]: dump_lines = [] if site_property in frame.site_properties: charges = np.array(frame.site_properties[site_property]) @@ -1014,9 +1017,18 @@ def _log_mulliken_charges(frame: Structure | Molecule, site_property: str = "cha " Mulliken charges:", " 1", ] + spins: None | NDArray = None + if spin_site_property in frame.site_properties: + spins = np.array(frame.site_properties[spin_site_property]) + dump_lines[-2].replace(":", " and spin densities:") + dump_lines[-1] += " 2" for i in range(nAtoms): dump_lines.append(f"{int(i + 1)} {symbols[i]} {charges[i]} ") + if spins is not None: + dump_lines[-1] += f"{spins[i]} " dump_lines.append(f" Sum of Mulliken charges = {np.sum(charges)}") + if spins is not None: + dump_lines[-1] += f" {np.sum(spins)}" return dump_lines @@ -1045,28 +1057,79 @@ def _log_esp_charges(frame: Structure | Molecule, site_property: str = "charges" return dump_lines -def _log_nbo_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: +_gview_nbo_header_lines = [ + " Summary of Natural Population Analysis: ", + "", + " Natural Population ", + " Natural ----------------------------------------------- ", + " Atom No Charge Core Valence Rydberg Total", + " -----------------------------------------------------------------------", +] + + +def _gview_nbo_spin_title_lines(spin: str): + return [ + "***************************************************", + f"******* {spin} spin orbitals *******", + "***************************************************", + ] + + +def _gview_log_nbo_fields( + frame: Structure | Molecule, + charges: NDArray, + cores: NDArray | None = None, + valences: NDArray | None = None, + rydbergs: NDArray | None = None, + totals: NDArray | None = None, +) -> list[str]: + # As far as I can tell, cores, valences, and rydbergs don't affect anything in gview. + if cores is None: + cores = np.zeros_like(charges) + if valences is None: + valences = np.zeros_like(charges) + if rydbergs is None: + rydbergs = np.zeros_like(charges) + if totals is None: + totals = cores + valences + rydbergs + dump_lines = _gview_nbo_header_lines.copy() + symbols = [site.specie.symbol for site in frame.sites] + for i, symbol in enumerate(symbols): + dump_lines.append( + f"{symbol} {i + 1} {charges[i]:.6f} {cores[i]:.6f} {valences[i]:.6f} {rydbergs[i]:.6f} {totals[i]:.6f}" + ) + dump_lines.append(" =======================================================================") + dump_lines.append( + " * Total * " + " ".join(f"{np.sum(arr):.6f}" for arr in (charges, cores, valences, rydbergs, totals)) + ) + return dump_lines + + +def _log_nbo_charges( + frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" +) -> list[str]: dump_lines = [] if site_property in frame.site_properties: - dump_lines = [ - " Summary of Natural Population Analysis: ", - "", - " Natural Population ", - " Natural ----------------------------------------------- ", - " Atom No Charge Core Valence Rydberg Total", - " -----------------------------------------------------------------------", - ] nbo_charges = np.array(frame.site_properties[site_property]) - nAtoms = len(nbo_charges) - symbols = [site.specie.symbol for site in frame.sites] - for i in range(nAtoms): - dump_lines.append(f"{symbols[i]} {i + 1} {nbo_charges[i]} {0.0} {0.0} {0.0} {0.0}") - dump_lines.append(" =======================================================================") - dump_lines.append(f" * Total * {np.sum(nbo_charges)} 0.00000 0.00000 0.00000 0.00000") + dump_lines += _gview_log_nbo_fields(frame, nbo_charges) + spins: None | NDArray = None + if spin_site_property in frame.site_properties: + spins = np.array(frame.site_properties[spin_site_property]) + # NBO spins are evaluated as the difference between the "total" nbo fields for alpha and beta densities. + 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 += _gview_nbo_spin_title_lines(spin_type) + dump_lines += _gview_log_nbo_fields( + frame, np.zeros_like(spin_arr), cores=None, valences=None, rydbergs=None, totals=spin_arr + ) return dump_lines def _log_forces(frame: Structure | Molecule) -> list[str]: + # Gaussian calculations with periodic boundary conditions will log forces on lattice vectors, but these forces + # do not appear to be read by Gaussview, and writing forces for lattice vectors does not seem to fix the issue + # of atomic forces not being read by Gaussview (for log files with lattice vectors). dump_lines = [] if "forces" in frame.site_properties: dump_lines = [ From a92dc3ff900f89e11fc2ffb58681abfac67d5ada Mon Sep 17 00:00:00 2001 From: benrich37 Date: Sat, 31 May 2025 17:05:35 -0600 Subject: [PATCH 07/18] Fixed improper use of `str.replace` --- src/pymatgen/core/trajectory.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 0347bd683be..e953e3c1ac4 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -1020,7 +1020,7 @@ def _log_mulliken_charges( spins: None | NDArray = None if spin_site_property in frame.site_properties: spins = np.array(frame.site_properties[spin_site_property]) - dump_lines[-2].replace(":", " and spin densities:") + dump_lines[-2] = dump_lines[-2].replace(":", " and spin densities:") dump_lines[-1] += " 2" for i in range(nAtoms): dump_lines.append(f"{int(i + 1)} {symbols[i]} {charges[i]} ") From bd131a767ddf7e5dd676521ac6f0e93888e0a4cd Mon Sep 17 00:00:00 2001 From: benrich37 Date: Sat, 31 May 2025 17:13:49 -0600 Subject: [PATCH 08/18] Removing unrequired lines for esp --- src/pymatgen/core/trajectory.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index e953e3c1ac4..d6a86f54347 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -1038,13 +1038,6 @@ def _log_esp_charges(frame: Structure | Molecule, site_property: str = "charges" posns = frame.cart_coords symbols = [site.specie.symbol for site in frame.sites] charges = np.array(frame.site_properties[site_property]) - dump_lines = [ - " **********************************************************************", - " Electrostatic Properties Using The SCF Density.", - " **********************************************************************", - ] - for i, posn in enumerate(posns): - dump_lines.append(f" Atomic Center {i + 1} is at {' '.join(f'{x:.6f}' for x in posn)}") dump_lines += ["", " Charges from ESP fit", " ESP charges:", " 1"] for i, symbol in enumerate(symbols): dump_lines.append(f" {i + 1} {symbol} {charges[i]:.6f}") From 9014e70755e1b9464c00041e2ce6a7b47ba09722 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Sat, 31 May 2025 17:17:32 -0600 Subject: [PATCH 09/18] Removing unrequired lines and whitespace for mulliken logging --- src/pymatgen/core/trajectory.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index d6a86f54347..6b1a7b43fda 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -1008,20 +1008,14 @@ def _log_mulliken_charges( nAtoms = len(charges) symbols = [site.specie.symbol for site in frame.sites] dump_lines = [ - " **********************************************************************", - "", - " Population analysis using the SCF Density.", - "", - " **********************************************************************", - "", - " Mulliken charges:", - " 1", + "Mulliken charges:", + "1", ] spins: None | NDArray = None if spin_site_property in frame.site_properties: spins = np.array(frame.site_properties[spin_site_property]) dump_lines[-2] = dump_lines[-2].replace(":", " and spin densities:") - dump_lines[-1] += " 2" + dump_lines[-1] += " 2" for i in range(nAtoms): dump_lines.append(f"{int(i + 1)} {symbols[i]} {charges[i]} ") if spins is not None: From 1051362c2c5a2e2e72f4843cad6a49ceee87821d Mon Sep 17 00:00:00 2001 From: benrich37 Date: Sat, 31 May 2025 17:26:25 -0600 Subject: [PATCH 10/18] Removing unneeded lines for nbo --- src/pymatgen/core/trajectory.py | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 6b1a7b43fda..8f13de38a61 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -1054,14 +1054,6 @@ def _log_esp_charges(frame: Structure | Molecule, site_property: str = "charges" ] -def _gview_nbo_spin_title_lines(spin: str): - return [ - "***************************************************", - f"******* {spin} spin orbitals *******", - "***************************************************", - ] - - def _gview_log_nbo_fields( frame: Structure | Molecule, charges: NDArray, @@ -1106,7 +1098,7 @@ def _log_nbo_charges( 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 += _gview_nbo_spin_title_lines(spin_type) + dump_lines += [f" {spin_type} spin orbitals "] dump_lines += _gview_log_nbo_fields( frame, np.zeros_like(spin_arr), cores=None, valences=None, rydbergs=None, totals=spin_arr ) From d4f5e65454ed0ef86f9e31e92d45e573bd720914 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 18 Jun 2025 18:52:03 -0600 Subject: [PATCH 11/18] Moving helper functions to `io.gaussian`, working to rewrite helper functions in as few lines as possible --- src/pymatgen/core/trajectory.py | 220 ++------------------------------ src/pymatgen/io/gaussian.py | 157 +++++++++++++++++++++++ 2 files changed, 165 insertions(+), 212 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 8f13de38a61..f1285890903 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -17,6 +17,7 @@ from pymatgen.core.structure import Composition, DummySpecies, Element, Lattice, Molecule, Species, Structure from pymatgen.io.ase import NO_ASE_ERR, AseAtomsAdaptor +from pymatgen.io.gaussian import _traj_to_gaussian_log_lines if NO_ASE_ERR is None: from ase.io.trajectory import Trajectory as AseTrajectory @@ -28,7 +29,6 @@ from collections.abc import Iterator, Sequence from typing import Any - from numpy.typing import NDArray from typing_extensions import Self from pymatgen.util.typing import PathLike, SitePropsType @@ -500,8 +500,8 @@ def write_Gaussian_log( in the form of `gview__key`, e.g. `gview_mulliken_charges_key`. The current implemented Gaussian properties and the data format they expect are: - `gview_mulliken_charges_key`: list[float] - - `gview_esp_charges_key`: list[float] - - `gview_nbo_charges_key`: list[float] + - ` "esp charges"`: list[float] + - `"nbo charges"`: list[float] If None (default), the mapping will be set to fill all currently implemented Gaussian properties with the "charges" property from the site properties. write_lattice: Whether to write the lattice information in the log file. This will enable the logged @@ -510,9 +510,9 @@ def write_Gaussian_log( if site_property_map is None: # Default mapping fills all currently implemented Gaussian properties site_property_map = { - gview_mulliken_charges_key: "charges", - gview_esp_charges_key: "charges", - gview_nbo_charges_key: "charges", + "mulliken charges": "charges", + "esp charges": "charges", + "nbo charges": "charges", } # Each frame must have an energy to be read properly by Gaussview. if hasattr(self, "frame_properties") and self.frame_properties is not None: @@ -526,9 +526,9 @@ def write_Gaussian_log( energies = [0] * len(self) # Ensure trajectory is in position form self.to_positions() - lines = _traj_to_logx_lines(self, write_lattice, energies, site_property_map, is_done=True) + lines = _traj_to_gaussian_log_lines(self, write_lattice, energies, site_property_map) with open(filename, "w") as file: - file.write("\n".join(lines)) + file.write(lines) file.close() def as_dict(self) -> dict: @@ -931,207 +931,3 @@ def to_ase( temp_file.close() return ase_traj - - -gview_mulliken_charges_key = "mulliken charges" -gview_esp_charges_key = "esp charges" -gview_nbo_charges_key = "nbo charges" - - -def _traj_to_logx_lines( - traj: Trajectory, do_cell: bool, energies: list[float], site_property_map: dict, is_done: bool = False -) -> list[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) - dump_lines += [f" SCF Done: E = {energies[i]}", "", ""] - if gview_mulliken_charges_key in site_property_map: - dump_lines += _log_mulliken_charges(traj[i], site_property_map[gview_mulliken_charges_key]) - if gview_esp_charges_key in site_property_map: - dump_lines += _log_esp_charges(traj[i], site_property_map[gview_esp_charges_key]) - if gview_nbo_charges_key in site_property_map: - dump_lines += _log_nbo_charges(traj[i], site_property_map[gview_nbo_charges_key]) - dump_lines += _log_forces(traj[i]) - dump_lines += _opt_spacer(i, len(traj)) - if is_done: - dump_lines += _log_input_orientation(traj[-1]) - dump_lines += [" Normal termination of Gaussian 16"] - return dump_lines - - -def _log_input_orientation(frame: Structure | Molecule, do_cell=False) -> list[str]: - dump_lines = [ - " Standard orientation: ", - " ---------------------------------------------------------------------", - " Center Atomic Atomic Coordinates (Angstroms)", - " Number Number Type X Y Z", - " ---------------------------------------------------------------------", - ] - at_ns = [site.specie.number for site in frame.sites] - at_posns = frame.cart_coords - nAtoms = len(at_ns) - for i in range(nAtoms): - dump_lines.append(f" {i + 1} {at_ns[i]} 0 ") - for j in range(3): - dump_lines[-1] += f"{at_posns[i][j]} " - # TODO: Gaussian also logs the a/b/c angles and lattice vector lengths, but I don't know - # if those are read by Gaussview. - if do_cell: - cell = frame.lattice.matrix - for i in range(3): - dump_lines.append(f"{i + nAtoms + 1} -2 0 ") - for j in range(3): - dump_lines[-1] += f"{cell[i][j]} " - dump_lines.append(" ---------------------------------------------------------------------") - return dump_lines - - -def _opt_spacer(i, nSteps) -> list[str]: - dump_lines = [ - "", - " GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad", - "", - f" Step number {i + 1}", - ] - if i + 1 == nSteps: - dump_lines += [" Optimization completed.", " -- Stationary point found."] - dump_lines += ["", " GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad"] - return dump_lines - - -def _log_mulliken_charges( - frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" -) -> list[str]: - dump_lines = [] - if site_property in frame.site_properties: - charges = np.array(frame.site_properties[site_property]) - nAtoms = len(charges) - symbols = [site.specie.symbol for site in frame.sites] - dump_lines = [ - "Mulliken charges:", - "1", - ] - spins: None | NDArray = None - if spin_site_property in frame.site_properties: - spins = np.array(frame.site_properties[spin_site_property]) - dump_lines[-2] = dump_lines[-2].replace(":", " and spin densities:") - dump_lines[-1] += " 2" - for i in range(nAtoms): - dump_lines.append(f"{int(i + 1)} {symbols[i]} {charges[i]} ") - if spins is not None: - dump_lines[-1] += f"{spins[i]} " - dump_lines.append(f" Sum of Mulliken charges = {np.sum(charges)}") - if spins is not None: - dump_lines[-1] += f" {np.sum(spins)}" - return dump_lines - - -def _log_esp_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: - dump_lines = [] - if site_property in frame.site_properties: - posns = frame.cart_coords - symbols = [site.specie.symbol for site in frame.sites] - charges = np.array(frame.site_properties[site_property]) - dump_lines += ["", " Charges from ESP fit", " ESP charges:", " 1"] - for i, symbol in enumerate(symbols): - dump_lines.append(f" {i + 1} {symbol} {charges[i]:.6f}") - charge = np.sum(charges) - dipole = np.zeros(3) - for i, posn in enumerate(posns): - dipole += posn * charges[i] - tot = np.linalg.norm(dipole) - dump_lines.append(f" Charge= {charge:.6f} Dipole= {' '.join(str(v) for v in dipole)} Tot= {tot}") - return dump_lines - - -_gview_nbo_header_lines = [ - " Summary of Natural Population Analysis: ", - "", - " Natural Population ", - " Natural ----------------------------------------------- ", - " Atom No Charge Core Valence Rydberg Total", - " -----------------------------------------------------------------------", -] - - -def _gview_log_nbo_fields( - frame: Structure | Molecule, - charges: NDArray, - cores: NDArray | None = None, - valences: NDArray | None = None, - rydbergs: NDArray | None = None, - totals: NDArray | None = None, -) -> list[str]: - # As far as I can tell, cores, valences, and rydbergs don't affect anything in gview. - if cores is None: - cores = np.zeros_like(charges) - if valences is None: - valences = np.zeros_like(charges) - if rydbergs is None: - rydbergs = np.zeros_like(charges) - if totals is None: - totals = cores + valences + rydbergs - dump_lines = _gview_nbo_header_lines.copy() - symbols = [site.specie.symbol for site in frame.sites] - for i, symbol in enumerate(symbols): - dump_lines.append( - f"{symbol} {i + 1} {charges[i]:.6f} {cores[i]:.6f} {valences[i]:.6f} {rydbergs[i]:.6f} {totals[i]:.6f}" - ) - dump_lines.append(" =======================================================================") - dump_lines.append( - " * Total * " + " ".join(f"{np.sum(arr):.6f}" for arr in (charges, cores, valences, rydbergs, totals)) - ) - return dump_lines - - -def _log_nbo_charges( - frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" -) -> list[str]: - dump_lines = [] - if site_property in frame.site_properties: - nbo_charges = np.array(frame.site_properties[site_property]) - dump_lines += _gview_log_nbo_fields(frame, nbo_charges) - spins: None | NDArray = None - if spin_site_property in frame.site_properties: - spins = np.array(frame.site_properties[spin_site_property]) - # NBO spins are evaluated as the difference between the "total" nbo fields for alpha and beta densities. - 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 += _gview_log_nbo_fields( - frame, np.zeros_like(spin_arr), cores=None, valences=None, rydbergs=None, totals=spin_arr - ) - return dump_lines - - -def _log_forces(frame: Structure | Molecule) -> list[str]: - # Gaussian calculations with periodic boundary conditions will log forces on lattice vectors, but these forces - # do not appear to be read by Gaussview, and writing forces for lattice vectors does not seem to fix the issue - # of atomic forces not being read by Gaussview (for log files with lattice vectors). - dump_lines = [] - if "forces" in frame.site_properties: - dump_lines = [ - "-------------------------------------------------------------------", - " Center Atomic Forces (Hartrees/Bohr)", - " Number Number X Y Z", - " -------------------------------------------------------------------", - ] - forces = frame.site_properties["forces"] - atomic_numbers = [site.specie.number for site in frame.sites] - for i, number in enumerate(atomic_numbers): - add_str = f" {i + 1} {number}" - force = forces[i] - for j in range(3): - add_str += f"\t{force[j]:.9f}" - dump_lines.append(add_str) - lattice_forces = np.zeros((3, 3)) - for ilat in range(3): - dump_lines.append( - " -2 " + " ".join(f"{lattice_forces[ilat][j]:.9f}" for j in range(3)) - ) - dump_lines.append(" -------------------------------------------------------------------") - aforces = np.array(forces) - nforces = np.linalg.norm(aforces, axis=1) - dump_lines.append(f" Cartesian Forces: Max {max(nforces):.9f} RMS {np.std(nforces):.9f}") - return dump_lines diff --git a/src/pymatgen/io/gaussian.py b/src/pymatgen/io/gaussian.py index 861033b6392..c8736193c91 100644 --- a/src/pymatgen/io/gaussian.py +++ b/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" @@ -1333,3 +1336,157 @@ def to_input( link0_parameters=link0_parameters, dieze_tag=dieze_tag, ) + + +def _traj_to_gaussian_log_lines(traj: Trajectory, do_cell: bool, energies: list[float], site_property_map: dict) -> str: + _site_property_map = { + "mulliken charges": "charges", + "mulliken spin": "magmom", + "nbo charges": "charges", + "nbo spin": "magmom", + "esp charges": "charges", + } + _site_property_map.update(site_property_map) + 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), + *_log_mulliken_fields(traj[i], _site_property_map["mulliken charges"], _site_property_map["mulliken spin"]), + *_log_esp_charges(traj[i], _site_property_map["esp charges"]), + *_log_nbo_fields(traj[i], _site_property_map["nbo charges"], _site_property_map["nbo 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=False) -> list[str]: + dump_lines = [ + " " * 24 + "Standard orientation:" + " " * 26, + " " + "-" * 69, + " Center Atomic Atomic Coordinates (Angstroms)", + " Number Number Type X Y Z", + " " + "-" * 69, + ] + at_ns = [site.specie.number for site in frame.sites] + at_posns = frame.cart_coords + nAtoms = len(at_ns) + for i in range(nAtoms): + dump_lines.append(f" {i + 1} {at_ns[i]} 0 " + " ".join(f"{at_posns[i][j]:.6f}" for j in range(3)) + " ") + if do_cell: + cell = frame.lattice.matrix + for i in range(3): + dump_lines.append(f"{i + nAtoms + 1} -2 0 {cell[i][0]:.6f} {cell[i][1]:.6f} {cell[i][2]:.6f} ") + dump_lines.append(" ---------------------------------------------------------------------") + return dump_lines + + +def _log_mulliken_fields( + frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" +) -> list[str]: + dump_lines = [] + if site_property in frame.site_properties: + charges = np.array(frame.site_properties[site_property]) + dump_lines += ["Mulliken charges:", "1"] + spins: None | NDArray = None + if spin_site_property in frame.site_properties: + spins = np.array(frame.site_properties[spin_site_property]) + dump_lines[-2] = dump_lines[-2].replace(":", " and spin densities:") + dump_lines[-1] += " 2" + for i, charge in enumerate(charges): + dump_lines.append(f"{int(i + 1)} {frame.sites[i].specie.symbol} {charge} ") + if spins is not None: + dump_lines[-1] += f"{spins[i]} " + dump_lines.append(f" Sum of Mulliken charges = {np.sum(charges)}") + if spins is not None: + dump_lines[-1] += f" {np.sum(spins)}" + return dump_lines + + +def _log_esp_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: + dump_lines = [] + if site_property in frame.site_properties: + charges = np.array(frame.site_properties[site_property]) + dump_lines += ["", " Charges from ESP fit", " ESP charges:", " 1"] + for i, symbol in enumerate([site.specie.symbol for site in frame.sites]): + dump_lines.append(f" {i + 1} {symbol} {charges[i]:.6f}") + dipole = np.sum(frame.cart_coords * charges[:, np.newaxis], axis=0) + tot = np.linalg.norm(dipole) + dump_lines.append(f" Charge= {np.sum(charges):.6f} Dipole= {' '.join(str(v) for v in dipole)} Tot= {tot}") + return dump_lines + + +def _gview_log_nbo_fields( + frame: Structure | Molecule, + charges: NDArray, + totals: NDArray | None = None, +) -> list[str]: + dump_lines = [ + " Summary of Natural Population Analysis: ", + "", + " " * 39 + "Natural Population ", + " " * 17 + "Natural " + "-" * 47 + " ", + " Atom No Charge Core Valence Rydberg Total", + " " + "-" * 71, + ] + totals = np.zeros_like(charges) if totals is None else totals + for i, symbol in enumerate([site.specie.symbol for site in frame.sites]): + dump_lines.append( + f"{symbol} {i + 1} {charges[i]:.6f} {0.0:.6f} {0.0:.6f} {0.0:.6f} {totals[i]:.6f}" + ) # Values for core, valence, and rydberg fields don't affect anything in gview. + 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_fields( + frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" +) -> list[str]: + dump_lines = [] + if site_property in frame.site_properties: + nbo_charges = np.array(frame.site_properties[site_property]) + dump_lines += _gview_log_nbo_fields(frame, nbo_charges) + spins: None | NDArray = None + if spin_site_property in frame.site_properties: + spins = np.array(frame.site_properties[spin_site_property]) + # NBO spins are evaluated as the difference between the "total" nbo fields for alpha and beta densities. + 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 += _gview_log_nbo_fields(frame, np.zeros_like(spin_arr), totals=spin_arr) + return dump_lines + + +def _log_forces(frame: Structure | Molecule) -> list[str]: + # Gaussian calculations with periodic boundary conditions will log forces on lattice vectors, but these forces + # do not appear to be read by Gaussview, and writing forces for lattice vectors does not seem to fix the issue + # of atomic forces not being read by Gaussview (for log files with lattice vectors). + dump_lines = [] + if "forces" in frame.site_properties: + dump_lines += [ + "-" * 67, + " Center Atomic" + " " * 19 + "Forces (Hartrees/Bohr)", + " Number " + (" " * 14).join(["Number", "X", "Y", "Z"]), + "-" * 67, + ] + forces = frame.site_properties["forces"] + atomic_numbers = [site.specie.number for site in frame.sites] + for i, number in enumerate(atomic_numbers): + dump_lines.append(f" {i + 1} {number}\t" + "\t".join(f"{forces[i][j]:.9f}" for j in range(3))) + lat_forces = np.zeros((3, 3)) + dump_lines += [ + " " * 25 + "-2" + " " * 10 + " ".join(f"{lat_forces[i][j]:.9f}" for j in range(3)) for i in range(3) + ] + nforces = np.linalg.norm(np.array(forces), axis=1) + dump_lines += [" " + "-" * 67, f" Cartesian Forces: Max {max(nforces):.9f} RMS {np.std(nforces):.9f}"] + return dump_lines From ff1638d9e9f38bad20d72b00b8c19257c7cc5a56 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 18 Jun 2025 18:53:18 -0600 Subject: [PATCH 12/18] working to rewrite helper functions in as few lines as possible --- src/pymatgen/io/gaussian.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/pymatgen/io/gaussian.py b/src/pymatgen/io/gaussian.py index c8736193c91..1e955aaaac9 100644 --- a/src/pymatgen/io/gaussian.py +++ b/src/pymatgen/io/gaussian.py @@ -1355,12 +1355,7 @@ def _traj_to_gaussian_log_lines(traj: Trajectory, do_cell: bool, energies: list[ *_log_esp_charges(traj[i], _site_property_map["esp charges"]), *_log_nbo_fields(traj[i], _site_property_map["nbo charges"], _site_property_map["nbo spin"]), *_log_forces(traj[i]), - "", - " " + "Grad" * 18, - "", - f" Step number {i + 1}", - "", - " " + "Grad" * 18, + *["", " " + "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"] From bc492456b5542ef3e417fde3c75e7eb2d41e32ef Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 18 Jun 2025 19:00:41 -0600 Subject: [PATCH 13/18] working to rewrite helper functions in as few lines as possible --- src/pymatgen/io/gaussian.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/pymatgen/io/gaussian.py b/src/pymatgen/io/gaussian.py index 1e955aaaac9..7e4e56c60a4 100644 --- a/src/pymatgen/io/gaussian.py +++ b/src/pymatgen/io/gaussian.py @@ -1371,14 +1371,11 @@ def _log_input_orientation(frame: Structure | Molecule, do_cell=False) -> list[s " " + "-" * 69, ] at_ns = [site.specie.number for site in frame.sites] - at_posns = frame.cart_coords - nAtoms = len(at_ns) - for i in range(nAtoms): - dump_lines.append(f" {i + 1} {at_ns[i]} 0 " + " ".join(f"{at_posns[i][j]:.6f}" for j in range(3)) + " ") + dump_lines += [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 for i in range(3): - dump_lines.append(f"{i + nAtoms + 1} -2 0 {cell[i][0]:.6f} {cell[i][1]:.6f} {cell[i][2]:.6f} ") + dump_lines.append(f"{i + len(at_ns) + 1} -2 0 {cell[i][0]:.6f} {cell[i][1]:.6f} {cell[i][2]:.6f} ") dump_lines.append(" ---------------------------------------------------------------------") return dump_lines From 8ceb7ab52fe56dc11d2c9c75c262549b34d4f91e Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 18 Jun 2025 20:19:24 -0600 Subject: [PATCH 14/18] working to rewrite helper functions in as few lines as possible --- src/pymatgen/core/trajectory.py | 37 +++------- src/pymatgen/io/gaussian.py | 117 +++++++++++++------------------- 2 files changed, 56 insertions(+), 98 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index f1285890903..472bf00ea49 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -17,7 +17,7 @@ from pymatgen.core.structure import Composition, DummySpecies, Element, Lattice, Molecule, Species, Structure from pymatgen.io.ase import NO_ASE_ERR, AseAtomsAdaptor -from pymatgen.io.gaussian import _traj_to_gaussian_log_lines +from pymatgen.io.gaussian import _traj_to_gaussian_log if NO_ASE_ERR is None: from ase.io.trajectory import Trajectory as AseTrajectory @@ -494,39 +494,22 @@ def write_Gaussian_log( 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 {`gview_mulliken_charges_key`: `charges`} will write the data stored in + 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 keys of implemented Gaussian properties are stored in `pymatgen.core.trajectory` as constants - in the form of `gview__key`, e.g. `gview_mulliken_charges_key`. The current implemented - Gaussian properties and the data format they expect are: - - `gview_mulliken_charges_key`: list[float] - - ` "esp charges"`: list[float] - - `"nbo charges"`: list[float] - If None (default), the mapping will be set to fill all currently implemented Gaussian properties - with the "charges" property from the site properties. - write_lattice: Whether to write the lattice information in the log file. This will enable the logged - forces to be readable by Gaussview. + 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. """ - if site_property_map is None: - # Default mapping fills all currently implemented Gaussian properties - site_property_map = { - "mulliken charges": "charges", - "esp charges": "charges", - "nbo charges": "charges", - } # 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 = [] - for fp in self.frame_properties: - if "energy" in fp: - energies.append(fp["energy"]) - else: - energies.append(0) + energies = [fp.get("energy", 0) for fp in self.frame_properties] else: energies = [0] * len(self) - # Ensure trajectory is in position form self.to_positions() - lines = _traj_to_gaussian_log_lines(self, write_lattice, energies, site_property_map) + lines = _traj_to_gaussian_log(self, write_lattice, energies, site_property_map) with open(filename, "w") as file: file.write(lines) file.close() diff --git a/src/pymatgen/io/gaussian.py b/src/pymatgen/io/gaussian.py index 7e4e56c60a4..064eafca95d 100644 --- a/src/pymatgen/io/gaussian.py +++ b/src/pymatgen/io/gaussian.py @@ -1338,22 +1338,23 @@ def to_input( ) -def _traj_to_gaussian_log_lines(traj: Trajectory, do_cell: bool, energies: list[float], site_property_map: dict) -> str: - _site_property_map = { +def _traj_to_gaussian_log(traj: Trajectory, do_cell: bool, energies: list[float], sp_map: dict | None) -> str: + _sp_map = { "mulliken charges": "charges", "mulliken spin": "magmom", "nbo charges": "charges", "nbo spin": "magmom", "esp charges": "charges", } - _site_property_map.update(site_property_map) + _sp_map.update(sp_map or {}) 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), - *_log_mulliken_fields(traj[i], _site_property_map["mulliken charges"], _site_property_map["mulliken spin"]), - *_log_esp_charges(traj[i], _site_property_map["esp charges"]), - *_log_nbo_fields(traj[i], _site_property_map["nbo charges"], _site_property_map["nbo spin"]), + f"SCF Done: E = {energies[i]:.6f} ", + *_log_mulliken(traj[i], _sp_map["mulliken charges"], _sp_map["mulliken spin"]), + *_log_esp(traj[i], _sp_map["esp charges"]), + *_log_nbo(traj[i], _sp_map["nbo charges"], _sp_map["nbo spin"]), *_log_forces(traj[i]), *["", " " + "Grad" * 18, "", f" Step number {i + 1}", "", " " + "Grad" * 18], ] @@ -1362,123 +1363,97 @@ def _traj_to_gaussian_log_lines(traj: Trajectory, do_cell: bool, energies: list[ return "\n".join(dump_lines) -def _log_input_orientation(frame: Structure | Molecule, do_cell=False) -> list[str]: +def _log_input_orientation(frame: Structure | Molecule, do_cell=True) -> list[str]: + at_ns = frame.atomic_numbers dump_lines = [ - " " * 24 + "Standard orientation:" + " " * 26, - " " + "-" * 69, + *[" " * 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)], ] - at_ns = [site.specie.number for site in frame.sites] - dump_lines += [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 - for i in range(3): - dump_lines.append(f"{i + len(at_ns) + 1} -2 0 {cell[i][0]:.6f} {cell[i][1]:.6f} {cell[i][2]:.6f} ") + 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_fields( - frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" -) -> list[str]: +def _log_mulliken(frame: Structure | Molecule, charge_key: str, spin_key: str) -> list[str]: dump_lines = [] - if site_property in frame.site_properties: - charges = np.array(frame.site_properties[site_property]) + 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_site_property in frame.site_properties: - spins = np.array(frame.site_properties[spin_site_property]) + 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"{int(i + 1)} {frame.sites[i].specie.symbol} {charge} ") - if spins is not None: - dump_lines[-1] += f"{spins[i]} " - dump_lines.append(f" Sum of Mulliken charges = {np.sum(charges)}") - if spins is not None: - dump_lines[-1] += f" {np.sum(spins)}" + 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_charges(frame: Structure | Molecule, site_property: str = "charges") -> list[str]: +def _log_esp(frame: Structure | Molecule, charge_key: str) -> list[str]: dump_lines = [] - if site_property in frame.site_properties: - charges = np.array(frame.site_properties[site_property]) + if charge_key in frame.site_properties: + charges = np.array(frame.site_properties[charge_key]) dump_lines += ["", " Charges from ESP fit", " ESP charges:", " 1"] - for i, symbol in enumerate([site.specie.symbol for site in frame.sites]): - dump_lines.append(f" {i + 1} {symbol} {charges[i]:.6f}") + 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):.6f} Dipole= {' '.join(str(v) for v in dipole)} Tot= {tot}") return dump_lines -def _gview_log_nbo_fields( - frame: Structure | Molecule, - charges: NDArray, - totals: NDArray | None = None, -) -> list[str]: +def _log_nbo_fields(frame: Structure | Molecule, charges: NDArray, totals: NDArray | None = None) -> list[str]: dump_lines = [ - " Summary of Natural Population Analysis: ", - "", - " " * 39 + "Natural Population ", + *[" Summary of Natural Population Analysis: ", "", " " * 39 + "Natural Population "], " " * 17 + "Natural " + "-" * 47 + " ", - " Atom No Charge Core Valence Rydberg Total", - " " + "-" * 71, + *[" Atom No Charge Core Valence Rydberg Total", " " + "-" * 71], ] - totals = np.zeros_like(charges) if totals is None else totals - for i, symbol in enumerate([site.specie.symbol for site in frame.sites]): - dump_lines.append( - f"{symbol} {i + 1} {charges[i]:.6f} {0.0:.6f} {0.0:.6f} {0.0:.6f} {totals[i]:.6f}" - ) # Values for core, valence, and rydberg fields don't affect anything in gview. + _totals = totals if isinstance(totals, np.ndarray) else np.zeros_like(charges) + for i, symbol in enumerate(frame.symbol_set): + # Values for core, valence, and rydberg fields don't affect anything in gview. + dump_lines.append(f"{symbol} {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}", + *[" " + "-" * 71, f" * Total * {np.sum(charges):.6f} " + f"{0.0:.6f} " * 3 + f"{np.sum(_totals):.6f}"], ] return dump_lines -def _log_nbo_fields( - frame: Structure | Molecule, site_property: str = "charges", spin_site_property: str = "magmom" -) -> list[str]: +def _log_nbo(frame: Structure | Molecule, charge_key: str, spin_key: str) -> list[str]: dump_lines = [] - if site_property in frame.site_properties: - nbo_charges = np.array(frame.site_properties[site_property]) - dump_lines += _gview_log_nbo_fields(frame, nbo_charges) - spins: None | NDArray = None - if spin_site_property in frame.site_properties: - spins = np.array(frame.site_properties[spin_site_property]) + 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: + spins = np.array(frame.site_properties[spin_key]) # NBO spins are evaluated as the difference between the "total" nbo fields for alpha and beta densities. 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 += _gview_log_nbo_fields(frame, np.zeros_like(spin_arr), totals=spin_arr) + dump_lines += _log_nbo_fields(frame, np.zeros_like(spin_arr), totals=spin_arr) return dump_lines -def _log_forces(frame: Structure | Molecule) -> list[str]: +def _log_forces(frame: Structure | Molecule, do_cell: bool = True) -> list[str]: # Gaussian calculations with periodic boundary conditions will log forces on lattice vectors, but these forces # do not appear to be read by Gaussview, and writing forces for lattice vectors does not seem to fix the issue # of atomic forces not being read by Gaussview (for log files with lattice vectors). dump_lines = [] if "forces" in frame.site_properties: - dump_lines += [ - "-" * 67, - " Center Atomic" + " " * 19 + "Forces (Hartrees/Bohr)", - " Number " + (" " * 14).join(["Number", "X", "Y", "Z"]), - "-" * 67, + dump_lines = [ + *["-" * 67, " Center Atomic" + " " * 19 + "Forces (Hartrees/Bohr)"], + *[" Number " + (" " * 14).join(["Number", "X", "Y", "Z"]), "-" * 67], ] forces = frame.site_properties["forces"] - atomic_numbers = [site.specie.number for site in frame.sites] - for i, number in enumerate(atomic_numbers): - dump_lines.append(f" {i + 1} {number}\t" + "\t".join(f"{forces[i][j]:.9f}" for j in range(3))) - lat_forces = np.zeros((3, 3)) - dump_lines += [ - " " * 25 + "-2" + " " * 10 + " ".join(f"{lat_forces[i][j]:.9f}" for j in range(3)) for i in range(3) - ] + 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(np.array(forces), axis=1) dump_lines += [" " + "-" * 67, f" Cartesian Forces: Max {max(nforces):.9f} RMS {np.std(nforces):.9f}"] return dump_lines From d6c125ab7bc3c8ea6e01473903b4562ff97d1cc2 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 18 Jun 2025 20:36:26 -0600 Subject: [PATCH 15/18] working to rewrite helper functions in as few lines as possible --- src/pymatgen/io/gaussian.py | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/src/pymatgen/io/gaussian.py b/src/pymatgen/io/gaussian.py index 064eafca95d..9c8b4f0932b 100644 --- a/src/pymatgen/io/gaussian.py +++ b/src/pymatgen/io/gaussian.py @@ -1350,13 +1350,11 @@ def _traj_to_gaussian_log(traj: Trajectory, do_cell: bool, energies: list[float] 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_input_orientation(traj[i], do_cell=do_cell), f"SCF Done: E = {energies[i]:.6f} "], *_log_mulliken(traj[i], _sp_map["mulliken charges"], _sp_map["mulliken spin"]), *_log_esp(traj[i], _sp_map["esp charges"]), *_log_nbo(traj[i], _sp_map["nbo charges"], _sp_map["nbo spin"]), - *_log_forces(traj[i]), - *["", " " + "Grad" * 18, "", f" Step number {i + 1}", "", " " + "Grad" * 18], + *[*_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"] @@ -1368,8 +1366,7 @@ def _log_input_orientation(frame: Structure | Molecule, do_cell=True) -> list[st dump_lines = [ *[" " * 24 + "Standard orientation:" + " " * 26, " " + "-" * 69], " Center Atomic Atomic Coordinates (Angstroms)", - " Number Number Type X Y Z", - " " + "-" * 69, + *[" 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: @@ -1403,20 +1400,20 @@ def _log_esp(frame: Structure | Molecule, charge_key: str) -> list[str]: 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):.6f} Dipole= {' '.join(str(v) for v in dipole)} Tot= {tot}") + 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], ] - _totals = totals if isinstance(totals, np.ndarray) else np.zeros_like(charges) - for i, symbol in enumerate(frame.symbol_set): + 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"{symbol} {i + 1} {charges[i]:.6f} {0.0:.6f} {0.0:.6f} {0.0:.6f} {_totals[i]:.6f}") + 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}"], ] @@ -1428,9 +1425,8 @@ def _log_nbo(frame: Structure | Molecule, charge_key: str, spin_key: str) -> lis 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: + if spin_key in frame.site_properties: # NBO spins read as difference between alpha/beta totals spins = np.array(frame.site_properties[spin_key]) - # NBO spins are evaluated as the difference between the "total" nbo fields for alpha and beta densities. 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): @@ -1445,15 +1441,15 @@ def _log_forces(frame: Structure | Molecule, do_cell: bool = True) -> list[str]: # of atomic forces not being read by Gaussview (for log files with lattice vectors). 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], ] - forces = frame.site_properties["forces"] 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(np.array(forces), axis=1) + 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 From c5a38db89d1fa63c2f0607fec32a5467a833c082 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 18 Jun 2025 20:41:05 -0600 Subject: [PATCH 16/18] TODO --- src/pymatgen/io/gaussian.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/pymatgen/io/gaussian.py b/src/pymatgen/io/gaussian.py index 9c8b4f0932b..6cf644c69c3 100644 --- a/src/pymatgen/io/gaussian.py +++ b/src/pymatgen/io/gaussian.py @@ -1436,9 +1436,7 @@ def _log_nbo(frame: Structure | Molecule, charge_key: str, spin_key: str) -> lis def _log_forces(frame: Structure | Molecule, do_cell: bool = True) -> list[str]: - # Gaussian calculations with periodic boundary conditions will log forces on lattice vectors, but these forces - # do not appear to be read by Gaussview, and writing forces for lattice vectors does not seem to fix the issue - # of atomic forces not being read by Gaussview (for log files with lattice vectors). + # 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"]) From 6540945e68c991fd713412f9f7528cc7c28788f4 Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 18 Jun 2025 20:46:35 -0600 Subject: [PATCH 17/18] using `or` functionality to enforce defaults --- src/pymatgen/core/trajectory.py | 3 ++- src/pymatgen/io/gaussian.py | 14 +++----------- 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 472bf00ea49..25265a1fb49 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -509,7 +509,8 @@ def write_Gaussian_log( else: energies = [0] * len(self) self.to_positions() - lines = _traj_to_gaussian_log(self, write_lattice, energies, site_property_map) + _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() diff --git a/src/pymatgen/io/gaussian.py b/src/pymatgen/io/gaussian.py index 6cf644c69c3..977b07112d2 100644 --- a/src/pymatgen/io/gaussian.py +++ b/src/pymatgen/io/gaussian.py @@ -1339,21 +1339,13 @@ def to_input( def _traj_to_gaussian_log(traj: Trajectory, do_cell: bool, energies: list[float], sp_map: dict | None) -> str: - _sp_map = { - "mulliken charges": "charges", - "mulliken spin": "magmom", - "nbo charges": "charges", - "nbo spin": "magmom", - "esp charges": "charges", - } - _sp_map.update(sp_map or {}) 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], _sp_map["mulliken charges"], _sp_map["mulliken spin"]), - *_log_esp(traj[i], _sp_map["esp charges"]), - *_log_nbo(traj[i], _sp_map["nbo charges"], _sp_map["nbo spin"]), + *_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."] From 6f5a2a86b26fa6a7e8d10cfdda0b3894b5556dac Mon Sep 17 00:00:00 2001 From: benrich37 Date: Wed, 3 Sep 2025 15:21:37 -0600 Subject: [PATCH 18/18] Moving import into function to reduce general Trajectory import time --- src/pymatgen/core/trajectory.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/pymatgen/core/trajectory.py b/src/pymatgen/core/trajectory.py index 4c9ce41185d..4fc0695cef4 100644 --- a/src/pymatgen/core/trajectory.py +++ b/src/pymatgen/core/trajectory.py @@ -17,7 +17,6 @@ from pymatgen.core.structure import Composition, DummySpecies, Element, Lattice, Molecule, Species, Structure from pymatgen.io.ase import NO_ASE_ERR, AseAtomsAdaptor -from pymatgen.io.gaussian import _traj_to_gaussian_log if NO_ASE_ERR is None: from ase.io.trajectory import Trajectory as AseTrajectory @@ -501,6 +500,10 @@ def write_Gaussian_log( 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]