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

Hi,

I am quite familiar with that error message myself: it pops up when the fort.9 unit provided upon the restart is the wrong one (i.e. does not match with the current calculation).

Let me make a general comment on how to best approach these frequency calculations with CRYSTAL based on my experience.

I noticed that you tend to use the PREOPTGEOM option within FREQCALC. Now, this is a very nice feature of CRYSTAL that allows you to run a single job where everything is fully automated: the structure gets optimized, the numerical Hessian computed and diagonalized. Not many programs can do that, to the best of my knowledge. At the same time, while this is very convenient if you can indeed run everything in a single job, it complicates things if you then need to do a restart from a previous incomplete calculation.

If you envisage that this could be the case (maybe because of a wall clock limit on the cluster) then it is preferable (to put it mildly) to do things step by step:

You first run a geometry optimization with something like: [initial geometry] OPTGEOM END You prepare a new input file where you insert the optimized geometry from the previous run as a starting point, and perform a frequency calculation, with something like: [optimized geometry] FREQCALC END If this calculation stops, you can now restart it easily with: [optimized geometry] FREQCALC RESTART END

by providing the required restart files (FREQINFO.DAT from the first frequency calculation, fort.13 unit, and fort.9, to be renamed fort.20 in the new scratch folder).

Personally, I always do things separately, running first the geometry optimization and then the frequency calculation in two separate jobs.

Hope this clarifies things a little,

Whatever works, I guess 😉

OK, here we go. It just is stuck, always the same position in the output

(ceres20-compute-46:0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95)
(ceres24-compute-18:96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191)

export TMPDIR=/local/bgfs/jonas.baltrusaitis/15383115
export TMOUT=5400
export SINGULARITY_TMPDIR=/local/bgfs/jonas.baltrusaitis/15383115

MAX NUMBER OF SCF CYCLES 200 CONVERGENCE ON DELTAP 10**-20
WEIGHT OF F(I) IN F(I+1) 30% CONVERGENCE ON ENERGY 10**-10
SHRINK. FACT.(MONKH.) 6 6 6 NUMBER OF K POINTS IN THE IBZ 64
SHRINKING FACTOR(GILAT NET) 6 NUMBER OF K POINTS(GILAT NET) 64

*** K POINTS COORDINATES (OBLIQUE COORDINATES IN UNITS OF IS = 6)
1-R( 0 0 0) 2-C( 1 0 0) 3-C( 2 0 0) 4-R( 3 0 0)
5-C( 0 1 0) 6-C( 1 1 0) 7-C( 2 1 0) 8-C( 3 1 0)
9-C( 0 2 0) 10-C( 1 2 0) 11-C( 2 2 0) 12-C( 3 2 0)
13-R( 0 3 0) 14-C( 1 3 0) 15-C( 2 3 0) 16-R( 3 3 0)
17-C( 0 0 1) 18-C( 1 0 1) 19-C( 2 0 1) 20-C( 3 0 1)
21-C( 0 1 1) 22-C( 1 1 1) 23-C( 2 1 1) 24-C( 3 1 1)
25-C( 0 2 1) 26-C( 1 2 1) 27-C( 2 2 1) 28-C( 3 2 1)
29-C( 0 3 1) 30-C( 1 3 1) 31-C( 2 3 1) 32-C( 3 3 1)
33-C( 0 0 2) 34-C( 1 0 2) 35-C( 2 0 2) 36-C( 3 0 2)
37-C( 0 1 2) 38-C( 1 1 2) 39-C( 2 1 2) 40-C( 3 1 2)
41-C( 0 2 2) 42-C( 1 2 2) 43-C( 2 2 2) 44-C( 3 2 2)
45-C( 0 3 2) 46-C( 1 3 2) 47-C( 2 3 2) 48-C( 3 3 2)
49-R( 0 0 3) 50-C( 1 0 3) 51-C( 2 0 3) 52-R( 3 0 3)
53-C( 0 1 3) 54-C( 1 1 3) 55-C( 2 1 3) 56-C( 3 1 3)
57-C( 0 2 3) 58-C( 1 2 3) 59-C( 2 2 3) 60-C( 3 2 3)
61-R( 0 3 3) 62-C( 1 3 3) 63-C( 2 3 3) 64-R( 3 3 3)

DIRECT LATTICE VECTORS COMPON. (A.U.) RECIP. LATTICE VECTORS COMPON. (A.U.)
X Y Z X Y Z
13.1430453 0.0000000 0.0000000 0.4780616 0.0000000 0.0000000
0.0000000 11.6066979 0.0000000 0.0000000 0.5413413 0.0000000
0.0000000 0.0000000 21.1989478 0.0000000 0.0000000 0.2963914

DISK SPACE FOR EIGENVECTORS (FTN 10) 53868000 REALS

SYMMETRY ADAPTION OF THE BLOCH FUNCTIONS ENABLED
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT gordsh1 TELAPSE 186.18 TCPU 45.44

Good! The clean way is:

CRYSTAL
1 0 0
P 21/A

I guess the extra "1"s you put are safe as they correspond to the identity operator.

For the purpose of finding the minimum energy structure to then do Raman calculations, it is.

EOS gives you much more than that of course: the p(V) or, equivalently, V(p) relation (i.e. structure as a function of pressure), the bulk modulus K(p), and allows to compute the enthalpy H(p).

Thank you very much, that worked (but not restart)

you're the best, thank you for going extra mile

thank you

very grateful, I will run it again but not sure what to do here if it aborts again

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).

job314 Yes, it should print the geometry of the last optimization step. Please, be aware that the CIFPRT and CIFPRTSYM options are not as general as one would like them to be. For systems where the primitive and crystallographic cells differ, they may result in an incomplete list of atoms. In those cases, I recommend switching symmetry off with the SYMMREMO option in combination with TESTGEOM, just to generate the .cif files correctly.

thank you

job314 Good point! About the shrinking factor: the anisotropic shrinking factor in CRYSTAL does not work properly for those calculations where the symmetry of the system may change (for instance in frequency calculations, FREQCALC, where displaced nuclear configurations are explored, or elastic calculations, ELASTCON, where the lattice is strained, etc.). So in general, I personally tend to avoid using an anisotropic shrinking factor.

However, for symmetry-preserving calculations (such as SCF, OPTGEOM, EOS) the use of an anisotropic shrinking factor should be fine.

thank you

Sorry, I have been working fiercely on these cases, must have overwritten many times. Several EOS cases run out of cycles such as this after 102 cycles abive, including this one when I start range at 0.92 lattice constant (according to the online manual examples). Some of these problems went away when I narrowed the range, e.g. start at 0.96

Hi,

In CRYSTAL23, the anharmonicity of a given vibration mode - or of a set thereof - can be computed via: I) the evaluation of cubic and quartic interatomic force constants (in the basis of the normal modes) followed by II) the solution of the nuclear Schroedinger equation with either the vibrational self-consistent field (VSCF) or vibrational configuration interaction (VCI) method.

Details on the actual implementation of steps I) and II) can be found here:

I) Anharmonic force constants

II) VSCF and VCI for solids

Step-by-step, the procedure is as follows:

Geometry optimization to fully relax atomic positions within the cell (OPTGEOM keyword); Harmonic frequency calculation (FREQCALC keyword); Selection of the normal modes for which the anharmonic correction is to be computed (for instance, in your case, those corresponding to C-H stretching vibrations); Calculation of cubic and quartic interatomic force constants for the selected modes + VSCF (or VCI) calculation of anharmonic states.

As an example, let' s assume that C-H stretching vibrations correspond to modes 15-20 in the list generated from the harmonic calculation (that is, there are 6 different C-H stretching modes). The input for steps 3. and 4. above would read:

FREQCALC RESTART ANHAPES 6 15 16 17 18 19 20 3 0.9 VSCF END

that is, we restart the harmonic frequency calculation, and where 6 is the number of modes, 15 16 17 18 19 20 are the selected modes, and 3 0.9 are two parameters specifying the numerical approach used for the evaluation of cubic and quartic force constants.

Please, note that with such a calculation not only the "intrinsic" anharmonicity of each selected mode is evaluated but also the couplings among all selected modes. If, instead, one just wants to compute the "intrinsic" anharmonicity with no couplings, independent calculations can be run, one per each selected mode.

Visit this page for a tutorial on anharmonic calculations in CRYSTAL