Skip to content

Jake update#51

Open
j77king wants to merge 3 commits into
mainfrom
jake_update
Open

Jake update#51
j77king wants to merge 3 commits into
mainfrom
jake_update

Conversation

@j77king

@j77king j77king commented Jul 7, 2026

Copy link
Copy Markdown
Collaborator

adding mlip support and spin density

Summary by CodeRabbit

  • New Features

    • Added MLIP support to process neutral, reduced, and oxidized structure files and generate molecular, atomic, and bond-level outputs.
    • Added command-line options for supplying MLIP input directories.
    • Included sample structure files for the new workflow.
  • Bug Fixes

    • Improved handling of charge and spin data so missing values are managed more consistently.
    • Added support for parsing Mulliken spin information when direct spin data is unavailable.
  • Tests

    • Expanded automated coverage for charge parsing and MLIP output generation.

@coderabbitai

coderabbitai Bot commented Jul 7, 2026

Copy link
Copy Markdown

Review Change Stack

📝 Walkthrough

Walkthrough

Adds a new moldscript/MLIP.py module implementing an mlip class that parses .extxyz files across neutral/reduced/oxidized states into a structured dictionary with derived molecular/atomic properties, wires it into CLI arguments and the main pipeline, adds example fixture files and tests, and separately reformats/extends Mulliken spin parsing in charges.py.

Changes

MLIP Module and Pipeline Integration

Layer / File(s) Summary
mlip class construction and file discovery
moldscript/MLIP.py
Constructor accepts neutral/reduced/oxidized directories, counts and reports progress reading *.extxyz files.
extxyz parsing and header/property extraction
moldscript/MLIP.py
Parses extxyz headers, Properties schema, atom symbols/coordinates, and configurable per-atom charge property.
Derived molecular/atomic feature computation
moldscript/MLIP.py
Computes SCF energy, vertical IP/EA, dipole vectors, Fukui indices, and bond-length matrices with NaN fallback handling.
CLI flags and pipeline wiring
moldscript/argument_parser.py, moldscript/moldscript.py
Adds mlip_neutral/mlip_reduced/mlip_oxidized options and invokes the mlip stage in main(), replacing data_dicts with returned file_data.
Example extxyz fixtures and MLIP tests
moldscript/examples/arbr/mlip/.../*.extxyz, tests/test_moldscript.py
Adds sample neutral/reduced/oxidized extxyz files and tests validating parsed shapes, numeric values, and CSV outputs via get_df.

Charges Module Spin Parsing and Formatting

Layer / File(s) Summary
Mulliken spin parsing and formatting cleanup
moldscript/charges.py, tests/test_moldscript.py
Reformats imports and data_dict writes, adds get_atom_spins/parse_gaussian_mulliken_spins for Mulliken spin density fallback parsing, with a new test.

Estimated code review effort: 3 (Moderate) | ~25 minutes

Sequence Diagram(s)

sequenceDiagram
  participant CLI
  participant moldscript_main as moldscript.main
  participant mlip
  participant Filesystem

  CLI->>moldscript_main: mlip_neutral/reduced/oxidized args
  moldscript_main->>mlip: mlip(neutral, reduced, oxidized, data_dicts)
  mlip->>Filesystem: scan and read *.extxyz files
  Filesystem-->>mlip: header, symbols, coords, atom_props
  mlip->>mlip: compute SCF energy, IP/EA, dipole, Fukui indices, bond matrix
  mlip-->>moldscript_main: file_data dictionary
  moldscript_main->>moldscript_main: update data_dicts
Loading

Related Issues: None specified.

Related PRs: None specified.

Suggested labels: enhancement, feature, tests

Suggested reviewers: patonlab maintainers familiar with charge/spin parsing and MLIP data pipelines

🐰 A rabbit hopped through extxyz files galore,
Counting dipoles, spins, and charges by the score.
Fukui indices bloom from every state,
Neutral, oxidized, reduced — all cooperate.
New tests confirm the numbers ring true —
A tidy burrow of data, built anew!

🚥 Pre-merge checks | ✅ 4 | ❌ 1

❌ Failed checks (1 inconclusive)

Check name Status Explanation Resolution
Title check ❓ Inconclusive The title is too vague and does not describe the MLIP and spin-density changes in the pull request. Use a descriptive title such as "Add MLIP support and spin density parsing".
✅ Passed checks (4 passed)
Check name Status Explanation
Description Check ✅ Passed Check skipped - CodeRabbit’s high-level summary is enabled.
Docstring Coverage ✅ Passed No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.
Linked Issues check ✅ Passed Check skipped because no linked issues were found for this pull request.
Out of Scope Changes check ✅ Passed Check skipped because no linked issues were found for this pull request.
✨ Finishing Touches
📝 Generate docstrings
  • Create stacked PR
  • Commit on current branch
🧪 Generate unit tests (beta)
  • Create PR with unit tests
  • Commit unit tests in branch jake_update

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands.

@chatgpt-codex-connector chatgpt-codex-connector Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: aec37f9578

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

Comment thread moldscript/MLIP.py
entry["atom"]["frad"] = frad
entry["atom"]["charges_neutral"] = q_neutral

self.data_dict[key] = entry

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P1 Badge Preserve existing entries when adding MLIP descriptors

When MLIP is combined with any earlier module in the same CLI run, main() passes the accumulated data_dicts into mlip, but this assignment replaces the whole molecule entry for matching keys. For extxyz files whose source_file stem matches the Gaussian-derived key, descriptors already parsed from opt/SPC/NMR/NBO/charges are silently dropped from the final CSVs and only the MLIP fields remain; merge into the existing nested mol/atom/bond dicts instead of overwriting the entry.

Useful? React with 👍 / 👎.

Comment thread moldscript/MLIP.py
Comment on lines +261 to +263
source = header.get("source_file", "")
if source:
return Path(source).stem

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2 Badge Normalize source_file state suffixes before grouping

When extxyz headers contain state-specific source_file values such as foo_neutral_vertical.com, foo_anion_Nplus1_vertical.com, and foo_cation_Nminus1_vertical.com, this early return bypasses the suffix stripping below. The three supplied state files then land under separate molecule keys, so each states dict is missing the other charge states and the vertical IE/EA/Fukui outputs become NaN; apply the same suffix normalization to the source stem before returning it.

Useful? React with 👍 / 👎.

@coderabbitai coderabbitai Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actionable comments posted: 4

🧹 Nitpick comments (1)
tests/test_moldscript.py (1)

184-242: 📐 Maintainability & Code Quality | 🔵 Trivial | ⚡ Quick win

Consider adding a merge-conflict regression test.

These tests validate the standalone mlip parser output well, but none exercise mlip(..., data_dict=<pre-populated dict>) where the same key already has mol/atom/bond data from another stage (e.g. opt/charges). Given the overwrite behavior flagged in MLIP.py (Line 192), a regression test asserting pre-existing keys are preserved after merging would help catch this class of bug going forward.

🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

In `@tests/test_moldscript.py` around lines 184 - 242, Add a regression test for
the merge path in mlip when data_dict already contains an entry for the same
filename. In tests/test_moldscript.py, exercise mlip(...,
data_dict=<pre-populated dict>) with an existing key that already has
mol/atom/bond data from another stage, then assert the parser merge preserves
the pre-existing values instead of overwriting them. Use build_arbr141_parser,
mlip, and the existing arbr141_wb97xd fixture shape as the reference points.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

Inline comments:
In `@moldscript/charges.py`:
- Around line 62-73: Guard chg_data before accessing atomcharges in the charge
parsing flow, because parse_cc_data() can return None and the current
unconditional use in the charges handling block will raise AttributeError and
stop the batch. Update the logic around the atom charge assignment and the
subsequent get_atom_spins call in charges.py to first verify chg_data is present
(following the existing pattern used for cpu_time), and only then read
atomcharges and populate self.data_dict[filename]['atom'].

In `@moldscript/MLIP.py`:
- Around line 291-303: The _state_charges helper in MLIP.py silently pads or
truncates charge arrays when len(arr) does not match natoms, which can hide
state/geometry mismatches. Update _state_charges to emit a warning whenever the
atom charge array length differs from natoms, while keeping the existing
NaN-padding behavior; use the _state_charges method and
self.atom_charge_property to locate the mismatch path and make the warning
include the expected and actual lengths.
- Around line 196-206: The _read_extxyz method currently slices atom_lines
without checking that the file actually contains the declared number of atoms,
so truncated or malformed input can slip through silently. In _read_extxyz,
after building atom_lines from lines[2:2 + natoms], validate that its length
matches natoms and raise a clear ValueError when it does not; keep the check
close to the existing natoms/header parsing logic so the failure is reported
before _parse_properties_schema or downstream atomnos/coords construction.
- Around line 117-194: The per-key record in the MLIP aggregation loop is being
overwritten by assigning a new entry to self.data_dict[key], which drops any
existing mol, atom, bond, or CPU_time data already stored for that key. Update
the logic in the state-processing block inside the MLIP class so it merges into
the existing self.data_dict[key] entry in place, preserving previously populated
nested fields while adding the new MLIP values.

---

Nitpick comments:
In `@tests/test_moldscript.py`:
- Around line 184-242: Add a regression test for the merge path in mlip when
data_dict already contains an entry for the same filename. In
tests/test_moldscript.py, exercise mlip(..., data_dict=<pre-populated dict>)
with an existing key that already has mol/atom/bond data from another stage,
then assert the parser merge preserves the pre-existing values instead of
overwriting them. Use build_arbr141_parser, mlip, and the existing
arbr141_wb97xd fixture shape as the reference points.
🪄 Autofix (Beta)

Fix all unresolved CodeRabbit comments on this PR:

  • Push a commit to this branch (recommended)
  • Create a new PR with the fixes

ℹ️ Review info
⚙️ Run configuration

Configuration used: defaults

Review profile: CHILL

Plan: Pro

Run ID: fdea8659-82ff-4828-952c-a048c85a5471

📥 Commits

Reviewing files that changed from the base of the PR and between b6afac9 and aec37f9.

⛔ Files ignored due to path filters (3)
  • moldscript/examples/spin_examples/A1a_cat_rad_opt.log is excluded by !**/*.log
  • moldscript/examples/spin_examples/A1b_cat_rad_opt.log is excluded by !**/*.log
  • moldscript/examples/spin_examples/A1c_cat_rad_opt.log is excluded by !**/*.log
📒 Files selected for processing (8)
  • moldscript/MLIP.py
  • moldscript/argument_parser.py
  • moldscript/charges.py
  • moldscript/examples/arbr/mlip/neutral/arbr141_wb97xd_neutral_vertical.extxyz
  • moldscript/examples/arbr/mlip/oxidized/arbr141_wb97xd_cation_Nminus1_vertical.extxyz
  • moldscript/examples/arbr/mlip/reduced/arbr141_wb97xd_anion_Nplus1_vertical.extxyz
  • moldscript/moldscript.py
  • tests/test_moldscript.py

Comment thread moldscript/charges.py
Comment on lines +62 to +73
if len(chg_data.atomcharges.keys()) == 1 and 'mulliken' in chg_data.atomcharges:
self.data_dict[filename]['atom']['mulliken_charge'] = chg_data.atomcharges['mulliken']
else:
for i in chg_data.atomcharges.keys():
if 'mulliken' not in i and 'sum' not in i:
self.data_dict[filename]['atom'][str(i)+'_charge'] = chg_data.atomcharges[i]

for spin_type, spins in self.get_atom_spins(chg_data, self.data[file_name]).items():
self.data_dict[filename]['atom'][str(spin_type)+'_spin'] = spins

cpu_times = chg_data.metadata.get("cpu_time") if chg_data and hasattr(chg_data, "metadata") else None
self.module_cpu_seconds += record_cpu_time(self.data_dict, file_name, self.data[file_name], cpu_times)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🩺 Stability & Availability | 🟠 Major | ⚡ Quick win

🧩 Analysis chain

🏁 Script executed:

#!/bin/bash
set -euo pipefail

# Map the target file and inspect the relevant region with line numbers.
git ls-files moldscript/charges.py tests/test_moldscript.py
echo '--- moldscript/charges.py ---'
cat -n moldscript/charges.py | sed -n '1,220p'
echo '--- tests/test_moldscript.py (relevant excerpts) ---'
cat -n tests/test_moldscript.py | sed -n '1,240p'

Repository: patonlab/molDscript

Length of output: 20734


Guard chg_data before reading atomcharges. parse_cc_data() returns None on parse errors, so this unconditional access raises AttributeError and aborts the batch when one file fails; the new spin/cpu-time guards already show the intended pattern.

🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

In `@moldscript/charges.py` around lines 62 - 73, Guard chg_data before accessing
atomcharges in the charge parsing flow, because parse_cc_data() can return None
and the current unconditional use in the charges handling block will raise
AttributeError and stop the batch. Update the logic around the atom charge
assignment and the subsequent get_atom_spins call in charges.py to first verify
chg_data is present (following the existing pattern used for cpu_time), and only
then read atomcharges and populate self.data_dict[filename]['atom'].

Comment thread moldscript/MLIP.py
Comment on lines +117 to +194
for key, states in self.state_records.items():
base_state = self._pick_base_state(states)
base_record = states[base_state]

symbols = base_record["symbols"]
coords = np.array(base_record["coords"], dtype=float)
atomnos = np.array([self.ptable.GetAtomicNumber(sym) for sym in symbols], dtype=int)

entry = {
"mol": {},
"atom": {},
"bond": {},
"CPU_time": datetime.timedelta(0),
}

entry["atom"]["atomnos"] = atomnos
entry["bond"]["bond_length"] = self._bond_length_matrix(coords)

entry["mol"].setdefault("smiles", "")

# State energies
e_neutral_ev = self._state_energy_ev(states.get("neutral"))
e_reduced_ev = self._state_energy_ev(states.get("reduced"))
e_oxidized_ev = self._state_energy_ev(states.get("oxidized"))

if np.isfinite(e_neutral_ev):
entry["mol"]["scfenergy"] = e_neutral_ev * eV_to_hartree

else:
entry["mol"]["scfenergy"] = np.nan

if np.isfinite(e_oxidized_ev) and np.isfinite(e_neutral_ev):
ie_h = (e_oxidized_ev - e_neutral_ev) * eV_to_hartree
entry["mol"]["vertical_ie"] = ie_h
entry["mol"]["ionization_potential"] = ie_h
else:
entry["mol"]["vertical_ie"] = np.nan
entry["mol"]["ionization_potential"] = np.nan

if np.isfinite(e_reduced_ev) and np.isfinite(e_neutral_ev):
ea_h = (e_reduced_ev - e_neutral_ev) * eV_to_hartree
entry["mol"]["vertical_ea"] = ea_h
entry["mol"]["electron_affinity"] = ea_h
else:
entry["mol"]["vertical_ea"] = np.nan
entry["mol"]["electron_affinity"] = np.nan

# Dipole from the first state that provides it (neutral preferred).
dip = self._dipole_vector(states)
entry["mol"]["dipole_x"] = dip[0]
entry["mol"]["dipole_y"] = dip[1]
entry["mol"]["dipole_z"] = dip[2]
entry["mol"]["dipole"] = float(np.linalg.norm(dip)) if np.all(np.isfinite(dip)) else np.nan

q_neutral = self._state_charges(states.get("neutral"), len(atomnos))
q_reduced = self._state_charges(states.get("reduced"), len(atomnos))
q_oxidized = self._state_charges(states.get("oxidized"), len(atomnos))

# Follow sign convention already used in fukui.py
if np.all(np.isfinite(q_neutral)) and np.all(np.isfinite(q_reduced)):
fplus = -1.0 * (q_reduced - q_neutral)
else:
fplus = np.full(len(atomnos), np.nan)

if np.all(np.isfinite(q_neutral)) and np.all(np.isfinite(q_oxidized)):
fminus = -1.0 * (q_neutral - q_oxidized)
else:
fminus = np.full(len(atomnos), np.nan)

frad = 0.5 * (fplus + fminus)
entry["atom"]["fplus"] = fplus
entry["atom"]["fminus"] = fminus
entry["atom"]["frad"] = frad
entry["atom"]["charges_neutral"] = q_neutral

self.data_dict[key] = entry

return self.data_dict

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🗄️ Data Integrity & Integration | 🟠 Major | ⚡ Quick win

🧩 Analysis chain

🏁 Script executed:

#!/bin/bash
set -euo pipefail

# Map the relevant file first, then inspect the surrounding code paths.
printf '\n== moldscript/MLIP.py outline ==\n'
ast-grep outline moldscript/MLIP.py --view expanded || true

printf '\n== moldscript/MLIP.py around the assignment ==\n'
sed -n '1,260p' moldscript/MLIP.py | cat -n

printf '\n== search for data_dict construction / reuse ==\n'
rg -n "data_dict|mlip\(" moldscript.py moldscript -g '!**/__pycache__/**' || true

printf '\n== inspect charges.py pattern mentioned in the comment ==\n'
if [ -f moldscript/charges.py ]; then
  sed -n '1,240p' moldscript/charges.py | cat -n
fi

Repository: patonlab/molDscript

Length of output: 33898


🏁 Script executed:

#!/bin/bash
set -euo pipefail

printf '\n== moldscript.py around the MLIP call ==\n'
sed -n '120,150p' moldscript/moldscript.py | cat -n

printf '\n== MLIP key normalization and base-state selection ==\n'
sed -n '260,340p' moldscript/MLIP.py | cat -n

printf '\n== utils.initiate_data_dict and record_cpu_time ==\n'
sed -n '130,210p' moldscript/utils.py | cat -n
sed -n '80,115p' moldscript/utils.py | cat -n

Repository: patonlab/molDscript

Length of output: 9986


Merge MLIP fields into the existing entry
self.data_dict[key] = entry replaces the whole per-file record, so any earlier mol/atom/bond/CPU_time fields for the same key are dropped. Update the nested dicts in place instead of overwriting the entry.

🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

In `@moldscript/MLIP.py` around lines 117 - 194, The per-key record in the MLIP
aggregation loop is being overwritten by assigning a new entry to
self.data_dict[key], which drops any existing mol, atom, bond, or CPU_time data
already stored for that key. Update the logic in the state-processing block
inside the MLIP class so it merges into the existing self.data_dict[key] entry
in place, preserving previously populated nested fields while adding the new
MLIP values.

Comment thread moldscript/MLIP.py
Comment on lines +196 to +206
def _read_extxyz(self, file_path):
lines = file_path.read_text(encoding="utf-8").splitlines()
if len(lines) < 2:
raise ValueError(f"Invalid extxyz file: {file_path}")

natoms = int(lines[0].strip())
header_line = lines[1].strip()
header = self._parse_header_line(header_line)

schema = self._parse_properties_schema(header.get("Properties", "species:S:1:pos:R:3"))
atom_lines = lines[2 : 2 + natoms]

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🎯 Functional Correctness | 🟡 Minor | ⚡ Quick win

No validation that the file actually contains natoms atom lines.

atom_lines = lines[2 : 2 + natoms] (Line 206) will silently yield fewer rows than natoms if the file is truncated/malformed, producing an atomnos/coords array shorter than the declared count with no error raised.

Consider asserting len(atom_lines) == natoms and raising a clear ValueError otherwise.

🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

In `@moldscript/MLIP.py` around lines 196 - 206, The _read_extxyz method currently
slices atom_lines without checking that the file actually contains the declared
number of atoms, so truncated or malformed input can slip through silently. In
_read_extxyz, after building atom_lines from lines[2:2 + natoms], validate that
its length matches natoms and raise a clear ValueError when it does not; keep
the check close to the existing natoms/header parsing logic so the failure is
reported before _parse_properties_schema or downstream atomnos/coords
construction.

Comment thread moldscript/MLIP.py
Comment on lines +291 to +303
def _state_charges(self, record, natoms):
if record is None:
return np.full(natoms, np.nan)
values = record["atom_props"].get(self.atom_charge_property)
if values is None:
return np.full(natoms, np.nan)
arr = np.array(values, dtype=float)
if len(arr) == natoms:
return arr
padded = np.full(natoms, np.nan)
ncopy = min(natoms, len(arr))
padded[:ncopy] = arr[:ncopy]
return padded

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🗄️ Data Integrity & Integration | 🟡 Minor | ⚡ Quick win

Silent NaN-padding on atom-count mismatch.

If a state's charge array length differs from natoms, values are silently truncated/padded with NaN with no warning (Lines 300-303). This could mask a real geometry/state mismatch and quietly produce partially-NaN Fukui indices downstream.

Consider logging a warning when len(arr) != natoms so mismatches are visible rather than silently absorbed.

🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

In `@moldscript/MLIP.py` around lines 291 - 303, The _state_charges helper in
MLIP.py silently pads or truncates charge arrays when len(arr) does not match
natoms, which can hide state/geometry mismatches. Update _state_charges to emit
a warning whenever the atom charge array length differs from natoms, while
keeping the existing NaN-padding behavior; use the _state_charges method and
self.atom_charge_property to locate the mismatch path and make the warning
include the expected and actual lengths.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant