Vibrational Spectroscopies: IR, Raman, INS

I see, thank you

​Dear Prof. Erba,

Thank you for your clear explanation.

Thank you, this is too clear and very helpful.

Hi,

job314 said in SCANMODE io error Read_int_1d:

So it turns out having fort.13 and fort.20 is not optional for restart.

Yes, the lack of these two files was the origin of the I/O error.

job314 said in SCANMODE io error Read_int_1d:

PS. I would also like to ask developers to allow uploading compressed files to this forum

Now also .zip, .tar, .tgz, .tar.gz files can be uploaded.

Cheers,

Either way, seems you got your geometry for further optimization 🙂
Cheers,
A

aerba I've noticed for my systems, these SCF convergence issues appear when I try to compute for intensities in a restart job. Then when removing intensities, they converge fine and progress forward.

Is it possible to finish the frequency calculations first and then compute for the Raman intensities afterwards? Would that be a work around? Maybe with RAMANREA - because I have the TENS_RAMAN.DAT, but it just gets stuck at some point during HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
FORCE CONSTANT MATRIX - NUMERICAL ESTIMATE
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
and doesn't progress with restarts

Hi,

I understand your confusion, which originates from the fact that PRESSURE is an old option of FREQCALC that indeed does not apply any pressure. With this option, the structure is unchanged, the forces are unchanged, and ultimately, the harmonic vibration frequencies are unchanged. The only bit of information that is affected by this keyword is the printed value of the PV term in the thermodynamic analysis.

So, this PRESSURE keyword within FREQCALC should be used in just one way: to set the value of pressure corresponding to the structure used to run the harmonic frequency calculation, so as to have the PV term entering the thermodynamic functions (enthalpy and Gibbs free energy) right. Let me make a couple of examples:

You perform a full structural relaxation (atomic positions + lattice parameters). Thus, you have a zero pressure (p=0) structure. Then you compute the harmonic frequencies. If you are interested in thermodynamic potentials (i.e. enthalpy, Gibbs, Helmholtz, etc.) you should use the PRESSURE keyword and set the pressure to zero to get the PV term right in the output (otherwise by default it would be computed at p = 1 atm - do not ask me why).

You perform a pressure-constrained geometry optimization with the EXTPRESS keyword of OPTGEOM (let's say at 29.457 GPa as in your example). Then you compute harmonic frequencies with FREQCALC on this compressed structure. Now, again, if you are interested in thermodynamic potentials (i.e. enthalpy, Gibbs, Helmholtz, etc.) you should use the PRESSURE keyword and set the pressure to 29.457 GPa to get the PV term right in the output.

So, as I said at the beginning, the PRESSURE keyword here does not change anything on the structure/forces. It only allows you to tell the program what is the pressure of the structure you are working on, so that the printed PV term is right.

A more formally-sound way to combine temperature and pressure is provided by the quasi-harmonic approximation: See Eq. (13) of the tutorial page on the Quasi-Harmonic Approximation (QHA).

Hope this helps clarifying things a little,

Hi Jonas,

Would you mind adding the INPUT file as well, for easier checks?

It generally doesn't look too bad, I would wait until the calculations finishes and see how the results look like. There is nothing obviously "wrong", but you can always try tightening the convergence criteria. See earlier post for a great overview of options (especially if you notice in the end that you have negative frequencies appearing!):
Earlier post on "Imaginary frequencies", see comment by Prof. Erba
I personally wouldn't run a geometry optimization together with a frequency calculation, I would split them in too, but I guess that is a matter of taste as the code can handle it : )

Cheers,
Aleks

It makes sense, thank you.

Hi Erba,
Honestly, I have ANDERSON acceleration on because I inherited a file with these keywords written in when I was taught how to use CRYSTAL for frequency calculations - and worked pretty well so far with our other systems so I thought to leave it alone (though it was more for pharmacuetical systems). Thanks for the suggestion, I'll try it out if with the original TOLINTEG 10 10 10 10 20 and see what happens.

Thanks,
Peter

Hi,

In CRYSTAL, the most general syntax to specify exchange-correlation (xc) functionals is, within the DFT input block through the EXCHANGE and CORRELAT keywords as:

DFT EXCHANGE label of exchange functional (e.g. VBH, BECKE, PBE, ...) CORRELAT label of correlation functional (e.g. VWN, LYP, PBE, ...) ENDDFT

For certain standard combinations of exchange and correlation functionals, we have implemented single keywords. For instance:

DFT BLYP ENDDFT

is equivalent to

DFT EXCHANGE BECKE CORRELAT LYP ENDDFT

Another example (here the "XC" letters are appended at the end of the name to reflect that the functional is used for both the exchange and correlation part):

DFT PBESOLXC ENDDFT

is equivalent to

DFT EXCHANGE PBESOL CORRELAT PBESOL ENDDFT

For PBE, there are three equivalent ways to define it in CRYSTAL:

DFT PBEXC ENDDFT

or

DFT PBE ENDDFT

or

DFT EXCHANGE PBE CORRELAT PBE ENDDFT

There isn't a general rule, I'm afraid. This syntax aspects are discussed at page 134 of CRYSTAL23 User's Manual.

Hope this clarifies things a little,

I tried rerunning it with fewer nodes - thought it is some parallel issue. A problem again

ANGULAR INTEGRATION - INTERVALS (ACCURACY LEVEL [N. POINTS] UPPER LIMIT):
1( 4[ 86] 0.2) 2( 8[ 194] 0.5) 3( 12[ 350] 0.9) 4( 16[ 974] 3.5)
5( 12[ 350]9999.0)
CYCLE 0 ALPHA 227.814788 EPSILON 1.894274 DELTA 2.2781E+02
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 1902.88 TCPU 1885.42
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 1904.27 TCPU 1886.79
DIIS TEST: 0.61205E+01 AT CPHF CYCLE 1 - MIX 60 %
CYCLE 1 ALPHA 257.133404 EPSILON 2.009363 DELTA 2.9319E+01
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2002.77 TCPU 1984.78
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2004.16 TCPU 1986.16
DIIS TEST: 0.71887E+01 AT CPHF CYCLE 2 - DIIS ACTIVE - HISTORY: 2 CYCLES
CYCLE 2 ALPHA 268.265588 EPSILON 2.053062 DELTA 1.1132E+01
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2102.04 TCPU 2083.54
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2103.42 TCPU 2084.92
DIIS TEST: 0.36370E+00 AT CPHF CYCLE 3 - DIIS ACTIVE - HISTORY: 3 CYCLES
CYCLE 3 ALPHA 276.769385 EPSILON 2.086443 DELTA 8.5038E+00
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2202.03 TCPU 2183.04
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2203.42 TCPU 2184.42
DIIS TEST: 0.54051E-01 AT CPHF CYCLE 4 - DIIS ACTIVE - HISTORY: 4 CYCLES
CYCLE 4 ALPHA 278.095061 EPSILON 2.091647 DELTA 1.3257E+00
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2302.12 TCPU 2282.64
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2303.51 TCPU 2284.02
DIIS TEST: 0.85023E-02 AT CPHF CYCLE 5 - DIIS ACTIVE - HISTORY: 5 CYCLES
CYCLE 5 ALPHA 278.435921 EPSILON 2.092985 DELTA 3.4086E-01
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2402.16 TCPU 2382.20
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2403.54 TCPU 2383.57
DIIS TEST: 0.38480E-03 AT CPHF CYCLE 6 - DIIS ACTIVE - HISTORY: 6 CYCLES
CYCLE 6 ALPHA 278.461661 EPSILON 2.093086 DELTA 2.5739E-02
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2502.06 TCPU 2481.62
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2503.45 TCPU 2482.99
DIIS TEST: 0.44991E-03 AT CPHF CYCLE 7 - DIIS ACTIVE - HISTORY: 7 CYCLES
CYCLE 7 ALPHA 278.460154 EPSILON 2.093080 DELTA -1.5071E-03
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2601.70 TCPU 2580.74
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2603.08 TCPU 2582.11
DIIS TEST: 0.36243E-03 AT CPHF CYCLE 8 - DIIS ACTIVE - HISTORY: 8 CYCLES
CYCLE 8 ALPHA 278.473843 EPSILON 2.093134 DELTA 1.3689E-02
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 2701.77 TCPU 2680.26
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT CP_MONMON TELAPSE 2703.15 TCPU 2681.62
DIIS TEST: 0.85073E-04 AT CPHF CYCLE 9 - DIIS ACTIVE - HISTORY: 9 CYCLES
CYCLE 9 ALPHA 278.474328 EPSILON 2.093136 DELTA 4.8487E-04
forrtl: severe (256): unformatted I/O to unit open for formatted transfers, unit 85, file /dev/null
Image PC Routine Line Source
Pcrystal 0000000007374206 Unknown Unknown Unknown
Pcrystal 0000000001BA179E Unknown Unknown Unknown
Pcrystal 0000000000A8038B Unknown Unknown Unknown
Pcrystal 0000000000A63D97 Unknown Unknown Unknown
Pcrystal 0000000000D4DAD1 Unknown Unknown Unknown
Pcrystal 000000000074B942 Unknown Unknown Unknown
Pcrystal 000000000040591E Unknown Unknown Unknown
Pcrystal 00000000004053FD Unknown Unknown Unknown
libc.so.6 000014B7C14295D0 Unknown Unknown Unknown
libc.so.6 000014B7C1429680 __libc_start_main Unknown Unknown
Pcrystal 0000000000405315 Unknown Unknown Unknown

Yes sure. No problem.
I'm copying it here below.
Thank you for the help

Daria

Furfural CRYSTAL 0 0 0 19 7.83542762 10.92181175 14.94182734 90.000000 90.000000 90.000000 33 8 1.907561246186E-01 3.035506248522E-01 4.408415126674E-01 8 -4.953260939445E-01 3.221701934644E-01 3.435788389827E-01 6 2.036695655476E-01 3.000622663369E-01 3.482645168957E-01 6 4.363374159111E-02 2.908170991680E-01 3.122160878064E-01 1 1.581732742072E-02 2.880586476013E-01 2.419150131741E-01 6 -7.402527715589E-02 2.893444053747E-01 3.842913887226E-01 1 -2.109486401803E-01 2.822752418693E-01 3.798189138779E-01 6 2.099722211847E-02 2.968688203007E-01 4.606664256414E-01 1 -9.310483322231E-03 2.937832972654E-01 -4.691825366176E-01 6 3.647591498990E-01 3.132132214593E-01 3.052219591009E-01 1 3.550467304924E-01 3.168499490771E-01 2.320173417360E-01 8 2.803460822693E-01 -3.091160241316E-01 4.442088720927E-01 8 -4.110526913885E-01 -3.641439806750E-01 3.516630126321E-01 6 2.872037985115E-01 -3.695128622796E-01 3.626185432797E-01 6 1.245767077225E-01 -3.991101749711E-01 3.362778399038E-01 1 9.163169834882E-02 -4.467733390550E-01 2.755816936910E-01 6 1.226510841197E-02 -3.550438914407E-01 4.034527830796E-01 1 -1.250093786328E-01 -3.609739865930E-01 4.033488316714E-01 6 1.122081912874E-01 -3.009269785992E-01 4.671425614917E-01 1 8.596288473646E-02 -2.539158941058E-01 -4.709004162095E-01 6 4.474098992917E-01 -3.907868742186E-01 3.203804132925E-01 1 4.349587971852E-01 -4.336723578587E-01 2.540887571847E-01 8 1.068738713332E-01 -3.851604875541E-02 3.194441316968E-01 8 -2.354778627755E-01 -1.170305051641E-02 3.709929081706E-01 6 5.921416878240E-02 -5.109842318475E-03 4.052174692072E-01 6 2.027905991547E-01 1.001107546720E-02 4.563966208009E-01 1 2.038099671676E-01 3.508046669402E-02 -4.738476409789E-01 6 3.453540637796E-01 -1.523126326522E-02 4.003413693017E-01 1 4.790458182320E-01 -1.335338648893E-02 4.173410231626E-01 6 2.808235437569E-01 -4.461452496729E-02 3.184022003331E-01 1 3.367732131084E-01 -7.306049141011E-02 2.558072164743E-01 6 -1.187612609562E-01 6.696688015001E-03 4.250681762815E-01 1 -1.454796120848E-01 3.516201574611E-02 4.940534466886E-01 FREQCALC NUMDERIV 2 INTENS INTRAMAN INTCPHF END END BASISSET POB-TZVP-REV2 DFT B3LYP-D3 XLGRID ENDdft SHRINK 4 4 TOLDEE 10 SCFDIR LEVSHIFT 6 1 FMIXING 30 TOLINTEG 7 7 7 7 14 EXCHSIZE 13220300 BIPOSIZE 13220300 ENDscf

Dear Jonas,

Thanks for reaching out and being one of the most active users of these early days of the forum. Your question gives us the chance to clarify some aspects of the output file that might not be obvious to non expert users. Below, I will refer to your output file.

Harmonic Frequencies and IR intensities

To compute harmonic frequencies and IR intensities (with the default approach of the Berry phase) the input looks like:

FREQCALC INTENS ENDFREQ

In the output file, the following table is printed:

HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH EIGENVALUES (EIGV) OF THE MASS WEIGHTED HESSIAN MATRIX AND HARMONIC TRANSVERSE OPTICAL (TO) FREQUENCIES. IRREP LABELS REFER TO SYMMETRY REPRESENTATION ANALYSIS; A AND I INDICATE WHETHER THE MODE IS ACTIVE OR INACTIVE, RESPECTIVELY, FOR IR AND RAMAN; INTEGRATED IR INTENSITIES IN BRACKETS. CONVERSION FACTORS FOR FREQUENCIES: 1 CM**(-1) = 0.4556335E-05 HARTREE 1 THZ = 0.3335641E+02 CM**(-1) HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH MODES EIGV FREQUENCIES IRREP IR INTENS RAMAN (HARTREE**2) (CM**-1) (THZ) (KM/MOL) 1- 1 0.3488E-07 40.9894 1.2288 (A ) A ( 0.40) A 2- 2 0.6077E-07 54.1037 1.6220 (A ) A ( 0.96) A 3- 3 0.6801E-07 57.2344 1.7158 (A ) A ( 2.45) A 4- 4 0.2238E-06 103.8371 3.1130 (A ) A ( 11.47) A [...]

For each mode (or set of degenerate modes) its eigenvalue (in Ha\(^2\)), harmonic frequency (in cm\(^{-1}\) and THz) and irreducible representation get printed. In addition, labels specifying whether the mode is IR/Raman active are also displayed (A and I indicate whether the mode is active or inactive, respectively).

Raman intensities

Raman intensities can be computed via a coupled-perturbed approach by inserting the INTRAMAN keyword followed by the INTCPHF block in the input deck:

FREQCALC INTRAMAN INTCPHF END ENDFREQ

Raman intensities are computed for each independent component of the polarizability tensor (xx, xy, xz, yy, yz, zz, labeled as "Single Crystal" in the output file) and are also averaged to mimic polycrystalline powder samples (total, parallel polarisation, perpendicular polarisation averages are printed in the output).

POLYCRYSTALLINE ISOTROPIC INTENSITIES (ARBITRARY UNITS) MODES FREQUENCIES I_tot I_par I_perp ---------------------------------------------------------------- 1- 1 40.9894 (A ) 0.46 0.27 0.19 2- 2 54.1037 (A ) 7.35 4.23 3.12 3- 3 57.2344 (A ) 12.79 8.82 3.96 4- 4 103.8371 (A ) 13.66 7.89 5.77 SINGLE CRYSTAL DIRECTIONAL INTENSITIES (ARBITRARY UNITS) MODES FREQUENCIES I_xx I_xy I_xz I_yy I_yz I_zz ---------------------------------------------------------------------------- 1- 1 40.9894 (A ) 0.00 0.37 0.02 0.63 0.00 0.21 2- 2 54.1037 (A ) 3.17 0.69 0.00 4.35 3.66 10.05 3- 3 57.2344 (A ) 3.82 3.54 0.02 3.50 0.03 27.53 4- 4 103.8371 (A ) 2.57 1.81 0.01 16.25 3.34 19.62

For more details on such polycrystalline averages, please refer to sections 8.4 and 8.7 of the CRYSTAL23 manual.

Raman spectrum

A continuous Raman spectrum can be simulated by use of the RAMSPEC block, as in:

FREQCALC INTRAMAN INTCPHF END RAMSPEC END ENDFREQ

The simulated spectrum is printed in an external file named RAMSPEC.DAT that contains several columns: column 1 with frequencies in cm\(^{-1}\), columns 2-4 with polycrystalline intensities (total, parallel, perpendicular), columns 5-10 with single crystal intensities (xx, xy, xz, yy, yz, zz).

Effect of Temperature and Laser wavelength

The effect of temperature and laser wavelength on computed Raman intensities can be accounted for by use of the RAMANEXP keyword, as in:

FREQCALC INTRAMAN INTCPHF END RAMANEXP 298 532 RAMSPEC END ENDFREQ

Here we set 298 K for the temperature and 532 nm for the laser wavelength. This option modifies the values of all computed Raman intensities (in the output and in the RAMSPEC.DAT file accordingly).

Please, note that other properties (harmonic frequencies and IR intensities) are not affected by this option and thus remain unchanged in the output.

Plots

When CRYSPLOT reads the CRYSTAL output file it only plots the total intensity of the polycrystalline powder model.

When CRYSPLOT reads the RAMSPEC.DAT file it plots all components:

Screenshot 2025-03-13 at 12.47.37.png

Other plotting tools can be used to plot specific columns of the RAMSPEC.DAT file (e.g., CRYSTALClear, gnuplot).

Thank you very much Jacques and Alessandro.
Now it works :).
Best regards
Xavier

Hi,

The presence of imaginary frequencies is a sign that the geometry is not a minimum of the potential energy surface (PES). In general, this may due to two main factors:

1) A somewhat loose overall numerical precision in the geometry optimization + harmonic frequencies calculation. Here, it seems that you have already explored a few parameters. We can distinguish between parameters governing the overall numerical precision of the SCF + forces calculations, and those that are specific to the evaluation of the Hessian:

1.1) Precision of SCF+forces

One may increase a bit the thresholds for the screening of two-electron integrals (TOLINTEG keyword), switch-off the bipolar approximation (NOBIPOLA keyword), increase the shrinking factor (SHRINK keyword), use a denser grid for numerical integration of the exchange-correlation term (see XLGRID, XXLGRID keywords), tighten the convergence criteria for the SCF (setting TOLDEE to 10 or 11 for instance), tighten the convergence criteria for the geometry optimization step (see TOLDEG and TOLDEX keywords).

1.2) Numerical evaluation of the Hessian

In many cases, a more numerically stable evaluation of the Hessian is achieved by use of a two-sided finite difference approach (rather than the default one-sided approach). This can be activated with the NUMDERIV keyword within the FREQCALC input block as follows:

FREQCALC NUMDERIV 2 ENDFREQ

2) The presence of symmetry-constraints that prevent the optimizer to get to the minimum of the PES. If this is the case, removing symmetry constraints may be key to reach the minimum. This can be done by use of the SYMMREMO keyword (to be inserted in the geometry input block).

Hi,

By default, thermodynamic properties on top of harmonic frequencies are computed at room temperature. Different values of temperature can be explored by use of the TEMPERAT keyword. To restart the harmonic frequency calculation, use the RESTART keyword. An example is given below:

FREQCALC RESTART TEMPERAT 5 200 600 END

With the input above, thermodynamic properties will be computed and printed at 5 temperatures, equally spaced in the range 200-600 K.