How to obtain the radial probability distribution function from a quantum chemical calculation?

In the book "Quantum Chemistry" by Ira Levine I found a plot of the radial probability distribution function of argon:

The figure appears in the context of the following section:

Electron densities calculated from Hartree–Fock wave functions are quite accurate. Figure 11.1 compares the radial distribution function of argon (found by integrating the electron density over the angles θ and φ and multiplying the result by r2) calculated by the Hartree–Fock method with the experimental radial distribution function found by electron diffraction. (Recall from Section 6.6 that the radial distribution function is proportional to the probability of finding an electron in a thin spherical shell at a distance r from the nucleus.) Note the electronic shell structure in Fig. 11.1. The high nuclear charge in 18Ar makes the average distance of the 1s electrons from the nucleus far less than in H or He. Thus there is only a moderate increase in atomic size as we go down a given group in the periodic table. Calculations show that the radius of a sphere containing 98% of the Hartree–Fock electron probability density gives an atomic radius in good agreement with the empirically determined van der Waals radius. [See C. W. Kammeyer and D. R. Whitman, J. Chem. Phys., 56, 4419 (1972).]

Source: Ira N. Levine. Quantum Chemistry 7th ed.; Pearson, 2014; page 294.

What kind of chemical software allows one to do this?

It is possible to do with Gaussian, and, of course, how is it possible?

• The electron density can be calculated with maybe all the common quantum chemical software, however it is generally exported as a bunch of value on a grid (cube files), and this histogram is not directly calculated. To make this particular plot you most probably have to do the calculation yourself from the the calculated electron distribution.
– Greg
Mar 8, 2017 at 18:26
• @Greg Ok, I get .cube file for Argon like this rgho.st/8MSlxwyMZ How can I process this file for relevant data? Mar 8, 2017 at 19:43

This is gladly not a big problem as MultiWFN (open source) can do that for you. I'll guide you through the process from start to finish. My current setup is Gaussian 09 Rev D.01 and MultiWFN 3.3.8; there are newer versions out there with which this should work, too.

Run the calculation, for example ar.com at the HF/STO-3G level of theory, but anything else will also be fine:

%chk=ar.chk
#P HF/STO-3G
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
pop=full pop=nbo6

title card

0 1
Ar    0.0  0.0  0.0



This will generate a checkpoint file ar.chk which needs to be formatted:

formchk -3 ar.chk ar.fchk


Now fire up your MultiWFN program, you should get a prompt similar to the following:

 Multiwfn -- A Multifunctional Wavefunction Analyzer (for Linux)
Version 3.3.8, release date: 2015-Dec-1
Project leader: Tian Lu  (Beijing Kein Research Center for Natural Sciences)
Citation of Multiwfn: Tian Lu, Feiwu Chen, J. Comp. Chem. 33, 580-592 (2012)
Multiwfn official website: http://Multiwfn.codeplex.com
Multiwfn official forum (in Chinese): http://bbs.keinsci.com
( The number of threads:  2   Current date: 2017-03-09   Time: 14:46:57 )

Input file path, for example E:\Haiyore!Nyaruko-san\Nyaruko.wfn
(Supported types: .wfn/.wfx/.fch/.31/.chg/.pdb/.xyz/.cub/.grd/.molden, etc.)
Hint: To reload the file last time used, simply input the letter "o". Input such as ?miku.fch can open miku.fch in the same folder of the file last time used


If you loaded MultiWFN from the directory you're files are located, simply put in ar.fchk. It will load the file and output some basic information and take you to the "Main function menu":

 Please wait...
Converting basis function information to GTF information...
Generating overlap matrix...
Generating SCF density matrix...

Total/Alpha/Beta electrons:     18.0000      9.0000      9.0000
Net charge:     0.00000      Expected multiplicity:    1
Atoms:      1,  Basis functions:      9,  GTFs:     27
Total energy:  -521.222880802727 Hartree,   Virial ratio:  2.01704078
This is a restricted single-determinant wavefunction
Orbitals from 1 to     9 are occupied
Title line of this file: title card
Formula: Ar1
Molecule weight:        39.94800



For the sake of sanity I will not copy every menu screen here, just simply what you have to put in. On this screen the radial distribution function is hidden in the "Other functions (Part 2)" and there it is option 5, also see the manual chapter 3.200.5 on page 175 (v.3.3.8). So choose

200 // Other functions (2)


and then

5 // RDF


You should get the following menu:

   ====== Plot radial distribution function for a real space function ======
-1 Exit
0 Calculate radial distribution function and its integration curve
1 Select real space function, current:   1
2 Set sphere center, current    0.0000    0.0000    0.0000 Ang
3 Set lower and upper limit of radial distance, current:   0.0000   5.0000 Ang
4 Set the number of integration point in each shell, current:  2030
5 Set the number of radial points, current:   500


Here you can select a couple of options. The electron density is already preselected, but you can change it in option 1. You'll get a table with all available functions.
The only thing you might want to adjust are the limits, so choose

3 // Limits
0.0,2.0


and calculate the RDF:

0 // Calculate RDF


Choose 1 to plot and view the function. Or choose

3 // Save graph


or choose

5 // Export to RDF.txt


for post-procesing (with gnuplot).

You can shorten the journey if you prepare a file like this mwfn.in

ar.fchk
200 // Other functions p2
1   // Select real space fucntion
1   // Electron density
2   // Set sphere centre
0.0,0.0,0.0
3   // Set limits
0.0,2.0
0   // Calculate RDF
3   // Save graph
5   // export to text file


and pipe it into the program:

multiwfn < mwfn.in


The native graphing capabilities of MultiWFN don't allow for customisation, so you might want to run it through gnuplot. Here is a file which gives you the graph below, I am using gnuplot 4.4, but later versions are more feature rich.

# General Settings
set terminal svg enhanced font "Arial,12" size 800,600
set tics font "Arial, 9"
set output "ar-rdf.svg"
unset key
#or set key top right Left # etc
set title 'Argon Radial Distribution Function HF/STO-3G'
#
set xlabel 'r / Å'
set xtics in scale 0.6, 0.4 nomirror
set xtics 0,1
set mxtics 2
set xrange [0:2]
# Y-Axis: RDF/a.u.
set ylabel "RDF / a.u." rotate left
set ytics  in  scale 0.6, 0.4 nomirror
set ytics 0,10
set mytics 2
set yrange [0:25]
#
#
# PLOT!
plot \
"RDF.txt"   using ($2):($3) title "HF/STO-3G"   with lines lt 1 axes x1y1


This is the resulting graph:

• This is incredibly useful. Mar 9, 2017 at 18:40
• Tested with Gaussian 16 B.01 and Multiwfn 3.6.0 and it worked the same. Same for gnuplot 4.6 patchlevel 6. Jan 9, 2020 at 17:45
• hi, Gaussian is a must here? may I know if I dont have Gaussian, is there any way (open source software) combined with multiwfn to get the RDF ? @Martin-マーチン Oct 5, 2022 at 6:14
• @Jiadong probably not. Multiwfn understands quite a few packages, essentially you'll need a wave function. Please consult their manual, I'm not active in the field anymore, so I can't really help. Oct 5, 2022 at 17:40