​Dear Prof. Erba,

Thank you for your clear explanation.

Hi,

With CRYSTAL (from the PROPERTIES module actually) one can output the Hartree+EN potential in the all-electron case in 2D or 3D grids with POTM and POT3 keywords, respectively. In the latter case, the output is in .cube format. But I am afraid that currently there is no keyword to plot the XC part of the potential on a grid.

Dear CRYSTAL users,

I'm working with a molecular crystal, and I need to calculate the forces acting on the atoms in different configurations. We want to verify if the BSSE correction is important in the determination of forces. In addition, is the calculation of forces with MOLEBSSE comparable with the calculation without this correction, since we don't employ the reciprocal space in the first one?

We need some light on the right approach for the calculation with MOLEBSSE. We've tested a few configurations, and the resulting forces varie considerable. For instance, do we need to "isolate" all the molecules of the cell, or just one? Also, how big should the distance R and the number of neighbours to be searched be? In our system, the conventional cell has 4 diatomic molecules.

Thank you so much! I'll give this a try; it helps a lot.

Hi Alessandro, thanks for the info! Much appreciated, Chris.

It was a great week! Let me share a group picture from the event:

QMMC2026_Volta_Redonda.jpeg

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,

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)