Allow changing molecular weights to compute spectroscopic data for isotopologues
zpalmer618 opened this issue · 6 comments
Since the normal coordinates contain information specific to the mass of the atom, there is, at this moment, no way to update that mass for studying isotope shifts in the rovibrational spectroscopic data. Please consider adding this feature in your next update of spectro.
I've updated the title a bit because this isn't related to pbqff, per se. Especially because I don't expect to add any isotopic information to pbqff itself. I'm picturing this being something you re-run spectro itself on after a pbqff run. Does that sound reasonable?
Also can you attach an example input and output file (spectro.in and spectro2.out) with isotopic substitutions for testing?
Thanks, and yes that does sound reasonable. Attached are the input and output files as requested. The system is H2SC that has been doubly deuterated.
spectro.in.txt
spectro2.out.txt
The first part of this request (read weights from the input file) should be pretty trivial to implement. Instead of accessing the default weights via self.geom.weights() here, it should combine those defaults with any overrides in self.weights. I would tentatively call that method self.weights and replace all calls to self.geom.weights() with just self.weights().
The second part (converting force constants from normal coordinates to Cartesian coordinates) is much more involved. I think the procedure for doing so is outlined in this paper. I'm not entirely sure what interface I would want for doing these transformations in spectro itself (a command line flag? read it from the input file? and is this part of a normal run, or is the flag/input option just for doing the transformation?) so I think a good first step would be writing some kind of script to figure out the transformations. This script would read any information it needs from finish.spectro or spectro.json (presumably the LXM matrix, harmonic frequencies, and geometry/masses) and output files (spectro.in, fort.15,30,40) that the updated version of spectro from step 1 can read.
Just like they do in the paper, we'd also want some round-trip tests on this. A normal spectro run from the Fortran inputs reads Cartesian force constants and converts them to arrays called f3qcm and f4qcm. This request involves turning those back into Cartesian force constants, so the transformation function should take f3qcm and output whatever f3qcm is made of in the first place. I guess this kind of test is actually an argument against writing this as an external script because you need access to the Rust API that actually makes f*qcm. We'll see how it goes.
Here are the fort.15,30,40 files for your test purposes. Thanks, Dr. Westbrook.
fort.15.txt
fort.30.txt
fort.40.txt
As usual, the first part is not as easy as expected. I made the changes I mentioned above, but the provided D2SC test case is failing. I suspect now that the input weights need to be used in process_geom as well. This in turn means that I need to handle isotopic masses in symm::normalize, and probably just in symm::Molecule more generally. com and moi both refer directly to WEIGHTS, but I should probably move weights to a field on Molecule and allow it to be modified from the defaults by spectro.