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,

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,

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,

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

Thankyou so much. This will help me

Thankyou so much. It worked perfectly. This is a big help

wonderful, thank you both

GiacomoAmbrogio Great! It works now! Thanks a lot!

Hi Giacomo,

thank you for your answer.
You are right, D4 is mentioned in the JCTC article as "further developments to the code", which I misunderstood.
The article clearly mentions it as future development, not as list of further developments not laid out in detail in the article.

Thank you for your help in clarifying this.

Kind regards,
Georg

Dear Jonas,

A good option just for visualization is using the CRYSPLOT webpage, which includes a tool for this based on JSMol and only requires your output file. Another option would be to use JMol, which requires a local installation of Java. Using JSMol locally is not an easy task, as its developed for its usage on web servers. If you really need to do so, you can find some guidance in the the following sites (it will depend on which browser you want to use, and probably also on your OS):

Recent thread on usage of JSMol locally JSMol Wiki post

Hope this helps.

Marcos

Hi Jonas,

Would you mind adding the INPUT file as well, for easier checks?

It generally doesn't look too bad, I would wait until the calculations finishes and see how the results look like. There is nothing obviously "wrong", but you can always try tightening the convergence criteria. See earlier post for a great overview of options (especially if you notice in the end that you have negative frequencies appearing!):
Earlier post on "Imaginary frequencies", see comment by Prof. Erba
I personally wouldn't run a geometry optimization together with a frequency calculation, I would split them in too, but I guess that is a matter of taste as the code can handle it : )

Cheers,
Aleks

Dear Aparajita,

Gryffindor said in Problem with restarting CPKS calculation:

I’m reaching out again regarding the CPKS step. After several attempts (and a fair bit of persistence!), I was finally able to complete the SCF calculations.

Good!

Gryffindor said in Problem with restarting CPKS calculation:

As per your previous suggestions, I ensured the SCF was converged beforehand, but unfortunately, the CPKS has never been able to finish successfully on my side.

Is there a reason for switching DIIS off in your CPKS calculation?

Also, you can change the convergence criterion for the CPKS with the TOLALPHA keyword. In your case, it may be enough to run it with:

CPKS TOLALPHA 2 END

Gryffindor said in Problem with restarting CPKS calculation:

Since I need these results quite urgently, I was wondering if you could try running it on your end to see if there's anything I'm missing?

I am afraid I can not. It is a huge calculation that you are running on 800+ atoms/cell and I do not have the computing power for this.

Thanks a lot! I will take your advice.

Dear Jonas,

first of all, 450 atoms is impressive indeed!
Now moving on to a few comments/thoughts that might be of use to you.

i) From your input file, it appears you still have a few symmetry operators in your simulation cell. Worth checking whether/if that applies to the atoms in the Fe-centred octahedra as that might force a particular solution to the electronic configurations.

ii) Given that your octahedra are not aligned with any particular crystallographic axis in your system, it is not always straightforward to extract the (d-)orbital occupations. You can have a look at the ROTREF keyword in the manual by which you can rotate the eigenvectors in the properties calculation to orient the frame along a principal axis.

iii) In Scenario 1 you are right - you removed 1 H atom, from which indeed the symmetry of the octahedron has been reduced/broken. You could test the case without "neutralization", to see whether the extra electron would stay on the site or delocalize...formation energies might be a good guide...?

iv) Scenario 2 is a bit more complicated. Given that you removed a whole NH4+ tetrahedron, I would expect significant structural distortions occurring in the surrounding of the defect. From there, a conducting state might be not so unrealistic - did you optimize the atomic positions and/or cell parameters following the inclusion of that defect? If so, could be worth seeing where the metallic states originate from (e.g., in the DOS or a charge density difference plot would be sufficient...don't forget SMEAR). The question is where is the extra electron from the anticipated Fe3+ ending if you remove the whole tetrahedron?

v) Finally, I guess it really depends on what you see experimentally and what would be the most representative simulation cell matching the measurements to be able to extract meaningful data for the oxidation state. What could be potentially useful is to take simple bulk structures for which you know the oxidation state of Fe (e.g., Fe2O3, Fe3O4, ...), analyse the corresponding charge densities (e.g., Mulliken) and have a reference state for comparing the charge on the Fe-ion in your Struvite structure.

Hope any of this is useful.

Cheers,
Aleks