[Q-e-developers] gamma_only inconsistency in xspectra
Guido Fratesi
fratesi at mater.unimib.it
Wed Dec 18 12:23:06 CET 2013
Dear QE developers,
let me report about xspectra an inconsistency between the flag
gamma_only (read from pw save file) and the kpoint mesh, which is read
from input.
That would produce wrong results for NEXAFS calculations of, e.g.,
molecules in the gas phase, where pw.x is commonly run with gamma_only
algorithms, unless the user is especially aware and runs an intermediate
pw.x nscf with k-points.
Comments, results, and scripts for a quick test for C2H2 molecule are
attached for your consideration. The problem is also subtle because in
small molecules the spectrum may immediately appear as "wrong", but only
careful analysis shows that the spectrum is indeed incorrect for larger
molecules including several atoms of the same species.
Best wishes,
Guido Fratesi
PS the archive (1.5MB) is available from
http://www.sendspace.com/file/cy3sr3
--
Guido Fratesi
Dipartimento di Scienza dei Materiali
Universita` degli Studi di Milano-Bicocca
via Cozzi 55, 20125 Milano, Italy
Phone: +39 02 6448 5183
email: fratesi at mater.unimib.it
web: https://sites.google.com/site/guidofratesi/
-------------- next part --------------
Comments on the effect of gamma_only flag in xspectra.
The K-point mesh in xspectra.x (XSPECTRAKP) is in principle not
related to the one (PWKP) used for the previous pw.x execution, since
one only needs the charge density to build the Hamiltonian, and wave
functions are neglected. Hence, one can define the mesh XSPECTRAKP
differently than PWKP. However, the gamma_only flag in xspectra.x is
currently set according to PWKP and not to XSPECTRAKP, linking
unnecessarily the two meshes.
There are indeed cases in which the PWKP mesh can be coarser than
XSPECTRAKP. Consider a molecule-in-a-box, for which PWKP=gamma is
fine, but the description of unbound states (as relevant for xspectra)
is improved by kpoint sampling. So one could be in the position of
running xspectra, with k points different than (0,0,0), after a
gamma_only pw.x calculation, and that will erroneously set the
gamma_only flag to true also in xspectra.x (with consequences
described below).
One could avoid the problem by running an intermediate nscf
calculation with pw.x without gamma_only algorithms, but the code
should rather set the variable correctly. Alternatively, but this
seems an unnecessary limitation, issue an error if gamma_only flag
read from PWKP is inconsistent with XSPECTRAKP.
Consequences, i.e. how the error shows up:
Transitions should be at energies of unoccupied states, so we can
compare xspectra results to the PDOS with the proper symmetry (2p_x,
if E||x for C1s to valence excitations).
Here I take C2H2 molecule (aligned along x axis, E||x) and compute the
following cases:
1) PWKP=gamma_only; XSPECTRAKP = (0,0,0)
2) PWKP=gamma_only; XSPECTRAKP not (0,0,0)
3) PWKP=(0,0,0), not gamma; XSPECTRAKP = (0,0,0)
4) PWKP=(0,0,0), not gamma; XSPECTRAKP not (0,0,0)
The case (2), blue curve in the upper panel of the figure, gives
grossly wrong results, while the others are consistent with the DOS.
Notice that the first transitions (at 6-7eV) correspond to bound
states and do not depend on K, so we can compare results at different
K's. At higher energies, spectra do depend on K.
In the attached archive, I enclose my results as well as the scripts
to generate them:
1-run-scf-nscf.bash runs the scf at gamma and the nscf with gamma_only
or without gamma_only
2-run-postproc.bash projwfc to get pdos_x at the ionized atom
3-run-NEXAFS.bash xspectra calculation starting from different PWKP,
and adopting a kpoint=(0,0,0) or different than (0,0,0) in xspectra
4-make-plots.bash self-evident
The problem arises in 5.0.2, as well as in earlier versions (e.g.,
4.3.1).
Guido Fratesi
December 18th 2013
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plot-dos-vs-nexafs.gif
Type: image/gif
Size: 55890 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20131218/825e5e75/attachment.gif>
More information about the developers
mailing list