Bug: ABACUS cannot support orbitals with components whose angular momentums are higher than 4

Describe the bug

I am working on orbitals with angular momentum much higher than all present versions, trying to achieve the chemical accuracy as much as possible. For some elements, the H orbital (l = 5) is needed, however, in ABACUS what is implemented is:

void UnitCell::read_orb_file(int it, std::string &orb_file, std::ofstream &ofs_running, Atom* atom)
        // omit
        if (strcmp("Sorbital-->", word) == 0)
        {
            ModuleBase::GlobalFunc::READ_VALUE(ifs, atom->l_nchi[L]);
            atom->nw += (2*L + 1) * atom->l_nchi[L];
            std::stringstream ss;
            ss << "L=" << L << ", number of zeta";
            ModuleBase::GlobalFunc::OUT(ofs_running,ss.str(),atom->l_nchi[L]);
            L++;
        }
        if (strcmp("Porbital-->", word) == 0)
        {
            ModuleBase::GlobalFunc::READ_VALUE(ifs, atom->l_nchi[L]);
            atom->nw += (2*L + 1) * atom->l_nchi[L];
            std::stringstream ss;
            ss << "L=" << L << ", number of zeta";
            ModuleBase::GlobalFunc::OUT(ofs_running,ss.str(),atom->l_nchi[L]);
            L++;
        }
        if (strcmp("Dorbital-->", word) == 0)
        {
            ModuleBase::GlobalFunc::READ_VALUE(ifs, atom->l_nchi[L]);
            atom->nw += (2*L + 1) * atom->l_nchi[L];
            std::stringstream ss;
            ss << "L=" << L << ", number of zeta";
            ModuleBase::GlobalFunc::OUT(ofs_running,ss.str(),atom->l_nchi[L]);
            L++;
        }
        if (strcmp("Forbital-->", word) == 0)
        {
            ModuleBase::GlobalFunc::READ_VALUE(ifs, atom->l_nchi[L]);
            atom->nw += (2*L + 1) * atom->l_nchi[L];
            std::stringstream ss;
            ss << "L=" << L << ", number of zeta";
            ModuleBase::GlobalFunc::OUT(ofs_running,ss.str(),atom->l_nchi[L]);
            L++;
        }
        if (strcmp("Gorbital-->", word) == 0)
        {
            ModuleBase::GlobalFunc::READ_VALUE(ifs, atom->l_nchi[L]);
            atom->nw += (2*L + 1) * atom->l_nchi[L];
            std::stringstream ss;
            ss << "L=" << L << ", number of zeta";
            ModuleBase::GlobalFunc::OUT(ofs_running,ss.str(),atom->l_nchi[L]);
            L++;
        }
    }
    ifs.close();
    if(!atom->nw)
    {
        ModuleBase::WARNING_QUIT("read_orb_file","get nw = 0");
    }

The need of high angular momentum's orbital is not quite intuitive, it is not the full story that only orbital with plus one angular momentum is enough, instead, take GTO as example, one of the most famous GTO basis, Dunning Correlation-Consistent basis set even includes I-orbitals for cc-pv6Z. See: https://www.basissetexchange.org/

The Complete Basis Sets (CBS) limit is very important for precise quantum chemistry calculations, in which both a systematic construction of basis and some extrapolation technique are required.

Expected behavior

No response

To Reproduce

No response

Environment

No response

Additional Context

No response

Task list for Issue attackers (only for developers)

[ ] Verify the issue is not a duplicate.
[ ] Describe the bug.
[ ] Steps to reproduce.
[ ] Expected behavior.
[ ] Error message.
[ ] Environment details.
[ ] Additional context.
[ ] Assign a priority level (low, medium, high, urgent).
[ ] Assign the issue to a team member.
[ ] Label the issue with relevant tags.
[ ] Identify possible related issues.
[ ] Create a unit test or automated test to reproduce the bug (if applicable).
[ ] Fix the bug.
[ ] Test the fix.
[ ] Update documentation (if necessary).
[ ] Close the issue and inform the reporter (if applicable).

deepmodeling / abacus-develop