Multiwfn forum

Dear Prof Tian Lu,

Thank you for your response and for your suggestions. That protocol is definitely an option, thank you for that. However, for my problem case I would need to generate numerous input/output files (90*7*19) and run many calculations (doable via a script as per your suggestion), but it would be easier to analyze a .wfx directly using Multiwfn and this also avoid the generation of many input/output files to keep track of.

The request: I was wondering if a specific function "Calculate density curvature perpendicular to a specific plane at a point" could be (easily?) modified in the Multiwfn code to accept cartesian coordinates in addition to atom indices, that will allow for the evaluation of electron density curvature along user-defined axis, in an roundabout/indirect way:
Function 2>21>0.0,0.0,0.0>1><three sets of x,y,z coordinate data (in Bohr) (see the comment below).

The rest of the following is about the protocol that I would like to implement (on my own):
For the system that I'm trying to analyze:
1) the Gaussian input geometries were oriented such that the BCP bond path was approximately colinear to the z-axis, such that I would only need to evaluate the electron density (ED) curvature along the z-axis once and then evaluate the ED curvature along several axis/lines on the xy plane in order to identify the axis of maximum absolute ED curvature (coplanar).
2) The Gaussian input geometries were also translated such that the coordinate origin would be - approximately - in the same location at the BCP of interest.

My suggested Protocol/Workflow:
1) Identify the BCP coordinate via 'Topology analysis'
2) Translate the geometry such that the BCP (or any critical point) is located at the origin (0.0, 0.0, 0.0), via:  300>7>1>Enter>{x,y,z-translation vector, in Angstrom} <- Here the correct curvature information would be obtained after a translation operation was applied, but not for the rotation operation (I think, please confirm this)
3) Use the modified "Calculate density curvature perpendicular to a specific plane at a point" Multiwfn executable function [2>21>0.0,0.0,0.0>1><three sets of x,y,z coordinate data (in Bohr) (see the comment below)] to determine the ED curvature along a user-defined axis, using the following approach:
a) The "point" coordinate would be 0.0,0.0,0.0
b) the three coordinates defining a plane that would have the "point" located on the plane:
1 = 0.0, 0.0, 1.0 <- Comment: 1st point is on the z-axis (the axis of rotation)
2 = 0.0, 0.0, -1.0 <- Comment: 2nd point is on the z-axis (the axis of rotation)
3 = -cos(alpha), sin(alpha), 0.0 <- Comment: 3rd point is located on the xy plane, that defines a plane that is perpendicular to the axis/line of interest, where "alpha" would be an angle of clockwise rotation about the z-axis afterwards the y-axis will now be located 'alpha' degrees relative to its original orientations/direction.
With this approach a 'Z-rotation' could be simulated from 0 to 180 degrees for example and the appropriate x & y axes could be assigned based on 'alpha' and the curvatures, without needing to resubmit any jobs to Gaussian (this also includes the ED gradient information as well).

21 Calculate density curvature perpendicular to a specific plane at a point
21
Input the coordinate, e.g. 2.0,2.4,1.1   or input indices of a CP, e.g. 4
0.0,0.0,0.0
The coordinate you inputted is in which unit?  1=Bohr  2=Angstrom
1
Input indices of three atoms to define a plane, e.g. 3,4,9 <- Comment: Add the option of providing three coordinates to define a plane(eg.: 0.0,0.0,0.0; 0.0,0.0,-1.0; <a>,<b>,0.0), where <a> = -cos(alpha) & <b> = sin(alpha) (these values can be entered as floating point numbers after result has been determined via a script/'by hand' )
1,3,5
The unit normal vector is    0.09744032   -0.99275034   -0.07037154
Electron density is                            0.0454599448
Electron density gradient is                   0.0000001537
Electron density curvature is                 -0.0344065605

BTW: The X,Y,Z coordinate (row) of current point, the points below and above 1 Angstrom of the plane from current point, respectively (in Angstrom).
    0.0000000000    0.0000000000    0.0000000000
   -0.0974403249    0.9927503359    0.0703715396
    0.0974403249   -0.9927503359   -0.0703715396

Please let me know if you would be willing to modify Multiwfn to allow the above request, and if you think this protocol will work.

PS: Please keep in mind, I know a bit Python programming (enough to automate/script the above protocol) but not any Fortran.

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

Dear Prof 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

Dear Prof Tian Lu,

I have been trying to investigate the natural orbital approach that you suggested for determining the orbitally optimised CCD density (No singles excitation contributions are available for the density calculation, so I did not determine OO-CCSD instead I determined OO-CCD), as well as trying to understand the composition of the various ORCA output files (.gbw, .optorb, .nat, .mdcip) in terms of their generated electron density .cube files.

I have received a response from Stanpa (Stan, a PhD Student at Max-Planck-Institut für Kohlenforschung) on the ORCA forums that is a bit concerning (see https://orcaforum.kofo.mpg.de/viewtopic … 16#p38419). Stan says that the only way to obtain the OO-CCD density is via generating the .cube file by using the ELDENSMDCI keyword in the %plots block in ORCA (point 1 below). From your approach it seems as if you have employed point 4 below, which he says will produce the "HF electron density in the optimised OO-CCD basis (orbitals are rotated to make the OO-CCD Lagrangian stationary)." I am quite confused, as there appear to be multiple types of 'HF electron densities'.

So it seems I have to use approach (1) below to generate the OO-CCD electron density .cube files, since I cannot use the generated .molden files (via orca_2mkl).

Here is an excerpt from the Stanpa's response (his responses are in bold, while my questions are in plain text):

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

Dear Prof Tian Lu,

My apologies for the delay. Thank you for introducing this new feature in Multiwfn (24/03/2022). This is long post, so feel free to focus your replies to question 1 and 3. The errors shown in point 3 are the most concerning to me at the moment, and I need to decide on which protocol I should use for interaction density .cube file generation :1) Multwfn's 'calc. grid data>set custom operation' route where is specify the .molden files or 2) ORCA's plot generation utility to generate each component .cube file (for each .mdcip file) followed by Multiwfn's 'Grid data. calculation').

1) I have tested the new special feature and it produces output. I would like to know how the atomic space - used in the option "2 Calculate atomic and molecular multipole moments and <r^2>" - and corresponding atomic dipole moments are determined.

Does it use the fuzzy atomic space (i.e. the atomic weighing function) generated via the selected fuzzy partitioning method (for example, the default Becke's method), where it evaluates the trilinearly interpolated .cube data (charge density data) instead of the electron density derived from GTFs - using the aforementioned weighing function- to yield the atomic dipole moments? Does the code take into consideration for charge density contamination originating at atoms other than the selected atom for which the atomic dipole moment is to be determined?

I thank you for your efforts. If I'm not able to use the aforementioned method in my work, I should be able to write a Python script to generate approximate atomic dipole moments via the domain-based method that I proposed, and with the help of the domain export utility that you provided in the domain analysis module. An initial cursory investigation has revealed some potential concerns though, see below.

---------------------------------------------------------------------------------------------------------------------------------------------

2) For the interpyridyl interaction charge density data: I have evaluated several thresholds of positive and negative domains and I have determined the net charge for each set of domains, via integration of all positive domains and all negative domains, as defined by the domain boundary conditions/thresholds shown in the figure below. It can be seen from the figure that the relative polarisation error [(integral of + density) + (integral of - density )/(integral of + density)*100] contribute significantly for certain domain thresholds. The assumption being that the net charge should be zero for each set of domains, i.e. polarisation of x electrons, should lead to a electron peak of x electrons (+ domain) and a electron hole of x electrons (- domain).

There are two domain thresholds where the relative polarisation errors are ideal or close to 0: -0.002 > ρ > 0.002 and -0.00005 > ρ > 0.00005 a.u. As expected the polarisation error for the full .cube (all data) is close to zero, however significant errors approaching -12.44% and >14.82% can be observed for other domain thresholds. This leads me to conclude that the charge density modelled for the non-ideal domain threshold cases do not show the correct polarisation behaviour. Determination of approximate atomic dipole moments from the | 0.00005 a.u.| data will clearly be problematic, since the atomic charge density is significantly delocalised and merges/mixes significantly with nearby atomic charge density, whereas the atomic dipole moments determined from the |0.002 a.u.| show clearly unmerged/defined atom-centred polarision lobes, which would be ideal for atomic dipole moment determination. The drawback of the latter data being that only a fraction of the total charge density data will be considered, i.e. the |0.002 a.u.| data with a total polarisation of 29.68 milielectrons, compared to the full |0.0 a.u.| data with a total polarisation of 96.06 milielectrons.

Perhaps by applying a scaling factor to the sum of the individual atomic dipoles (as determined from the below equations) calculated from the |0.002 a.u.| data to the produce the correct total polarisation  for the the full .cube data (i.e. |0.0 (a.u.)| data) can be approximated?

I have written the equations and definitions for the interaction-specific atomic dipole moments that I would like to determine (See below). Please note that I do not have a strong mathematical background, so I may have made several mistakes along the way. From the equations it can be seen that symmetrical polarisation of electron density would produce a norm of the total atomic dipole moment of 0. The atom dipole moment as provided here, would give an indication of the deviation from symmetrical distribution of charge density. Do you have any concerns or suggestions here? 

Note: The charge density threshold of + or - 0.00004 a.u. given here is only a placeholder, until I decide how I would like to approach this problem then I will update it.

---------------------------------------------------------------------------------------------------------------------------------------------

3) I would like to know if the following discrepancy in electron densities between the generated .cube data given by ORCA and Multiwfn, was due to some setting that I did not adjust correctly in MultiWFN (I did update some values in the settings.ini file, see below for more details)? If not, then was this error due to the conversion of the .optorb file to .molden file format, by the orca_2mkl program? I assume that the ORCA %plots utility correctly determines the electron density that it outputs in the .cube file, when it reads the CCSD MDCI density from the .mdcip file.

I have checked whether the ORCA gaussian .cube files generated match the corresponding Multiwfn-generated .cube files. These data are for the molecular electron densities and not the interaction densities, as I wanted to test/evaluate potential sources of error. I have found the following deviation/error in the number of electrons across all grid points, between the two files (Multwfn .cube relative to ORCA .cube):

The procedure used to generate the two Gaussian .cube files (and the difference .cube) was as follows:

3.a) a .cube was generated, via the %plots module within ORCA 4.2.0, using the CCSD orbitally-optimised .mdcip density.
3.b) a .molden file generated from an orbitally-optimised coupled-cluster wave function file (.mdci.optorb ) was produced after the CCSD calculation finished, via the orca_2mkl program of ORCA 4.2.0., where the .molden file was loaded in Multiwfn (with the following settings.ini options set: 'expcutoff= 1' and 'idelvirorb= 0', in order to generate high accuracy electron density data from the wave-function file). The following sequence of commands were entered after loading the .molden file: '5 Output and plot specific property within a spatial region (calc. grid data)' > '1 Electron density' > '8 Use grid setting of another cube file' > 'here the .cube file generated in (3.a) was supplied' > '2 Export data to a Gaussian-type cube file in current folder'.

A new session of Multiwfn was started, where cube 3.a was loaded and cube 3.b was subtracted from cube 3.a, which was then saved as a new difference .cube file (cube 3.c):
'13 Process grid data' > '11 Grid data calculation' > '4 Subtract a grid file' > 'cube 2.b file path' > '0 Export present grid data to Gaussian-type cube file (.cub)' > 'cube 3.c'.

Next, a new session of Multiwfn was started, where the difference cube 'cube 3.c' was loaded ('2_DPA_41ac_mdcip_orca_min_molden_diff.cube'), which produced the following summary:

I have uploaded the files to a shared Google drive directory (I have compressed all files larger than 300MB, via 7zip):
1) the '2_DPA_41ac_DLPNO_MP2_gas_opt_AVDZ_mdcip_900' file is the .cube generated by ORCA 4.2.0 (using the ELDENSMDCI keyword) (3.a)
2) the '2_DPA_41ac_DLPNO_MP2_gas_opt_CCSD_AVDZ.mdci.optorb_molden' file is the .cube generated by Multiwfn (3.b)
3) the '2_DPA_41ac_mdcip_orca_min_molden_diff' file is the difference cube (3.c)
4) I've also added the .optorb file, the .molden file, .mdcip file and the orca_2mkl binary from ORCA 4.2.0

Here is the link: https://drive.google.com/drive/folders/ … sp=sharing

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

Dear Prof Tian Lu,

Thank you so much for updating Multiwfn and for adding option 12, it is much appreciated. It will be much easier for me to extract the required data now.

Please note that I am still familiarising myself with the theory behind determining (atomic) dipole moments. I may not be explaining myself properly.

As I understand it, an atomic dipole moment is determined from the (presence of) electron density data in all regions that are associated with a specific atom's, via discrete or fuzzy partitioning methods, as well as the the nuclear coordinate and nuclear charge. Whereas, the 'intramolecular charge density' consists of regions of space associated with the presence of electron density (+ values/lobes) and the absence of electron density (- values/lobes), relative to the electron density distribution of the conformer/molecule of interest (eg. DMA) and are obtained for a specific set of fragments, i.e. the CH3...CH3 interaction in DMA. For larger molecules (see the .cube file, where >6 sets can be defined), there may be several interfragment sets/pairs that may be defined, where the interaction energy and interaction charge density may be determined for each set. As you can see from the .cube data the dipoles are centred on each atom (astride each atom), but are not spherically distributed around the atoms, so the scenario is different from standard atom-centred GTO data. I'm not certain if these type of interaction-specific atom-centred dipoles are still considered to be a type of atomic dipole.

From the Gaussian .cube file, I suppose I could determine the following contributing components to the (atomic) dipole moment for each interaction set:[+(charge density),+(nuclear)] and [-(charge density),+(nuclear)], where +(nuclear) would be the nuclear charge and where unmerged charge density domains could be analysed for each atom (this would define the atomic space), and then all component dipoles could be combined to generate the total dipole associated with each atom. I don't know if the method/code in the 'Fuzzy atomic space'>'Calculate atomic and molecular multipole moments and <r^2>' module in Multiwfn will give the desired output, these methods utilise generated atomic densities and I do not believe this will model the interaction-specific atom-centred polarised density correctly.

I would like to visualise and evaluate the atom-centred dipole vectors for a chosen interaction for several conformers of a molecule, and I would like to establish what type of correlation exists between the calculated interfragment interaction energies and the 'charge density polarisation'/dipole data.

Thank you, I was unaware of the the 'main function 5 > sub-function 0' method of calculating the 'interaction charge density' .cube files. I could generate .molden files, via ORCA, for the molecular fragments and load them using this method. This should speed up the generation of these Guassian .cube files considerably, assuming that 'interaction only' wave function files cannot be generated for analytical evaluation instead.

Please let me know if you have any other questions or suggestions.

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

Dear Prof Tian Lu,

Thank you for your prompt response.

The Gaussian .cube file is quite large (it's about 10 GB, for 900^3 points). I've made the uncompressed file available via a shared Google drive link (as I don't think compression will be very effective here):
https://drive.google.com/file/d/1ty0eDf … sp=sharing
(Please let me know if you can access this link)

On a separate note:
1) is it possible to export the data for each individual domain as a .pdb or a type of coordinate file (x,y,z,value), instead of only exporting the domain boundary points, in multiwfn. As I would like to 1) integrate the domains, 2) determine the minima/maxima of the charge density and 3) calculate the (atomic?) dipole moments for all the dipole pairs associated with each atom in the system from the Gaussian .cube file.

2) I'm not quite certain how I would go about calculating the aforementioned (atomic?) dipole moments, for the dipoles generated by the molecular tailoring approach.

From section '3.18.3 Calculate atomic and molecular multipole moments and <r2> (2)', in the manual, and from the requirements of an electric dipole moment: 'When two electrical charges, of opposite sign and equal magnitude, are separated by a distance, an electric dipole is established' I have interpreted the information as follows:

The atomic dipole moment is calculated based on the assumption that the 'number of electrons associated with an atom' = 'nuclear charge of an atom' and are defined according to the atomic dipole vector components (involving the relative spatial positions only) that are centred on the nuclear position and terminate at all regions in space, where the atomic dipole vector components  may be separated into x,y,z components and mulitplied by the number of electrons - associated with a specific atom- at all points in space, thereafter the sum of all dipole components produce in the x-, y- and z-dimensions may then be converted into the total atomic dipole moment vector, which has a magnitude and direction, by the formula supplied. The molecular dipole moment would then be the sum of all atomic dipole vectors.

For the case of interaction charge given by the molecular tailoring approach, the positive charge would no longer be situated on the nuclear position, but would be delocalised into the + component of the dipole (i.e. the region of electron depletion or negative charge density). How would one compute the dipole moment for this case, perhaps by approximating the dipole moment by determining the vector from the centroid/barycenter of each lobe (from - to + charge density) multiplied by the total number of electrons (minimum?/average?/maximum?) in the lobe(s)?

3) (Please read below for context) Is it possible to generate an 'interaction only' wave function file (eg. molden) (and ultimately obtaining the 'interaction only' electron density data) by applying the below operations to the individual fragment .molden files, i.e. by modifying the MO coefficients, or is it only possible by applying the below operations to the individual fragment Gaussian .cube files? Since the cube files contain Cartesian/spatial data, whereas the .molden files contain atom-centred data.

Here is an example of the molecular tailoring approach for determination of the intramolecular interaction energy/density:
For example, to calculate the intramolecular interaction between two methyl groups in dimethylamine (DMA) in a specific conformation, one would determine the electron density data of
1) dimethylamine (DMA) - of that specific conformation given in Cartessian coordinates- ,
2) the left-hand side(LHS), neutral methylamine (generated by substituting the RHS methyl group with H and performing a cartesian-based constrained geometry optimisation where only the non-physical H is optimised) -> fragment MA1,
3) the right-hand side (RHS) methylamine (generated by substituting the LHS methyl group with H and performing a cartesian-based constrained geometry optimisation where only the non-physical H is optimised)-> fragment MA2,
4) the central NH3 (generated by substituting the both methyl groups with the coordinates of the H atoms as calculated in (2) and (3)) -> fragment NH3

The approximate 'intramolecular charge density' = ρ(DMA) -  ρ(MA1) - ρ(MA2) + ρ(NH3), where the scalar field or Gaussian .cube grids - for all components - are determined according to the grid assigned to ρ(DMA).

Kind regards
Cameron