    
    export generate


    """
    `Basics.generate(representation::AtomicState.Representation)`  
        ... to generate an atomic representation as specified by the representation.repType::AbstractRepresentationType.
            All relevant intermediate and final results are printed to screen (stdout). Nothing is returned.

    `Basics.perform(representation::AtomicState.Representation; output=true)`  
        ... to generate the same but to return the complete output in a dictionary; the particular output depends on the type and 
            specifications of the representation but can easily accessed by the keys of this dictionary.
    """
    function Basics.generate(representation::AtomicState.Representation; output::Bool=false)
        results = Basics.generate(representation.repType, representation; output=output)
        
        Defaults.warn(PrintWarnings())
        Defaults.warn(ResetWarnings())
        println(" ")
        return( results )
    end



    """
    `Basics.generate(repType::AtomicState.MeanFieldBasis, representation::AtomicState.Representation)`  
        ... to generate a mean-field basis (representation) for a set of reference configurations; all relevant intermediate 
            and final results are printed to screen (stdout). Nothing is returned.

    `Basics.generate(repType::AtomicState.MeanFieldBasis, representation::AtomicState.Representation; output=true)`  
        ... to generate the same but to return the complete output in a dictionary; the particular output depends on the type and 
            specifications of the representation but can easily accessed by the keys of this dictionary.
    """
    function Basics.generate(repType::AtomicState.MeanFieldBasis, rep::AtomicState.Representation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel    = rep.nuclearModel

        # The asfSettings only define the SCF part and are partly derived from the MeanFieldSettings
        asfSettings = AsfSettings(AsfSettings(); scField=repType.settings.scField) 
        
        basis      = Basics.performSCF(rep.refConfigs, nModel, rep.grid, asfSettings; printout=true)
        if output    results = Base.merge( results, Dict("mean-field basis" => basis) )          end
        
        return( results )
    end



    """
    `Basics.generate(repType::AtomicState.CiExpansion, representation::AtomicState.Representation)`  
        ... to generate a configuration-interaction expansion for a single level symmetry and based on a set of reference configurations
            and a number of pre-specified steps. All relevant intermediate and final results are printed to screen (stdout). 
            Nothing is returned.

    `Basics.generate(repType::AtomicState.CiExpansion, representation::AtomicState.Representation; output=true)`  
        ... to generate the same but to return the complete output in a dictionary; the particular output depends on the type and 
            specifications of the computations but can easily accessed by the keys of this dictionary.
    """
    function Basics.generate(repType::AtomicState.CiExpansion, rep::AtomicState.Representation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel    = rep.nuclearModel
        orbitals  = repType.applyOrbitals
        
        # Generate a list of relativistic configurations and  CSF's for the given subshell list
        relconfList = ConfigurationR[]
        for  conf in rep.refConfigs
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        ##x if  false    for  i = 1:length(relconfList)    println(">> include ", relconfList[i])    end   end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        ##x Defaults.setDefaults("relativistic subshell list", subshellList; printout=printout)
        csfList = CsfR[]
        for  relconf in relconfList
            newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
            append!( csfList, newCsfs)
        end
        
        # Generate all level symmetries that need to be taken into account
        if     length(repType.settings.levelSelectionCI.symmetries) != 0    symmetries = repType.settings.levelSelectionCI.symmetries
        else   symmetries = LevelSymmetry[]
            for csf in csfList
                if      LevelSymmetry(csf.J, csf.parity)  in  symmetries
                else    push!( symmetries, LevelSymmetry(csf.J, csf.parity) )
                end
            end
        end
        println("*** Level symmetries = $symmetries ")

        # The asfSettings only define the CI part and are partly derived from the CiSettings
        asfSettings = AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(),    0, 0., Subshell[], Subshell[], 
                                  repType.settings.eeInteractionCI, NoneQed(), LSjjSettings(false), 
                                  repType.settings.levelSelectionCI) 
        
        basis      = Basics.generateBasis(rep.refConfigs, symmetries, repType.excitations)
        basis      = Basis( true, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals )  
        # Test that all orbitals are available 
        for  subshell in  basis.subshells   
            if   !(haskey(orbitals, subshell))  error("No orbital available for subshell $subshell ")     end
        end
        multiplet  = Basics.performCI(basis,  nModel, rep.grid, asfSettings; printout=true) 
        if output    results = Base.merge( results, Dict("CI multiplet" => Multiplet("CI multiplet:", multiplet.levels)) )              end
        
        return( results )
    end



    """
    `Basics.generate(repType::AtomicState.RasExpansion, representation::AtomicState.Representation)`  
        ... to generate a restricted active-space expansion for a single level symmetry and based on a set of reference configurations
            and a number of pre-specified steps. All relevant intermediate and final results are printed to screen (stdout). 
            Nothing is returned.

    `Basics.generate(repType::AtomicState.RasExpansion, representation::AtomicState.Representation; output=true)`  
        ... to generate the same but to return the complete output in a dictionary; the particular output depends on the type and 
            specifications of the computations but can easily accessed by the keys of this dictionary.
    """
    function Basics.generate(repType::AtomicState.RasExpansion, rep::AtomicState.Representation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel   = rep.nuclearModel
        # First perform a SCF+CI computations for the reference configurations below to generate a spectrum of start orbitals
        asfSettings    = AsfSettings()  ## Use default settings to define a first multiplet from the reference configurations;
                                        ## all further details are specified for each step
        priorBasis     = Basics.performSCF(rep.refConfigs, nModel, rep.grid, asfSettings; printout=true)
        priorMultiplet = Basics.performCI(priorBasis, nModel, rep.grid, asfSettings; printout=true)
        nuclearPot     = Nuclear.nuclearPotential(nModel, rep.grid)
        ## electronicPot  = Basics.compute("radial potential: Dirac-Fock-Slater", rep.grid, priorMultiplet.levels[1].basis)
        ## meanPot        = Basics.add(nuclearPot, electronicPot)
        subshellList   = Basics.extractRelativisticSubshellList(rep)             ## extract all subshells that occur in the RAS computation
        startOrbitals  = Basics.generateOrbitalsForPotential(rep.grid, nuclearPot, subshellList)  ## generate a spectrum of sufficient size
        if output    results = Base.merge( results, Dict("reference multiplet" => Multiplet("Reference multiplet:", priorMultiplet.levels) ) )  end

        # The asfSettings only define the CI part of the RAS steps and partly derived from the RasSettings
        asfSettings = AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(),    0, 0., Subshell[], Subshell[], 
                                  repType.settings.eeInteractionCI, NoneQed(), LSjjSettings(true), repType.settings.levelSelectionCI ) 
        
        # Now, cycle over all steps of the RasExpansion
        for (istep, step)  in  enumerate(repType.steps)
            println("")
            printstyled("++ Compute the orbitals, orbitals and multiplet for step $istep ... \n", color=:light_green)
            printstyled("--------------------------------------------------------------      \n", color=:light_green)
            basis      = Basics.generateBasis(rep.refConfigs, [repType.symmetry], step)
            orbitals   = Basics.generateOrbitalsForBasis(basis, step.frozenShells, priorMultiplet.levels[1].basis, startOrbitals)
            basis      = Basis( true, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals )  
            basis      = Basics.performSCF(basis, nModel, rep.grid, step.frozenShells, repType.settings; printout=true)
            multiplet  = Basics.performCI(basis,  nModel, rep.grid, asfSettings; printout=true) 
            if output    results = Base.merge( results, Dict("step"*string(istep) => Multiplet("Multiplet:", multiplet.levels)) )              end
            priorMultiplet = multiplet
        end
        
        return( results )
    end



    """
    `Basics.generate(repType::AtomicState.GreenExpansion, representation::AtomicState.Representation)`  
        ... to generate a Green (function) expansion for a given approach and excitation scheme of the electron,
            based on a set of reference configurations, a list of level symmetries as well as for given settings.
            All relevant intermediate and final results are printed to screen (stdout). Nothing is returned.

    `Basics.generate(repType::AtomicState.GreenExpansion, representation::AtomicState.Representation; output=true)`  
        ... to generate the same but to return the complete output in a dictionary; the particular output depends on the type and 
            specifications of the representation but can easily accessed by the keys of this dictionary.
    """
    function Basics.generate(repType::AtomicState.GreenExpansion, rep::AtomicState.Representation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel    = rep.nuclearModel
        settings  = repType.settings
        # First perform a SCF+CI computations for the reference configurations below to generate a spectrum of start orbitals
        asfSettings   = AsfSettings()  ## Use default settings to define a first multiplet from the reference configurations;
                                       ## all further details are specified for each step
        refBasis      = Basics.performSCF(rep.refConfigs, nModel, rep.grid, asfSettings; printout=true)
        refMultiplet  = Basics.performCI(refBasis, nModel, rep.grid, asfSettings; printout=true)
        nuclearPot    = Nuclear.nuclearPotential(nModel, rep.grid)
        electronicPot = Basics.compute("radial potential: Dirac-Fock-Slater", rep.grid, refBasis)
        meanPot       = Basics.add(nuclearPot, electronicPot)
        
        println("")
        printstyled("Compute an approximate Green function expansion ... \n", color=:light_green)
        printstyled("--------------------------------------------------- \n", color=:light_green)
        
        channels = AtomicState.GreenChannel[]
        
        # Generate all (non-relativistic) configurations from the bound configurations due to the given excitation scheme 
        confList = Basics.generateConfigurationsForExcitationScheme(rep.refConfigs, repType.excitationScheme, settings.nMax, settings.lValues)
        # Print (if required) information about the generated configuration list
        if  settings.printBefore    Basics.display(stdout, confList)    end
        
        # Generate shell list abd a full single-electron spectrum for this potential
        subshellList = Basics.extractRelativisticSubshellList(confList)                      ## extract all subshells that occur in confList
        orbitals     = Basics.generateOrbitalsForPotential(rep.grid, meanPot, subshellList)  ## generate a spectrum of sufficient size
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=true)
        Basics.display(stdout, orbitals, rep.grid)

        # The asfSettings only define the CI part of the Green channels and are partly derived from the GreenSettings
        asfSettings = AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(),    0, 0., Subshell[], Subshell[], 
                                  CoulombInteraction(), NoneQed(), LSjjSettings(false), settings.levelSelection ) 
        
        # Cycle over all selected level symmetries to generate the requested channels
        for  levelSymmetry  in  repType.levelSymmetries
            basis      = Basics.generateBasis(confList, [levelSymmetry])
            basis      = Basis( true, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals )  
            multiplet  = Basics.computeMultipletForGreenApproach(repType.approach, basis, nModel, rep.grid, asfSettings, settings; printout=true) 
            push!( channels, GreenChannel( levelSymmetry, multiplet) )
        end        
        
        # Print all results to screen
        Basics.display(stdout, channels)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    Basics.display(iostream, channels)         end
        
        if output    results = Base.merge( results, Dict("Green channels" => channels) )   end
        return( results )
    end
    
    

    """
    `Basics.generate("condensed multiplet: by single weight", multiplet::Multiplet)`  
        ... to condense/reduce the number of CSF in the basis of the given multiplet due to a single 'weight'; 
            a multiplet::Multiplet is returned.  **Not yet implemented !**
    """
    function Basics.generate(sa::String, multiplet::Multiplet)
        !(sa == "condensed multiplet: by single weight")   &&   error("Unsupported keystring = $sa")
        error("Not yet implemented !")
    end



    """
    `Basics.generate("configuration list: NR, from basis", basis::Basis)`  
        ... to (re-) generate the list of NR configurations from the given basis; a confList::Array{Configuration,1} is returned.
    """
    function Basics.generate(sa::String, basis::Basis)
        !(sa == "configuration list: NR, from basis")   &&   error("Unsupported keystring = $sa")
        confList    = Configuration[]
        NoElectrons = sum( basis.csfs[1].occupation )
        for  csf in basis.csfs
            shellList = Dict{Shell,Int64}[]
            for  i = 1:length(csf.occupation)
                n = csf.subshell[i].n;    l = Basics.subshell_l( csf.subshell[i] );    sh = Shell(n, l)
                if     haskey(shellList, sh)    shellList(sh) = shellList(sh) + csf.occupation[i]
                else   shellList = merge( shellList, Dict(sh => csf.occupation[i]) )
                end
            end
            confNew = Configuration( shellList, NoElectrons )
            push!( confList, confNew )
        end
        confListNew = ConfigurationsExcludeDoubles(confList)
        
        return( confListNew )
    end


    """
    `Basics.generate("configuration list: NR, single-configuration", refConf::Configuration, NoExcitations::Int64, fromShells::Array{Shell,1},`
                     toShells::Array{Shell,1})  
        ... to generate a non-relativistic configuration list, including the given reference configuration (refConf) and with 
            all configurations that differ by NoExcitations from the fromShells into the toShells; an Array{Configuration,1} 
            is returned.
    """
    function Basics.generate(sa::String, refConf::Configuration, NoExcitations::Int64, fromShells::Array{Shell,1}, toShells::Array{Shell,1})
        !(sa == "configuration list: NR, single-configuration")   &&   error("Unsupported keystring = $sa")
        confList = [refConf]
        # First prepare a proper reference configuration that also includes all fromShells and toShells with zero occupation
        shellDict = deepcopy(refConf.shells)
        allShells = deepcopy(fromShells)
        for  sh in toShells
            if     sh in allShells
            else   append!(allShells, [sh] )
            end
        end
        
        for  sh in allShells
        if  !haskey(shellDict, sh )    shellDict = merge( shellDict, Dict( sh => 0 ) )    end
        end

        rconf = Configuration(shellDict, refConf.NoElectrons)
        println("aa", rconf)

        # Add single excitations
        if      NoExcitations == 0
        elseif  NoExcitations == 1
            for  fromsh in fromShells
                for  tosh in toShells
                    confnew = deepcopy(rconf)
                    confnew.shells[fromsh] = confnew.shells[fromsh] - 1
                    confnew.shells[tosh]   = confnew.shells[tosh]   + 1
                    println("bb", confnew)

                    # Check that this is valid configuration and append, if appropriate
                    add = true
                    for (k,v) in  confnew.shells
                        if  v <= -1  ||  Parity(rconf) != Parity(confnew)    add = false    end
                    end
                    if  add 
                        for (k,v) in confnew.shells   
                        if  v == 0   delete!(confnew.shells, k)   end
                        end
                        push!(confList, confnew);   println("fromsh = $fromsh  tosh = $tosh")   
                    end
                end
            end

        # Add double excitations
        elseif  NoExcitations == 2
            for  fromsha in fromShells
            for  fromshb in fromShells
                for  tosha in toShells
                for  toshb in toShells
                    confnew = deepcopy(rconf)
                    confnew.shells[fromsha] = confnew.shells[fromsha] - 1
                    confnew.shells[fromshb] = confnew.shells[fromshb] - 1
                    confnew.shells[tosha]   = confnew.shells[tosha]   + 1
                    confnew.shells[toshb]   = confnew.shells[toshb]   + 1
                    println("bb", confnew)

                    # Check that this is valid configuration and append, if appropriate
                    add = true
                    for (k,v) in  confnew.shells
                        if  v <= -1  ||  Parity(rconf) != Parity(confnew)    add = false    end
                    end
                    if  add    
                        for (k,v) in confnew.shells   
                        if  v == 0   delete!(confnew.shells, k)   end
                        end
                        push!(confList, confnew)    
                    end
                end
                end
            end
            end

        else
            error("stop a")
        end

        # Exclude all configurations that appear twice (and more) in the list
        confList = Basics.ConfigurationsExcludeDoubles(confList)

        return( confList )
    end


    """
    `Basics.generate("configuration list: relativistic", conf::Configuration)`  
        ... to split/decompose a non-relativistic configuration into an list of relativistic ConfigurationR[]. The proper 
            occupuation of the relativistic subshells is taken into account.
    """
    function Basics.generate(sa::String, conf::Configuration)
        subshellList = Subshell[]
        confList     = ConfigurationR[]

        !(sa == "configuration list: relativistic")  &&   error("Unsupported keystring = $sa")

        initialize = true;    NoElectrons = 0
        for (k, occ)  in conf.shells
            NoElectrons     = NoElectrons + occ
            wa              = Basics.shellSplitOccupation(k, occ)
            subshellListNew = Dict{Subshell,Int64}[]

            if  initialize
                subshellListNew = wa;    initialize = false
            else
                for  s in 1:length(subshellList)
                    for  a in 1:length(wa)
                        wb = Base.merge( subshellList[s], wa[a] )
                        push!(subshellListNew, wb)
                    end
                end
            end
            subshellList = deepcopy(subshellListNew)
        end
        
        for subsh in subshellList
        wa = ConfigurationR(subsh, NoElectrons)
        push!(confList, wa)
        end

        return( confList )
    end


    """
    `Basics.generate("CSF list: from single ConfigurationR", conf::ConfigurationR, subshellList::Array{Subshell,1})` 
        ... to construct from a given (relativistic) configuration all possible CSF with regard to the subshell order as specified 
            by subshellList; a list::Array{CsfR,1} is returned.
    """
    function Basics.generate(sa::String, conf::ConfigurationR, subshellList::Array{Subshell,1})
        parity  = Basics.determineParity(conf)
        csfList = CsfR[];   useStandardSubshells = true;    first = true;    previousCsfs = CsfR[]
        # subhshellList = Subshell[];   
        for  subsh in subshellList
            if   subsh in keys(conf.subshells)    occ = conf.subshells[subsh]    else    occ = 0    end
            if   first
                stateList   = ManyElectron.provideSubshellStates(subsh, occ)
                currentCsfs = CsfR[]
                for  state in stateList
                    push!( currentCsfs, CsfR( true, AngularJ64(state.Jsub2//2), parity, [state.occ], [state.nu],
                                                [AngularJ64(state.Jsub2//2)], [AngularJ64(state.Jsub2//2)], Subshell[]) )
                end
                previousCsfs = copy(currentCsfs)
                first        = false
            else
                # Now support also all couplings of the subshell states with the CSFs that were built-up so far
                stateList   = ManyElectron.provideSubshellStates(subsh, occ)
                currentCsfs = CsfR[]
                for  csf in  previousCsfs
                    for  state in stateList
                        occupation = deepcopy(csf.occupation);    seniorityNr = deepcopy(csf.seniorityNr);    
                        subshellJ  = deepcopy(csf.subshellJ);     subshells = deepcopy(csf.subshells)
                        push!(occupation, state.occ);   push!(seniorityNr, state.nu);   push!(subshellJ, AngularJ64(state.Jsub2//2) ) 
                        push!(subshells, subsh)
                        newXList = oplus( csf.subshellX[end], AngularJ64(state.Jsub2//2) )
                        for  newX in newXList
                            subshellX = deepcopy(csf.subshellX);   push!(subshellX, newX) 
                            push!( currentCsfs, CsfR( true, subshellX[end], parity, occupation, seniorityNr, subshellJ, subshellX, Subshell[]) ) 
                        end
                    end
                end
                previousCsfs = copy(currentCsfs)
            end
        end
        
        return( previousCsfs )
    end

    
    """
    `Basics.generate("shells: ordered list for NR configurations", confs::Array{Configuration,1})`  
        ... to generate for confs, i.e. all the given (non-relativistic) configurations, a common and ordered shell list; 
            a list::Array{Shell,1} is returned.
    """
    function Basics.generate(sa::String, confs::Array{Configuration,1})
        shells = Shell[]   

        !(sa == "shells: ordered list for NR configurations")  &&   error("Unsupported keystring = $sa")

        wa = Defaults.getDefaults("ordered shell list: non-relativistic", 7)
        for  a in wa
            for  cf in 1:length(confs)
                ks = keys(confs[cf].shells)
                if  a in ks   push!(shells, a);    break    end
            end 
        end

        return( shells )
    end

    
    """
    `Basics.generate("subshells: ordered list for relativistic configurations", confs::Array{ConfigurationR,1})`  
        ... to generate for confs, i.e. all the given (relativistic) configurations, common and ordered subshell list; 
            a list::Array{Subshell,1} is returned.
    """
    function Basics.generate(sa::String, confs::Array{ConfigurationR,1})
        subshells = Subshell[]   

        !(sa == "subshells: ordered list for relativistic configurations")  &&   error("Unsupported keystring = $sa")
        
        for  conf in confs
            for  subsh in keys(conf.subshells)      push!(subshells, subsh)     end
        end
        subshells = Base.unique(subshells)
        subshells = Base.sort( subshells, lt=Base.isless)
        ## Do include 'empty' subshells into the subshell list if they are specified by the given configurations
        ## if  a in ks   &&   confs[cf].subshells[a]  !=  0     push!(subshells, a);    break    end
        

        #== wa = Defaults.getDefaults("ordered subshell list: relativistic", 7)
        for  a in wa
            for  cf in 1:length(confs)
                ks = keys(confs[cf].subshells)
                ## Do include 'empty' subshells into the subshell list if they are specified by the given configurations
                ## if  a in ks   &&   confs[cf].subshells[a]  !=  0     push!(subshells, a);    break    end
                if  a in ks     push!(subshells, a);    break    end
            end 
        end ==#

        return( subshells )
    end

    
    """
    `Basics.generate("subshells: ordered list for two bases", basisA::Basis,  basisB::Basis)`  
        ... to generate common and ordered subshell list for the two basis A and B; a list::Array{Subshell,1} is returned.
    """
    function Basics.generate(sa::String, basisA::Basis,  basisB::Basis)
        function areEqual(nn::Int64, sha::Array{Subshell,1}, shb::Array{Subshell,1})
            # Determines whether the first nn subshells are equal in sha and shb (true) or not (false)
            for  i = 1:nx   
                if    sha[i] != shb[i]    return( false )   end
            end
            return( true )
        end
            
        subshells = Subshell[]   

        !(sa == "subshells: ordered list for two bases")  &&   error("Unsupported keystring = $sa")

        nx = min(length(basisA.subshells), length(basisB.subshells))
        if  areEqual(nx, basisA.subshells, basisB.subshells)
            # If subshell order is equal, any subshell order is accepted for those subshells that occur in both basis
            for  i = 1:nx   
                if    basisA.subshells[i] != basisB.subshells[i]   error("Inconsistent subshells of two bases.")
                else  push!( subshells, basisA.subshells[i])
                end
            end
            #
            if       length(basisA.subshells) > nx   
                for  i = nx+1:length(basisA.subshells)    push!(subshells, basisA.subshells[i])    end
            elseif   length(basisB.subshells) > nx   
                for  i = nx+1:length(basisB.subshells)    push!(subshells, basisB.subshells[i])    end
            end
        else
            println("basisA.subshells = $(basisA.subshells)")
            println("basisB.subshells = $(basisB.subshells)")
            # If subshell order is NOT equal, all subshell in basisA and basisB must follow standard order
            standardList = Defaults.getDefaults("ordered subshell list: relativistic", 7)
            na = 0;  nb = 0
            for  sh in standardList   
                if  sh in basisA.subshells  ||   sh in basisB.subshells   push!( subshells, sh)    end    
            end
            nn = 0;   
            for sh in basisA.subshells   
                if    !(sh in subshells) push!( subshells, sh)    
                else  wb = findall(x->x==sh, subshells);    if wb[1] <= nn   error("stop a")   else   nn = wb[1]  end    
                end
            end
            nn = 0;   
            for sh in basisB.subshells   
                if    !(sh in subshells) push!( subshells, sh)    
                else  wb = findall(x->x==sh, subshells);    if wb[1] <= nn   error("stop a")   else   nn = wb[1]  end    
                end
            end
            println("*** extended subshells from two basis = $subshells")
        end

        return( subshells )
    end


    
    """
    `Basics.generate("single-electron spectrum: STO", N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0,
                     beta_0::Float64=1.1)` 
        ... to generate a complete one-electron spectrum with N positive and N negative states, and by using even-tempered Slater-type 
            orbitals (STO) with parameters ``\alpha_i = \alpha_0 \beta_0^i``; a spectrum::SingleElecSpectrum is returned where just 
            N0 positive and N_0 negative are kept for later use.  **Not yet implemented !**

    `Basics.generate("single-electron spectrum: STO, positive", N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0,
                     beta_0::Float64=1.1)`  
        ... to generate the same but to return only the N_0 positive states.  **Not yet implemented !**
    """
    function Basics.generate(sa::String, N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0, beta_0::Float64=1.1)
        !(sa == "single-electron spectrum: STO"  &&  sa == "single-electron spectrum: STO, positive")   &&   error("Unsupported keystring = $sa")
        error("Not yet implemented !")

        #= * define N normalized states in stoplus::Array{Vector{Float64},1} and stominus::Array{Vector{Float64},1} 
            on the given grid and with given parameters
        * calculate one-electron (NxN) Hamiltonian matrix with potential h_pq and overlap S_pq, p,q = 1..N
        * solve generalized eigenvalue problem h c_vec = epsilon S c_vec
        * set spectrum::SingleElecSpectrum with just N0 functions
        * check normalization and orthogonality of the functions   =#

        return( nothing )  
    end



    """
    `Basics.generateBasis(confList::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1})`  
        ... generates the CSF basis for the given configuration list and level symmetries; a basis::Basis is returned but 
            without a valid representation of the radial orbitals.
    """
    function Basics.generateBasis(confList::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1})
        #
        relconfList = ConfigurationR[]
        for  conf in confList
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=true)

        # Generate the relativistic CSF's for the given subshell list
        csfList = CsfR[]
        for  relconf in relconfList
            newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
            for  csf in newCsfs     if  LevelSymmetry(csf.J, csf.parity) in symmetries   push!(csfList, csf)   end   end
        end
        #
        if length(csfList) == 0     error("There are no CSF with $symmetries in the given configuration list.")      end

        # Determine the number of electrons and the list of coreSubshells
        NoElectrons      = sum( csfList[1].occupation )
        coreSubshellList = Subshell[]
        for  k in 1:length(subshellList)
            mocc = Basics.subshell_2j(subshellList[k]) + 1;    is_filled = true
            for  csf in csfList
                if  csf.occupation[k] != mocc    is_filled = false;    break   end
            end
            if   is_filled    push!( coreSubshellList, subshellList[k])    end
        end
        
        basis = Basis(true, NoElectrons, subshellList, csfList, coreSubshellList, Dict{Subshell, Orbital}())
        
        println("Construct a basis with $(length(basis.csfs)) CSF for J^P = $symmetries with $(length(basis.subshells)) subshells: " *
                "$(basis.subshells[1])  $(basis.subshells[2]) ...  $(basis.subshells[end-1])  $(basis.subshells[end])")
        return( basis )
    end


    """
    `Basics.generateBasis(refConfigs::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1}, step::AtomicState.RasStep)`  
        ... generates the CSF basis for the given reference configurations and single, double, ... excitations of
            electrons from the corresponding fromShells --> toShells. A basis::Basis is returned but without a valid
            representation of the radial orbitals
    """
    function Basics.generateBasis(refConfigs::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1}, step::AtomicState.RasStep)
        # Single excitations
        if      step.seFrom == Shell[]  ||  step.seTo == Shell[]    confSingles = Configuration[]
        else    confSingles = Basics.generateConfigurations(refConfigs, step.seFrom, step.seTo)
        end
        # Double excitations
        if      step.deFrom == Shell[]  ||  step.deTo == Shell[]    confDoubles = Configuration[]
        else    confDoubles = Basics.generateConfigurations(refConfigs, step.deFrom, step.deTo)
                confDoubles = Basics.generateConfigurations(confDoubles,     step.deFrom, step.deTo)
        end
        # Triple excitations
        if      step.teFrom == Shell[]  ||  step.teTo == Shell[]    confTriples = Configuration[]
        else    confTriples = Basics.generateConfigurations(refConfigs,      step.teFrom, step.teTo)
                confTriples = Basics.generateConfigurations(confTriples,     step.teFrom, step.teTo)
                confTriples = Basics.generateConfigurations(confTriples,     step.teFrom, step.teTo)
        end
        # Quadruple excitations
        if      step.teFrom == Shell[]  ||  step.teTo == Shell[]    confQuadruples = Configuration[]
        else    confQuadruples = Basics.generateConfigurations(refConfigs, step.qeFrom, step.qeTo)
                confQuadruples = Basics.generateConfigurations(confQuadruples,  step.qeFrom, step.qeTo)
                confQuadruples = Basics.generateConfigurations(confQuadruples,  step.qeFrom, step.qeTo)
                confQuadruples = Basics.generateConfigurations(confQuadruples,  step.qeFrom, step.qeTo)
        end
        
        # Now get a unique set of configurations and generate the relativistic configurations
        # configurations = Base.unique((==), configurations) does not work here
        confList       = Configuration[]  
        configurations = deepcopy(refConfigs);      append!(configurations, confSingles);    append!(configurations, confDoubles)
                                                    append!(configurations, confTriples);    append!(configurations, confQuadruples)
        for  confa in configurations
            addTo = true
            for confb in confList    if   confa == confb    addTo = false;    break     end     end
            if  addTo    push!(confList, confa)     end
        end
        for  i = 1:length(confList)    println(">> include ", confList[i])    end
        #
        relconfList = ConfigurationR[]
        for  conf in confList
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=true)

        # Generate the relativistic CSF's for the given subshell list
        csfList = CsfR[]
        for  relconf in relconfList
            newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
            for  csf in newCsfs     if  LevelSymmetry(csf.J, csf.parity) in symmetries   push!(csfList, csf)   end   end
        end

        # Determine the number of electrons and the list of coreSubshells
        NoElectrons      = sum( csfList[1].occupation )
        coreSubshellList = Subshell[]
        for  k in 1:length(subshellList)
            mocc = Basics.subshell_2j(subshellList[k]) + 1;    is_filled = true
            for  csf in csfList
                if  csf.occupation[k] != mocc    is_filled = false;    break   end
            end
            if   is_filled    push!( coreSubshellList, subshellList[k])    end
        end
        
        basis = Basis(true, NoElectrons, subshellList, csfList, coreSubshellList, Dict{Subshell, Orbital}())
        
        println("Construct a basis with $(length(basis.csfs)) CSF for J^P = $symmetries with $(length(basis.subshells)) subshells: " *
                "$(basis.subshells[1])  $(basis.subshells[2]) ...  $(basis.subshells[end-1])  $(basis.subshells[end])")
        return( basis )
    end



    """
    `Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1})`  
        ... generates all nonrelativistic configurations due to the (single) excitation of an electron fromShells --> toShells
            for all given reference configurations. A confList::Array{Configuration,1} is returned.
    """
    function Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1})
        confList = Configuration[]
        for  config in refConfigs
            for  fromShell in fromShells
                for  toShell in toShells
                    # Cycle through all shells of config to built-up a new configuration with a single excitation if fromShell occurs
                    if  haskey(config.shells, fromShell )
                        newShells = Dict{Shell,Int64}();   addShell = true;   addConfiguration = true
                        for (k,v) in config.shells
                            if        k == fromShell == toShell       newShells = Base.merge( newShells, Dict( k => v))
                                      addShell = false   
                            elseif    k == fromShell    if  v-1 > 0   newShells = Base.merge( newShells, Dict( fromShell => v-1))  end
                            elseif    k == toShell      
                                if    v+1 <= 2*(2*toShell.l + 1)      newShells = Base.merge( newShells, Dict( toShell => v+1))
                                      addShell = false   
                                else  addConfiguration = false
                                end
                            else      newShells = Base.merge( newShells, Dict( k => v))
                            end
                            #
                            if  addShell                              newShells = Base.merge( newShells, Dict( toShell => 1))      end
                        end
                        if  addConfiguration   push!(confList, Configuration(newShells, config.NoElectrons ))          end
                    else    # No configuration to add if fromShell does not occur in the reference configuration
                    end
                end
            end
        end
        
        confList = Base.unique(confList)
        
        return( confList )
    end



    """
    `Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1}, 
                                   noex::Int64; restrictions::Array{AbstractConfigurationRestriction,1}=AbstractConfigurationRestriction[])`  
        ... generates all nonrelativistic configurations with excitation of up to noex electrons fromShells --> toShells
            for all given reference configurations. Moreover, the list of restrictions (if any is given) is finally applied 
            to restrict the configurations to a given set of limitations. It remains the reponsibility of the user to make sure that
            the given restrictions do not contradict each other and are consistent with what is to be achieved. The given set of 
            restrictions can be easily extended if this need arises by the users. A confList::Array{Configuration,1} is returned.
    """
    function Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1}, 
                                           noex::Int64; restrictions::Array{AbstractConfigurationRestriction,1}=AbstractConfigurationRestriction[])
        if  noex > 0
            newConfigs = Basics.generateConfigurations(refConfigs, fromShells, toShells)
            newConfigs = Base.unique(newConfigs)
            for ne = 2:noex    
                newConfigs = Basics.generateConfigurations(newConfigs, fromShells, toShells)    
                newConfigs = Base.unique(newConfigs)
            end
        end
        
        # Now apply in turn all given restrictions, if any, and append if no restriction is violated
        confList = Configuration[]
        for  conf  in  newConfigs
            addConf = true
            for res in restrictions
                if  Basics.isViolated(conf,res)     addConf = false;    break   end
            end
            if  addConf     push!(confList, conf)   end
        end
        
        return( confList )
    end


    """
    `Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.NoExcitationScheme, 
                                                      nMax::Int64, lValues::Array{Int64,1})`  
        ... generates a list of non-relativistic configurations for the given (reference) confs and the excitation scheme. 
            All orbitals in standard order are considered up to maxShell = (n_max, l_max).
    """
    function Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.NoExcitationScheme, 
                                                            nMax::Int64, lValues::Array{Int64,1})
        @warn(  "No excitations are included if exScheme::NoExcitationScheme. ")
        println("No excitations are included if exScheme::NoExcitationScheme. ")
        #
        newConfList = confs
        return( newConfList )
    end



    """
    `Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.DeExciteSingleElectron, 
                                                      nMax::Int64, lValues::Array{Int64,1})`  
        ... generates a list of non-relativistic configurations for the given (reference) confs and the excitation scheme. 
            All orbitals in standard order are considered up to maxShell = (n_max, l_max).
    """
    function Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.DeExciteSingleElectron, 
                                                              nMax::Int64, lValues::Array{Int64,1})
        confList  = Configuration[];    NoElectrons = confs[1].NoElectrons
        shellList = Basics.generateShellList(confs, nMax, lValues)
        
        # Create all configurations with single excitations from the given list of configurations
        for  conf in confs
            # Take one electron fromShell and add toShell
            for  (fromShell,occ)  in  conf.shells
                for  toShell  in  shellList
                    newShells = deepcopy( conf.shells )
                    if      fromShell == toShell    
                    elseif  haskey(conf.shells, toShell )   fromOcc = occ;   toOcc = conf.shells[toShell]
                            if  fromOcc - 1 < 0                        println("..");    continue    end
                            if  toOcc   + 1 > 2*(2*toShell.l + 1)      println("..");    continue    end
                            newShells[fromShell] = newShells[fromShell] - 1
                            newShells[toShell]   = newShells[toShell] + 1
                    else    fromOcc = occ
                            if  fromOcc - 1 < 0                        println("..");    continue    end
                            newShells[fromShell] = newShells[fromShell] - 1
                            newShells = Base.merge( newShells, Dict( toShell => 1))
                    end
                    # Add a new configuration
                    push!( confList, Configuration( newShells, NoElectrons))
                    if  true  println(">> Generate $(Configuration( newShells, NoElectrons))")   end
                end
            end
        end

        nbefore     = length(confList)
        newConfList = Configuration[]
        for  conf  in  confList
            if  conf in newConfList    continue;    else    push!( newConfList,  conf)    end
        end
        nafter      = length(newConfList)
        println(">> Number of generated configurations for $exScheme is: $nbefore (before) and $nafter (after).")

        return( newConfList )
    end



    """
    `Basics.generateConfigurationsWithElectronCapture(confs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1}, noex::Int64)`  
        ... generates a list of non-relativistic configurations for the given (reference) confs and with one additional (cpatured) 
            electron. All (doubly) excited configurations with upto :NoExcitations displacements of electrons fromShells into toShells 
            and 'one' additional electron in the toShells are taken into account. The may result in large configuration lists even for 
            a moderate number of fromShell and/or toShells.
    """
    function Basics.generateConfigurationsWithElectronCapture(confs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1},
                                                              noex::Int64)
        newConfList = Configuration[];     NoElectrons = confs[1].NoElectrons + 1
        confList    = Basics.generateConfigurations(confs, fromShells, toShells, noex)
        # Now add one (captured) electron from the toShells to all configurations in confList
        for  conf in confList
            # Take one electron toShells and `add' it to conf
            for  toShell  in  toShells
                newShells = deepcopy( conf.shells )
                if      haskey(conf.shells, toShell )  &&  conf.shells[toShell]  + 1 > 2*(2*toShell.l + 1)     continue    
                elseif  haskey(conf.shells, toShell )  newShells[toShell] = newShells[toShell] + 1
                        push!( newConfList, Configuration( newShells, NoElectrons))
                else    newShells = Base.merge( newShells, Dict( toShell => 1))
                        push!( newConfList, Configuration( newShells, NoElectrons))
                end
                if  false  println(">> Generate $(Configuration( newShells, NoElectrons)) with electron capture.")   end
            end
        end
        newConfList = unique(newConfList)
        
        return( newConfList )
    end


    """
    `Basics.generateLevelWithExtraElectron(newOrbital, symt::LevelSymmetry, level::Level)`  
        ... generates a (new) level with one extra electron in subshell newOrbital.subshell and with overall symmetry symt. The function 
            assumes that all CSF in the basis of level have the same symmetry as the level itself, that the new subshell is not yet part 
            of the basis and that CSF with the total symmetry can be constructed; it terminates with an error if one of this assumptions
            is violated. A newLevel::Level with the same number of CSF and the same representation (eigenvector) as the given level is 
            returned. From a physics viewpoint, newLevel refers to a antisymmetrized product state of [level x orbital(sh, energy)] J^P 
            (symt) and with an total energy = level.energy + newOrbital.energy.
    """
    function Basics.generateLevelWithExtraElectron(newOrbital, symt::LevelSymmetry, level::Level)
        basis = level.basis;    newSubshells = copy(basis.subshells);    newCsfs = CsfR[];   J = level.J;   parity = level.parity
        
        push!(newSubshells, newOrbital.subshell)
        Defaults.setDefaults("relativistic subshell list", newSubshells; printout=false)
        newOrbitals = Base.merge( copy(basis.orbitals), Dict( newOrbital.subshell => newOrbital ))
    
        for  i = 1:length(basis.csfs)
            stateList   = ManyElectron.provideSubshellStates(newOrbital.subshell, 1)
            if  basis.csfs[i].J != J   ||   basis.csfs[i].parity != parity    error("Improper symmetry of CSF.")             end
            if  length(stateList) != 1                                        error("Improper number of subshell states.")   end
            if  AngularMomentum.triangularDelta(J, Basics.subshell_j(newOrbital.subshell), symt.J) != 1    
                                                                            error("Improper coupling of (new) subshell & total J.")   end
            substate   = stateList[1]
            occupation = copy(basis.csfs[i].occupation);    push!(occupation, substate.occ)
            seniorityNr  = copy(basis.csfs[i].seniorityNr);     push!(seniorityNr,  substate.nu)
            subshellJ  = copy(basis.csfs[i].subshellJ);     push!(subshellJ,  AngularJ64(substate.Jsub2//2) )
            subshellX  = copy(basis.csfs[i].subshellX);     push!(subshellX,  symt.J)
            push!(newCsfs, CsfR(true, symt.J, symt.parity, occupation, seniorityNr, subshellJ, subshellX, Subshell[] ) )
        end

        newBasis = Basis(true, basis.NoElectrons+1, newSubshells, newCsfs, copy(basis.coreSubshells), newOrbitals)
        newLevel = Level(symt.J, AngularM64(symt.J), symt.parity, level.index, level.energy+newOrbital.energy, level.relativeOcc, 
                        true, newBasis, level.mc)
        return( newLevel )
    end


    """
    `Basics.generateLevelWithExtraSubshell(sh::Subshell, level::Level)`  
        ... generates a (new) level with one extra subshell sh but with the same overall symmetry as before. The function assumes that the 
            new subshell is not yet part of the basis; it terminates with an error if one of this assumptions is violated. 
            A newLevel::Level with the same symmetry, number of CSF, energy and the same representation (eigenvector) as the given 
            level is returned. 
    """
    function Basics.generateLevelWithExtraSubshell(sh::Subshell, level::Level)
        basis = level.basis;    newSubshells = copy(basis.subshells);    newCsfs = CsfR[];   J = level.J;   parity = level.parity
        newOrbitals = copy(basis.orbitals)

        push!(newSubshells, sh)
        newOrbitals = Base.merge( newOrbitals, Dict( sh => Orbital(sh, -10000.) ))
        Defaults.setDefaults("relativistic subshell list", newSubshells; printout=false)
    
        for  i = 1:length(basis.csfs)
            stateList   = ManyElectron.provideSubshellStates(sh, 0)
            if  length(stateList) != 1                                        error("Improper number of subshell states.")   end
            substate   = stateList[1]
            occupation = copy(basis.csfs[i].occupation);    push!(occupation, substate.occ)
            seniorityNr  = copy(basis.csfs[i].seniorityNr);     push!(seniorityNr,  substate.nu)
            subshellJ  = copy(basis.csfs[i].subshellJ);     push!(subshellJ,  AngularJ64(substate.Jsub2//2) )
            subshellX  = copy(basis.csfs[i].subshellX);     push!(subshellX,  subshellX[end])
            push!(newCsfs, CsfR(true, basis.csfs[i].J, basis.csfs[i].parity, occupation, seniorityNr, subshellJ, subshellX, Subshell[] ) )
        end

        newBasis = Basis(true, basis.NoElectrons, newSubshells, newCsfs, copy(basis.coreSubshells), newOrbitals)
        newLevel = Level(level.J, level.M, level.parity, level.index, level.energy, level.relativeOcc, true, newBasis, level.mc)
        return( newLevel )
    end


    """
    `Basics.generateLevelWithSymmetryReducedBasis(level::Level, subshells::Array{Subshell,1})`  
        ... generates a level with a new basis and representation that only includes the CSF with the same symmetry as the level 
            itself. It also ensures that the new basis is based on the given subshells so that it can be later used with levels
            which have one or several subshells more (at it often appears in ionization processes). Both, the underlying CSF basis 
            and the eigenvectors are adopted accordingly. The procedures only supports the 'addition' of subshells but it stops 
            if (1) the individual sequence differs or if the basis of level contains more subshells than the given subshell-list.
            A (new) level::levelR is returned.
    """
    function Basics.generateLevelWithSymmetryReducedBasis(level::Level, subshells::Array{Subshell,1})
        testBasis = level.basis;    newCsfs = CsfR[];    newMc = Float64[]

        # Decide of whether the basis need to be extended by some subshells
        nt = length(testBasis.subshells)
        if      testBasis.subshells == subshells      basis = testBasis;    nLevel = level
        elseif  nt >= length(subshells)               error("stop a")
        elseif  nt + 1 == length(subshells)  &&  testBasis.subshells == subshells[1:nt]
                nLevel = Basics.generateLevelWithExtraSubshell(subshells[end], level)
                basis  = nLevel.basis
        else
            nLevel  = level
            for  ns = nt+1:length(subshells)
                nLevel = Basics.generateLevelWithExtraSubshell(subshells[ns], nLevel)
            end
            basis  = nLevel.basis
        end
            
        
        for  i = 1:length(basis.csfs)
            if  basis.csfs[i].J == nLevel.J  &&   basis.csfs[i].parity == nLevel.parity        &&    abs(nLevel.mc[i]) > 1.0e-7
                push!(newCsfs, deepcopy(basis.csfs[i]));    push!(newMc, deepcopy(nLevel.mc[i]))
            end
        end
        newBasis = Basis(true, basis.NoElectrons, deepcopy(basis.subshells), newCsfs, deepcopy(basis.coreSubshells), deepcopy(basis.orbitals))
        newLevel = Level(nLevel.J, nLevel.M, nLevel.parity, nLevel.index, nLevel.energy, nLevel.relativeOcc, true, newBasis, newMc)
        return( newLevel )
    end
    

    """
    `Basics.generateOrbitalsForPotential(grid::Radial.Grid, meanPot::Radial.Potential, subshellList::Array{Subshell,1})`  
        ... generates a set of (start) orbitals from the given potential and for all the subshells in subshellList. 
            A set of orbitals::Dict{Subshell, Orbital} is returned.
    """
    function  Basics.generateOrbitalsForPotential(grid::Radial.Grid, meanPot::Radial.Potential, subshellList::Array{Subshell,1})
        orbitals = Dict{Subshell, Orbital}()
        
        # Determine kappaMin and kappaMax
        kappaMin = 0;   kappaMax = 0
        for  subsh in subshellList   
            if       subsh.kappa < kappaMin    kappaMin = subsh.kappa
            elseif   subsh.kappa > kappaMax    kappaMax = subsh.kappa
            end
        end
        
        # Generate the primitives for a B-spline basis
        wa = Bsplines.generatePrimitives(grid)
        
        # Now cycle through all kappa symmetries
        for  kappa = kappaMin:kappaMax
            # Determine all requested subshells of symmetry kappa ... to compute them together
            shList     = Subshell[];    for  subsh in subshellList    if kappa == subsh.kappa    push!( shList, subsh)    end    end
            shOrbitals = Bsplines.generateOrbitalsForPotential(wa, kappa, shList, meanPot; printout=false)
            for  sh in shList
                orbitals = Base.merge( orbitals, Dict( sh => shOrbitals[sh]))
            end
        end

        return( orbitals )
    end



    """
    `Basics.generateOrbitalsForBasis(basis::Basis, frozenShells::Array{Shell,1}, priorBasis::Basis, startOrbitals::Dict{Subshell, Orbital})`  
        ... generates a dict of (relativistic) orbitals as specificed by subshells of the basis. These orbitals are taken from priorBasis
            if contained in frozen-shells, and from spectrum (start orbitals) otherwises. An error message is issued if some
            requested orbital is not found. A orbList::Dict{Subshell, Orbital} is returned.
    """
    function Basics.generateOrbitalsForBasis(basis::Basis, frozenShells::Array{Shell,1}, priorBasis::Basis, 
                                             startOrbitals::Dict{Subshell, Orbital})
        orbitals = Dict{Subshell, Orbital}()
        for  subsh in basis.subshells
            if  subsh in frozenShells    &&   haskey(priorBasis.orbitals, subsh)
                orbitals = Base.merge( orbitals, Dict( subsh => priorBasis.orbitals[subsh]) )
            elseif  subsh in frozenShells     error("Frozen orbital $subsh not found in prior basis.")
            elseif  haskey(priorBasis.orbitals, subsh)
                println(">> Start orbital $subsh is taken from prior basis")
                orbitals = Base.merge( orbitals, Dict( subsh => priorBasis.orbitals[subsh]) )
            else
                println(">> Start orbital $subsh is taken from hydrogenic orbitals")
                orbitals = Base.merge( orbitals, Dict( subsh => startOrbitals[subsh]) )
            end
        end
        return( orbitals )
    end



    """
    `Basics.generateOrbitalSuperposition(a::Orbital, b::Orbital, cx::Float64, grid::Radial.Grid)`  
        ... generates a superposition  a + cx * b of two given orbitals; the function just takes the linear combination
            of the large and small components as well as the energy. The function assumes that both orbitals are defined on the same
            grid. A re-normalized newOrbital::Orbital is returned for which the derivatives are not defined. 
            This function is used to accelerate the convergence (hopefully).
    """
    function Basics.generateOrbitalSuperposition(a::Orbital, b::Orbital, cx::Float64, grid::Radial.Grid)
        if  a.subshell != b.subshell  ||  !(a.isBound)  ||  !(b.isBound)   error("stop a")     end
        mtp  = max(size(a.P, 1), size(b.P, 1));     newP = zeros(mtp);  newQ = zeros(mtp);   newEnergy = (a.energy + cx * b.energy) / (1 + cx)
        mtpa = size(a.P, 1);    newP[1:mtpa] = a.P[1:mtpa];    newQ[1:mtpa] = a.Q[1:mtpa]
        mtpb = size(b.P, 1);    newP[1:mtpb] = newP[1:mtpb] + cx * b.P[1:mtpb]    
                                newQ[1:mtpb] = newQ[1:mtpb] + cx * b.Q[1:mtpb]
        newOrbital = Orbital(a.subshell, a.isBound, a.useStandardGrid, newEnergy, newP, newQ, Float64[], Float64[], a.grid)
        norm       = RadialIntegrals.overlap(newOrbital, newOrbital, grid)
        newP       = newP / sqrt(norm);      newQ = newQ / sqrt(norm)
        newOrbital = Orbital(a.subshell, a.isBound, a.useStandardGrid, newEnergy, newP, newQ, Float64[], Float64[], a.grid)

        return( newOrbital )
    end

    

    """
    `Basics.generateShellList(confs::Array{Configuration,1}, nMax::Int64, lValues::Array{Int64,1})`  
        ... generates a list of (non-relativistic) shells in standard order that contains all shells from the given list of
            configurations as well as the shells up to the principal quantum number nMax and the list of orbital quantum
            numbers lValues. A shellList::Array{Shell,1} is returned.
    """
    function Basics.generateShellList(confs::Array{Configuration,1}, nMax::Int64, lValues::Array{Int64,1})
        shellList = Shell[]
        # Add all shells from the given configurations
        for  conf in confs
            for (k,v) in conf.shells    push!( shellList, k)    end
        end
        # Add all shells for the given quantum numbers
        for  n = 1:nMax
            for  ll in  lValues  
                if  n >= ll + 1     push!( shellList, Shell(n,ll))      end    
            end
        end
        
        # Now bring the shells in standard order
        newShellList = Shell[]
        for  n = 1:nMax
            for  l = 0:30  
                if  n >= l + 1  && Shell(n,l) in shellList    push!( newShellList, Shell(n,l))     end    
            end
        end
        println(">> From configurations generated shell list $newShellList ")
        
        return( newShellList )
    end

    
