# How to calculate ionic radius in Gaussian?

I am working with lithium sulfate and am wondering how to figure out the ionic radius of the Li atoms. I did a DFT/B3LYP/G-311G ++ (d,p) calculation.

• Ionic radius is not a well-defined quantum-mechanical quantity, so estimating it properly in an unbiased can be quite involved. Jul 10, 2023 at 21:44

## Does it really need to be done in Gaussian?

I found a paper(1) titled "Atomic and Ionic Radii of Elements 1-96" where they calculated it for a large fraction of the periodic table:

Abstract: Atomic and cationic radii have been calculated for the first 96 elements, together with selected anionic radii. The metric adopted is the average distance from the nucleus where the electron density falls to 0.001 electrons per bohr3, following earlier work by Boyd. Our radii are derived using relativistic all-electron density functional theory calculations, close to the basis set limit. They offer a systematic quantitative measure of the sizes of non-interacting atoms, commonly invoked in the rationalization of chemical bonding, structure, and different properties. Remarkably, the atomic radii as defined in this way correlate well with van der Waals radii derived from crystal structures. A rationalization for trends and exceptions in those correlations is provided.

To do it, they used the non-empirical hybrid-exchange correlation functional PBE0, together with a very large and uncontracted atomic natural orbital-relativistic correlation consistent(ANO-RCC) basis set and the Douglas–Kroll–Hess second-order scalar relativistic Hamiltonian. Since I dropped out of my master degree back in the pandemic I don't have access to a Gaussian anymore, but I posted here previously a question regarding volumes of surfaces with uniform electronic density with some examples. Now I adapted the examples to try reproduce these calculations of Rahm et al., at least for Li you are interested in, and also for the next three elements in the alkali metal group, Na, K and Rb.

My attempt was done using NWChem(2), a major free software computational chemistry package(3) available in the repositories from many linux distributions. The particular version I used was the one currently available from Fedora 38 repositories:

[user@fedora ~]$$dnf list installed | grep 'nwchem' nwchem.x86_64 7.0.2-12.fc38 @fedora nwchem-common.noarch 7.0.2-12.fc38 @fedora nwchem-openmpi.x86_64 7.0.2-12.fc38 @fedora [user@fedora ~]$$


So I created input files for single point energy calculations, using PBE0 and ANO-RCC, where each ion is centered at x=0, y=0, z=0 coordinates, as shown below. The largest value for radius calculated in the paper we are trying reproduce is a little under 3Å, so the input files also direct NWChem to output density cubes measuring 6Åx6Åx6Å, that should be enough to avoid any part of the volume with electronic density above 0.001 electrons per bohr³ to be left outside the cube. The grid spacings were also chosen so our voxel volume (21.4μÅ³) is close to the one used in the paper (18.5μÅ³):

[user@fedora ionic_radii2]$more *.nw | cat :::::::::::::: input_lithium.nw :::::::::::::: echo start molecule title "Single Point | pbe0/ANO-RCC" charge 1 memory total 2 gb geometry units angstroms print xyz autosym Li 0.000000 0.000000 0.000000 end basis Li library ANO-RCC rel end relativistic douglas-kroll end dft xc pbe0 mult 1 end task dft energy dplot TITLE lithium LimitXYZ -3.0 3.0 216 -3.0 3.0 216 -3.0 3.0 216 spin total gaussian output lithium_density.cube end task dplot :::::::::::::: input.nw :::::::::::::: echo start molecule title "Single Point | pbe0/aug-pcS-1/pcS-0" charge 0 memory total 2 gb geometry units angstroms print xyz autosym C 0.913139 0.913139 0.000000 H 0.676070 1.981751 0.000000 C 0.000000 0.000000 0.780217 C 0.000000 0.000000 -0.780217 H 1.378212 -1.576370 0.000000 C 0.334232 -1.247371 0.000000 H 1.981751 0.676070 0.000000 C -1.247371 0.334232 0.000000 H -0.405381 -2.054282 0.000000 H -2.054282 -0.405381 0.000000 H -1.576370 1.378212 0.000000 end basis H library pcS-0 C library aug-pcS-1 end dft xc pbe0 mult 1 end task dft energy dplot TITLE propellane LimitXYZ -4.0 4.0 10 -4.0 4.0 10 -2.0 2.0 10 spin total gaussian output propel_density.cube end task dplot:::::::::::::: input_potassium.nw :::::::::::::: echo start molecule title "Single Point | pbe0/ANO-RCC" charge 1 memory total 2 gb geometry units angstroms print xyz autosym K 0.000000 0.000000 0.000000 end basis K library ANO-RCC rel end relativistic douglas-kroll end dft xc pbe0 mult 1 end task dft energy dplot TITLE Potassium LimitXYZ -3.0 3.0 216 -3.0 3.0 216 -3.0 3.0 216 spin total gaussian output potassium_density.cube end task dplot :::::::::::::: input_rubidium.nw :::::::::::::: echo start molecule title "Single Point | pbe0/ANO-RCC" charge 1 memory total 2 gb geometry units angstroms print xyz autosym Rb 0.000000 0.000000 0.000000 end basis Rb library ANO-RCC rel end relativistic douglas-kroll end dft xc pbe0 mult 1 end task dft energy dplot TITLE Rubidium LimitXYZ -3.0 3.0 216 -3.0 3.0 216 -3.0 3.0 216 spin total gaussian output rubidium_density.cube end task dplot :::::::::::::: input_sodium.nw :::::::::::::: echo start molecule title "Single Point | pbe0/ANO-RCC" charge 1 memory total 2 gb geometry units angstroms print xyz autosym Na 0.000000 0.000000 0.000000 end basis Na library ANO-RCC rel end relativistic douglas-kroll end dft xc pbe0 mult 1 end task dft energy dplot TITLE Sodium LimitXYZ -3.0 3.0 216 -3.0 3.0 216 -3.0 3.0 216 spin total gaussian output sodium_density.cube end task dplot [user@fedora ionic_radii2]$


This page has instructions on how to install and run NWChem under Fedora. To install it there, you can run this command in the terminal:

dnf install nwchem-openmpi


Once installed, I did a parallel run on each input, using 8 threads on a AMD Ryzen 7 1700 machine, as follows:

module load mpi/openmpi
mpirun -np 8 nwchem_openmpi input_lithium.nw > output_Li.out
mpirun -np 8 nwchem_openmpi input_sodium.nw > output_Na.out
mpirun -np 8 nwchem_openmpi input_potassium.nw > output_K.out
mpirun -np 8 nwchem_openmpi input_rubidium.nw > output_Rb.out


Once finished, I checked the outputs produced no error and the cube files were properly recorded. The runtimes are also shown, by printing the last line of each output file:

[user@fedora ionic_radii2]$$grep -i 'error' *.out [user@fedora ionic_radii2]$$ tail -n 1 *.out
==> output_K.out <==
Total times  cpu:     2009.8s     wall:     2116.0s

==> output_Li.out <==
Total times  cpu:       49.0s     wall:      102.6s

==> output_Na.out <==
Total times  cpu:      394.4s     wall:      474.6s

==> output_Rb.out <==
Total times  cpu:    11640.9s     wall:    11949.9s
[user@fedora ionic_radii2]$$ls -l -s *.cube 131428 -rw-r--r--. 1 user user 134580624 jul 16 12:09 lithium_density.cube 131428 -rw-r--r--. 1 user user 134580626 jul 16 13:42 potassium_density.cube 131428 -rw-r--r--. 1 user user 134580625 jul 16 17:02 rubidium_density.cube 131428 -rw-r--r--. 1 user user 134580623 jul 16 13:04 sodium_density.cube [user@fedora ionic_radii2]$$


The timings show computational cost scales quickly with the size of the system. While the calculation for Li+ takes a bit over 100s wall time, the one for Rb+ takes over 3h. My initial idea was to try reproduce the results for all the alkali metals, but when I saw how quick processing time increased I decided to stop at rubidium.

Now the remaining problem is to do volume integration on the electron densities data present in the cube files. In the original study they used the partitioning method decribed in this paper(4). I sticked to the less sophisticated method I used in my past question, of just counting points above the given threshold of 0.001 electrons per bohr³, and then calculating the ion volume as a ratio of the box volume. In the original study, that covers 96 elements, they raise the issue of most of them except those in groups 1,2 and 18 have non-spherical grids, what complicates estimating average radii. But as all elements we are dealing with are from group 1, we can estimate the radii from the respective volumes by just assuming them spherical. The following Python 3 code implements it (must be saved to a .py file and run in the same folder the cube files are located):

from math import pi
from os import listdir
import re

def main():
cubes = [file for file in listdir() if 'density.cube' in file]
for file in cubes:
with open(file, 'r') as f:

##Matches data in format the cube file uses to store
##density values, like 0.77410E-07 or 0.38588E+01.
p = re.compile(r"0\.\d{5}E[-+]\d{2}")

## From the LimitXYZ section in the input files. If you change the cube volume, this must be updated.
cube_volume = 6*6*6
isosurface_threshold = 0.001

data = p.findall(text)
numeric_data = [float(value) for value in data]
fraction_inside = len([density
for density in numeric_data
if density >= isosurface_threshold]) / len(numeric_data)

bound_volume = cube_volume * fraction_inside
## From the sphere volume formula.
print('For',
file,
'angstroms')

main()


It then outputs:

Python 3.9.16 (15ea147a6c65, May 30 2023, 13:28:11)
[PyPy 7.3.11 with GCC 13.1.1 20230511 (Red Hat 13.1.1-2)] on linux
>>>
For lithium_density.cube the calculated average radius is 0.97 angstroms
For potassium_density.cube the calculated average radius is 1.74 angstroms
For rubidium_density.cube the calculated average radius is 1.91 angstroms
For sodium_density.cube the calculated average radius is 1.32 angstroms
>>>


The results I got were very close to those in the original work:

Ion Results from Rahm et al. (Å) My results (Å)
Li+ 0.98 0.97
Na+ 1.33 1.32
K+ 1.75 1.74
Rb+ 1.91 1.91

On how these results compare with values derived from empirical data, they don't discuss the case for monopositive cations as in depth as is done with the neutral atoms. For the neutral atoms, they plot the difference between their calculated radii against crystallographic radii, as shown bellow. On average their radii are smaller by 0.01 Å, and the standard deviation corresponds to 0.211Å, with most alkali metals acting as outliers:

The determination of radii also depend on other factors, like coordination numbers. Regarding the ionic species, bellow I put some values compiled by Shannon(5) from empirical data, for comparison:

Li+ IV 0.730 0.590
Li+ VI 0.90 0.76
Na+ IV 1.13 0.99
Na+ V 1.14 1.00
Na+ VI 1.16 1.02
Na+ VII 1.26 1.12
Na+ VIII 1.32 1.18
Na+ XII 1.53 1.39
K+ IV 1.51 1.37
K+ VI 1.52 1.38
K+ VII 1.60 1.46
K+ VIII 1.65 1.51
K+ IX 1.69 1.55
K+ X 1.73 1.59
K+ XII 1.78 1.64
Rb+ VI 1.66 1.52
Rb+ VII 1.70 1.56
Rb+ VIII 1.75 1.61
Rb+ X 1.80 1.66
Rb+ XI 1.83 1.69
Rb+ XII 1.86 1.72
Rb+ XIV 1.97 1.83

References:

(1). Rahm, Martin, et al. “Atomic and Ionic Radii of Elements 1 –96.” Chemistry – A European Journal, vol. 22, no. 41, Oct. 2016, pp. 14625–32. DOI.org (Crossref), https://doi.org/10.1002/chem.201602949.

(2). Aprà, E., et al. “NWChem: Past, Present, and Future.” The Journal of Chemical Physics, vol. 152, no. 18, May 2020, p. 184102. DOI.org (Crossref), https://doi.org/10.1063/5.0004997.

(3). Lehtola, Susi, and Antti J. Karttunen. “Free and Open Source Software for Computational Chemistry Education.” WIREs Computational Molecular Science, vol. 12, no. 5, Sept. 2022. DOI.org (Crossref), https://doi.org/10.1002/wcms.1610.

(4). Tang, W., et al. “A Grid-Based Bader Analysis Algorithm without Lattice Bias.” Journal of Physics: Condensed Matter, vol. 21, no. 8, Feb. 2009, p. 084204. DOI.org (Crossref), https://doi.org/10.1088/0953-8984/21/8/084204.

(5). Shannon, R. D. “Revised Effective Ionic Radii and Systematic Studies of Interatomic Distances in Halides and Chalcogenides.” Acta Crystallographica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography, vol. 32, no. 5, Sept. 1976, pp. 751–67. scripts.iucr.org, https://doi.org/10.1107/S0567739476001551.

• I gave my +1 long ago, but I want to say that I'm very surprised that this only has 2 upvotes so far! It was a very high-quality and high-effort answer! Sep 30, 2023 at 21:00