Thank you so much for adding the requested feature to Multiwfn, as well as the - even more convenient - vector-based method. I've managed to code my protocol, and everything seems to be checking out so far...

21 Calculate gradient and curvature of electron density along a direction
21
Input XYZ coordinate of the point to be calculated, e.g. 2.0,2.4,1.1
or input index 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
Calculate gradient and curvature of electron density in which direction?
1 Along the direction by specifying a vector
2 Along the direction perpendicular to a given plane
1
Input X,Y,Z of the direction vector, e.g. 3.2,1.0,5
1.0,1.0,1.0
X,Y,Z of unit direction vector is    0.57735027    0.57735027    0.57735027
Electron density is                            0.0431233950 a.u.
Electron density gradient is                  -0.0000005835 a.u.
Electron density curvature is                  0.0388861044 a.u.

BTW: The X,Y,Z coordinates (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.5773502692   -0.5773502692   -0.5773502692
    0.5773502692    0.5773502692    0.5773502692

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

Your represention is a bit complicated to me... I just updated Multiwfn on its website, now option 21 of topology analysis is able to calculate curvature of electron density for arbitrary direction, user just needs to specify the direction vector. Hope this helps. As for the scan, you need to write a Bash or Python script to generate a series direction vectors, which are used in input stream for invoking Multiwfn to calculate the curvatures.

Best,

Tian

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

Unfortunately, to realize this aim, you have to generate wavefunction files by Gaussian for a series of rotated geometries. This is not difficult, it is easy to write a Bash shell script to automatically finish all operations. In each calculation, you can use guess=read keyword to read converged wavefunction in chk file of previous calculation, the convergence will only need one cycle, and thus performing single point calculations for all the geometries will not take much time.

Best regards,

Tian

I'm trying to evaluate the curvature of electron density at a point (obtained from a .wfx file), along a specific user-defined axis/vector. And not the Laplacian eigenvalues/vectors themselves, i.e. the three orthogonal axes of maximum absolute electron density curvature. I've tried to evaluate the data generated by 'Output all properties at a point', at the origin (0.0, 0.0, 0.0) of my molecular system, before and after applying a rotation operation about the z-axis.  I expect the curvature along the z-axis to remain the same after applying the z-rotation (eg. of 30 degrees), however the values of curvature along the x- and y-axes should change - for systems with an asymmetrical electron density distribution about the origin. However the results I'm getting are not consistent (see below). The value of the Laplacian is also changing, which shouldn't be the case, as the same point is being evaluated both times: x=0.0, y=0.0, z=0.0.

I believe this is probably related to the following warning in the manual: 'Note that transformation of coefficient matrix of orbitals due to rotation and inversion operations is not taken into account by this function.' I've also used the translation operation (300>7>1>Enter>) in order to set a critical point coordinate to the origin of the coordinate system so that a z-operation could be employed, as described above.

I would like to be able to extract the 'component of Laplacian' (or curvatures) along any axes (x/y) that can be generated after rotation of the z-axis from 0 to 90 degrees (in 2 degree increments, for example). Or perhaps by specifying an axis/vector along which the electron density curvature will be evaluated at a given point in space? Is there a way to do this in Multiwfn, without rotating the input geometry in  the Gaussian 16 input file and generating new .wfx files of the desired orientations?

Components of Laplacian in x/y/z are:
-0.3422890839E-01 -0.3548794863E-01  0.1890675051E+00
Total:  0.1193506481E+00
After rotation about the z-axis of 180 degrees (Functions: 300>7>3>Enter>3>180.0):
Components of Laplacian in x/y/z are:
-0.3547969583E-01 -0.4232451008E-01  0.1901760811E+00
Total:  0.1123718752E+00

PS: I'm using the latest Windows version of Multiwfn [Version 3.8(dev), release date: 2023-May-18]

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