Is there a way to print a list of computed force constants?

[phi-l-l-ip-thomas]

Hi Brent,

I just built and ran my first test run of PBQFF, and I am really impressed at how much effort you put into automating a process that can consume many-an-hour for a computational spectroscopist!

I have a question (and perhaps a feature request?): I would like to print out the list of quadratic / cubic / quartic force constants in normal mode representation produced in a PBQFF run. From the documentation of the ANPASS and SPECTRO packages, I think that ANPASS generates these in some form in the fort.15, fort.30, and fort.40 files which SPECTRO then reads to perform the VPT2 analysis, but I am struggling to figure the structure and units. Is there a straightforward way (perhaps via input flags of one of the sub-packages) to parse these into an ordered list? Especially helpful would be something like that generated in Gaussian outputs during 'freq=anharm' runs:

 ........................................................
 :      QUADRATIC FORCE CONSTANTS IN NORMAL MODES       :
 :                                                      :
 : FI =  Frequency [cm-1]                               :
 : k  =  Force Const.[ attoJ * amu(-1) * ang(-2) ]      :
 : K  =  Force Const.[ Hartrees * amu(-1) * bohr(-2) ]  :
 :......................................................:

     I  J      FI(I,J)      k(I,J)       K(I,J)

     1  1       3185.49899     5.97868     0.38401
     2  2       3185.34846     5.97811     0.38398
     ...
 ........................................................
 :        CUBIC FORCE CONSTANTS IN NORMAL MODES         :
 :                                                      :
 : FI =  Reduced values [cm-1]  (default input)         :
 : k  =  Cubic Force Const.[AttoJ*amu(-3/2)*Ang(-3)]    :
 : K  =  Cubic Force Const.[Hartree*amu(-3/2)*Bohr(-3)] :
 :......................................................:

     I  J  K     FI(I,J,K)    k(I,J,K)    K(I,J,K)

     1  1  1     805.84435    14.70125     0.49969
     2  2  1     808.62858    14.75135     0.50139
       ...
 ........................................................
 :                                                      :
 :       QUARTIC FORCE CONSTANTS IN NORMAL MODES        :
 :                                                      :
 : FI =  Reduced values [cm-1]  (default input)         :
 : k  =  Quartic Force Const.[AttoJ*amu(-2)*Ang(-4)]    :
 : K  =  Quartic Force Const.[Hartree*amu(-2)*Bohr(-4)] :
 :......................................................:

  I  J  K  L   FI(I,J,K,L)  k(I,J,K,L)  K(I,J,K,L)

  1  1  1  1     198.22710    35.15132     0.63224
  2  2  1  1     199.12180    35.30831     0.63507
     ...

Also, I tried to check the original SPECTRO manual but the link appears to no longer load the manual. If you know of a way to generate this list (or if you feel like adding this in a future release), then I would be very grateful!

[ntBre]

Hi Phillip, thanks for the kind words and for opening an issue!

I thought there was an existing way to do this, but I'm not seeing it from poking around in the code a bit. I did manage to find a working link to the spectro manual though: https://r410berry.com/static/media/spectro_manual.8f984ff89e638ac74583.pdf

If you have the fort.15,30,40 files, the formulas for the indices are given on pages 12-13 in the Cartesian Fields section. Unfortunately it's still a bit inscrutable (and in Cartesian coordinates rather than normal coordinates), and the manual has this to say about the units:

The energy units are mdyne*A = attojoules

which based on the Gaussian output is maybe not as unusual as I always thought.

Anyway, I think it should be pretty easy to add this, and I might even have some time to look at it this week!