Posts by pbjo

pbjo · Multiwfn and wavefunction analysis

When I look at the sections plots (figure 4c inhttps://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

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

Board footer