# How to obtain the radial probability distribution function of a given orbital from a quantum chemical calculation?

I know that a very similar question has already been asked here and a very competent and useful answer has been provided (see here). However, that previous question was related to the total electron density whereas I am interested in plotting the radial probability distribution function of a given single orbital. Can I also use MultiWFN tool and Gaussian output to do this?

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

>>

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