Posts

• ---

  Unexpectedly cannot get the geometry after optimization: KEYWORD EXTPRT NOT ALLOWED
  undefined GiacomoAmbrogio

  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.

• ---

  SCF fails spinlock with POB-DZVP-REV2
  undefined GiacomoAmbrogio

  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.

• ---

  Which keywords can help the geometry optimization to what you expected ?
  undefined aerba

  Hi,

  The FINALRUN option works fine in P-CRYSTAL: that sentence in the Manual must be a leftover from the past, sorry for that.

  Can you elaborate on how the optimized structures are unreasonable? Is it just for the underestimation of the computed band gap or also for some structural aspects? In this case, it would help if you could show us the "expected" structure, as well as the one you obtain from the optimization.

• ---

  Negative value in the SCF of an anion vacancy
  undefined aerba

  Hi,

  Jefferson Maul has run a few tests and indeed we were able to reproduce the odd behavior you've been experiencing. Those negative values for the electron population of the vacancy are indeed very strange and may be suggestive of a basis set unbalance. It seems that whenever basis functions are left on the vacant site (either with your initial ATOMSUBS approach or with the more "canonical" GHOSTS approach) those functions are much "needed" by surrounding atoms.

  We have found a possible way forward for your system through the use of the ATOMREMO option. With this option, you create the vacancy by removing the selected atom alongside its basis functions. For instance, we have used this approach to study Oxygen vacancies in CaSnO3 here.

  We have used this approach on your system (on a smaller supercell to be able to run the tests more efficiently) and the geometry optimization went well. See the attached output file as a guide.

  Hope this helps,

• ---

  Question about basic principle of pressure and frequency calculation
  undefined aerba

  Hi,

  I understand your confusion, which originates from the fact that PRESSURE is an old option of FREQCALC that indeed does not apply any pressure. With this option, the structure is unchanged, the forces are unchanged, and ultimately, the harmonic vibration frequencies are unchanged. The only bit of information that is affected by this keyword is the printed value of the PV term in the thermodynamic analysis.

  So, this PRESSURE keyword within FREQCALC should be used in just one way: to set the value of pressure corresponding to the structure used to run the harmonic frequency calculation, so as to have the PV term entering the thermodynamic functions (enthalpy and Gibbs free energy) right. Let me make a couple of examples:

  1. You perform a full structural relaxation (atomic positions + lattice parameters). Thus, you have a zero pressure (p=0) structure. Then you compute the harmonic frequencies. If you are interested in thermodynamic potentials (i.e. enthalpy, Gibbs, Helmholtz, etc.) you should use the PRESSURE keyword and set the pressure to zero to get the PV term right in the output (otherwise by default it would be computed at p = 1 atm - do not ask me why).

  2. You perform a pressure-constrained geometry optimization with the EXTPRESS keyword of OPTGEOM (let's say at 29.457 GPa as in your example). Then you compute harmonic frequencies with FREQCALC on this compressed structure. Now, again, if you are interested in thermodynamic potentials (i.e. enthalpy, Gibbs, Helmholtz, etc.) you should use the PRESSURE keyword and set the pressure to 29.457 GPa to get the PV term right in the output.

  So, as I said at the beginning, the PRESSURE keyword here does not change anything on the structure/forces. It only allows you to tell the program what is the pressure of the structure you are working on, so that the printed PV term is right.

  A more formally-sound way to combine temperature and pressure is provided by the quasi-harmonic approximation: See Eq. (13) of the tutorial page on the Quasi-Harmonic Approximation (QHA).

  Hope this helps clarifying things a little,

• ---

  CRYSTAL23 Installation Linux one node
  undefined GiacomoAmbrogio

  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.

• ---

  OpenMP problem
  undefined GiacomoAmbrogio

  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

• ---

  OpenMP problem
  undefined GiacomoAmbrogio

  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
  
  
  

• ---

  Generatinng INPUT file for difference map of electron charge density and electrostatic potential plot
  undefined aerba

  Hi,

  When computing (and then plotting) 2D maps of the electron density, or spin density (by use of the ECHG keyword), or electrostatic potential (by use of the POTM keyword), a plane must be specified. Now, in CRYSTAL a plane is specified by three points A, B, and C (i.e. if you specify the coordinates of three -- non rectilinear -- points, you identify a unique plane).

  In CRYSTAL, you can specify the three points A, B and C in two ways:

  1. You can directly input the Cartesian coordinates of the three points by use of the COORDINA keyword;

  2. You can select three atoms as your three points, with the ATOMS keyword. In this case, as the system is periodic and you may need to refer to atoms outside of the reference cell, each atom is identified by four integer numbers (the sequential index of the atom within the reference cell and three indices identifying the cell where your selected atom is centered -- i.e. 4 0 0 1 for instance means atom number 4 in the list, found in cell 0 0 1, that is a cell translated by one lattice vector along the c crystallographic axis).

  Now, if the three points A, B and C that you have defined are not "orthogonal", i.e. the angle \( \widehat{ABC} \neq 90^\circ\) this would produce an oblique map. To make it rectangular, you can insert the RECTANGU keyword, as shown below:

  

  You can also change the explored range in the selected plane (relative to that identified by the three original points A, B and C) by use of the keyword MARGINS, as shown below:

  

  Input examples can be found at this tutorial page.

  Hope this helps,

• ---

  OpenMP problem
  undefined GiacomoAmbrogio

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