1 You don't need to fit an analytic function for all data in the .rad. You can simply use linear interpolation to determine rho at a given radial distance based on rho at its adjacent two points. Using more advanced interpolation technique such as Lagrangian interpolation would be slightly better.

2 Please search "subroutine proatmgrad" in function.f90 of Multiwfn source code package, you will find answer. Note that this subroutine has a bug, which has been fixed very recently, so please download latest version of source code package.

Best,

Tian Lu

Thank you for your response, and apologies for my severely delayed response. I would greatly appreciate it if you would assist me with the following:

I've recently considered whether the Iterative-Hirshfeld electron density values may be determined via a Python script, for any point in space, based on the following approach:
1) Calculating HI atomic charges for a molecular system using symmetrized atomic radial densities .rad files generated either Manually or via Gaussian (via Multiwfn), at the appropriate level of theory. For example if the HI charge for 1(N) was -0.15 then the N-1.rad and N_0.rad files would be used in the next step,
2) a function could then be fitted to the the corresponding x-data[radial distance w.r.t Nucleus (bohr)] and y-data[electron density, a.u.] for each of the .rad files associated with an atom (N-1.rad and N_0.rad),
3) subsequently a linear interpolation could be performed, at a specific point in space or radial distance, and the atomic charges could be used as a weighing function, for example, if a) the electron density 1.01575 bohr from the nucleus is 0.269600 au in the N_0.rad file and b) the electron density 1.01575 bohr from the nucleus is 0.275784 au in the N-1.rad file, then a HI atomic charge of -0.15 would indicate a 85% weight of (a) and 15% weight of (b) yielding an HI-based electron density of 0.270528 a.u. at 1.01575 bohr from the nucleus [0.269600+(0.275784-0.269600)*0.15 OR (1+(-0.15)*0.269600-(-0.15)*0.275784],
4) calculating the total HI electron density, at a specific point in space, would involve using a similar approach to determine the sum of the HI-(promolecular) electron density contributions from all atoms in the molecular system.
5) The deformation density could then be determined by subtracting the HI-promolecular electron density (4) from the molecular electron density, at the corresponding point in space, i.e p_def = p_mol - p_pro.

Question 1: If this approach is sound, what exact type/form of function would need to be used to fit the .rad data, to determine the electron density at radial distances that are not reported in the .rad files? Some type of exponential function, perhaps? Which values would be constant and which coefficients would need to be optimised via a least-squares regression procedure?

Question 2: For evaluation of the delta-g parameter of the Independent Gradient Model (IGM), how would the gradient norm for a point in 3D space, associated with a specific atom (eg. 1(N)), be determined analytically for the HI-radial electron densities? Furthermore, how would this be determined numerically, perhaps, by calculating the electron density at the point of interest in addition to 6 orthogonal points (+x,-x,+y,-y,+z,-z) around the point of interest (using a small step size, what distance should be used?), then determining, via finite/central difference in the x/y/z direction the gradient vector and the gradient norm?

Thank you for your assistance with this matter.
Kind regards
Cameron

Currently, it is not formally supported in Multiwfn to calculate Hirsheld-I promolecular- or deformation-electron densities.

However, you can use Multiwfn to calculate H-I charges first, and then try to manually perform fractional occupation calculation for every atom by specifying net charge as corresponding H-I charge. Some programs such as NWChem supports this feature. Via this calculation, it is possible to obtain atomic wavefunction file corresponding to final Hirshfeld-I states. I am not quite sure if this idea works, but it is worth to consider and seemingly feasible.

About the last question, Multiwfn doesn't know which .rad files will be actually involved. For example, if at a step charge of N is -0.15, then in this step only .rad files of N(0) and N(-1) are needed; if after that the charge changes to e.g. 0.05, then next step N(0) and N(+1) will used. If needed .rad file(s) cannot be found during H-I iteration, the calculation will abnormally stopped.

Best regards,

Tian

I would like to know if it is possible to calculate the promolecular and deformation properties with Multiwfn, via the iterative Hirshfeld method (Hirshfeld-I), perhaps by using the mainfunction '5 Output and plot specific property within a spatial region (calc. grid data)' module? Specifically, my aim is to produce the corresponding Gaussian .cube files of the Hirsheld-I promolecular- or deformation-electron densities.

As described in section 3.9.13, the symmetrized atomic radial densities (.rad) can be computed from the AIM .wfn files for the required atomic oxidation states - obtained for the theoretical method of interest, eg. orbital-optimised CCSD/aug-cc-pVDZ data, for the N(0), N(+1) and N(-1) species - can be utilized instead of the default atomrad/atomwfn data.

I have investigated the following methods for determining the H-I promolecular electron densities in Multiwfn, without success:

1) As I understand it, the mainfunction '5 Output and plot specific property within a spatial region (calc. grid data)' >  '-1 Obtain promolecule property' > '1 Electron density' route determines the Hirshfeld promolecular atomic densities for the neutral species, and not the more accurate or 'physically meaningful' Hirshfeld-I promolecular density. Presumably, the Special functions (1000) > '17 Generate promolecular wavefunction by calculating and combining atomic ones' route also produces the Hirshfled promolecular density.

2) The mainfunction '7 Population analysis and calculation of atomic charges' > '15 Hirshfeld-I atomic charge', calculates the HI atomic charges only.

3) Additionally, the '15 Fuzzy atomic space analysis' module, provides the function to '1 Perform integration in fuzzy atomic spaces for a real space function' using the Hirshfeld-I partitioning space (via the subfunction '-1 Select method for partitioning atomic space, current:'). This allow one to determine the integral of the Hirshfeld-I proatom electron density to yield the atomic population value, similar to point (2).

Finally, is it possible to identify which .rad files are being utilised for the iterative Hirshfeld (H-I) method from the Multiwfn output? Does Multiwfn firstly evaluate the initial Hirshfeld atomic charges and then determine the required .rad files to be used from the atomic charge integers adjacent to that value, for example, if the Hirshfeld atomic charge is -0.150251 for 1(N), then Multiwfn will always use the N(0) and N(-1) .rad files for linear interpolation?

Kind regards
Cameron