How to obtain the radial probability distribution function of a given orbital from a quantum chemical calculation?
This is actually a lot simpler than I initially thought. I'll be using the same example as previously explained in How to obtain the radial probability distribution function from a quantum chemical calculation? My test setup has slightly changed since then, i.e. I've tried it with Gaussian 16 B.01 and Multiwfn 3.6.0 and it worked the same.

We start the same way as in the linked question, until you have your ar.fchk and load that into Multiwfn, you should see something like this:


 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.6, release date: 2019-May-21
 Project leader: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Below paper *MUST BE CITED* if Multiwfn is utilized in your work:
          Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580-592 (2012)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn

>>

 Loaded ar.fchk successfully!
 Formula: Ar1 
 Molecule weight:        39.94800


We need to modify the wave function now:


                    ************ Main function menu ************
 0 Show molecular structure and view orbitals
 1 Output all properties at a point
 2 Topology analysis
 3 Output and plot specific property in a line
 4 Output and plot specific property in a plane
 5 Output and plot specific property within a spatial region (calc. grid data)
 6 Check & modify wavefunction
 7 Population analysis and atomic charges
 8 Orbital composition analysis
 9 Bond order analysis
 10 Plot total DOS, partial DOS, OPDOS, local DOS and photoelectron spectrum
 11 Plot IR/Raman/UV-Vis/ECD/VCD/ROA spectrum
 12 Quantitative analysis of molecular surface
 13 Process grid data (No grid data is presented currently)
 14 Adaptive natural density partitioning (AdNDP) analysis
 15 Fuzzy atomic space analysis
 16 Charge decomposition analysis (CDA) and extended CDA (ECDA)
 17 Basin analysis
 18 Electron excitation analysis
 19 Orbital localization analysis
 20 Visual study of weak interaction
 21 Energy decomposition analysis
 100 Other functions (Part1)        200 Other functions (Part2)


Chose 6:


           ============ Modify & Check wavefunction ============ 
 Number of GTFs:    27, Orb:     9, Atoms:    1, A/B elec:   9.000   9.000
 -4 Exclude contribution of some atoms to real space functions
 -3 Only retain contribution of some atoms to real space functions
 -1 Return
 0 Save the modified wavefunction to a new .wfn file
 1 List all primitive function      2 List all basis function
 3 List all orbitals                4 Print detail information of an orbital
 5 Print coefficient matrix in basis function
 6 Print density matrix in basis function
 7 Print various kinds of integral matrix between basis functions
 11 Swap some information of two primitive functions
 21 Set center of a primitive       22 Set type of a primitive
 23 Set exponent of a primitive     24 Set coefficient of a primitive
 25 Set coefficients of GTFs/basis functions that satisfied certain conditions
 26 Set occupation number of some orbitals
 27 Set type of some orbitals       28 Set energy of some orbitals
 31 Translate the system            32 Translate and duplicate the system
 33 Rotate wavefunction, namely X->Y, Y->Z, Z->X
 34 Set occupation number of inner orbitals to zero
 36 Invert phase of some orbitals


As you can see there is plenty you can do. For this example I will reset the occupation of some orbitals (26). More precisely, I'll transform the wave function so that only the 2s orbital is doubly occupied. Here is how I did it:


           ============ Modify & Check wavefunction ============ 
 Number of GTFs:    27, Orb:     9, Atoms:    1, A/B elec:   9.000   9.000
 -4 Exclude contribution of some atoms to real space functions
 -3 Only retain contribution of some atoms to real space functions
 -1 Return
 0 Save the modified wavefunction to a new .wfn file
 1 List all primitive function      2 List all basis function
 3 List all orbitals                4 Print detail information of an orbital
 5 Print coefficient matrix in basis function
 6 Print density matrix in basis function
 7 Print various kinds of integral matrix between basis functions
 11 Swap some information of two primitive functions
 21 Set center of a primitive       22 Set type of a primitive
 23 Set exponent of a primitive     24 Set coefficient of a primitive
 25 Set coefficients of GTFs/basis functions that satisfied certain conditions
 26 Set occupation number of some orbitals
 27 Set type of some orbitals       28 Set energy of some orbitals
 31 Translate the system            32 Translate and duplicate the system
 33 Rotate wavefunction, namely X->Y, Y->Z, Z->X
 34 Set occupation number of inner orbitals to zero
 36 Invert phase of some orbitals
26

 Select the orbitals for which the occupation numbers are needed to be changed
 e.g. 2,4,13-16,20 means selecting orbitals 2,4,13,14,15,16,20
 Input 0 can select all orbitals, input q or 00 can return
0
 Set to which value? e.g. 1.2
 Note:
 You can also input for example "+1.1" "-1.1" "*1.1" "/1.1" to add, minus, multiply and divide the occupation numbers by 1.1
 To recover the initial occupation numbers, input "i"
 To generate occupation state for calculating odd electron density, input "odd"
0
 Done!

 Select the orbitals for which the occupation numbers are needed to be changed
 e.g. 2,4,13-16,20 means selecting orbitals 2,4,13,14,15,16,20
 Input 0 can select all orbitals, input q or 00 can return
2
 Set to which value? e.g. 1.2
 Note:
 You can also input for example "+1.1" "-1.1" "*1.1" "/1.1" to add, minus, multiply and divide the occupation numbers by 1.1
 To recover the initial occupation numbers, input "i"
 To generate occupation state for calculating odd electron density, input "odd"
2
 Done!

 Select the orbitals for which the occupation numbers are needed to be changed
 e.g. 2,4,13-16,20 means selecting orbitals 2,4,13,14,15,16,20
 Input 0 can select all orbitals, input q or 00 can return
q

           ============ Modify & Check wavefunction ============ 
 Number of GTFs:    27, Orb:     9, Atoms:    1, A/B elec:   1.000   1.000
 -4 Exclude contribution of some atoms to real space functions
 -3 Only retain contribution of some atoms to real space functions
 -1 Return
 0 Save the modified wavefunction to a new .wfn file
 1 List all primitive function      2 List all basis function
 3 List all orbitals                4 Print detail information of an orbital
 5 Print coefficient matrix in basis function
 6 Print density matrix in basis function
 7 Print various kinds of integral matrix between basis functions
 11 Swap some information of two primitive functions
 21 Set center of a primitive       22 Set type of a primitive
 23 Set exponent of a primitive     24 Set coefficient of a primitive
 25 Set coefficients of GTFs/basis functions that satisfied certain conditions
 26 Set occupation number of some orbitals
 27 Set type of some orbitals       28 Set energy of some orbitals
 31 Translate the system            32 Translate and duplicate the system
 33 Rotate wavefunction, namely X->Y, Y->Z, Z->X
 36 Invert phase of some orbitals
3

 Basic information of all orbitals:
 Orb:     1 Ene(au/eV):  -117.263218   -3190.8944 Occ: 0.000000 Type:A+B       
 Orb:     2 Ene(au/eV):   -12.054866    -328.0296 Occ: 2.000000 Type:A+B       
 Orb:     3 Ene(au/eV):    -9.308431    -253.2953 Occ: 0.000000 Type:A+B       
 Orb:     4 Ene(au/eV):    -9.308431    -253.2953 Occ: 0.000000 Type:A+B       
 Orb:     5 Ene(au/eV):    -9.308431    -253.2953 Occ: 0.000000 Type:A+B       
 Orb:     6 Ene(au/eV):    -1.154432     -31.4137 Occ: 0.000000 Type:A+B       
 Orb:     7 Ene(au/eV):    -0.495941     -13.4953 Occ: 0.000000 Type:A+B       
 Orb:     8 Ene(au/eV):    -0.495941     -13.4953 Occ: 0.000000 Type:A+B       
 Orb:     9 Ene(au/eV):    -0.495941     -13.4953 Occ: 0.000000 Type:A+B       

           ============ Modify & Check wavefunction ============ 
 Number of GTFs:    27, Orb:     9, Atoms:    1, A/B elec:   1.000   1.000
 -4 Exclude contribution of some atoms to real space functions
 -3 Only retain contribution of some atoms to real space functions
 -1 Return
 0 Save the modified wavefunction to a new .wfn file
 1 List all primitive function      2 List all basis function
 3 List all orbitals                4 Print detail information of an orbital
 5 Print coefficient matrix in basis function
 6 Print density matrix in basis function
 7 Print various kinds of integral matrix between basis functions
 11 Swap some information of two primitive functions
 21 Set center of a primitive       22 Set type of a primitive
 23 Set exponent of a primitive     24 Set coefficient of a primitive
 25 Set coefficients of GTFs/basis functions that satisfied certain conditions
 26 Set occupation number of some orbitals
 27 Set type of some orbitals       28 Set energy of some orbitals
 31 Translate the system            32 Translate and duplicate the system
 33 Rotate wavefunction, namely X->Y, Y->Z, Z->X
 36 Invert phase of some orbitals


(The values I have typed in all start in the first column.)

Now leave this screen. You may want to safe with 0 the resulting wave function as new.wfn; but it is not necessary. To get back to the main menu simply type -1.

Now that you have modified the wave function, you can again follow the steps from the linked post.

TL;DR: 


200 // other functions
5   // RDF
0   // calc
1   // plot
5   // Export to RDF.txt


I'm sure there  is a lot of fancy stuff you can do with that, I used the same gnuplot file to plot the result:





Shameless self-plug:

I have written some startup scripts for Multiwfn, which make some things easier or extend the functionality to run calculations within a queueing system (e.g. QTAIM charges). Find it on GitHub: https://github.com/polyluxus/runMultiwfn.bash