undefined
Hi Ilias,
Even though it’s not documented in the manual, you can use the SMALLDIST keyword (in gemetry section of the input) to convert that error into a warning and allow the calculation to proceed. However, please note that atoms placed very close together can lead to numerical instabilities and linear dependencies in the basis set. Therefore, it’s important to carefully check your geometry and computational setup (i.e., basis set, TOLINTEG, etc.) before continuing.
Hi Chris,
The reason you only see a few k-points printed is indeed that the printing options are not fully supported in parallel execution. In a parallel run, each MPI process handles a subset of the k-points, and only process 0 writes to the output file. As a result, you’ll only see the k-points assigned to process 0, which can be just a few, or even none, depending on how they are distributed.
By default, the code distributes k-points in a round-robin fashion starting from the last process, so the gamma point is typically assigned to the last process and therefore never printed in a parallel execution.
This behavior applies to both options 66 and 67.
Try running the following input:
NEWK
4 4
1 2
66 999
67 999
END
However, you’ll need to run it on a single process for the eigenvectors to be printed correctly.
If the calculation is too expensive to run on a single process, we may need to find an alternative approach to extract those values (though that won’t be straightforward).
Hi,
We provide a distribution of precompiled object files that can be linked on this architecture. You should be able to build and run CRYSTAL23 successfully using those.
We do not have specific performance optimization data or configuration recommendations for the KNL platform. Performance is highly dependent on the system setup and workload characteristics, so we recommend running a few short benchmarks to identify the best configuration for your case.
As regards OpenMP, the optimal number of threads depends on your memory limitations; however, we recommend not exceeding 8 threads per MPI process.
Let me know if you have any further questions.
The AVX (or AVX2/AVX512) option can indeed cause compatibility issues with old processors. Since the object files you are using for compilation were built with AVX2 enabled, there is no way to remove that option and make those files work on processors without AVX support.
However, for this exact reason, we also provide object files compiled without AVX enabled. You should download the appropriate version of the distribution from the CRYSTAL Solutions website (should be named CRYSTAL23 - Precompiled object files for parallel version no avx2) and compile it using the same .inc file you provided.
Let me know if you manage to do it.
Hi,
The intel oneAPI toolkit should in principle be compatible with any Xeon processor, have you compiled the code from the object files or have you used directly the executables distributed?
In the first case, could you provide the .inc file that you used for the compilation?
Dear CRYSTAL community,
From Monday 15 to Thursday 18 December 2025, the 39th Winter School in Theoretical Chemistry will take place in Helsinki, Finland.
This year’s focus is on Electronic Structure Theory! The School is organised by the Department of Chemistry, University of Helsinki, and is open to everyone in academia: students, postdocs, and professors. Best of all, it’s free of charge!
Location:
Helsinki - Finland
City Center Campus of the University of Helsinki
Lecturers:
Jacques Desmarais [Jacques] (University of Turin)
Janus Juul Eriksen (Technical University of Denmark)
Arno Förster (VU Amsterdam)
Christof Holzer (Karlsruhe Institute of Technology)
Ida-Marie Høyvik (Norwegian University of Science and Technology)
Stanislav Komorovsky (Slovak Academy of Sciences)
Frank Neese (Max Planck Institut für Kohlenforschung)
More info & Registration:
For full details, registration form, and updates, please check the official Winter School page.
Don’t miss this opportunity to meet leading experts, learn the latest in electronic structure theory, and connect with fellow scientists from around the world! 
Hi Chris,
The error comes from the missing END line that should close the NEWK block.
As for print option 67: it should work in properties as well as in crystal. That said, these print flags are not very well tested, and some of them may have broken in newer versions. I took a look in the code, and I can confirm that the -505 trick will not work here. Instead, you should set the value to 999 to get the maximum possible output (ie all the virtual states, or 499 for just the occupied states).
One more note: in the NEWK section this option was never extended to parallel execution, so you’ll need to run it on a single CPU core if you want it to print.
Hope this helps.
Hi Chris,
If you want to switch on printing options 66 and 67, your input should be:
NEWK
6 6
1 2
66 xxx
67 yyy
where xxx and yyy are the numbers "N" specified in the manual (see page 443, column input)
Let me know if you manage to get the correct printings
The block
TWOCOMPON
SOC
END
can be inserted anywhere in the third section of the INPUT file. It is common practice for us to place it at the top of this section, but this is not strictly required.
Hi Chris,
Not quite, the pob-DZVP and pob-TZVP basis sets are all-electron valence basis sets, thus they are not paired with relativistic pseudopotentials. This means they cannot be used directly for spin–orbit coupling calculations in CRYSTAL.
For Nickel and Cobalt I suggest you to use the COLUSC ECP paired with the basis set that you can find here for Ni (the one under "Nickel Relativistic Effective Potential and Basis set including the 3s and 3p subshells...") and here for Co (the one under "Cobalt Relativistic Effective Potential and Basis set including the 3s and 3p subshells...").
You just need to convert the valence basis set in the CRYSTAL format.
Hi Rams,
The need for relatively high values of FMIXING is not unusual, so what you are observing is expected. That said, you are correct that if FMIXING is set extremely high (close to 100%), there is a risk that the SCF procedure may reach converge before reaching the "true" ground state. This is not unique to your system, it is a general trade-off with strong damping.
To improve robustness while still guiding the calculation toward the correct solution, you can try:
TOLDEE: this tightens the SCF energy convergence criterion and will allow the SCF to continue reaching the minimum even with high FMIXING. Keep in mind, however, that this will typically require many more SCF cycles.LEVSHIFT) or electronic smearing (SMEAR) (if applicable to your system), which often reduce the need for such extreme damping.Hi Chris,
For SOC calculations in CRYSTAL, you cannot use all-electron basis sets directly (such as your cc-pVDZ-DK). SOC in CRYSTAL is implemented through the use of relativistic effective core potentials (ECPs), which already incorporate scalar-relativistic and spin-orbit effects for the core electrons. Thus, at least one pseudopotential is required in your system (usually on transition metals and heavier elements).
The INPSOC keyword can be used to manually insert an ECP, but the associated basis set must be designed to describe only the valence electrons that remain after introducing the pseudopotential.
I recommend checking the CRYSTAL basis set library, where you’ll find pseudopotentials and corresponding valence basis sets that are already formatted for SOC calculations in CRYSTAL. These can be used directly without modification and will ensure consistency between the basis and the relativistic treatment. However, we don't have a SOC basis set for evry element.
In the case you are interested in a specific element for wich we don't have the SOC basis, you can try using one of the internal ECP options (see page 88-92 of the manual) and combine it with a molecular basi set.
For example, you can find the basis set for the Stuttgart ECPs (STUTSC and STUTLC keyword in CRYSTAL) here. Or here for the COLUSC and COLULC ECPs.
Hi Chris,
You’re right that the code uses the symmetry of the crystal to reduce the number of k-points it actually computes. That’s why you don’t see all 216 k-points from the 6×6×6 Monkhorst-Pack grid: only the irreducible k-points are listed, since the energies at the symmetry-equivalent points can be obtained by applying the crystal’s symmetry operations.
That said, one important distinction: the program does not literally “reuse” the eigenvalues from one representative k-point for its equivalents, it only knows that due to symmetry, those results will be the same, so it avoids recalculating them.
For band structure plots along a high-symmetry path (e.g. Gamma-X-K-Gamma), the reduced k-point mesh from the SCF run is not sufficient. Instead, you should perform a separate band structure calculation (BAND keyword in PROPERTIES, page 309 of the user manual). In this type of calculation, you explicitly define the path through reciprocal space, and PROPERTIES computes eigenvalues along that path directly, without relying on the Monkhorst-Pack grid.
Let me know if you have any truble running the properties calculation.
Hi Rams,
I’ll look into this issue further, but in the meantime, could you please try running the MPP version with the same input, just without the DIIS keyword, and let me know if that converge?
Hi esmuigors,
When you run a geometry optimization (within the OPTGEOM block), CRYSTAL will always produce a fort.34 file (that contains the optimized structure and can be used in a subsequent calculation with the EXTERNAL keyword).
If you don’t see it after the run, a common cause is that it wasn’t copied back by the script used to manage scratch directories for the execution as the runCRY script provided with the code.
As for the error that you get, the issue is that EXTPRT should be placed outside the OPTGEOM block, as in this input.d12 example.
If you prefer, you can also extract the optimized geometry manually from the .out file. However, note that you must take both the final cell parameters (and angles, if present) and the final fractional atomic coordinates from the “FINAL OPTIMIZED GEOMETRY” section (those in the primitive cell), using only updated cell parameters without matching coordinates will lead to inconsistencies and errors like the ones you saw.
Let me know if you manage to get the optimized structure.
Hi jquertin,
As Aleks correctly pointed out, it’s better not to use DIIS (which is active by default) when using SPINLOCK.
A good way to combine the two is like this:
SPINLOCK
1
-2
THREDIIS
0.005
This way, DIIS is only activated when SPINLOCK is turned off.
However, I tried running your input, and it’s likely that this is actually a bug. Unfortunately, I don’t have a better solution at the moment other than changing the functional.
Hi Emmanuel,
For an updated guide on the installation of CRYSTAL23, please follow the instructions provided here.
Specifically, to install the executable, refer to Section 1. If you need to recompile the code from object files, see Section 3.
To test the parallel executable, you can refer either to Section 2, or to this tutorial.
Let me know if you need further help.
Hi Fabio,
A few clarifications:
i) I'm aware that setting KMP_DETERMINISTIC_REDUCTION=true can help, but in practice, it doesn't guarantee reproducibility on its own. Even building crystalOMP with stricter floating-point flags like -fp-model precise (or equivalent) doesn't always lead to fully consistent results, at least not in my experience with CRYSTAL.
ii) Yes, MPI reductions are typically deterministic, as most implementations use pairwise summation or other stable schemes. OpenMP, on the other hand, can still introduce variability due to threading and compiler optimizations.
iii) If you're seeing discrepancies around 1e-5, it's plausible that small numerical differences from OpenMP reductions are enough to drive the SCF, especially in extremely delicate cases like metallic systems treated at the HF level, toward slightly different solutions.
At this point, I’ve investigated what I could from our end. If maximum reproducibility is essential, I strongly recommend sticking with MPI-only parallelism.
Have a great day,
Giacomo
Hi Fabio,
After an in-depth investigation, I can confirm that the behavior you're observing is definitely not due to uninitialized variables. The most likely explanation is the use of non-deterministic summation in reduction operations, which is a known behavior when running floating-point operations with OpenMP parallelization.
In OpenMP, when large arrays or matrices are reduced (e.g., summing elements of the density or Fock matrices), the order of operations can change depending on how the threads are scheduled at runtime. Since floating-point addition is not associative, this can lead to slightly different results depending on the number of threads, their distribution, or the memory layout at runtime, even if the input is exactly the same.
To confirm this, I ran your script using the public OpenMP build, with both 16 threads and a single thread. Results are reported below.
As expected, the more threads you use, the larger these differences may become, although they typically remain very small.
Increasing the number of threads likely increases the occurence for these small numerical differences, although they typically remain insignificant in practice. For this reason, we recommend limiting the number of OpenMP threads to 4, and primarily use MPI for more extended parallelism, which generally offers a better balance between performance and numerical stability.
Additionally, it's good practice in our workflow to run each job in a clean, isolated SCRATCH folder, to avoid any unintended interference from leftover files or cached data. Your current script doesn’t do this, so we recommend updating it to clean or define a separate working directory for each run.
Energies at each SCF Cycle - Copper bulk
-----------------------------------------------------------------------------
OpenMP 16 th
-----------------------------------------------------------------------------
run 1 run 2 diff
CYC 0 -1.959786174741E+02 -1.959786174741E+02 0.0
CYC 1 -1.959830545289E+02 -1.959830545289E+02 0.0
CYC 2 -1.959842257372E+02 -1.959842257372E+02 0.0
CYC 3 -1.959843727892E+02 -1.959843727892E+02 0.0
CYC 4 -1.959844091796E+02 -1.959844091796E+02 0.0
CYC 5 -1.959844090374E+02 -1.959844090374E+02 0.0
CYC 6 -1.959844090603E+02 -1.959844077969E+02 -1.2634000086109154e-06
CYC 7 -1.959844090705E+02 -1.959844090714E+02 9.000018508231733e-10
CYC 8 -1.959844090772E+02 6.700020094285719e-09
CYC 9 -1.959844090706E+02 1.000159954855917e-10
-----------------------------------------------------------------------------
Open MP 1 th
-----------------------------------------------------------------------------
run 1 run 2 diff
CYC 0 -1.959786174741E+02 -1.959786174741E+02 0.0
CYC 1 -1.959830545289E+02 -1.959830545289E+02 0.0
CYC 2 -1.959842257372E+02 -1.959842257372E+02 0.0
CYC 3 -1.959843727892E+02 -1.959843727892E+02 0.0
CYC 4 -1.959844091796E+02 -1.959844091796E+02 0.0
CYC 5 -1.959844090374E+02 -1.959844090374E+02 0.0
CYC 6 -1.959844090603E+02 -1.959844090603E+02 0.0
CYC 7 -1.959844090705E+02 -1.959844090705E+02 0.0
Hi Fabio,
Regarding your initial question, I ran the public OpenMP version of the code (using 4 threads) multiple times with your input, and I consistently obtain the same results across runs. The SCF procedure converges smoothly in 6 cycles, as shown below:
$ grep DETOT input_4th.out
CYC 0 ETOT(AU) -1.959786174741E+02 DETOT -1.96E+02 tst 0.00E+00 PX 1.00E+00
CYC 1 ETOT(AU) -1.959830545289E+02 DETOT -4.44E-03 tst 0.00E+00 PX 1.00E+00
CYC 2 ETOT(AU) -1.959842257372E+02 DETOT -1.17E-03 tst 5.19E-04 PX 2.32E-02
CYC 3 ETOT(AU) -1.959843829641E+02 DETOT -1.57E-04 tst 6.44E-04 PX 5.14E-02
CYC 4 ETOT(AU) -1.959844090345E+02 DETOT -2.61E-05 tst 4.59E-06 PX 4.09E-03
CYC 5 ETOT(AU) -1.959844090249E+02 DETOT 9.60E-09 tst 1.92E-07 PX 1.01E-03
CYC 6 ETOT(AU) -1.959844090704E+02 DETOT -4.56E-08 tst 8.72E-09 PX 1.01E-03
$ grep DETOT input_4th_2.out
CYC 0 ETOT(AU) -1.959786174741E+02 DETOT -1.96E+02 tst 0.00E+00 PX 1.00E+00
CYC 1 ETOT(AU) -1.959830545289E+02 DETOT -4.44E-03 tst 0.00E+00 PX 1.00E+00
CYC 2 ETOT(AU) -1.959842257372E+02 DETOT -1.17E-03 tst 5.19E-04 PX 2.32E-02
CYC 3 ETOT(AU) -1.959843829641E+02 DETOT -1.57E-04 tst 6.44E-04 PX 5.14E-02
CYC 4 ETOT(AU) -1.959844090345E+02 DETOT -2.61E-05 tst 4.59E-06 PX 4.09E-03
CYC 5 ETOT(AU) -1.959844090249E+02 DETOT 9.60E-09 tst 1.92E-07 PX 1.01E-03
CYC 6 ETOT(AU) -1.959844090704E+02 DETOT -4.56E-08 tst 8.72E-09 PX 1.01E-03
$ grep DETOT input_4th_3.out
CYC 0 ETOT(AU) -1.959786174741E+02 DETOT -1.96E+02 tst 0.00E+00 PX 1.00E+00
CYC 1 ETOT(AU) -1.959830545289E+02 DETOT -4.44E-03 tst 0.00E+00 PX 1.00E+00
CYC 2 ETOT(AU) -1.959842257372E+02 DETOT -1.17E-03 tst 5.19E-04 PX 2.32E-02
CYC 3 ETOT(AU) -1.959843829641E+02 DETOT -1.57E-04 tst 6.44E-04 PX 5.14E-02
CYC 4 ETOT(AU) -1.959844090345E+02 DETOT -2.61E-05 tst 4.59E-06 PX 4.09E-03
CYC 5 ETOT(AU) -1.959844090249E+02 DETOT 9.60E-09 tst 1.92E-07 PX 1.01E-03
CYC 6 ETOT(AU) -1.959844090704E+02 DETOT -4.56E-08 tst 8.72E-09 PX 1.01E-03
$ grep DETOT input_4th_4.out
CYC 0 ETOT(AU) -1.959786174741E+02 DETOT -1.96E+02 tst 0.00E+00 PX 1.00E+00
CYC 1 ETOT(AU) -1.959830545289E+02 DETOT -4.44E-03 tst 0.00E+00 PX 1.00E+00
CYC 2 ETOT(AU) -1.959842257372E+02 DETOT -1.17E-03 tst 5.19E-04 PX 2.32E-02
CYC 3 ETOT(AU) -1.959843829641E+02 DETOT -1.57E-04 tst 6.44E-04 PX 5.14E-02
CYC 4 ETOT(AU) -1.959844090345E+02 DETOT -2.61E-05 tst 4.59E-06 PX 4.09E-03
CYC 5 ETOT(AU) -1.959844090249E+02 DETOT 9.60E-09 tst 1.92E-07 PX 1.01E-03
CYC 6 ETOT(AU) -1.959844090704E+02 DETOT -4.56E-08 tst 8.72E-09 PX 1.01E-03
I suspect that the behavior you're observing might be related to differences in the OpenMP and/or MKL library versions used in your environment. In my experience, such differences can influence numerical results to some extent (especially in delicate systems like metals treated with HF, where convergence can be sensitive to small perturbations or numerical noise).