gen_f_orbital.sh

#!/bin/bash
##################################
# Author : cndaqiang             #
# Update : 2024-12-15            #
# Build  : 2024-12-15            #
# What   : Generate ONCVPSP with #
#          additional f orbital  #
#          for band projection   #
##################################

# ONCVPSP directory (version 3.3.1)
ONCVDIR=/home/chendq/soft/mvapich/source/oncvpsp-3.3.1

# Functional
FUN=PBE

# Elements with lmax=3 that need f orbital for projection
# These elements have empty 4f in ground state, suitable for adding as valence
ELEMENTS="La"
# Other lmax=3 elements (Au Hf Hg Ir Os Pb Pt Re Ta Tl W) have filled 4f in core
# They can be added if needed, but may require different treatment

# Input directory
INDIR=../sg15

echo "=========================================="
echo "Generating ONCVPSP with f orbital"
echo "ONCVPSP version: 3.3.1"
echo "Functional: $FUN"
echo "=========================================="

for element in $ELEMENTS; do
    infile="${INDIR}/${element}.in"

    if [ ! -f "$infile" ]; then
        echo "ERROR: $infile not found"
        continue
    fi

    echo ""
    echo "Processing: $element"
    echo "----------------------------------------"

    # Read original nc and nv
    nc=$(sed -n '3p' "$infile" | awk '{print $3}')
    nv=$(sed -n '3p' "$infile" | awk '{print $4}')
    new_nv=$((nv + 1))

    echo "  Original: nc=$nc, nv=$nv"
    echo "  Modified: nc=$nc, nv=$new_nv (added 4f)"

    # Create modified input file
    # 1. Change nv from $nv to $new_nv in line 3
    # 2. Add "4 3 0.0000000001" after the last valence line (line nc+nv+5)
    #    and adjust the previous valence occupation to maintain charge neutrality

    # For La: 5d has 2 electrons, we change to 1.9999999999 and add 4f with 0.0000000001

    outfile="${element}_f.in"

    # Create the modified input
    awk -v nc="$nc" -v nv="$nv" -v new_nv="$new_nv" '
    BEGIN {
        valence_start = nc + 5 + 1  # Line number where valence starts (after header and core)
        valence_end = nc + nv + 5   # Line number of last valence
    }
    NR == 3 {
        # Modify nv in header line
        $4 = new_nv
        print $1, $2, $3, $4, $5, $6
        next
    }
    NR == valence_end {
        # Last valence line - reduce occupation by 1e-10
        n = $1
        l = $2
        f = $3 - 0.0000000001
        printf "%d %d %.10f\n", n, l, f
        next
    }
    NR == valence_end + 1 {
        # Insert 4f line before this line
        print "4 3 0.0000000001"
        print $0
        next
    }
    { print }
    ' "$infile" > "$outfile"

    echo "  Created: $outfile"

    # Generate scalar relativistic UPF
    echo "  Generating scalar relativistic..."
    ${ONCVDIR}/src/oncvpsp.x < "$outfile" > "${element}_f.out" 2>&1

    if grep -q "ERROR" "${element}_f.out"; then
        echo "  ERROR in generation! Check ${element}_f.out"
        mv "${element}_f.out" "${element}_ONCV_${FUN}_sr.fail"
    else
        # Extract UPF
        awk 'BEGIN{out=0};/END_PSP/{out=0}; {if(out == 1) {print}}; /Begin PSP_UPF/{out=1}' \
            "${element}_f.out" > "${element}_ONCV_${FUN}_sr.upf"

        # Verify 4F label exists
        if grep -q 'label="4F"' "${element}_ONCV_${FUN}_sr.upf"; then
            echo "  SUCCESS: ${element}_ONCV_${FUN}_sr.upf (with 4F orbital)"
        else
            echo "  WARNING: 4F orbital not found in output!"
        fi
        rm -f "${element}_f.out"
    fi

    # Generate fully relativistic UPF
    echo "  Generating fully relativistic..."
    ${ONCVDIR}/src/oncvpspr.x < "$outfile" > "${element}_f.out" 2>&1

    if grep -q "ERROR" "${element}_f.out"; then
        echo "  ERROR in generation! Check ${element}_f.out"
        mv "${element}_f.out" "${element}_ONCV_${FUN}_fr.fail"
    else
        # Extract UPF
        awk 'BEGIN{out=0};/END_PSP/{out=0}; {if(out == 1) {print}}; /Begin PSP_UPF/{out=1}' \
            "${element}_f.out" > "${element}_ONCV_${FUN}_fr.upf"

        # Verify 4F label exists
        if grep -q 'label="4F"' "${element}_ONCV_${FUN}_fr.upf"; then
            echo "  SUCCESS: ${element}_ONCV_${FUN}_fr.upf (with 4F orbital)"
        else
            echo "  WARNING: 4F orbital not found in output!"
        fi
        rm -f "${element}_f.out"
    fi
done

echo ""
echo "=========================================="
echo "Generation complete!"
echo "=========================================="
echo ""
echo "Generated files:"
ls -la *.upf 2>/dev/null
echo ""
echo "Failed files:"
ls -la *.fail 2>/dev/null || echo "  None"