Multiwfn forum

pbjo

Member

Registered: 2021-04-28

Posts: 1

IRI values are inverted compared to IRI paper?

When I look at the sections plots (figure 4c in https://chemrxiv.org/articles/preprint/ … /13591142)
the IRI value is above 2 when going far away from the molecule and in the atom cores, but when I use multiwfn tool to calculate IRI and visualize the IRI with VESTA, it looks like the IRI value goes to 0 far away from the molecule and in the atom cores.
How is the IRI calculated in multiwfn? Do you use the default value of a=1.1?

I have tried to also calculate IRI using numerical differentation from a file calculated with Gaussian cubegen and here it follows the convention in the paper (large values in cores and outside molecule and small values for the bonds).

Below is the Python script I use to compute IRI numerically after converting using cubegen:
cubegen 10  Density=SCF  qm9_000214_PBE1PBE_pcS-3.fchk qm9_000214_PBE1PBE_pcS-3_scf.cube -10

For some reason, the result I get from Cubegen->cube_to_iri.py does not look the same as when I use multiwfn directly.
Here are the two examples illustrated with VESTA. In both cases the scale is 0->2 where 0 is blue and red is 2.

Any idea what the problem here is?

import sys
import os
import io

import numpy as np
import ase
import ase.io
import ase.io.cube
from ase.calculators.vasp import VaspChargeDensity

def main():
    if len(sys.argv) > 1:
        cubefilepath = sys.argv[1]
    else:
        print("Usage: python cube_to_iri.py cube_file_path")
        sys.exit(1)

    if cubefilepath.endswith(".cube"):
        density, atoms, origin = read_cube(cubefilepath)
    else:
        density, atoms, origin = read_vasp(cubefilepath)

    cell = atoms.get_cell()

    if not cell.orthorhombic:
        print("Only works with orthorombic cell")
        sys.exit(1)

    step_sizes = cell.lengths()/density.shape[0:3]
    print(step_sizes)

    grad = np.gradient(density, *step_sizes)

    ##
    ## IRI(r) = ||grad(p(r))||/(p(r)**a)
    ##
    norm_grad = np.linalg.norm(grad, axis=0)
    a = 1.1 # Default value for a in IRI
    iri = norm_grad / (np.abs(density) ** 1.1)

    with open("iri.cube", "w") as f:
        ase.io.cube.write_cube(f, atoms, iri, origin, "IRI computed by finite difference")

def read_vasp(filepath):
    vasp_charge = VaspChargeDensity(filename=filepath)
    density = vasp_charge.chg[-1]  # separate density
    atoms = vasp_charge.atoms[-1]  # separate atom positions

    return density, atoms, np.zeros(3)  # TODO: Can we always assume origin at 0,0,0?


def read_cube(filepath):
    with open(filepath, "r") as f:
        cube = ase.io.cube.read_cube(f)
    # sometimes there is an entry at index 3
    # denoting the number of values for each grid position
    origin = cube["origin"][0:3]
    # by convention the cube electron density is given in electrons/Bohr^3,
    # and ase read_cube does not convert to electrons/Å^3, so we do the conversion here
    cube["data"] *= 1.0 / ase.units.Bohr ** 3
    return cube["data"], cube["atoms"], origin

if __name__ == "__main__":
    main()

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 1,830

Website

Re: IRI values are inverted compared to IRI paper?

Thank you for your attention to IRI.

Please note that definition of IRI has been changed in Multiwfn 3.8(dev) compared to Multiwfn 3.7.
IRI paper has been published very recently on Chemistry-Methods journal, see https://doi.org/10.1002/cmtd.202100007. The definition of the IRI in this paper is fully in line with latest Multiwfn version on Multiwfn website.

The code for evaluating IRI is "function IRIfunc" in function.f90 in Multiwfn source code package.