Fix: report error when fr-uspp used in linear calculations

[YuLiu98]

Linked Issue

Fix #4342

Unit Tests and/or Case Tests for my changes

• A unit test is updated.

What's changed?

• add an error report in function `Pseudopot_upf::average_p(const double& lambda).

[kirk0830]

There is still one thing I cannot understand. You can run with the following two pseudopotentials respectively:
has_so_test.zip
Why still the full-relativistic version of SG15 pseudopotential can run but not for the RRKJUS? Actually they both have flags like has_so as true or T. Or say, why for ultrasoft pseudopotential one must turn on lspinorb?
I have tested QE, yes QE acts the same, for SG15 full-relativistic version, the calculation did not crash:

     PseudoPot. # 1 for Si read from file:
     ./Si_ONCV_PBE_FR-1.1.upf
     MD5 check sum: 9ba06b18bf99a35e84ea732b067bf0db
     Pseudo is Norm-conserving, Zval =  4.0
     Generated using ONCVPSP code by D. R. Hamann
     Using radial grid of  602 points,  4 beta functions with: 
                l(1) =   0
                l(2) =   0
                l(3) =   1
                l(4) =   1

     atomic species   valence    mass     pseudopotential
     Si                4.00    28.08550     Si( 1.00)

     16 Sym. Ops., with inversion, found (12 have fractional translation)



   Cartesian axes

     site n.     atom                  positions (alat units)
         1        Si     tau(   1) = (   0.7500000   0.0000000   0.4356497  )

     number of k points=    30  Gaussian smearing, width (Ry)=  0.0020
                       cart. coord. in units 2pi/alat
        k(    1) = (   0.0000000   0.0000000   0.0000000), wk =   0.0092593
        k(    2) = (   0.0000000   0.0000000   0.1912852), wk =   0.0370370
        k(    3) = (   0.0000000   0.0000000   0.3825704), wk =   0.0370370
        k(    4) = (   0.0000000   0.0000000  -0.5738556), wk =   0.0185185
        k(    5) = (   0.0000000   0.1912852   0.1912852), wk =   0.0370370
        k(    6) = (   0.0000000   0.1912852   0.3825704), wk =   0.0740741
        k(    7) = (   0.0000000   0.1912852  -0.5738556), wk =   0.0370370
        k(    8) = (   0.0000000   0.3825704   0.3825704), wk =   0.0370370
        k(    9) = (   0.0000000   0.3825704  -0.5738556), wk =   0.0370370
        k(   10) = (   0.0000000  -0.5738556  -0.5738556), wk =   0.0092593
        k(   11) = (   0.1666667  -0.0956426  -0.0956426), wk =   0.0740741
        k(   12) = (   0.1666667  -0.0956426   0.2869278), wk =   0.1481481
        k(   13) = (   0.1666667  -0.0956426  -0.6694982), wk =   0.1481481
        k(   14) = (   0.1666667   0.2869278   0.2869278), wk =   0.0740741
        k(   15) = (   0.1666667   0.2869278  -0.6694982), wk =   0.1481481
        k(   16) = (   0.1666667  -0.6694982  -0.6694982), wk =   0.0740741
        k(   17) = (   0.3333333  -0.1912852  -0.1912852), wk =   0.0740741
        k(   18) = (   0.3333333  -0.1912852   0.0000000), wk =   0.0740741
        k(   19) = (   0.3333333  -0.1912852  -0.7651408), wk =   0.1481481
        k(   20) = (   0.3333333  -0.1912852  -0.5738556), wk =   0.0740741
        k(   21) = (   0.3333333   0.0000000   0.0000000), wk =   0.0185185
        k(   22) = (   0.3333333   0.0000000  -0.7651408), wk =   0.0740741
        k(   23) = (   0.3333333   0.0000000  -0.5738556), wk =   0.0370370
        k(   24) = (   0.3333333  -0.7651408  -0.7651408), wk =   0.0740741
        k(   25) = (   0.3333333  -0.7651408  -0.5738556), wk =   0.0740741
        k(   26) = (   0.3333333  -0.5738556  -0.5738556), wk =   0.0185185
        k(   27) = (  -0.5000000   0.2869278   0.2869278), wk =   0.0370370
        k(   28) = (  -0.5000000   0.2869278   0.4782130), wk =   0.1481481
        k(   29) = (  -0.5000000   0.4782130   0.4782130), wk =   0.0740741
        k(   30) = (  -0.5000000   0.4782130  -0.0956426), wk =   0.0740741

     Dense  grid:     2113 G-vectors     FFT dimensions: (  18,  18,  18)

     Estimated max dynamical RAM per process >       2.35 MB

     Initial potential from superposition of free atoms

     starting charge       3.9538, renormalised to       4.0000
     Starting wfcs are random

     total cpu time spent up to now is        0.2 secs

     Self-consistent Calculation

     iteration #  1     ecut=    30.00 Ry     beta= 0.50
     CG style diagonalization
     ethr =  1.00E-02,  avg # of iterations =  7.4

     Threshold (ethr) on eigenvalues was too large:
     Diagonalizing with lowered threshold

     CG style diagonalization
     ethr =  1.05E-04,  avg # of iterations =  5.1

     total cpu time spent up to now is        0.3 secs

     total energy              =      -7.83982561 Ry
     estimated scf accuracy    <       0.00434072 Ry

[YuLiu98]

There is still one thing I cannot understand. You can run with the following two pseudopotentials respectively: has_so_test.zip Why still the full-relativistic version of SG15 pseudopotential can run but not for the RRKJUS? Actually they both have flags like has_so as true or T. Or say, why for ultrasoft pseudopotential one must turn on lspinorb? I have tested QE, yes QE acts the same, for SG15 full-relativistic version, the calculation did not crash:

     PseudoPot. # 1 for Si read from file:
     ./Si_ONCV_PBE_FR-1.1.upf
     MD5 check sum: 9ba06b18bf99a35e84ea732b067bf0db
     Pseudo is Norm-conserving, Zval =  4.0
     Generated using ONCVPSP code by D. R. Hamann
     Using radial grid of  602 points,  4 beta functions with: 
                l(1) =   0
                l(2) =   0
                l(3) =   1
                l(4) =   1

     atomic species   valence    mass     pseudopotential
     Si                4.00    28.08550     Si( 1.00)

     16 Sym. Ops., with inversion, found (12 have fractional translation)



   Cartesian axes

     site n.     atom                  positions (alat units)
         1        Si     tau(   1) = (   0.7500000   0.0000000   0.4356497  )

     number of k points=    30  Gaussian smearing, width (Ry)=  0.0020
                       cart. coord. in units 2pi/alat
        k(    1) = (   0.0000000   0.0000000   0.0000000), wk =   0.0092593
        k(    2) = (   0.0000000   0.0000000   0.1912852), wk =   0.0370370
        k(    3) = (   0.0000000   0.0000000   0.3825704), wk =   0.0370370
        k(    4) = (   0.0000000   0.0000000  -0.5738556), wk =   0.0185185
        k(    5) = (   0.0000000   0.1912852   0.1912852), wk =   0.0370370
        k(    6) = (   0.0000000   0.1912852   0.3825704), wk =   0.0740741
        k(    7) = (   0.0000000   0.1912852  -0.5738556), wk =   0.0370370
        k(    8) = (   0.0000000   0.3825704   0.3825704), wk =   0.0370370
        k(    9) = (   0.0000000   0.3825704  -0.5738556), wk =   0.0370370
        k(   10) = (   0.0000000  -0.5738556  -0.5738556), wk =   0.0092593
        k(   11) = (   0.1666667  -0.0956426  -0.0956426), wk =   0.0740741
        k(   12) = (   0.1666667  -0.0956426   0.2869278), wk =   0.1481481
        k(   13) = (   0.1666667  -0.0956426  -0.6694982), wk =   0.1481481
        k(   14) = (   0.1666667   0.2869278   0.2869278), wk =   0.0740741
        k(   15) = (   0.1666667   0.2869278  -0.6694982), wk =   0.1481481
        k(   16) = (   0.1666667  -0.6694982  -0.6694982), wk =   0.0740741
        k(   17) = (   0.3333333  -0.1912852  -0.1912852), wk =   0.0740741
        k(   18) = (   0.3333333  -0.1912852   0.0000000), wk =   0.0740741
        k(   19) = (   0.3333333  -0.1912852  -0.7651408), wk =   0.1481481
        k(   20) = (   0.3333333  -0.1912852  -0.5738556), wk =   0.0740741
        k(   21) = (   0.3333333   0.0000000   0.0000000), wk =   0.0185185
        k(   22) = (   0.3333333   0.0000000  -0.7651408), wk =   0.0740741
        k(   23) = (   0.3333333   0.0000000  -0.5738556), wk =   0.0370370
        k(   24) = (   0.3333333  -0.7651408  -0.7651408), wk =   0.0740741
        k(   25) = (   0.3333333  -0.7651408  -0.5738556), wk =   0.0740741
        k(   26) = (   0.3333333  -0.5738556  -0.5738556), wk =   0.0185185
        k(   27) = (  -0.5000000   0.2869278   0.2869278), wk =   0.0370370
        k(   28) = (  -0.5000000   0.2869278   0.4782130), wk =   0.1481481
        k(   29) = (  -0.5000000   0.4782130   0.4782130), wk =   0.0740741
        k(   30) = (  -0.5000000   0.4782130  -0.0956426), wk =   0.0740741

     Dense  grid:     2113 G-vectors     FFT dimensions: (  18,  18,  18)

     Estimated max dynamical RAM per process >       2.35 MB

     Initial potential from superposition of free atoms

     starting charge       3.9538, renormalised to       4.0000
     Starting wfcs are random

     total cpu time spent up to now is        0.2 secs

     Self-consistent Calculation

     iteration #  1     ecut=    30.00 Ry     beta= 0.50
     CG style diagonalization
     ethr =  1.00E-02,  avg # of iterations =  7.4

     Threshold (ethr) on eigenvalues was too large:
     Diagonalizing with lowered threshold

     CG style diagonalization
     ethr =  1.05E-04,  avg # of iterations =  5.1

     total cpu time spent up to now is        0.3 secs

     total energy              =      -7.83982561 Ry
     estimated scf accuracy    <       0.00434072 Ry

Actually, I'm also confused about it. I speculate that it is related to the overlap matrix S.

QE calculations about fr-uspp report the following error:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine cdiaghg (33):
     S matrix not positive definite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S matrix in ncpp is a unit matrix. Thus, it will not suffer the risk that the S matrix is not positive definite. However, fr-uspp may be at risk of instability in non-SOC computing, so QE directly prohibits this use in the code. However, as you mentioned in the issue, the fr uspp of the H element can be calculated normally, so this instability is not absolute.

As for the theory, maybe @dyzheng knows something about it.

[kirk0830]

There is still one thing I cannot understand. You can run with the following two pseudopotentials respectively: has_so_test.zip Why still the full-relativistic version of SG15 pseudopotential can run but not for the RRKJUS? Actually they both have flags like has_so as true or T. Or say, why for ultrasoft pseudopotential one must turn on lspinorb? I have tested QE, yes QE acts the same, for SG15 full-relativistic version, the calculation did not crash:

     PseudoPot. # 1 for Si read from file:
     ./Si_ONCV_PBE_FR-1.1.upf
     MD5 check sum: 9ba06b18bf99a35e84ea732b067bf0db
     Pseudo is Norm-conserving, Zval =  4.0
     Generated using ONCVPSP code by D. R. Hamann
     Using radial grid of  602 points,  4 beta functions with: 
                l(1) =   0
                l(2) =   0
                l(3) =   1
                l(4) =   1

     atomic species   valence    mass     pseudopotential
     Si                4.00    28.08550     Si( 1.00)

     16 Sym. Ops., with inversion, found (12 have fractional translation)



   Cartesian axes

     site n.     atom                  positions (alat units)
         1        Si     tau(   1) = (   0.7500000   0.0000000   0.4356497  )

     number of k points=    30  Gaussian smearing, width (Ry)=  0.0020
                       cart. coord. in units 2pi/alat
        k(    1) = (   0.0000000   0.0000000   0.0000000), wk =   0.0092593
        k(    2) = (   0.0000000   0.0000000   0.1912852), wk =   0.0370370
        k(    3) = (   0.0000000   0.0000000   0.3825704), wk =   0.0370370
        k(    4) = (   0.0000000   0.0000000  -0.5738556), wk =   0.0185185
        k(    5) = (   0.0000000   0.1912852   0.1912852), wk =   0.0370370
        k(    6) = (   0.0000000   0.1912852   0.3825704), wk =   0.0740741
        k(    7) = (   0.0000000   0.1912852  -0.5738556), wk =   0.0370370
        k(    8) = (   0.0000000   0.3825704   0.3825704), wk =   0.0370370
        k(    9) = (   0.0000000   0.3825704  -0.5738556), wk =   0.0370370
        k(   10) = (   0.0000000  -0.5738556  -0.5738556), wk =   0.0092593
        k(   11) = (   0.1666667  -0.0956426  -0.0956426), wk =   0.0740741
        k(   12) = (   0.1666667  -0.0956426   0.2869278), wk =   0.1481481
        k(   13) = (   0.1666667  -0.0956426  -0.6694982), wk =   0.1481481
        k(   14) = (   0.1666667   0.2869278   0.2869278), wk =   0.0740741
        k(   15) = (   0.1666667   0.2869278  -0.6694982), wk =   0.1481481
        k(   16) = (   0.1666667  -0.6694982  -0.6694982), wk =   0.0740741
        k(   17) = (   0.3333333  -0.1912852  -0.1912852), wk =   0.0740741
        k(   18) = (   0.3333333  -0.1912852   0.0000000), wk =   0.0740741
        k(   19) = (   0.3333333  -0.1912852  -0.7651408), wk =   0.1481481
        k(   20) = (   0.3333333  -0.1912852  -0.5738556), wk =   0.0740741
        k(   21) = (   0.3333333   0.0000000   0.0000000), wk =   0.0185185
        k(   22) = (   0.3333333   0.0000000  -0.7651408), wk =   0.0740741
        k(   23) = (   0.3333333   0.0000000  -0.5738556), wk =   0.0370370
        k(   24) = (   0.3333333  -0.7651408  -0.7651408), wk =   0.0740741
        k(   25) = (   0.3333333  -0.7651408  -0.5738556), wk =   0.0740741
        k(   26) = (   0.3333333  -0.5738556  -0.5738556), wk =   0.0185185
        k(   27) = (  -0.5000000   0.2869278   0.2869278), wk =   0.0370370
        k(   28) = (  -0.5000000   0.2869278   0.4782130), wk =   0.1481481
        k(   29) = (  -0.5000000   0.4782130   0.4782130), wk =   0.0740741
        k(   30) = (  -0.5000000   0.4782130  -0.0956426), wk =   0.0740741

     Dense  grid:     2113 G-vectors     FFT dimensions: (  18,  18,  18)

     Estimated max dynamical RAM per process >       2.35 MB

     Initial potential from superposition of free atoms

     starting charge       3.9538, renormalised to       4.0000
     Starting wfcs are random

     total cpu time spent up to now is        0.2 secs

     Self-consistent Calculation

     iteration #  1     ecut=    30.00 Ry     beta= 0.50
     CG style diagonalization
     ethr =  1.00E-02,  avg # of iterations =  7.4

     Threshold (ethr) on eigenvalues was too large:
     Diagonalizing with lowered threshold

     CG style diagonalization
     ethr =  1.05E-04,  avg # of iterations =  5.1

     total cpu time spent up to now is        0.3 secs

     total energy              =      -7.83982561 Ry
     estimated scf accuracy    <       0.00434072 Ry

Actually, I'm also confused about it. I speculate that it is related to the overlap matrix S.

QE calculations about fr-uspp report the following error:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine cdiaghg (33):
     S matrix not positive definite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S matrix in ncpp is a unit matrix. Thus, it will not suffer the risk that the S matrix is not positive definite. However, fr-uspp may be at risk of instability in non-SOC computing, so QE directly prohibits this use in the code. However, as you mentioned in the issue, the fr uspp of the H element can be calculated normally, so this instability is not absolute.

As for the theory, maybe @dyzheng knows something about it.

I got some information from @dyzheng , it is because the average_pp (sum over the spin degree of freedom) is not implemented for USPP. In principle if lspinorb = .false., the FR PP should get result identical with SR.