Hi,

After running a few tests, we found a bug towards the end of the EOS module due to the size mismatch of different arrays (specifically, in the full list of direct lattice vectors of the equilibrium configuration and maximally compressed configuration). While this requires an actual fix to the code, we also found a workaround solution for you. Indeed, if the maximally compressed configuration has the same list as for the equilibrium configuration, everything is fine. Thus, in your case, it was enough to explore a maximum compression of 5% (rather than 8%) to be able to run smoothly the job.

Please, find here a working input file (where we have also reduced the number of explored volumes to 4+1 to speed up the calculation without resulting in any significant loss of accuracy).

job314 I wouldn't be THAT optimistic

Dear CRYSTALlers,

The Theoretical Chemistry Group of the University of Torino and the Department of Chemistry and the Thomas Young Centre at Imperial College London are organising the

MSSC2025 Summer School on Ab initio Modelling in Solid State Chemistry

to be held in Torino (Italy), 8-12 September 2025. The School is designed for Master and Ph.D. students, as well as for post-docs and researchers who have an interest in getting or strengthening a background in Computational Solid State Chemistry, Physics, Materials Science, Surface- and Nano-Science.

The week-long School consists of morning lectures and afternoon hands-on tutorial sessions, where the formal framework and functionalities of the CRYSTAL electronic structure package will be explored.

Registration opens on March 15 2025. See the School website for further details.

The School takes place in the elegant city of Torino (Italy):

­
While we strongly encourage in-person participation, we also offer the possibility to attend remotely through streaming of morning lectures and afternoon hands-on tutorials.

Participants will have the opportunity to present their research at a poster session.

See you in Torino!

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.

Hi,

Elastic constants of a 3D lattice in CRYSTAL are defined and computed as:

$$
C_{ij}^\textup{3D} = \frac{1}{V} \left( \frac{\partial^2 E}{\partial \eta_i \partial \eta_j}\right)
$$

where the strain \( \eta \) is dimensionless and thus the elastic constants have units of \( \textup{energy/length}^3 \), which corresponds to force/surface (i.e. a pressure).

For 1D and 2D periodic lattices, as the volume \(V\) is not defined, one may divide by the length \(l \) and area \( A\) of the lattice cell instead:

$$
C_{ij}^\textup{1D} = \frac{1}{l} \left( \frac{\partial^2 E}{\partial \eta_i \partial \eta_j}\right) \qquad \textup{and} \qquad C_{ij}^\textup{2D} = \frac{1}{A} \left( \frac{\partial^2 E}{\partial \eta_i \partial \eta_j}\right)
$$

However, in CRYSTAL for 1D and 2D lattices we do not divide by \(l \) or \( A\) , and just define and compute the elastic constants as:

$$
C_{ij}^\textup{1D and 2D} = \left( \frac{\partial^2 E}{\partial \eta_i \partial \eta_j}\right)
$$

with units of energy. So, starting from the constants printed in the CRYSTAL output all you have to do is divide by the area of the 2D cell if you need them expressed in \( \textup{energy/length}^2 \).

We are running some tests and I think we are close to identifying the problem. We’ll probably have a solution by the beginning of next week.

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.

job314 I have checked your attached output files. The OPTGEOM calculation converges to the optimized structure in 29 steps. The PREOPTGEOM calculation within the EOS option converges to the same optimized structure in 25 steps. A lower number of steps in EOS is likely due to the higher accuracy of the computed forces because of tighter default settings in EOS than in OPTGEOM (for instance, TOLDEE is set to 8 in EOS and to 7 in OPTGEOM).

The EOS.out file you shared looks perfectly fine. I can not find any failed geometry optimizations there. After the pre-optimization, the actual constant-volume optimizations are performed (converged in 40, 19, 16, 16, 8 steps, respectively). The output is not complete because I guess the calculation was still running.

Hi,

Find below a modified input file for which the SCF converges in 25 iterations:

Na2Si2O5, Rakic Physics and Chemistry of Minerals 29 (2002) 477-484 - Ale mod
CRYSTAL
0 0 0
14
4.725582 23.296569 7.936255 90.15
18
11 .2299 .26138 .5247
11 .2743 .51408 .3548
11 .2534 .47281 -.1019
11 .7548 .28108 .2783
14 .2961 .36325 .1966
14 .6843 .34077 .6429
14 .1831 .40821 .5373
14 .7948 .38941 -.0185
8 .1210 .3844 .0352
8 .2339 .3015 .2523
8 .6209 .3690 .1429
8 .2440 .4105 .3392
8 .7431 .3414 -.1590
8 .7456 .2819 .5655
8 .3573 .3561 .6162
8 -.1423 .3920 .5597
8 .2452 .4667 .6166
8 .7336 .4511 -.0767
EOS
RANGE
0.95 1.05 8
PREOPTGEOM
MAXCYCLE
500
END
BASISSET
POB-TZVP-REV2
DFT
cam-B3LYP
XLGRID
END
TOLINTEG
10 10 10 10 20
SHRINK
6 6
BIPOSIZE
11202400
EXCHSIZE
11202400
MAXCYCLE
200
END

We have increased the values of TOLINTEG, used an isotropic shrinking factor, and re-activated the DIIS accelerator.

job314 I see, interesting! I don't seem to be able to find the input/output files

Yes, the pre-optimization of the EOS is indeed a full structural relaxation where cell/volume/atomic positions are all relaxed at once. If you want me to have a closer look at this case, please send me the input/output files.

Hi,

Yes, it looks strange, and 102 is a funny number.

Could you try to run a regular optimization (i.e. with OPTGEOM and not within the EOS) and let me know if the same happens?

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:

1. Geometry optimization to fully relax atomic positions within the cell (OPTGEOM keyword);
2. Harmonic frequency calculation (FREQCALC keyword);
3. 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);
4. 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

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.