Posts by pbjo

pbjo · Multiwfn and wavefunction analysis

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()

Multiwfn forum

#1 Multiwfn and wavefunction analysis » IRI values are inverted compared to IRI paper? » 2021-04-28 12:41:17

Board footer