Thank you!

Hi Alessandro,

Thank you for your reply. I will use an isotropic shrink. I find it strange, though, to have the warning message as the structure was optimised. Anyway, I will rerun that.

Kind regards,
Pierre

Either way, seems you got your geometry for further optimization 🙂
Cheers,
A

Hi,

There may not be an option to "plot" orbital-resolved bands directly, but there is an option to compute and print principal atomic orbitals contributions to selected eigenvectors. See the ANBD keyword of the PROPERTIES module (page 307 of the CRYSTAL23 User's Manual).

Yes, if your system is 2D, you can use the usual BAND keyword of the PROPERTIES module to define a 2D path to compute and plot the band structure (page 309 of the CRYSTAL23 User's Manual).

If you need help on a specific system, just let me know

Thanks Giacomo, much appreciated,

Chris

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:

Increasing 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. Combining mixing with other stabilization strategies, such as level shifting (LEVSHIFT) or electronic smearing (SMEAR) (if applicable to your system), which often reduce the need for such extreme damping.

Hi Giacomo,

That's a very useful answer, thank you!

Chris

Dear esmuigors
esmuigors said in Unexpectedly cannot get the geometry after optimization: KEYWORD EXTPRT NOT ALLOWED:

Should I only take the atoms labelled by "T" if the lattice is tetragonal, and which ones are to be taken in case of an orthorombic one? Or does "T"'and "F" stand just for "TRUE"'and "FALSE"?

You are correct T and F stand for 'TRUE' and 'FALSE'.

esmuigors said in Unexpectedly cannot get the geometry after optimization: KEYWORD EXTPRT NOT ALLOWED:

I am actually talking about the standard runPcry23 script. Perhaps it should not be like that? Or should I modify this script?

That script already contain the few line of codes I suggested you, if not feel free to modify it and add those lines.

esmuigors said in Unexpectedly cannot get the geometry after optimization: KEYWORD EXTPRT NOT ALLOWED:

Actually, I have now tried using only the "T"-labeled atoms and got this error:

ERROR **** geometry **** FORMAT ERROR IN INPUT DECK

I checked the input on top of the CS2_B1WC.pob_tzvp_rev2_gamma_onlyT.out file it seems like you are defining six atoms in the asymmetric unit, but only two atomic positions are specified in the input.

CRYSTAL 0 0 0 64 6.34350113 5.54788046 9.71085545 6. ! Definitions of the number of atoms in the asymmetric unit 6 0.000000000000e+00 0.000000000000e+00 -5.000000000000e-01 16 2.303037287581e-17 3.119080311720e-01 1.207617820468e-01 ATOMSYMM ...

The number of atoms in the asymmetric unit and the position specified should always match.

I hope this helps you,

Best,

eascrizzi Thank you for your help.

Sure. I will try it. Thankyou so much.

Thanks for your response, aerba.
ATOMREMO works but when you calculate the defect formation energies you end with values (for anion vacancies) that are higher than with using GHOSTS or the ATOMSUBS, with the latter values agreeing more with published defect formation energies.

aerba I've noticed for my systems, these SCF convergence issues appear when I try to compute for intensities in a restart job. Then when removing intensities, they converge fine and progress forward.

Is it possible to finish the frequency calculations first and then compute for the Raman intensities afterwards? Would that be a work around? Maybe with RAMANREA - because I have the TENS_RAMAN.DAT, but it just gets stuck at some point during HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
FORCE CONSTANT MATRIX - NUMERICAL ESTIMATE
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
and doesn't progress with restarts

Hi Giu, it worked!

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.

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:

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

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,

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.

Dear Rams,

From your input/output files, it doesn't seem you are doing anything particularly wrong. It is more that the system is perhaps a bit more challenging. Here are a few ideas that might be of use regarding speed up and convergence:

i) You are having quite large lattice parameters (15 Angstrom in each direction). In this case having a 4x4x4 mesh of k-points is perhaps even more than needed to converge the total energy, so you can try reducing it to a 3x3x3 or 2x2x2 mesh to reduce the computational cost and get some speed up (without compromising quality). You could probably go to 1 k-point just as well, but then you might not get a correct Fermi level sampled. Which brings me to the second point...

ii) Do you expect your system to be metallic? If so, you can test a few values of the smearing parameter to see whether it makes a difference to the convergence behaviour. If you know you need to get a semiconductor, then you might want to try shifting the occupied states down in energy (keyword LEVSHIFT to prevent the system from getting stuck in a conducting state and causing the SCF to get stuck...

iii) You have the DIISALLK keyword in your input file, have you tested that as necessary for the convergence of the undoped cell? It does increase the memory/storage requirements, so could be worth trying with just the default (gamma-point) DIIS...

iv) One alternative is to turn off the DIIS convergence accelerator and revert to "classical" Fock mixing. Albeit slower in convergence than DIIS, it can stabilize the SCF procedure avoiding strong, sudden oscillation in the solution (while DIIS is turned on, the FMIXING is used only in the very first SCF step)...

v) Coming back to the charge oscillations you observed, they do indeed start appearing very quickly in your second run.
Here are the corresponding total energies from your OUTPUT_Run1.d12:

CYC 0 ETOT(AU) -1.128999989618E+05 DETOT -1.13E+05 tst 0.00E+00 PX 1.00E+00 CYC 1 ETOT(AU) -1.129553711002E+05 DETOT -5.54E+01 tst 0.00E+00 PX 1.00E+00 CYC 2 ETOT(AU) -1.129699205883E+05 DETOT -1.45E+01 tst 7.60E-05 PX 9.78E-02 CYC 3 ETOT(AU) -1.123277877060E+05 DETOT 6.42E+02 tst 4.05E+00 PX 8.33E+00 CYC 4 ETOT(AU) -1.123934069238E+05 DETOT -6.56E+01 tst 1.68E-01 PX 6.66E+00 CYC 5 ETOT(AU) -1.128505608381E+05 DETOT -4.57E+02 tst 2.53E-01 PX 1.17E+00 CYC 6 ETOT(AU) -1.126570396187E+05 DETOT 1.94E+02 tst 1.41E-01 PX 1.13E+00 CYC 7 ETOT(AU) -1.129397986733E+05 DETOT -2.83E+02 tst 1.36E-01 PX 1.16E+00 CYC 8 ETOT(AU) -1.128184205493E+05 DETOT 1.21E+02 tst 9.81E-02 PX 6.95E-01 CYC 9 ETOT(AU) -1.129501102455E+05 DETOT -1.32E+02 tst 9.49E-02 PX 6.74E-01 CYC 10 ETOT(AU) -1.129223250516E+05 DETOT 2.78E+01 tst 5.80E-01 PX 6.62E-01

And then from OUTPUT_Run2.d12

CYC 0 ETOT(AU) -1.129223250516E+05 DETOT -1.13E+05 tst 0.00E+00 PX 1.00E+00 CYC 1 ETOT(AU) -8.145313982597E+04 DETOT 3.15E+04 tst 0.00E+00 PX 1.71E+02 CYC 2 ETOT(AU) -1.058417907360E+05 DETOT -2.44E+04 tst 2.27E+00 PX 1.71E+02 CYC 3 ETOT(AU) -7.201760490303E+04 DETOT 3.38E+04 tst 9.37E+00 PX 3.62E+02 CYC 4 ETOT(AU) -3.830339502397E+04 DETOT 3.37E+04 tst 1.36E+02 PX 7.83E+02 CYC 5 ETOT(AU) -5.612117026766E+04 DETOT -1.78E+04 tst 2.27E+01 PX 5.54E+02 CYC 6 ETOT(AU) -1.006357459536E+05 DETOT -4.45E+04 tst 3.23E+01 PX 5.67E+02 CYC 7 ETOT(AU) -1.131945182713E+05 DETOT -1.26E+04 tst 9.68E+00 PX 2.13E+02 CYC 8 ETOT(AU) -1.203304661632E+05 DETOT -7.14E+03 tst 8.51E+00 PX 1.06E+02 CYC 9 ETOT(AU) -1.148704842219E+05 DETOT 5.46E+03 tst 1.77E+01 PX 6.53E+01 CYC 10 ETOT(AU) -1.235795679474E+05 DETOT -8.71E+03 tst 2.63E+00 PX 4.40E+01

In the first run the total energy does not oscillate, hence why the atomic charges are still somewhat as expected. In the second run, the SCF solutions starts fluctuating much more (see difference in total energies) and as a result the charges get destabilized. The SCF might not recover from that at all (even though it seems it could), so I would suggest first trying speeding up the simulation to try and get the job done in one run. Alternatively you might have to go with slower mixing and more restarts, but that is not a guarantee for a solution either. Further options include increasing the DFT grid (XXLGRID, XXXLGRID) or further increasing the last number of the integral tolerance (even though yours is pretty high already).

Hope any of this is helpful in the end!

Cheers,
Aleks

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,

If I may, I'll add a few more lines on Alessandro's comment that might be of use (or not...)

Let's say you want to change the charge state of one of the Li atoms in your cell. Assuming that it is atom number XY with atomic number 3, your input might look something like this:

Title [geometry] END [basis set] ... 3 3 -> Li example basis set 0 0 6 2. 1. 700.0 0.001421 220.0 0.003973 70.0 0.01639 20.0 0.089954 5.0 0.315646 1.5 0.494595 0 0 1 1. 1. 0.5 1.0 0 2 1 0. 1. 0.6 1.0 ... 99 0 CHEMOD 1 XY -> atom number whose charge you change for the initial guess 2.0 2.0 0.0 -> electronic charge of all shells in the basis set, here you alter the initial charge guess (note that the initial configuration is 2.0 1.0 0.0, so you will end up in the "-1" charge state) CHARGED END [SCF parameters] END

As usual, supercell size should be converged, localization of the defect confirmed, etc. Good luck with the defect formation energy corrections...one way to go is to use the multipole correction (aka Makov-Payne, see e.g., 2010 J. Phys.: Conf. Ser. 242 012004 or PRB 81, 205214 (2010))...or you can write your own piece of code to interface CRYSTAL's output to for example the scheme of Kumagai and Oba (Phys. Rev. B 89, 195205) !

Hope it helps.

Cheers,
Aleks