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

Hi,

Hmmm...even though the SI file of the paper you have shared does point to using HSE06, it remains not entirely clear whether its parameters were kept default or altered (for all systems together or individually for instance).

There are also differences between what different codes consider as "default" settings. For example, within the code used in the paper (VASP, nice code, no doubt), the defaults for HSE06 read (taken from https://www.vasp.at/wiki/index.php/List_of_hybrid_functionals) :
$$ \omega= 0.2\ \mathring{A} , \quad c = 0.25, \quad \text{correlation}=\text{PBE}, $$ with the first number reading the range separation parameter (omega) and the second the fraction of exact exchange used (c).

Within CRYSTAL (also nice code, no doubt), these read (taken from the manual, page 138):
$$ \omega= 0.11\ a_0^{-1}, \quad c = 0.25, \quad \text{correlation}=\text{PBE}, $$ adopting the same labels.

Not sure about the exact definition of units (perhaps a developer can comment if this is indeed Bohr radius as assumed?), but you can already see the subtle differences having to be taken into account when comparing between codes.

A few other thoughts worth considering:

• In the paper, the structure was optimized with PBEsol and on top of that geometry HSE06 was applied as a single-point calculation. Not sure about the exact composition of those ZIFs, but the structural differences could play a significant role as well (planewave codes are very costly when optimizing a structure with hybrid functionals). Here is also a good read on this topic: doi.org/10.1088/2516-1075/aafc4b

• One final small comment. Within the PAW formalism implemented in VASP, scalar relativistic effects are included in the pseudopotentials by default. No problem, cool feature, but should be taken into account when comparing results, especially for heavier elements (longer discussion found here https://blog.vasp.at/forum/viewtopic.php?t=902)

Hope this helps!

Cheers,
Aleks

Hi,

from a quick look at your input/output files, it doesn't seem that the intel error is the final one, as the code started struggling with converging the SCF solution. You can monitor the behaviour with each cycle, where at cycle 8 you end up with:

...
 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT NUMDFT      TELAPSE      156.79 TCPU      151.49
 CYC   8 ETOT(AU)                 NaN DETOT       NaN tst  7.97E+01 PX  1.03E+02
 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT MPP_TRAF    TELAPSE      157.35 TCPU      152.05
...

after which the code stops. Without unnecessarily repetition, I would rather point you to a great source itself:
https://tutorials.crystalsolutions.eu/tutorial.html?td=SCF_conv&tf=conv_tut

Give a shout if you still struggle after trying those things.
And perhaps this part could be moved to the geometry optimization sub-topic.

Cheers,
A

Dear Gryffindor,

I had a quick look at your input/output files. You have increased the number of cycles for which the number of (up-down) electrons will be fixed, without increased the actual number of SCF cycles. This can be done with the MAXCYCLE keyword. For example, in your input file:

...
SHRINK
1 1
MAXCYCLE
200
SPINLOCK
0 100
...

You can always verify that by checking the exact value in the output file, which in your case reads:

...
*******************************************************************************
 MAX NUMBER OF SCF CYCLES      50  CONVERGENCE ON DELTAP        10**-16               <---
 WEIGHT OF F(I) IN F(I+1)      30% CONVERGENCE ON ENERGY        10**- 6
 SHRINK. FACT.(MONKH.)    1  1  1  NUMBER OF K POINTS IN THE IBZ      1
 SHRINKING FACTOR(GILAT NET)    1  NUMBER OF K POINTS(GILAT NET)      1
 *******************************************************************************
...

Also, a final note, care is advised to not have the number of SPINLOCK cycles for too long as you might converge while locking is active. Ideally, the SPINLOCK directive should remain active until the SCF process has found a numerically stable path along for which the system converges to the required configuration.

Hope this helps,
Cheers,
Aleks

Hi,
if you could add your input/output files or something like screenshots of the same, that could be a good start for debugging!
Cheers,
A