Anharmonic intensities

[isolated-matrix]

Hi again,

I just had a question about anharmonic intensities.

Assuming calculation of the anharmonic intensities was ported into the rust version of Spectro, would they require any additional calculations to be run by molpro, four, etc.? Looking at the molpro output, it does seem to provide some kind of dipole moments... In the .xml files, it gives:

<property name="Dipole moment" method="RHF" principal="true" stateSymmetry="1" stateNumber="1" stateMultiplicity="1"
     value="0.0 0.459105794240303E-02 0.0"/>

While the normal CCSD(T)-F12 doesn't seem to provide a dipole moment, using density fitted CCSD(T)-F12 can give you unrelaxed CCSD(T)-F12 dipole moments (albeit with a very demanding calculation):

   <property name="Dipole moment" method="CCSD-F12(T)" principal="true" stateSymmetry="1" stateNumber="1" stateMultiplicity="1"
     value="0.194305562308815E-13 -0.368269338899864E-07 0.341455176247557E-02"/>

I'm pretty certain CFOUR will produce CCSD(T) dipole moments, too.

Would this be enough for pbqff to generate the 65, 66 and 67 files for the original Spectro (assuming the intensity code isn't lying dormant in the rust version of Spectro)?

[ntBre]

I'm honestly not too familiar with intensity calculations and not working in chemistry anymore, so I consider it pretty unlikely that this feature will get added. The intensity code is definitely not lying dormant in the Rust version 😄

From the little I remember, I think a separate dipole moment surface would be required to compute intensities. In principle I think that could reuse the same points from the potential energy surface, but, again in my limited recollection, I think you typically want to use a different displacement size.

[isolated-matrix]

In that case, was the rust version of Spectro basically a direct port of the fortran version? If it is, maybe I could try and convert the code in the original version with the f2rust program and, if that works, maybe try to work out how to get pbqff to read the dipole moment data from CCSD(T) calculations to generate the fort.65, 66 and 67 files...

One unrelated question, by the way: I'm a little confused about the F12-XZ-cCR method. In the articles about it, it's not clear exactly how the HFF is supposed to be calculated (i.e., whether it's also calculated using the compound method, or just the standard CCSD(T)-F12/cc-pCVXZ-F12 level of theory). When I use qffbuddy, it seems as though the HFF should be calculated using the compound method (with all the extra CCSD(T) corrections). Is that true?

[isolated-matrix]

I took a look at your spectro code and it definitely wasn't a straightforward port...

So I have another idea. How difficult would it be to basically replicate what you did for the force constants (in terms of extracting the values from a molpro output, calculating derivatives and formatting them into the fort.* files) and applying them to the dipole moments? As a non-expert, I am assuming the dipole derivatives are calculated in the same way as the forcefields, and, again as a non-expert, I am assuming that the force constants are determined by fitting the calculated energies and then differentiating them...

Are those assumptions correct? If they are, I'd be happy to give doing it myself a go. If so, could you point out the relevant files where you generate the force constants? If I work from your code, hopefully I can wrangle something together that extracts dipole moments, differentiates them and prints them into fort files. Then I can try putting them through the original Spectro to see if I get anharmonic intensities... The chances of me being successful are very low, but I'm willing to give it a shot :)

[ntBre]

Thanks for the suggestions! I don't want you to think I'm ignoring them :)

These were all interesting ideas, but as you noted, the port of the spectro code (and the spectro code itself really) are not straightforward. I spent a lot of time in gdb stepping through the old Fortran code and exporting intermediate values to debug my implementation, or at least to match the bugs in the Fortran version...

As I alluded to on #515, it's been quite a while since I worked on this project and my familiarity with it has faded, so some of these big features would require a lot of investment, first just to catch up, and then to actually implement the feature. And intensity/dipole calculations are actually beyond even my peak chemistry expertise from my PhD. So basically I'm happy to try to support your efforts when I have time, and I wish you luck, this would be awesome, but I'm not sure how helpful I will be, and I think the scale of this feature is such that I wouldn't really feel confident reviewing/merging it here.

That said, I threw your questions at Codex, and it produced an answer in case it's helpful:

Details

On the specific force-constant question: your assumptions are broadly correct. Force
constants are derivatives of the potential energy with respect to the selected
coordinates. In pbqff’s finite-difference workflows they are calculated directly as
weighted combinations of energies at displaced geometries—there is no fitting step.
In the fitted SIC and normal-coordinate workflows, the energies are fitted to a
Taylor expansion, and the fitted coefficients are converted into force constants.

The most relevant files are:

• findiff.rs

  (https://github.com/ntBre/pbqff/blob/master/crates/pbqff/src/coord_type/findiff.rs#L81-L260)
  generates the displaced geometries and accumulates their energies into force
  constants. The actual second-, third-, and fourth-derivative formulas start here

  (https://github.com/ntBre/pbqff/blob/master/crates/pbqff/src/coord_type/findiff.rs#L278).

• molpro/mod.rs

  (https://github.com/ntBre/pbqff/blob/master/crates/psqs/src/program/molpro/mod.rs#L145-L215)
  extracts the specially printed PBQFF = ... energy from each Molpro output.

• The fitted normal-coordinate route is in normal.rs

  (https://github.com/ntBre/pbqff/blob/master/crates/pbqff/src/coord_type/normal.rs#L648-L708),
  while the SIC route is in sic.rs

  (https://github.com/ntBre/pbqff/blob/master/crates/pbqff/src/coord_type/sic.rs#L245-L330).

• intder/src/lib.rs
  (https://github.com/ntBre/pbqff/blob/master/crates/intder/src/lib.rs#L1326-L1358)
  performs the final coordinate conversion and writes fort.15, fort.30, and fort.40.

Conceptually, the same approach could be applied independently to the three Cartesian
components of the dipole moment. However, the existing finite-difference machinery
operates on scalar energies and calculates derivatives of orders two through four,
whereas a cubic dipole field requires derivatives of orders one through three. It
would therefore need a first-derivative formula, storage for three values per
geometry, and corresponding parsing and output code.

The SPECTRO manual says that fort.65 contains the equilibrium dipole followed by its
first derivatives, fort.66 contains the second derivatives, and fort.67 contains the
third derivatives. These must be in atomic units. The molecular orientation must also
remain consistent across all displaced calculations, since unlike an energy, the
Cartesian components of a dipole change when the molecule is rotated. SPECTRO
additionally requires the dipole and force fields to use compatible coordinate
systems and specific indexing conventions.

For an initial experiment, a standalone postprocessor that reads the existing job.*
outputs, calculates the dipole derivatives, and writes fort.65–fort.67 would probably
be simpler than modifying the job and result machinery inside pbqff.

Its response seems broadly reasonable, although I interpreted your question about force constant calculations more literally and was following the usage of f3qcm and f4qcm (the cubic and quartic force constants in normal coordinates and wavenumbers, if I understand the Fortran variable names) backwards in this order:

pbqff/crates/spectro/src/run.rs

Lines 218 to 224 in 221be4e

	// start of cubic analysis
	let f3x = self.rot3rd(f3x);
	let f3qcm = force3(self.n3n, f3x, &lx, self.nvib, &freq);
	// start of quartic analysis
	let f4x = self.rot4th(f4x);
	let f4qcm = force4(self.n3n, &f4x, &lx, self.nvib, &freq);

pbqff/crates/pbqff/src/coord_type/cart.rs

Lines 96 to 99 in 221be4e

	let fc3 = spectro::new_fc3(spectro.n3n, f3);
	let fc4 = spectro::new_fc4(spectro.n3n, f4);
	let output = spectro.run(spectro::Derivative::Quartic(fc2, fc3, fc4));

pbqff/crates/pbqff/src/coord_type/cart.rs

Lines 135 to 142 in 221be4e

	let freq_dir = &dir.as_ref().join("freqs");
	let (fc2, f3, f4) = self.make_fcs(
	targets,
	&energies,
	&mut fcs,
	n,
	Derivative::quartic(n),

although maybe we ended up in the same place.

[ntBre]

I threw the question about F12-XZ-cCR to Codex as well because I'm again a bit fuzzy on the details, and it said:

yes, qffbuddy is doing what the paper describes. The
HFF is calculated using the composite cCR energy, not merely the uncorrected
CCSD(T)-F12/cc-pCVXZ-F12 energy. The reference geometry is optimized at the core-
correlated CCSD(T)-F12/cc-pCVXZ-F12 level, and the HFF single-point energies add
the scalar-relativistic correction obtained from the difference between
relativistic and nonrelativistic canonical CCSD(T)/cc-pVTZ-DK calculations. For a
non-hybrid F12-XZ-cCR QFF, that composite energy is used throughout the force
field. For the hybrid F12-TcCR+XZ method, it is used for the harmonic terms, while
the cubic and quartic terms are taken from the cheaper, uncorrected F12-XZ QFF.
This is the split represented by template and hybrid_template in qffbuddy.

I think only the methods with an additional +XZ indicate a different level of theory for the HFF. Codex and I were referencing this paper.