Proper method to generate f-metals' numerical atomic orbital

kirk0830 commented 9 months ago

Details

The need to calculate f-metals like Ce, La, Er is not non-negligible. However, it seems the orbital generation is not a trivial task because even if I have included both sides of the equilibrium bond length of dimer, the yielded orbital can still not reproduce correct lattice parameter of Er crystal, the deviation is as large as 0.5 Angstrom, which is totally unacceptable. My SIAB_INPUT is:

#--------------------------------------------------------------------------------
#1. CMD & ENV
EXE_mpi    mpirun -np 16
EXE_pw     /usr/local/bin/abacus
EXE_opt    /root/ABACUS-orbitals/SIAB/opt_orb_pytorch_dpsi/main.py
#-------------------------------------------------------------------------------- 
#2. Electronic calculatation
 element     Er          # Element Name
 Ecut        100         # in Ry
 Rcut        7 8 9 10 # in Bohr
 Pseudo_dir  /root/pseudopotentials/NCPP-PD04-PBE
 Pseudo_name Er-sp.PD04.PBE.UPF
 sigma       0.05        # energy range for gauss smearing (in Ry) 

#--------------------------------------------------------------------------------
#3. Reference structure related parameters for PW calculation
#For the built-in structure types (including 'dimer', 'trimer' and 'tetramer'):
#STRU Name   #STRU Type  #nbands #MaxL   #nspin  #Bond Length list
 STRU1       dimer       66      4       1      3.2 3.4 3.6 4.0 4.4 5.1
 STRU2       trimer      88      4       1      3.2 3.4 3.6 4.0 4.4 5.1

#-------------------------------------------------------------------------------- 
#4. SIAB calculatation
 max_steps    9000
# Orbital configure and reference target for each level
#LevelIndex   #Ref STRU Name    #Ref Bands   #InputOrb    #OrbitalConf
 Level1      STRU1           auto        none        2s1p1d1f
 Level2      STRU1           auto        fix         4s2p2d2f1g
 Level3      STRU2           auto        fix         6s3p3d3f2g

#--------------------------------------------------------------------------------
#5. Save Orbitals
#Index       #LevelNum    #OrbitalType
 Save1       Level1       SZ           
 Save2       Level2       DZP         
 Save3       Level3       TZDP

What is interesting is that I find both dimer and trimer have the lowest energy bond length at near 4.0 Angstrom, this coincidences with the lattice parameter (Er crystal has lattice parameters a=b=c, alpha=beta=gamma=60 deg). I guess if there is any over-fitting happens.

Additional notes

For ABACUS pw, the lattice parameters are calculated correctly.

Have you read FAQ on the online manual http://abacus.deepmodeling.com/en/latest/community/faq.html

[X] Yes, I have read the FAQ part on online manual.

Task list for Issue attackers (only for developers)

[ ] Understand the problem or question described by the user.
[ ] Check if the issue is a known problem or has been addressed in the documentation.
[ ] Test the issue or problem on a similar system or environment, if possible.
[ ] Identify the root cause or provide clarification on the user's question.
[ ] Provide a step-by-step guide, including any necessary resources, to resolve the issue or answer the question.
[ ] If the issue is related to documentation, update the documentation to prevent future confusion (optional).
[ ] If the issue is related to code, consider implementing a fix or improvement (optional).
[ ] Review and incorporate any relevant feedback from users or developers.
[ ] Ensure the user's issue is resolved or their question is answered and close the ticket.

kirk0830 commented 9 months ago

see work log here: https://ucoyxk075n.feishu.cn/docx/YgSgd9EH8orUdWx9EBgcoITxnIc

kirk0830 commented 9 months ago

After orbital generation

Finite difference orbital test script (PW):

# ABACUS run command
nproc=1
omp=1
# Test setting
cell_parameter_factor=0.0
stepsize=0.1
nsteps=5
# System description
element="Er"
a1=3.6999610056
a2=0.0
a3=0.0
b1=-1.2333207420
b2=3.4883578692
b3=0.0
c1=-1.2333208824
c2=-1.7441786700
c3=3.0210047064
# Pseudopotential
pseudo_dir="./"
pseudopotential="Er-sp.PD04.PBE.UPF"

# ------------------------------ #
cat > INPUT << EOF
INPUT_PARAMETERS
calculation scf

pseudo_dir .

basis_type pw
ecutwfc 200

nbands 44
nspin 1

ks_solver cg
scf_nmax 200

symmetry 0

smearing_method gauss
smearing_sigma 0.05

EOF

cat > KPT << EOF
K_POINTS
0
Gamma
8 8 8 0 0 0

EOF

cat > STRU << EOF
ATOMIC_SPECIES
_element_ 1.0 _pseudopotential_

LATTICE_CONSTANT
1.889726877

LATTICE_VECTORS
_a1_ _a2_ _a3_
_b1_ _b2_ _b3_
_c1_ _c2_ _c3_

ATOMIC_POSITIONS
Cartesian

_element_
0.00
1
0.0 0.0 0.0 m 0 0 0

EOF

abacus_command="OMP_NUM_THREADS=$omp mpirun -np $nproc abacus"
cell_parameter_factor=$(echo "$cell_parameter_factor - $stepsize * $nsteps" | bc)
for ((i=0; i<=2*nsteps; i++))
do
    echo "Cell parameter factor: $cell_parameter_factor, Step: $i"
    mkdir "cell_parameter_scan_$cell_parameter_factor"
    cd "cell_parameter_scan_$cell_parameter_factor"
    cp ../INPUT .
    cp ../KPT .
    cp ../$pseudopotential .
    # remember to avoid 0.1 write as .1
    a1_str=$(echo "scale=10; $a1 * (1.0 + $cell_parameter_factor)" | bc)
    a2_str=$(echo "scale=10; $a2 * (1.0 + $cell_parameter_factor)" | bc)
    a3_str=$(echo "scale=10; $a3 * (1.0 + $cell_parameter_factor)" | bc)
    b1_str=$(echo "scale=10; $b1 * (1.0 + $cell_parameter_factor)" | bc)
    b2_str=$(echo "scale=10; $b2 * (1.0 + $cell_parameter_factor)" | bc)
    b3_str=$(echo "scale=10; $b3 * (1.0 + $cell_parameter_factor)" | bc)
    c1_str=$(echo "scale=10; $c1 * (1.0 + $cell_parameter_factor)" | bc)
    c2_str=$(echo "scale=10; $c2 * (1.0 + $cell_parameter_factor)" | bc)
    c3_str=$(echo "scale=10; $c3 * (1.0 + $cell_parameter_factor)" | bc)
    sed -e "s/_element_/$element/g" -e "s/_pseudopotential_/$pseudopotential/g" -e "s/_a1_/$a1_str/g" -e "s/_a2_/$a2_str/g" -e "s/_a3_/$a3_str/g" -e "s/_b1_/$b1_str/g" -e "s/_b2_/$b2_str/g" -e "s/_b3_/$b3_str/g" -e "s/_c1_/$c1_str/g" -e "s/_c2_/$c2_str/g" -e "s/_c3_/$c3_str/g" ../STRU > STRU
    #nohup $abacus_command > abacus.out &
    cd ..
    cell_parameter_factor=$(echo "$cell_parameter_factor + $stepsize" | bc)
done

tar -czvf cell_parameter_scan.tar.gz cell_parameter_scan_*
rm -r cell_parameter_scan_*
rm INPUT KPT STRU

Finite difference orbital test script (LCAO):

# ABACUS run command
nproc=1
omp=1
# Test setting
cell_parameter_factor=0.0
stepsize=0.1
nsteps=5
# System description
element="Er"
a1=3.6999610056
a2=0.0
a3=0.0
b1=-1.2333207420
b2=3.4883578692
b3=0.0
c1=-1.2333208824
c2=-1.7441786700
c3=3.0210047064
# Pseudopotential
pseudo_dir="./"
pseudopotential="Er-sp.PD04.PBE.UPF"
# Numerical orbital
orbital_dir="./"
rcuts=(7 8 9 10)
zeta=(2 3)

# ------------------------------ #
cat > INPUT << EOF
INPUT_PARAMETERS
calculation scf

pseudo_dir .
orbital_dir .

basis_type lcao
ecutwfc 200

nbands 44
nspin 1

ks_solver genelpa
scf_nmax 200

symmetry 0

smearing_method gauss
smearing_sigma 0.05

EOF

cat > KPT << EOF
K_POINTS
0
Gamma
8 8 8 0 0 0

EOF

cat > STRU << EOF
ATOMIC_SPECIES
_element_ 1.0 _pseudopotential_

NUMERICAL_ORBITAL
_numerical_orbital_

LATTICE_CONSTANT
1.889726877

LATTICE_VECTORS
_a1_ _a2_ _a3_
_b1_ _b2_ _b3_
_c1_ _c2_ _c3_

ATOMIC_POSITIONS
Cartesian

_element_
0.00
1
0.0 0.0 0.0 m 0 0 0

EOF

abacus_command="OMP_NUM_THREADS=$omp mpirun -np $nproc abacus"
cell_parameter_factor=$(echo "$cell_parameter_factor - $stepsize * $nsteps" | bc)
for ((i=0; i<=2*nsteps; i++))
do
    for rcut in ${rcuts[@]}
    do
        for z in ${zeta[@]}
        do
            echo "Cell parameter factor: $cell_parameter_factor, Step: $i, Zeta: $z, Rcut: $rcut"
            if [ $z -eq 2 ]
            then
                orbital_configuration="DZP"
            elif [ $z -eq 3 ]
            then
                orbital_configuration="TZDP"
            fi
            folder=$orbital_configuration"_"$rcut"_cell_parameter_scan_"$cell_parameter_factor
            mkdir $folder
            cd $folder
            cp ../INPUT .
            cp ../KPT .
            cp ../$pseudopotential .
            # remember to avoid 0.1 write as .1
            a1_str=$(echo "scale=10; $a1 * (1.0 + $cell_parameter_factor)" | bc)
            a2_str=$(echo "scale=10; $a2 * (1.0 + $cell_parameter_factor)" | bc)
            a3_str=$(echo "scale=10; $a3 * (1.0 + $cell_parameter_factor)" | bc)
            b1_str=$(echo "scale=10; $b1 * (1.0 + $cell_parameter_factor)" | bc)
            b2_str=$(echo "scale=10; $b2 * (1.0 + $cell_parameter_factor)" | bc)
            b3_str=$(echo "scale=10; $b3 * (1.0 + $cell_parameter_factor)" | bc)
            c1_str=$(echo "scale=10; $c1 * (1.0 + $cell_parameter_factor)" | bc)
            c2_str=$(echo "scale=10; $c2 * (1.0 + $cell_parameter_factor)" | bc)
            c3_str=$(echo "scale=10; $c3 * (1.0 + $cell_parameter_factor)" | bc)
            cp ../STRU .
            sed -i "s/_element_/$element/g" STRU
            sed -i "s/_pseudopotential_/$pseudopotential/g" STRU
            sed -i "s/_a1_/$a1_str/g" STRU
            sed -i "s/_a2_/$a2_str/g" STRU
            sed -i "s/_a3_/$a3_str/g" STRU
            sed -i "s/_b1_/$b1_str/g" STRU
            sed -i "s/_b2_/$b2_str/g" STRU
            sed -i "s/_b3_/$b3_str/g" STRU
            sed -i "s/_c1_/$c1_str/g" STRU
            sed -i "s/_c2_/$c2_str/g" STRU
            sed -i "s/_c3_/$c3_str/g" STRU
            if [ $z -eq 2 ]
            then
                numerical_orbital=$element"_gga_"$rcut"au_200Ry_4s2p2d2f1g.orb"
            elif [ $z -eq 3 ]
            then
                numerical_orbital=$element"_gga_"$rcut"au_200Ry_6s3p3d3f2g.orb"
            fi
            sed -i "s/_numerical_orbital_/$numerical_orbital/g" STRU
            cp "../Er_sigma_0.05_200Ry_dimer/Orbital_"$element"_"$orbital_configuration"/"$numerical_orbital .

            cd ..
        done
    done

    cell_parameter_factor=$(echo "$cell_parameter_factor + $stepsize" | bc)
done

for rcut in ${rcuts[@]}
do
    for z in ${zeta[@]}
    do
        if [ $z -eq 2 ]
        then
            orbital_configuration="DZP"
        elif [ $z -eq 3 ]
        then
            orbital_configuration="TZDP"
        fi
        tar -czvf $orbital_configuration"_"$rcut"_cell_parameter_scan.tar.gz" $orbital_configuration"_"$rcut"_cell_parameter_scan_"*
        rm -r $orbital_configuration"_"$rcut"_cell_parameter_scan_"*
    done
done

# tar -czvf cell_parameter_scan.tar.gz *_cell_parameter_scan_*
# rm -r *_cell_parameter_scan_*
rm INPUT KPT STRU

kirk0830 commented 9 months ago

PW cell-relax convergence test

# ABACUS run command
nproc=1
omp=1
# Test setting
ecutwfc=20
stepsize=20
nsteps=14
# System description
element="Er"
# Pseudopotential
pseudo_dir="./"
pseudopotential="Er-sp.PD04.PBE.UPF"

cat > INPUT << EOF
INPUT_PARAMETERS
calculation scf

pseudo_dir .

basis_type pw
ecutwfc _ecutwfc_

nbands 44
nspin 1

ks_solver cg
scf_nmax 200

symmetry 0

smearing_method gauss
smearing_sigma 0.05

EOF

cat > KPT << EOF
K_POINTS
0
Gamma
8 8 8 0 0 0

EOF

cat > STRU << EOF
ATOMIC_SPECIES
_element_ 1.0 _pseudopotential_

LATTICE_CONSTANT
1.889726877

LATTICE_VECTORS
  3.42588982   0.00000000   0.00000000
 -1.14196365   3.22996099   0.00000000
 -1.14196378  -1.61498025   2.79722658

ATOMIC_POSITIONS
Cartesian

_element_
0.00
1
0.0 0.0 0.0 m 0 0 0

EOF

abacus_command="OMP_NUM_THREADS=$omp mpirun -np $nproc abacus"
sed -i "s/_element_/$element/g" STRU
sed -i "s/_pseudopotential_/$pseudopotential/g" STRU
for ((i=0; i<=nsteps; i++))
do
    echo "ecutwfc: $ecutwfc, Step: $i"
    mkdir "cell_energy_ecutwfc_scan_$ecutwfc"
    cd "cell_energy_ecutwfc_scan_$ecutwfc"
    cp ../INPUT .
    cp ../KPT .
    cp ../$pseudopotential .
    sed -i "s/_ecutwfc_/$ecutwfc/g" INPUT
    cp ../STRU .
    #nohup $abacus_command > abacus.out &
    cd ..
    ecutwfc=$(echo "$ecutwfc + $stepsize" | bc)
done

tar -czvf cell_energy_ecutwfc_scan.tar.gz cell_energy_ecutwfc_scan_*
rm -r cell_energy_ecutwfc_scan_*
rm INPUT KPT STRU

kirk0830 commented 9 months ago

Dimer (reference structure) energy convergence scan

nproc=1
omp=1
bond_length=3.6
stepsize=0.2
nsteps=10
element="Er"
pseudopotential="Er-sp.PD04.PBE.UPF"

cat > INPUT << EOF
INPUT_PARAMETERS
calculation scf

pseudo_dir .

basis_type pw
ecutwfc 200

nbands 44
nspin 1

ks_solver cg
scf_nmax 200

symmetry 0

smearing_method gauss
smearing_sigma 0.05

EOF

cat > KPT << EOF
K_POINTS
0
Gamma
1 1 1 0 0 0

EOF

cat > STRU << EOF
ATOMIC_SPECIES
_element_ 1.0 _pseudopotential_

LATTICE_CONSTANT
1.889726877

LATTICE_VECTORS
20.0 0.0 0.0
0.0 20.0 0.0
0.0 0.0 20.0

ATOMIC_POSITIONS
Cartesian

_element_
0.00
2
0.0 0.0 0.0 m 0 0 0
0.0 0.0 _bond_length_ m 0 0 1

EOF

abacus_command="OMP_NUM_THREADS=$omp mpirun -np $nproc abacus"
bond_length=$(echo "$bond_length - $stepsize * $nsteps" | bc)
for ((i=0; i<=3*nsteps; i++))
do
    echo "Bond length: $bond_length, Step: $i"
    mkdir "bond_length_scan_$bond_length"
    cd "bond_length_scan_$bond_length"
    cp ../INPUT .
    cp ../KPT .
    cp ../$pseudopotential .
    sed -e "s/_element_/$element/g" -e "s/_pseudopotential_/$pseudopotential/g" -e "s/_bond_length_/$bond_length/g" ../STRU > STRU
    #nohup $abacus_command > abacus.out &
    cd ..
    bond_length=$(echo "$bond_length + $stepsize" | bc)
done

tar -czvf bond_length_scan.tar.gz bond_length_scan_*
rm -r bond_length_scan_*
rm INPUT KPT STRU

hongriTianqi commented 9 months ago

This issue provides insights on generating f-metals' numerical atomic orbital. It is a good practice issue.

deepmodeling / abacus-develop