Questions about printing out unoccupied molecular orbitals (LUMO, LUMO+x)

[zh]

Dear all,

When I printed the unoccupied molecular orbitals of C6H6 into the cube
files and analyzed them, I met the followings problem:
The shapes of unoccupied molecular orbitals obtained with cp2k (GAPW
calculation) are quite different from those obtained with other codes
(LCAO method and Plane-wave pseudopotential method).

The keywords specified in the input file for generating the cube files
of unoccupied molecular orbitals:
....
MAX_ITER_LUMO 5999
EPS_LUMO 3.E-7
ADDED_MOS 10
...
&PRINT
&MO_CUBES
NHOMO 3
NLUMO 10
WRITE_CUBE .TRUE.
&END MO_CUBES
&END PRINT
....

I also examined the log file and found the following information:
---
Eigenvalues of the occupied subspace spin 1
and 10 added MO eigenvalues
---------------------------------------------
-9.87289173 -9.87276638 -9.87276248 -9.87242760
-9.87242383 -9.87226307 -0.76294011 -0.66142985
-0.66141836 -0.52895626 -0.52895605 -0.45746985
-0.39629304 -0.38443742 -0.36038167 -0.36037151
-0.31546968 -0.28755156 -0.28754954 -0.21567922
-0.21566949 -0.02617647 -0.02617516 0.03512159
0.07206783 0.07215809 0.10479005 0.10527223
0.11638929 0.13027147 0.24632554
Fermi Energy [eV] : -5.868665

Lowest Eigenvalues of the unoccupied subspace spin 1
-----------------------------------------------------
Reached convergence in 4483 iterations
0.24706383 0.25840280 0.25869375 0.28401710
0.32929682 0.37178993 0.37837600 0.37872848
0.39283865 0.40999697

HOMO - LUMO gap [eV] : 12.591614

****
The problem in the log file is that the 10 lowest eigenvalues of the
unoccupied subspace are quite larger than the 10 added MO eigenvalues.

I also read the source code related to these parts (i.e.,
qs_scf_post_gpw.F) and found that the unoccupied orbitals are
recalculated through the subroutines "make_lumo()", "
"ot_eigensolver()", and "calculate_subspace_eigenvalues()".

My questions are the followings:
i) Since the unoccupied molecular orbitals are already obtained by
other subroutine "get_mo_set()" before calling "make_lumo()", why are
the unoccupied molecular orbitals needed to be recalculated?
Actually, the convergence of the recalculation of unoccupied molecular
orbitals in "ot_eigensolver()" is quite slow and the eigenvalues of
unoccupied molecular orbitals become unreasonable.
ii) How to get the reasonable unoccupied molecular orbitals,
especially for the LUMO, LUMO+1, LUMO+2, and LUMO+3?

Thanks for the help.

[flo]

Hi,

at first, I agree, that there is something strange, as the HOMO-LUMO
gap is simply wrong. 0.007au will never be 12 eV, but as always,
without an input file, the it is almost impossible to reproduce the
problem and to solve it.

Related to your question i):
get_mo_set is only used to have access to the data and data structures
in the mo_set type. Nothing is computed there. The unoccupied orbitals
are indeed never computed before (for CP2K it is enough to know about
the occupied orbitals during SCF), that's why you need to compute them
when make_lumo() is called. The eigensolver is needed, as the orbitals
are the eigenvectors of the Kohn-Sham matrix.

Question ii):
see problem of missing input file

Flo

[zh]

Thank you very much for your reply.

The input and log files are included in the attached 'lumos-
problems.tar.gz'

My input file:

&FORCE_EVAL
METHOD Quickstep
&DFT
# LSD
&MGRID
CUTOFF 280
&END MGRID
&QS
METHOD GAPW
EPS_DEFAULT 1.0D-14
EPS_GVG 1.0E-8
EPS_PGF_ORB 1.0E-8
QUADRATURE GC_LOG
EPSFIT 1.E-4
EPSISO 1.0E-12
EPSRHO0 1.E-10
LMAXN0 2
LMAXN1 6
ALPHA0_H 10
ALPHA0_S 10
&END QS
&SCF
SCF_GUESS RESTART
EPS_SCF 3.0E-7
MAX_SCF 500
&DIAGONALIZATION
ALGORITHM STANDARD
&END DIAGONALIZATION

MAX_ITER_LUMO 5999
EPS_LUMO 3.E-7
ADDED_MOS 10

&MIXING
# ALPHA 0.05
METHOD BROYDEN_MIXING
N_SIMPLE_MIX 10
ALPHA 0.1
BETA 1.50
NBROYDEN 6
&END MIXING
&SMEAR F
METHOD FERMI_DIRAC
ELECTRONIC_TEMPERATURE [K] 600
&END SMEAR
&OT F
&END OT
&END SCF
BASIS_SET_FILE_NAME ../EMSL_BASIS_SETS
POTENTIAL_FILE_NAME ../POTENTIAL
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC XYZ
POISSON_SOLVER PERIODIC
&END POISSON
&PRINT
&MO
EIGENVALUES
OCCUPATION_NUMBERS
&END MO

&MO_CUBES
NHOMO 3
NLUMO 10
WRITE_CUBE .TRUE.
&END MO_CUBES
&END PRINT

&END DFT
&SUBSYS
&CELL
ABC 12 12 12
&END CELL
&COORD
C 0.0000000001 1.3998977949 6.000000000
C -1.2123471147 0.6999589170 6.000000000
C -1.2123471149 -0.6999589168 6.000000000
C -0.0000000001 -1.3998977949 6.000000000
C 1.2123471148 -0.6999589169 6.000000000
C 1.2123471149 0.6999589168 6.000000000
H 0.0000000000 2.4923734834 6.000000000
H -2.1584646684 1.2461924428 6.000000000
H -2.1584646683 -1.2461924428 6.000000000
H 0.0000000000 -2.4923734835 6.000000000
H 2.1584646684 -1.2461924428 6.000000000
H 2.1584646683 1.2461924428 6.000000000
&END COORD
&KIND H
BASIS_SET 6-311G**
POTENTIAL ALL
&END KIND
&KIND C
BASIS_SET 6-311G**
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT c6h6
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 50
OPTIMIZER BFGS
&END GEO_OPT
&END MOTION

[zh]

In my test calculations, I also found another problem:
If the keyword "ADDED_MOS" in the section &SCF is set as 0 or 1, the
LUMO and LUMO+1 orbitals of C6H6 molecule obtained with cp2k will be
similar to those obtained with other codes.

In the log file, there are two values printed out for the LUMO-HOMO
gap.

The first value of LUMO-HOMO gap is printed out due to the following
keyword:

"
&MO
EIGENVALUES
OCCUPATION_NUMBERS
&END MO
"

This one is reasonable. That is why I say that the unoccupied
molecular orbitals are already obtained before calling "make_lumo()".

The second value of LUMO-gap is printed out due to the following
keyword:

"
&MO_CUBES
NHOMO 3
NLUMO 10
WRITE_CUBE .TRUE.
&END MO_CUBES
"

While the second one is quite larger than the first one.

[marci]

Hi

Actually, I do not think this is really a problem of the
calculations,
but a misunderstanding in the definition of the MOs.
When you use the keyword ADDED_MOS together with the
diagonalization scheme, at each SCF step all the occupied MOs
plus a number of unoccupied MOs equal to the ADDED_MOS are
computed (eigenvectors, eigenvalues, and occupation numbers).
If you print through the section
&MO
...
&END
you can get the related eigenvalues and occupations, and from these
the gap.

Now, if on top of these you print through the section
&MO_CUBES
...
&END
and you require NLUMO different from zero,
NLUMO states are computed a posteriori, in addition to those
already computed along the SCF procedure. These are therefore higher
in energy, and do not correspond to the lowest LUMOS, but to
the lowest states laying above those already computed.
If you want the cubes of the lowest unoccupied states, and you have
ADDED_MOS
different from zero, you should consider that some of them are
obtained with
the keyword NHOMO.
I agree that in this situation the keyword names NHOMO and NLUMO in
the section MO_CUBES
are misleading.

kind regards
marcella

[Teodoro Laino]

Just a comment:

It is not the definition of the NHOMO and NLUMO to be not clear in
MO_CUBES rather it is the description of the keyword ADDED_MOS.
One should specify in that keyword that the ADDED MOS are considered
like HOMOS in post-processing calculations (although they are not filled).

Teo

[marci]

I do not fully agree. The ADDED_MOS keyword is mainly thought for
calculations
with smeared occupation numbers. In this case the number of partially
occupies MOs
is larger than the number of electrons, and definition homo and lumo
are not anymore a good
way to classify the MOs
True is that in the case illustrated by "zh"
there is no need to add MOs for the SCF and the unoccupied
states can be safely calculated a posteriori just for printing
purposes,
i.e. ADDED_MOS can be set to zero.

cheers
marcella

[Teodoro Laino]

Again: one should specify that if you use the ADDED_MOS keyword these
(in post-processing) will be considered like occupied orbitals
(irrespective of their occupation number) : see example provided by "zh".

I'm not discussing whether HOMO or LUMO definition holds when using
smearing or not.
I'm just saying that in MO_CUBES it is quite clear what a HOMO or a LUMO
is, in standard conditions. And the same definition still holds if one
classifies the ADDED_MOS as occupied orbitals (even if their occupation
number is zero).
Alternatively after one uses the ADDED_MOS, you should remove the
ADDED_MOS which have occupation number set to zero from the set of
occupied orbitals. This would definitely make all consistent..

Teo

[zh]

Thank all of your for your reply and discussion!

Now I understand the problem and know how to get the unoccupied
molecular orbitals in a reasonable way.

The problem can be summarized as follows:
i) ADDED_MOS = 0, NLUMO >=0
The cube files of unoccupied molecular orbitals will be generated in a
normal way (i.e., what we expected). The eigenvalues and eigevectors
of unoccupied molecular orbitals from LUMO to HOMO+NLUMO will be
calculated by calling "make_lumo()". But we should take care of
whether these unoccupied molecular orbitals are well converged or not.
ii) ADDED_MOS >0 , NLUMO =0
No cube files of unoccupied molecular orbitals will be generated. The
eigenvalues of unoccupied molecular orbitals from LUMO to HOMO
+ADDED_MOS will be calculated during the self-consistent calculation.
iii) ADDED_MOS >0, NLUMO >0
The eigenvalues and eigevectors of unoccupied molecular orbitals from
LUMO to HOMO+ADDED_MOS will be calculated during the self-consistent
calculation, while the eigenvalues and eigevectors of unoccupied
molecular orbitals from HOMO+ADD_MOS+1 to HOMO+ADD_MOS+NLUMO will be
calculated by calling "make_lumo()". The index for the cube files of
unoccupied molecular orbitals is wrong. In the current version of
source code, the index for the cube files of unoccupied molecular
orbitals starts from LUMO. Actually, it should start from HOMO
+ADDED_MOS+1.
iv) ADDED_MOS =0, NLUMO = -1
The cube files of unoccupied molecular orbitals will be generated in a
normal way (i.e., what we expected). The eigenvalues and eigevectors
of unoccupied molecular orbitals from LUMO to NAO-HOMO (NAO is the
number of atomic orbitals) will be calculated by calling
"make_lumo()". But we should take care of whether these unoccupied
molecular orbitals are well converged or not.
v)ADDED_MOS >0, NLUMO = -1
The eigenvalues and eigevectors of unoccupied molecular orbitals from
LUMO to HOMO+ADDED_MOS will be calculated during the self-consistent
calculation, while the eigenvalues and eigevectors of unoccupied
molecular orbitals from HOMO+ADD_MOS+1 to NAO-HOMO (NAO is the number
of atomic orbitals) will be calculated by calling "make_lumo()". The
index for the cube files of unoccupied molecular orbitals is wrong.
In the current version of source code, the index for the cube files of
unoccupied molecular orbitals starts from LUMO. Actually, it should
start from HOMO+ADDED_MOS+1.

[Leopold Talirz]

Dear all,

since members from our group also fall into this 'trap' every once in a while, I am wondering whether the situation could be improved or if at least one could print out a drastic warning if one uses &MO_CUBES together with &ADDED_MOS..

Just to sum up the "state of the art": Let's say I have a system of 10 occupied orbitals that I calculate with 5 ADDED_MOS.

Then, printing &MO_CUBES with NHOMO 2 and NLUMO 2 will print the cube files of orbitals number 9, 10 and 16, 17 (the cube files of the latter two will be incorrectly labeled with indices 11, 12).

I have no possibility to print the orbitals 11-15 unless I drop the ADDED_MOS keyword.

I think if the cube files were indexed correctly that would already help a lot (then, one would realize immediately that the orbitals 11-15 are missing).

Maybe it would also make sense to let NHOMO 2 print the orbitals 14,15 instead of 9,10 ... 

What do you think?

Best,

Leopold

[Leopold Talirz]

Hi,

First of all, a big thanks to Marcella for fixing this in rev. 14153.

We just noticed that there still seems to be a bug in the title line of the cube files.

E.g. for a system with 56 electrons (see attachment), we find

head  -n2 PROJ-WFN_00028_1-1_0.cube

-Quickstep-

  WAVEFUNCTION           28 spin            1 i.e. HOMO -            0

head  -n2 PROJ-WFN_00029_1-1_0.cube

-Quickstep-

  WAVEFUNCTION           29 spin            1 i.e. LUMO +           56

I.e. the number counting for the empty states is incorrect.

The problem probably lies around line 1394 of qs_scf_post_gpw.F

  WRITE(title,*) "WAVEFUNCTION ",index_mo," spin ",ispin," i.e. LUMO + ",ifirst+ivector-2

There are quite a few index variables in this function (ifirst, ivector, index_mo, homo and lumo) and I am not quite sure where the mistake lies, but I guess for the author it will be easy to spot.

Please also confirm that (as I suspect) the error only concerns the "LUMO + .." index and that the index_mo is printed correctly.

Best,

Leopold