Thank you, this is too clear and very helpful.

Hi QMQDCHEM ,
Could you share your BAND.DAT file so that I can do some quick test?
Thanks

Thank you! This is much more comfortable.

Dear Alexander,
We runned a few tests, and, indeed, we found the same behavior. This, though, is not due to an error in the definition in the basis set, but rather to some formatting issue, since the Cs goes up to P4 the code expects that also the Iodine pseudo goes up to P4 and fills the missing coefficients and exponents with zeros.

Luckily there is an easy workaround to this, it is sufficient to flip Cs and I definition in the geometry and the results are the same as the one you obtained by defining the basis set in the input as you can see from my test.

I will leave you here two output snippets hoping that they help clarifying the issue:

Cs defined before I in the geometry section INPUT COORDINATES ATOM AT. N. COORDINATES 1 38 5.000000000000E-01 5.000000000000E-01 5.000000000000E-01 2 55 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 3 53 0.000000000000E+00 5.000000000000E-01 5.000000000000E-01 [...] ******************************************************************************* *** PSEUDOPOTENTIAL INFORMATION *** ******************************************************************************* ATOMIC NUMBER 38, NUCLEAR CHARGE 10.000, PSEUDOPOTENTIAL TYPE EXPONENT COEFF. N EXPONENT COEFF. N P0 TMS 6.9334610 135.2710429 0 4.1140038 17.9440714 0 P1 TMS 7.2168166 29.4380813 0 7.1736962 58.8806749 0 3.0227988 4.9362827 0 2.8656990 9.7233521 0 P2 TMS 6.3215146 11.9072392 0 6.3914995 17.8595514 0 1.7697266 2.1991802 0 1.6367717 2.8935709 0 P3 TMS 4.2441984 -5.5093333 0 4.2291645 -7.3046417 0 ATOMIC NUMBER 55, NUCLEAR CHARGE 9.000, PSEUDOPOTENTIAL TYPE EXPONENT COEFF. N EXPONENT COEFF. N P0 TMS 4.0811192 84.5477223 0 2.4215224 16.6540350 0 P1 TMS 5.5339726 52.3496307 0 5.5067944 104.6994132 0 2.2809616 8.8065577 0 2.1034905 17.6166111 0 P2 TMS 1.8131494 5.2689855 0 1.8077217 7.9036419 0 0.8729040 1.3364313 0 0.8587203 2.0056513 0 P3 TMS 5.2170839 -16.4976543 0 5.1481965 -23.3081313 0 1.5805995 -2.2368273 0 1.3478959 -2.2269420 0 P4 TMS 1.8077398 -2.5041987 0 1.8050613 -3.1382445 0 ATOMIC NUMBER 53, NUCLEAR CHARGE 25.000, PSEUDOPOTENTIAL TYPE EXPONENT COEFF. N EXPONENT COEFF. N P0 TMS 40.0333760 49.9896490 0 17.3005760 281.0065560 0 8.8517200 61.4167390 0 P1 TMS 15.7201410 67.4162390 0 15.2082220 134.8076960 0 8.2941860 14.5665480 0 7.7539490 28.9684220 0 P2 TMS 13.8177510 35.5387560 0 13.5878050 53.3397590 0 6.9476300 9.7164660 0 6.9600990 14.9775000 0 P3 TMS 18.5229500 -20.1766180 0 18.2510350 -26.0880770 0 7.5579010 -0.2204340 0 7.5974040 -0.2216460 0 P4 TMS 0.0000000 0.0000000 0 0.0000000 0.0000000 0 I defined before Cs in the geometry section INPUT COORDINATES ATOM AT. N. COORDINATES 1 38 5.000000000000E-01 5.000000000000E-01 5.000000000000E-01 2 53 0.000000000000E+00 5.000000000000E-01 5.000000000000E-01 3 55 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 [...] ******************************************************************************* *** PSEUDOPOTENTIAL INFORMATION *** ******************************************************************************* ATOMIC NUMBER 38, NUCLEAR CHARGE 10.000, PSEUDOPOTENTIAL TYPE EXPONENT COEFF. N EXPONENT COEFF. N P0 TMS 6.9334610 135.2710429 0 4.1140038 17.9440714 0 P1 TMS 7.2168166 29.4380813 0 7.1736962 58.8806749 0 3.0227988 4.9362827 0 2.8656990 9.7233521 0 P2 TMS 6.3215146 11.9072392 0 6.3914995 17.8595514 0 1.7697266 2.1991802 0 1.6367717 2.8935709 0 P3 TMS 4.2441984 -5.5093333 0 4.2291645 -7.3046417 0 ATOMIC NUMBER 53, NUCLEAR CHARGE 25.000, PSEUDOPOTENTIAL TYPE EXPONENT COEFF. N EXPONENT COEFF. N P0 TMS 40.0333760 49.9896490 0 17.3005760 281.0065560 0 8.8517200 61.4167390 0 P1 TMS 15.7201410 67.4162390 0 15.2082220 134.8076960 0 8.2941860 14.5665480 0 7.7539490 28.9684220 0 P2 TMS 13.8177510 35.5387560 0 13.5878050 53.3397590 0 6.9476300 9.7164660 0 6.9600990 14.9775000 0 P3 TMS 18.5229500 -20.1766180 0 18.2510350 -26.0880770 0 7.5579010 -0.2204340 0 7.5974040 -0.2216460 0 ATOMIC NUMBER 55, NUCLEAR CHARGE 9.000, PSEUDOPOTENTIAL TYPE EXPONENT COEFF. N EXPONENT COEFF. N P0 TMS 4.0811192 84.5477223 0 2.4215224 16.6540350 0 P1 TMS 5.5339726 52.3496307 0 5.5067944 104.6994132 0 2.2809616 8.8065577 0 2.1034905 17.6166111 0 P2 TMS 1.8131494 5.2689855 0 1.8077217 7.9036419 0 0.8729040 1.3364313 0 0.8587203 2.0056513 0 P3 TMS 5.2170839 -16.4976543 0 5.1481965 -23.3081313 0 1.5805995 -2.2368273 0 1.3478959 -2.2269420 0 P4 TMS 1.8077398 -2.5041987 0 1.8050613 -3.1382445 0

I hope this helps

Hi,

We have run some tests and we have identified the origin of the problem. The calculation fails in the evaluation of the Fermi energy in the NEWK option (so before getting to the BOLTZTRA step) because of large memory requirements due to a very large number of k-points being asked and because of the replicated-memory parallel implementation of that bit of code.

In that part of the code, with Pproperties (parallel version), data are replicated in memory by each process.

We have run tests on this system in parallel with different number of processes (on a computing node with 128 CPU cores) and for different shrinking factor parameters of the NEWK keyword. Results are summarized in the table below:

analysis.png

"ok" marks combinations for which the calculation run without errors. The trend is clear and can be rationalized as follows:

reducing the number of k points reduces memory requirments reducing the number of MPI processes effectively increases the available memory/process

Hope this clarifies things and helps find a way forward,

Hi,

I used space group 186 for ZnO that corresponds to the one you mention. In the character table printed by CRYSTAL only those irreps that are actually used to build symmetry-adapted Bloch functions are shown. I have updated my original post above to show the irrep labels in the character tables, which match those found in the printing of the eigenvalues.

Hope this clarifies things,

Dear CRYSTAL community,

The CRYSTAL Team is heading to Brazil! 🇧🇷✨
Next week (26th Jan – 29th Jan 2026), we will be in Volta Redonda (Rio de Janeiro state) for the

logo.jpeg

QMMC 2026 will be hosted at the Universidade Federal Fluminense and it will be an exciting journey through quantum modelling of materials, covering a wide range of topics in computational chemistry and condensed matter physics.

We are truly excited to be in Volta Redonda and to share knowledge, experience, and, of course, to spread the CRYSTAL verb!

More information about the school can be found on the event website.

See you in Brazil!

Poking further, since the last error is that of MPI and not of CRYSTAL, I launched the same input with NEWK 8 8 as a single-core process (with Pproperties), and this time the error was different:

85-C( 3 2 4) 86-C( 4 2 4) 87-C( 5 2 4) 88-C( 3 3 4) 89-C( 4 3 4) 90-R( 4 4 4) ERROR **** PROJVR **** NULL COMPONENT 0.222045E-15 0.100000E-07

for

ROTREF ATOMS 10 # this is Manganese 0 0 0 # I also tried this with different cell indices 12 # P5 0 0 0 8 # P2 0 0 0

and then

85-C( 3 2 4) 86-C( 4 2 4) 87-C( 5 2 4) 88-C( 3 3 4) 89-C( 4 3 4) 90-R( 4 4 4) ERROR **** RHOLSK **** BASIS SET LINEARLY DEPENDENT

for the following input:

ROTREF ATOMS 10 1 0 0 12 1 0 1 8 1 1 0

As well as for other indices

JohnKendrick
I have had the same problem in a slightly different task. I have found a solution, but it may not work for your particular task: I have reduced the number of k-points.
I am modelling Mn5+ ion in a AlPO4 cell in which one P is replaced by Mn.
My input is

EXTERNAL FREQCALC ANALYSIS INTENS INTRAMAN INTCPHF END END BASISSET POB-TZVP-REV2 DFT WC1LYP SPIN END EXCHSIZE 44000000 BIPOSIZE 44000000 SHRINK 6 6 TOLINTEG 9 9 9 12 20 MAXCYCLE 400 TOLDEE 10 SPINLOCK 2 -6 ATOMSPIN 1 10 1 SLOSHING END

The same input with

SHRINK 8 8

Produced the same error as yours, regardless of convergence tools (DIIS / NODIIS (in CPHF block) / buffer sizes / etc):

ELECTRIC FIELD APPLIED ALONG CARTESIAN DIRECTIONS XX [some lines omitted] BECKE WEIGHT FUNCTION RADSAFE = 2.00 TOLERANCES - DENSITY:10**- 6; POTENTIAL:10**- 9; GRID WGT:10**-14 RADIAL INTEGRATION - INTERVALS (POINTS,UPPER LIMIT): 1( 75, 4.0*R) 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) TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MOQGAD TELAPSE 7695.33 TCPU 7664.63 forrtl: severe (67): input statement requires too much data, unit 81, file /scratch/tmp_p267436_student/fort.81.pe11

All running on a single machine: Dell EMC C6400 Server (2x20-Core Intel XEON Gold 6148 2.40GHz, 192GB RAM, 3x480GB SSD).

Hi Jack,
compiling from objects on Apple Silicon is possible, but there are two critical requirements:

You must use OpenMPI built with the same GNU Fortran version used to compile the object files (in particular gfortran 12.1)

You must use the MPI compiler wrappers (mpif90, mpicc, mpicxx) instead of the plain compilers for the final linking stage

The default include file is almost correct. The only necessary changes are the compiler definitions. Replace the first lines with:

F90 = mpif90 LD = $(F90) PLD = mpif90

Keep the rest unchanged.

Important notes

The OpenMPI you use must be built against gfortran 12.1. You can check with:

mpif90 --show

or

mpif90 --version

and verify that it points to gfortran-12.

Do not mix different GNU Fortran versions (e.g. gfortran 13 or Apple clang).
A mismatch here is the most common cause of runtime failures.

Hi Jack,

the parallel version of CRYSTAL23 shipped for Apple Silicon is built with OpenMPI 4.1.1, therefore it is essential that the code is executed using the same mpirun version (or at least the same major version, ie 4.x.x).

If a different OpenMPI installation is used (for example the Homebrew 5.x one), the program may start but fail internally, leading to errors such as the abnormal SCF termination you originally observed.

Concerning the message ls: No match. this is just a standard shell warning printed when the ls command does not find the files it is looking for.
It is probably produced by the run script when it tries to list some output or scratch files that may not exist (for example if the job stops before all files are written).

It is not an error of CRYSTAL itself.

You may try searching inside the run script to locate the line containing the ls command. From the path and filename it is trying to list, you can understand whether the file is genuinely not produced or if the script is looking in the wrong path.

Hope this helps.

dmitoli Thank you, it was a formatting problem with my notebook, now the issue has been solved.

Best regards
Wang

Let me just add that we do have a development version of the code for HSE+SOC, which we plan to include in the next release.

Hi,

Let me just add that of course in P1 there are no special high-symmetry points to guide you in the definition of the path, so you need to be a little creative. For instance, you can start from Gamma (0 0 0) and go to the edge of the FBZ along the b1 reciprocal lattice (1/2 0 0), to then go to (1/2 1/2 0), then to (0 1/2 0) then back to Gamma (0 0 0) and then to the edge along the b3 reciprocal lattice (0 0 1/2). Or something else! 🙂

Hi andrejsc,
There is no "hardcoded" switch in the TOLINTEG keyword, I suspect that the quote from the manual comes from the "old days", when a threshold of \( 10^{-20} \) was considered absurdly small. However, with advances in hardware and the increasingly complex structures we want to compute, such thresholds are sometimes necessary.

A small note: don't warry about going past machine precision with these small numbers. All evaluations of integral thresholds in the code are performed at the logaritmic level (ie only on the exponent), so you should be fine even with a value of 1 million in TOLINTEG (though maybe not fine in terms of computation time)

If you fear that a reduced symmetry might influence the thermodynamic results, you can simply run the calculation without imposing any symmetry at all. For a 3-atom system like CO₂ the computational cost is so minimal that you can run both symmetry-free and \(D_{4h}\) calculations. Comparing the two results should give you an indication of whether the imposed symmetry has any impact on the quantities you are interested in.

Hi Orion,
Looking at your input, I think you want to plot the orbital projected electronic DOS (in CRYSTAL this is done with the DOSS keyword, PDOS in CRYSTAL normally refers to phonon DOS).

To obtain separate px py pz projections, the correct approach is:
set the first parameter of DOSS to 2 (because you want two separate projections).
After the main DOSS line, you must then add one line per projection, specifying which atomic orbitals (AOs) belong to each projection.
Each line has the form:

n AOs index_1 index_2 ... index_n

The AO indices correspond to the ones printed in the CRYSTAL output. To find them, search in the SCF output for "LOCAL ATOMIC FUNCTIONS BASIS SET", you will see a table with all the basis set, that will look like this:

******************************************************************************* ATOM X(AU) Y(AU) Z(AU) N. TYPE EXPONENT S COEF P COEF D/F/G COEF ******************************************************************************* 1 O -2.793 -4.838 -4.110 1 S 8.589E+03 1.895E-03 0.000E+00 0.000E+00 1.297E+03 1.439E-02 0.000E+00 0.000E+00 2.993E+02 7.073E-02 0.000E+00 0.000E+00 8.738E+01 2.400E-01 0.000E+00 0.000E+00 2.568E+01 5.948E-01 0.000E+00 0.000E+00 3.740E+00 2.808E-01 0.000E+00 0.000E+00 2- 5 SP 4.212E+01 1.139E-01 3.651E-02 0.000E+00 9.628E+00 9.208E-01 2.372E-01 0.000E+00 2.853E+00-3.274E-03 8.197E-01 0.000E+00 6- 9 SP 9.057E-01 1.000E+00 1.000E+00 0.000E+00 10- 13 SP 2.556E-01 1.000E+00 1.000E+00 0.000E+00 14- 18 D 1.292E+00 0.000E+00 0.000E+00 1.000E+00

The indeces are the numbers befor each orbital shell (1 2-5 6-9 10-13 14-18).
The notation 2-5 means indices 2 through 5, because that shell contains 4 AOs.
From these blocks you must identify the p-type orbitals for each atom (Zn and Cl). The p shells will appear in groups of three basis functions (px, py, pz).

So, in your system:

Locate the Zn atom in this printed basis list Identify the block corresponding to its p orbitals Extract the AO indices for px, py, pz Do the same for Cl

Then you define one projection per line. For example, structurally something like:

NEWK 12 12 1 0 DOSS 2 100 3 6 1 12 0 3 <list of Zn p AOs> 3 <list of Cl p AOs> END

Keep in mind that, if your basis is double- or triple-Z, each p shell appears as a separate triplet of AOs (one triplet per Z). To correctly project all p you must include the corresponding AOs from every p-shell triplet for that atom.

Let me know if you manage to do this, or if you need further help.

Dear Jonas,

I tried using pymatgen to extract the point-symmetry information from your .xyz file (see the Python script below):

from pymatgen.core import Molecule from pymatgen.symmetry import analyzer bigstructure = Molecule.from_file("yourfile.xyz") PGstructure = analyzer.PointGroupAnalyzer(bigstructure) sym_mol = PGstructure.get_equivalent_atoms() print(sym_mol["eq_sets"])

This returns a Python data structure containing the symmetry-irreducible sets of atoms (only 6 for this system, which is indeed of Ih point group!).
When preparing the CRYSTAL input, be careful with the orientation of your asymmetric unit. In my case, for example, I had to change the sign of the x and y coordinates to make the symmetry consistent with CRYSTAL’s conventions.
Icosahedral point groups are available in CRYSTAL (Ih is point group number 47 in CRYSTAL), so the input fort this molecular cage reads:

Symm. structure MOLECULE 47 6 8 2.605032231 -11.914271806 11.762689798 6 4.344538598 15.664236912 4.366798884 6 3.427479683 14.428996321 8.326818862 6 -8.906580884 1.370529673 14.411150640 1 2.632960331 14.962480562 7.832453804 5 -8.776255231 -2.916256114 14.200279302 COORPRT TESTGEOM END

aerba Christmas is already in the air indeed!

Alright, it sounds like the difficulties originate from the system itself, rather than the particular CRYSTAL input.

Thanks again for all the help!

Best regards,
Wim

aerba, Thank you very much for your reply!

Unfortunately, I don't speak English very well, which is why I may have phrased my question incorrectly.

My question was that the function OPTGEOM does not work correctly with RUNCONFS. I had already tried the steps you suggested.

However, I was able to figure out the problem. For geometric optimization, I just needed to specify the keyword more correctly:

CRYSTAL 0 0 0 194 3.065 17.656 4 22 0 0 0.5 14 0 0 0.75 22 0.666666 0.333333 0.364919 6 0.333333 0.666666 0.427507 SCELCONF 1 0 0 0 1 0 0 0 1 RUNCONFS ATOMSUBS 22 273 OPTGEOM MAXCYCLE 500 ENDOPT END END

I'll leave it here, maybe it will be useful to someone in the future.

Nevertheless, thank you very much for your responsiveness