Skip to content

02_add_f_orbital_for_projection.md

Latest commit

 

History

History
321 lines (263 loc) · 8.84 KB

02_add_f_orbital_for_projection.md

File metadata and controls

321 lines (263 loc) · 8.84 KB
Raw

How to Add f Orbital for Band Structure Projection

Background

When performing DFT calculations with pseudopotentials, only valence electrons are explicitly treated. The atomic orbitals output by ONCVPSP correspond to valence states and are used for band structure projection analysis.

Important: If an orbital (e.g., 4f) is defined as a core state, it will NOT be available for projection analysis, even if lmax=3 is set.

Problem: Core vs Valence States

Understanding nc and nv Parameters

In ONCVPSP input files:

  • nc = number of core shells (first nc lines in the electron configuration)
  • nv = number of valence shells (last nv lines in the electron configuration)

Important: These are shell counts (line numbers), NOT electron counts!

The valence electron count (z_valence in UPF) = sum of occupation numbers in valence shells.

Example: Hf (Hafnium, Z=72)

In the Hf.in configuration:

# atsym  z   nc   nv     iexc    psfile
Hf   72.00   10   4   4   upf
#
#   n    l    f        energy (Ha)
1    0    2.00    # Core (line 1)
2    0    2.00    # Core (line 2)
2    1    6.00    # Core (line 3)
3    0    2.00    # Core (line 4)
3    1    6.00    # Core (line 5)
3    2   10.00    # Core (line 6)
4    0    2.00    # Core (line 7)
4    1    6.00    # Core (line 8)
4    2   10.00    # Core (line 9)
4    3   14.00    # Core (line 10) <-- 4f is in CORE because nc=10
5    0    2.00    # Valence (line 11)
5    1    6.00    # Valence (line 12)
5    2    2.00    # Valence (line 13)
6    0    2.00    # Valence (line 14) -- nv=4, last 4 lines are valence

Here nc=10 means first 10 shells are core, nv=4 means last 4 shells are valence. Valence electrons = 2+6+2+2 = 12.

Result: Output only shows 5S, 5P, 5D, 6S labels. No 4F orbital because 4f (line 10) is a core state.

Example: La (Lanthanum, Z=57)

La is a better example because its 4f orbital is empty in the ground state ([Xe]5d^1 6s^2), making it suitable for adding as a valence state for projection.

Original La.in:

# atsym  z   nc   nv     iexc    psfile
La 57.00    9    4       4      upf
#
#   n    l    f
1    0    2.00    # Core (line 1)
2    0    2.00    # Core (line 2)
2    1    6.00    # Core (line 3)
3    0    2.00    # Core (line 4)
3    1    6.00    # Core (line 5)
3    2   10.00    # Core (line 6)
4    0    2.00    # Core (line 7)
4    1    6.00    # Core (line 8)
4    2   10.00    # Core (line 9) -- nc=9, first 9 lines are core
5    0    2.00    # Valence (line 10)
5    1    6.00    # Valence (line 11)
5    2    1.00    # Valence (line 12)
6    0    2.00    # Valence (line 13) -- nv=4, last 4 lines are valence

Note:

  • nc=9 means the first 9 shells (lines) are core states
  • nv=4 means the last 4 shells (lines) are valence states
  • These are shell counts, NOT electron counts!
  • The valence electron count = sum of occupations in valence shells = 2+6+1+2 = 11 (matches z_valence in UPF)

Result: Output shows 5S, 5P, 5D, 6S. No 4F because it's not defined as a valence state.

Solution: Add f Orbital to Valence States

Method 1: Small Occupation Number (No Code Modification)

Modify the input file to add 4f as a valence state with a very small occupation:

# atsym  z   nc   nv     iexc    psfile
La 57.00    9    5       4      upf    # nv changed from 4 to 5
#
#   n    l    f
...
5    0    2.00
5    1    6.00
5    2    0.9999    # Reduced by 0.0001 to maintain charge neutrality
6    0    2.00
4    3    0.0001    # Added 4f with tiny occupation (nv=5)

Key points:

  1. Increase nv by 1 (from 4 to 5)
  2. Add 4 3 f line for 4f orbital
  3. Adjust other occupations to maintain total electron count = Z

Method 2: Modify Source Code to Allow f=0

Modify src/check_data.f90 line 91:

Original:

  if(fa(ii)<=0.0d0) then
   write(6,'(a,i4)') 'test_data: f <= 0.0, reference configuration line',ii
   ierr=ierr+1
  end if

Modified:

  if(fa(ii)<0.0d0) then
   write(6,'(a,i4)') 'test_data: f < 0.0, reference configuration line',ii
   ierr=ierr+1
  end if

Then recompile:

make clean && make

Now you can use f=0.00:

5    2    1.00    # Keep original 5d occupation
4    3    0.00    # 4f with zero occupation

Test Results: Effect of Occupation Number

We tested various occupation numbers for the 4f orbital:

f value total_psenergy (Ha) 4F pseudo_energy (Ha) Status
0.10 -63.4873446289 -0.2793709287 OK
1e-4 -63.4814332406 -0.3193059465 OK
1e-6 -63.4814258635 -0.3193459784 OK
1e-8 -63.4814257898 -0.3193463788 OK
1e-10 -63.4814257888 -0.3193463828 Converged
1e-12 -63.4814257888 -0.3193463828 Converged
1e-15 -63.4814257888 -0.3193463828 Converged
0.00 -63.4814257888 -0.3193463828 Converged (requires code mod)
Original (no 4f) -63.4814257887 - Reference

Key findings:

  1. Results converge from f=1e-10 onwards
  2. f=0 and f=1e-15 produce identical results
  3. Converged results match the original pseudopotential (without 4f valence)

Detailed Numerical Comparison: f=1e-15 vs f=0 vs Original

1. Total Pseudopotential Energy

f=1e-15:  -6.34814257888E+01 Ha
f=0:      -6.34814257888E+01 Ha
Original: -6.34814257887E+01 Ha

|f=1e-15 - f=0|:   0.00E+00 Ha (identical)
|f=0 - Original|:  1.00E-10 Ha (numerical precision)

2. Existing Orbitals (5S, 5P, 5D, 6S) - Are they affected?

| Orbital | Max |f=1e-15 - f=0| | Max |f=0 - Original| | Conclusion | |---------|---------------------|----------------------|------------| | 5S (PP_CHI.1) | 1.00E-13 | 1.00E-13 | Identical | | 5P (PP_CHI.2) | 8.67E-19 | 5.00E-14 | Identical | | 5D (PP_CHI.3) | 1.82E-18 | 1.00E-15 | Identical | | 6S (PP_CHI.4) | 1.17E-17 | 5.00E-13 | Identical |

3. l=3 Projector (PP_BETA.7)

f=1e-15:  [-1.4506469056e-09, 5.8127905831e-07, ...]
f=0:      [-1.4506471393e-09, 5.8127905811e-07, ...]
Original: [-1.4506454575e-09, 5.8127905831e-07, ...]

Max |f=1e-15 - f=0|:   4.00E-13 (identical)
Max |f=0 - Original|:  4.00E-13 (identical)

4. New 4F Orbital (PP_CHI.5)

f=1e-15:  [-1.6450954975e-10, 1.3107582116e-07, 2.0955881596e-06, ...]
f=0:      [-1.6450954932e-10, 1.3107582155e-07, 2.0955881659e-06, ...]
Original: N/A (no 4f valence state)

Max |f=1e-15 - f=0|:   9.30E-13 (identical)

Conclusions

  1. f=1e-15 and f=0 produce IDENTICAL results (differences < 1e-12)

    • No need to modify source code, just use f=1e-15

  2. Adding 4f valence state does NOT affect other orbitals:

    • 5S, 5P, 5D, 6S wavefunctions: IDENTICAL to original
    • l=3 projector: difference < 1e-12 (numerical precision)
    • Total energy: difference < 1e-10 Ha (numerical precision)

  3. The generated pseudopotential is RELIABLE and fully consistent with the original.

Recommended Approach

For Quick Implementation (No Code Change)

Use f=1e-10 or smaller:

La 57.00    9    5       4      upf
...
5    2    0.9999999999
4    3    0.0000000001

For Clean Implementation (With Code Change)

  1. Modify src/check_data.f90 as described above
  2. Recompile with make clean && make
  3. Use f=0.00:
La 57.00    9    5       4      upf
...
5    2    1.00
4    3    0.00

Complete Example: La.in with 4f Orbital

# ATOM AND REFERENCE CONFIGURATION
# atsym  z   nc   nv     iexc    psfile
La 57.00    9    5       4      upf
#
#   n    l    f        energy (Ha)
1    0    2.00
2    0    2.00
2    1    6.00
3    0    2.00
3    1    6.00
3    2   10.00
4    0    2.00
4    1    6.00
4    2   10.00
5    0    2.00
5    1    6.00
5    2    1.00
6    0    2.00
4    3    0.00
#
# PSEUDOPOTENTIAL AND OPTIMIZATION
# lmax
3
#
#   l,   rc,      ep,       ncon, nbas, qcut
0   2.10000  -1.43044    4    8   5.40000
1   1.60000  -0.82320    4    8   6.60000
2   1.85000  -0.11357    4    9   6.60000
3   1.65000  -0.15967    4    9  10.60000
#
# LOCAL POTENTIAL
# lloc, lpopt,  rc(5),   dvloc0
4    5   1.60000      2.00000
#
# VANDERBILT-KLEINMAN-BYLANDER PROJECTORs
# l, nproj, debl
0    2   1.29521
1    2   4.00000
2    2   4.00000
3    2   3.00000
#
# MODEL CORE CHARGE
# icmod, fcfact, rcfact
3   5.00000   1.35000
#
# LOG DERIVATIVE ANALYSIS
# epsh1, epsh2, depsh
-12.00   12.00    0.02
#
# OUTPUT GRID
# rlmax, drl
6.00    0.01
#
# TEST CONFIGURATIONS
# ncnf
0
# nvcnf
#   n    l    f

Verification

Run ONCVPSP and check for 4F label:

./src/oncvpsp.x < La.in | grep label=

Expected output:

       label="5S"
       label="5P"
       label="5D"
       label="6S"
       label="4F"

Summary

Approach Pros Cons
Small f (1e-10) No code modification needed Slightly less clean
f=0 (code mod) Clean, exact zero occupation Requires recompilation

Both approaches produce identical pseudopotential results and are suitable for band structure projection analysis.