Skip to content
  • Home
  • Recent
Collapse
Brand Logo
CRYSTAL23
Latest v1.0.1
Tutorials Try the Demo Get a License
Tutorials Try the Demo Get a License Instagram
GiacomoAmbrogio

Giacomo Ambrogio

@GiacomoAmbrogio
Developer
Unfollow Follow
About
Posts
74
Topics
5
Shares
0
Groups
1
Followers
8
Following
3

Posts

Recent Best Controversial

  • segmentation fault with double free or corruption during SCF
    GiacomoAmbrogio GiacomoAmbrogio
    Bug Reports

    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)


  • Input MOF geometry problem
    GiacomoAmbrogio GiacomoAmbrogio
    Geometry Editing

    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.


  • PCrystal job stuck when run between several nodes
    GiacomoAmbrogio GiacomoAmbrogio
    Running CRYSTAL in Parallel

    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.


  • OpenMP problem
    GiacomoAmbrogio GiacomoAmbrogio
    Running CRYSTAL in Parallel

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


  • OpenMP problem
    GiacomoAmbrogio GiacomoAmbrogio
    Running CRYSTAL in Parallel

    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
    
    
    

  • OpenMP problem
    GiacomoAmbrogio GiacomoAmbrogio
    Running CRYSTAL in Parallel

    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


  • How can I find the total number of irreducible representations corresponding to the k-points?
    GiacomoAmbrogio GiacomoAmbrogio
    Geometry Optimisations

    Hi Yachao_su,

    The relation you mention between MPI processes and efficiency is correct, but a bit of clarification might be necessary:limitations on efficiency apply only to the linear algebra steps and the density matrix construction (specifically, the FDIK and PDIG steps) during the SCF procedure.

    In the FDIK section we are essentially diagoalizing a set of Fock matrices (of size n atomic orbitals x n atomic orbitals), one for each k point. When using symmetry-adapted Bloch functions, these matrices become block-diagonal, and each block corresponds to an irreducible representation (irrep). These blocks can then be diagonalized independently (by different MPI processes).

    Diagonalization scales cubically with matrix size. However, for small systems, this time is usually negligible compared to other SCF steps. So in such cases, there's no need to worry too much about how many MPI processes you're using. 😊

    The other most time-consuming part is the computation of two-electron integrals (SHELLXN step), which scales almost linearly with system size and is not affected by the MPI-per-irrep limitation.

    The limitation becomes relevant when dealing with large systems, where the irreps (ie, the block sizes) are large. In that case, the time spent on FDIK and PDIG becomes significant. Since Pcrystal cannot assign more than one MPI process per irrep, any additional processes will stay idle during diagonalization, but they are still used in other parts of the calculation.

    68716be2-c92c-4dbf-a2db-f8f98e7bc3e6-image.png

    Here, I've reported an example. You can see the blue line (FDIK, i.e., diagonalization) is negligible on the left (small systems) but becomes dominant on the right (larger systems). So the MPI limitation becomes noticeable only in that region.

    If you're on the right side of the plot, where the time spent in FDIK exceeds that of other steps like SHELLXN or NUMDFT (DFT integration), you should limit your number of MPI processes using the rule:

    $$ \text{max MPI processes} = \sum_\mathbf{k} n_{irreps, \mathbf{k}} \times n_{spin}$$

    That is, the total number of processes should not exceed the sum of all irreps across all k points, multiplied by the number of spins (2 if you're doing a spin-polarized calculation, 1 for close shell).

    By default, the number of irreps is not printed in the output file. However, you can activate some advanced printing options to display this information. To do so, insert the following into the third block of your input file:

    KSYMMPRT
    SETPRINT
    1
    47 nk
    

    Here, nk is the maximum number of k points for which you want to print symmetry information. Likely this will be the number of k points in your calculation.

    You can find some reference about this option in the CRYSTAL User Manual, at pages 117-118.

    The output will look like this:

     +++ SYMMETRY ADAPTION OF THE BLOCH FUNCTIONS +++
    
     SYMMETRY INFORMATION:
     K-LITTLE GROUP: CLASS TABLE, CHARACTER TABLE.
     IRREP-(DIMENSION, NO. IRREDUCIBLE SETS)
     (P, D, RP, RD, STAND FOR PAIRING, DOUBLING, REAL PAIRING AND REAL DOUBLING
     OF THE IRREPS (SEE MANUAL))
    
     K[   1] (  0  0  0)                                                    <------ 1st k point
    
     CLASS | GROUP OPERATORS (COMPLEX PHASE IN DEGREES)
     ------------------------------------------------------------------------
       2   |   2(   0.0);
       3   |   3(   0.0);   4(   0.0);
       4   |   5(   0.0);   6(   0.0);
       5   |   7(   0.0);   8(   0.0);   9(   0.0);
       6   |  10(   0.0);  12(   0.0);  11(   0.0);
    
     IRREP/CLA     1     2     3     4     5     6
     ---------------------------------------------
      MULTIP |     1     1     2     2     3     3
     ---------------------------------------------
          1  |  1.00  1.00  1.00  1.00  1.00  1.00
          2  |  1.00 -1.00  1.00 -1.00 -1.00  1.00
          3  |  2.00  2.00 -1.00 -1.00  0.00  0.00
          4  |  1.00  1.00  1.00  1.00 -1.00 -1.00
          5  |  2.00 -2.00 -1.00  1.00  0.00  0.00
          6  |  1.00 -1.00  1.00 -1.00  1.00 -1.00
    
       1-(1,  27);   2-(1,  25);   3-(2,  37);   4-(1,  23);   5-(2,  37);  <------ This is the information about
       6-(1,  25);                                                          <------ the irreps
    
     K[   2] (  1  0  0)                                                    <------ 2nd k point.
    
     CLASS | GROUP OPERATORS (COMPLEX PHASE IN DEGREES)
     ------------------------------------------------------------------------
       2   |  11(   0.0);
    
     IRREP/CLA     1     2
     ---------------------
      MULTIP |     1     1
     ---------------------
          1  |  1.00  1.00
          2  |  1.00 -1.00
    
       1-(1, 126);   2-(1, 122);                                             <------ irreps of 2nd k point
    ....
    

    For example, here the first k point is adapted in 6 different irreps, while the second in 2

    Please note: this feature should be run using one single MPI process.


  • How use CRYSTALpytools?
    GiacomoAmbrogio GiacomoAmbrogio
    Band Structure

    Hi,
    I suggest you to use CRYSTALClear instead. Here you can find an installation guide, and on the GitHub page there is the documentation and some example notebooks that you can use as template.
    For any question related to CRYSTALClear there is a dedicated section here on the forum.


  • Compiliation of CRYSTAL 17 on Apple Silicon
    GiacomoAmbrogio GiacomoAmbrogio
    Installation and Configuration

    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.


  • running in parallel on OSX
    GiacomoAmbrogio GiacomoAmbrogio
    Running CRYSTAL in Parallel

    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.


  • Conversion of fort.25 file to .cube (Vesta file)
    GiacomoAmbrogio GiacomoAmbrogio
    Electron Charge Density Analysis

    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!


  • Unexpectedly cannot get the geometry after optimization: KEYWORD EXTPRT NOT ALLOWED
    GiacomoAmbrogio GiacomoAmbrogio
    Geometry Optimisations

    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.


  • Is it necessary to add the keyword NOSYMADA when using a user-defined anisotropic k-point grid?
    GiacomoAmbrogio GiacomoAmbrogio
    Single-Point Calculations

    We ran some tests regarding the inclusion of the keyword NOSYMADA in relation to anisotropic shrinking factors, and I can confirm that we did not find any significant difference. Therefore, anisotropic shrinking is fully compatible with symmetry-adapted block functions.

    Taking a closer look to the data and input/output files you submitted, the convergence with respect to the shrinking factors appears to be very similar. Comparisons of calculations using the same SHRINK value consistently show differences well within the default tolerance (10^-6 Ha/cell), except for one value (highlighted in red in the screenshot). In this case, the energy seems to be significantly off the trend. Could you please double-check that specific calculation?

    Screenshot 2025-04-16 112322.png


  • Input MOF geometry problem
    GiacomoAmbrogio GiacomoAmbrogio
    Geometry Editing

    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.

    image.png

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


  • R2SCAN not implemented with UHF or incompatible with EOS?
    GiacomoAmbrogio GiacomoAmbrogio
    Other Questions

    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
    [...]
    

  • Extracting eigenvalues for a certain high-symmetry path
    GiacomoAmbrogio GiacomoAmbrogio
    Single-Point Calculations

    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.


  • SMALLDIST keyword
    GiacomoAmbrogio GiacomoAmbrogio
    Single-Point Calculations

    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.

  • 1 / 1
  • Login

  • Don't have an account? Register

  • Login or register to search.
  • First post
    Last post
0
  • Home
  • Recent