The CRYSTAL code will have a dedicated full day at the upcoming Graduate School on Electronic Structure Theory at STFC Daresbury Laboratory, taking place 24–28 February 2025!

This event is a great opportunity to:

• Dive deep into the fundamentals of electronic structure theory.
• Gain hands-on experience with CRYSTAL, alongside other leading software packages like CASTEP and QUESTAAL.
• Learn both command-line and ASE (python based) interfaces.
• Connect with peers and experts in the field.

Whether you're a doctoral researcher, postdoc, or someone looking to enhance your expertise in electronic structure codes, the dedicated CRYSTAL day offers a valuable opportunity to develop a solid understanding of the code's features, gain practical experience, and understand how it fits into the broader landscape of computational tools.

Event Details
Dates: 24-28 February 2025
Application deadline: 19 January 2025.
Venue: Science and Technology Facilities Council, Daresbury Laboratory, UK.
Apply Now: https://ccp9.ac.uk/graduate_school_2025/

Take advantage of this opportunity to expand your knowledge and experience with the CRYSTAL code in a collaborative and dynamic environment. We look forward to seeing you there!

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 Antonio,
Unfortunately, I don't think there is any specific print statement in the .out file that indicates if the calculation was run using either SDFT or SCDFT formalism. Maybe in newer versions of the code, we will add more information about this!

I guess the only way to tell is by looking at the input file. A good practice could be to print the input file at the top of the output in your launch script, so that you always have a reference to your input for each output you generate.

In any case, you can take a look at this link. Here, you will find instructions on how to install the parallel version of CRYSTAL, both from executable and object files (see paragraphs 2 and 3). The instructions refer to the latest version of the code, which is CRYSTAL23, but I think the procedure is similar for CRYSTAL17.

Another useful reference is the tutorial on "How to run", see overview and parallel run sections.

Hope this helps!

Dear CRYSTAL community,

We’re excited to share our recent work on accelerating linear algebra operations in the CRYSTAL code using GPUs. Our implementation boosts the performance of self-consistent field (SCF) calculations by offloading key matrix operations like multiplication, diagonalization, inversion, and Cholesky decomposition to GPUs.

In the manuscript, we first analyze the performance and limitations of the standard parallel version of the code (Pcrystal) and then we evaluate the scalability of the new GPU-accelerated approach with 1 to 8 GPUs, observing remarkable scaling. To highlight these improvements, we present benchmark results on different systems, such as the example below.

We expected significant speedups for large systems due to the limited number of k points, each requiring substantial computational effort. To ensure a fair comparison, we ran calculations using the massively parallel version of CRYSTAL (MPPcrystal) on a large MOF structure with over 30000 basis functions. Surprisingly, a single GPU on one node performed comparably to 512–1024 CPU cores running across 4–8 nodes.

To find out more, read the full paper here.

We aim to make this GPU-accelerated version of CRYSTAL available in the upcoming release, allowing all users to benefit from its enhanced performance for large-scale simulations. We look forward to reading your thoughts and discussing potential applications or further improvements.

A big thanks to Lorenzo Donà, Chiara Ribaldone, and Filippo Spiga for their contributions to the development of this code!

Hi othmen1983,
What you need to do next depends on whether you have a Pcrystal executable or the OBJ files to compile.

If you have the executable, you can check if it runs correctly on your machine by downloading the input file INPUT and placing it in an empty folder (ensure the file is named INPUT, without any extension).
After that you can run this command in the same folder:

mpirun -np 2 path/to/crystal/exe/Pcrystal

If you see the output correctly on your screen, it means the executable is working as intended. The final step is to create launch scripts to manage input, output, and temporary files.

Example scripts for CRYSTAL23 are available here (they should also work with CRYSTAL17). A brief explanation on how to use them can be found in the How to run tutorial.

Unfortunately, these scripts only work on simple machines, such as workstations. If you need to run Pcrystal on a cluster with a queuing system like Slurm, you will need dedicated scripts.

Hi job314,
After looking into your issue further, I’m following up with more information.

The B973C functional is a composite method with built-in corrections specifically designed for the mTZVP basis set. Modifying the basis set can introduce errors and is not the right approach. This method and basis set were primarily developed for molecular systems and, at most, molecular crystals, not bulk materials like yours.

Explicit warnings about this functional can be found in the user manual on page 161.

Given this, I recommend choosing a different functional and basis set better suited for your system.

Hi Danny,

From a first look at your input file, I noticed that you are performing a combined PDOS and BANDS calculation. This has not been fully tested as we prefer to run two separate jobs, one for PDOS and one for BANDS. Once the harmonic frequencies are computed, they can both be run through restarts and are both very fast.

Indeed, the PDOS calculation produces a single set of data if run without the BANDS keyword.

I therefore suggest to split your calculation in two different ones: one for PDOS and one for BANDS. You can use the RESTART keyword in the FREQCALC block to avoid recomputing the Hessian matrix (see page 219 of the User manual).

I leave here a link to a tutorial webpage that we have recently updated about computing phonon-related properties with CRYSTAL.

Let me know if this helps!

Hi piquini,
You could try to add this to your mpirun line:

mpirun -np 4 $CRY23_EXEDIR/$VERSION/Pcrystal |& tee $PBS_O_WORKDIR/${PBS_JOBNAME}.out

Let me know if it works!

Hi heimurinn,

The error you're seeing is not actually related to the symmetry adaptation of the Bloch functions. What’s happening is that the error is triggered by another processor before the one responsible for writing the output reaches the actual failure point.

If you run the same calculation on one single process (serial mode), you should see the real error, which (for 81_from1.d12, I didn’t try 81_to1403845.d12, but I suspect it will be the same) is:
ERROR **** GROTA1 **** ERROR IN SYMMETRY EQUIVALENCE - CHECK INPUT COORDINATES

This issue is likely related to symmetry. One thing you can try is to follow the structure standardization procedure described in this thread.

If you have any questions or need further help, feel free to ask!

Hi Job,
It seems there might be a bit of confusion regarding how to report MPI bindings, maybe we will update the tutorial to make it more clear.

To display binding information correctly, you’ll need to include the appropriate flag in your usual mpirun command. For example:

• If you're using OpenMPI:

mpirun --report-bindings -np 192 Pcrystal

• If you're using Intel MPI:

mpirun -print-rank-map -np 192 Pcrystal

To check which MPI implementation you're using, you can inspect the loaded modules in your environment, or simply run:

mpirun --version

If you're using a different MPI implementation, feel free to let me know. I'd be happy to help you find the right way to print the bindings.

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

Hi R.Zosiamliana,

When using SLABCUT from a 3D system, or directly a SLAB input, CRYSTAL treats the resulting system as a lower-dimensional 2D structure. Specifically, the system retains periodicity only along two directions (conventionally in the x-y plane), while periodicity is inherently absent along the orthogonal direction (z).

Just to preserve the same output format as for a 3D calculation, an arbitrary value of 500 Å is printed in the output for the c lattice vector, even if it does not exist for such calculations. We understand this may be confusing and we are considering to remove it in future versions.

To summarise, in CRYSTAL when running a 2D (or 1D, or 0D) calculation, there is no need to repeat the structure along the non-periodic direction and thus there is no need to define a vacuum.

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.

• In the parallel run, there are slight differences in the final energy, but the magnitude of the difference is negligible, well within the SCF convergence tolerance.
• In the single-thread run (i.e., serial execution with OpenMP), the results are fully reproducible, with no differences at all.

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 ywang,
If you are running the parallel version of CRYSTAL, could you try run it on one single processor?
Sometimes, when running in parallel, an error can terminate the job before the output is printed, especially during the input reading section.

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 job314,

I tried to run your input (without the RESTART keyword), and it seems to work fine.
Can you double-check the file used for the restart? Or eventually, can you try run the code without RESTART?

Anyway, the proper way to perform a spin-polarized calculation in DFT is to use the SPIN keyword in the DFT block instead of UHF (but technically, both should work):

[...]
DFT
SPIN
R2SCAN
XLGRID
ENDdft
[...]

Hi ywang,

We managed to resolve the problem with your input file. The issue was that the system was not standardized according to the space group.

It should be possible to tell the program that the structure is not standard through a specific keyword. However, an alternative approach is to manually modify the CIF file. This can be done using VESTA, please refer to the image.

The new input geometry look like this:

Fe_complex
CRYSTAL
0 0 0
225
26.7506
6
8  0.05353    0.18260    0.24189
26 0.50000    0.21143    0.21143
6  0.20297    0.20297    0.07018
6  0.17805    0.17805    0.11438
6  0.13562    0.13562    0.19929
1  0.12115    0.12115    0.22832
ENDgeom

Fe.d12

This has also reduced the number of symmetry-irreducible atoms, so the calculation should be a bit faster.

Hi Chhatra,
This is because your system is no longer in 3D after you do a SLABCUT. In order to run a ECH3 properties calculation on a slab, you shoud add a RANGE for the non-periodic direction.

ECH3
100
RANGE
-2
+2
END

The two numbers after RANGE are the min and max distance from the 0 of the cell along the non periodic direction in atomic units.
The same spacing from the priodic dicrections is used to generate the points.

You can find more referneces on this at page 322 of the user manual.

Hope this helps!