Repository navigation
Parse Q-Chem custom SMD solvents and fragmented $molecule sections - #1537
rohithsrinivaas wants to merge 1 commit into
Conversation
Two fixes for Q-Chem task documents:
- Custom SMD solvents ("solvent other"): Calculation.from_qchem_files calls
lot_solvent_string without custom_smd, so TaskDoc.from_directory raised
"SMD calculation with solvent=other requires custom_smd!" for every calculation
in a custom SMD solvent, even though the solvent descriptors are in the input's
$smx section (atomate2 writes them there). solvent() now builds custom_smd from
$smx (epsilon, soln, sola, solb, solg, solc, solh) when it is not given, through
a new custom_smd_from_smx helper used by both solvent() implementations. An
explicit custom_smd still takes precedence.
- Fragments: a $molecule section split into fragments (e.g. fragment-based
Hirshfeld constrained DFT, EDA) is read by pymatgen as a list of Molecules, which
CalculationInput.from_qcinput rejected ("QCInput contains multiple molecules").
The fragments are now combined into one Molecule with the total charge and the
high-spin multiplicity that pymatgen's QCInput writes for them; "read" is still
rejected.
Test files: a wB97M-V/def2-SVPD single point in a custom SMD solvent and a
fragment-based Hirshfeld CDFT single point (Li+ / ethylene carbonate radical
anion, Q-Chem 6.2).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1537 +/- ##
==========================================
+ Coverage 85.37% 85.43% +0.06%
==========================================
Files 216 216
Lines 17793 17809 +16
==========================================
+ Hits 15190 15215 +25
+ Misses 2603 2594 -9 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
tsmathis
left a comment
tsmathis
left a comment
There was a problem hiding this comment.
A couple questions about the added test files. What are those? Just plain input/output files from these types of workflows?
I see you have a linked atomate2 PR, will this need an emmet release?
| #: $smx keys of a custom ("solvent other") SMD solvent, in the order of custom_smd: | ||
| #: dielectric constant, refractive index, Abraham acidity and basicity, surface | ||
| #: tension, aromaticity and electronegative halogenicity | ||
| CUSTOM_SMD_KEYS = ("epsilon", "soln", "sola", "solb", "solg", "solc", "solh") |
There was a problem hiding this comment.
Do these originate from some external docs?
Or are these just standard qchem terms/fields? Like how VASP INCARs have known fields, etc.
| else None | ||
| ), | ||
| initial_molecule=molecule, | ||
| charge=(int(molecule.as_dict()["charge"]) if molecule.as_dict() else None), |
There was a problem hiding this comment.
any reason not to do just int(molecule.charge)?
I know the as_dict() ternary existed before this PR, but do you know if it is actually necessary to guard on molecule.as_dict()?
There is already the early error above on if not isinstance(molecule, Molecule): ..., so the molecule should at least be valid here (?)
| if not isinstance(molecule, Molecule): | ||
| raise ValueError("QCInput contains multiple molecules.") |
There was a problem hiding this comment.
This error message doesn't really communicate to the caller what went wrong now that there are multiple validation checks.
There are at least two error paths: 1. molecule is a list and not every member of the list is a Molecule, 2. molecule is a str-literal, (and I guess 3. is that molecule from qcinput.molecule is not a Molecule?)
It would be good to disambiguate the reason from_qcinput failed.
| if isinstance(molecule, list) and all( | ||
| isinstance(m, Molecule) for m in molecule | ||
| ): | ||
| molecule = Molecule( | ||
| species=[site.species for fragment in molecule for site in fragment], | ||
| coords=[site.coords for fragment in molecule for site in fragment], | ||
| charge=sum(fragment.charge for fragment in molecule), | ||
| spin_multiplicity=1 | ||
| + sum(fragment.spin_multiplicity - 1 for fragment in molecule), | ||
| ) |
There was a problem hiding this comment.
This iterates over the sites twice and the fragments four times (and +1 for the initial iteration in the all(...)). A single-pass version should be possible, something like?
species, coords, charge, spin = [], [], 0, 1
for fragment in molecule:
for site in fragment:
species.append(site.species)
coords.append(site.coords)
charge += fragment.charge
spin += fragment.spin_multiplicity - 1
molecule = Molecule(
species=species,
coords=coords,
charge=charge,
spin_multiplicity=spin,
)could also technically skip the all(...), check if each fragment is a Molecule in the loop, and break early if any fragment is not a Molecule
Two fixes for Q-Chem task documents (
emmet.core.qchem):1. Custom SMD solvents (
solvent other).Calculation.from_qchem_filescallslot_solvent_string(input_doc, validate_lot=...)withoutcustom_smd, soTaskDoc.from_directoryraises "SMD calculation with solvent=other requirescustom_smd!" for every calculation in a custom SMD solvent, even though the seven solvent descriptors are in the input's
$smxsection (which is where atomate2'sQCInputGeneratorwrites them). In practice every atomate2 Q-Chem job withsmd_solvent="other"fizzles after Q-Chem has finished successfully.solvent()now buildscustom_smdfrom$smx(epsilon,soln,sola,solb,solg,solc,solh) when it is not given, through a newcustom_smd_from_smxhelper used by bothsolvent()implementations (calc_types/utils.pyandcalculation.py). An explicitcustom_smdstill takes precedence, andsolvent otherwithout descriptors still raises.2. Fragmented
$moleculesections. A$moleculesection split into fragments (e.g. fragment-based Hirshfeld constrained DFT,cdft_pop = fbh, or EDA) is read by pymatgen as a list ofMolecules, whichCalculationInput.from_qcinputrejected ("QCInput contains multiple molecules"). The fragments are now combined into one
Moleculewith the total charge and the high-spin multiplicity that pymatgen'sQCInput.molecule_templatewrites for them.molecule = "read"isstill rejected.
Test files (Q-Chem 6.2): a wB97M-V/def2-SVPD single point in a custom SMD solvent (
test_files/qchem/qchem_custom_smd_sp_test) and a fragment-based Hirshfeld CDFT single point of Li⁺ / ethylene carbonate radical anion with fragments and acustom SMD solvent (
test_files/qchem/qchem_cdft_fragments_sp_test).Contributor Checklist
custom_smdfrom the$smxsection in bothsolvent()implementations (newcustom_smd_from_smxhelper)$moleculesections into oneMoleculeinCalculationInput.from_qcinputtests/qchem,tests/test_qc_task.py,tests/molecules: 35 passed, 2 skipped — the skips are the existing "Pymatgen OBAlign needs fix" ones; ruff 0.12.0, black 25.1.0 and mypy clean)test_solvent; newtest_custom_smd_from_smx,test_custom_smd_solvation_lot,test_fragments_inputandtest_input_with_molecule_read. The new tests fail onmainwithout this change.)