Discrepancy between numerical and transformed derivatives

Question

This is related to a prior question of mine Derivatives with respect to user defined vibrational modes. While that one focuses on potential software to use for derivatives with respect to normal modes, I want to present the actual problem that led me to that.

I have the derivatives of some property $P$ with respect to the $3N-6=M$ normal vibrational modes $\big\{Q_i\big\}$ of a molecule. I wanted to convert these modes to a local mode basis to more directly relate these derivatives to functional groups of the molecule. Conversion of the modes to a local basis can done by a simple unitary transformation [1]: $$\mathbf{Q}'=\mathbf{QU}$$ Here, $\mathbf{Q}$ is a $3N\times M$ matrix of the normal modes, $\mathbf{U}$ is $M\times M$ unitary matrix defined via an iterative process described in the linked paper, and $\mathbf{Q'}$ is a matrix of the normal modes.

With the modes transformed, I also want derivatives in this local mode basis. I have two ways of doing this:

• Transform the original derivatives: $\frac{\partial P}{\partial Q_i'}=\sum_jU_{ji}\frac{\partial P}{\partial Q_i}$ where the derivatives are arranged as column vectors.
• Compute numerical derivatives along the new mode: $\frac{P(X+hQ_i')-P(X)}{h|Q_i'|}$ where $X$ is the initial molecule geometry.

However, the transformed derivatives and the numerical derivatives of the local modes do not seem to match. If I test my procedure on the normal modes, the numerical derivatives agree with the ones I get from Gaussian. I'm concerned that perhaps I have something mixed up with removing/keeping the mass weighting of the modes (Gaussian fiddles with the coordinate representation a lot during vibrational analysis). Is there something obviously wrong with the procedure I have outlined above? Can I transform mass-weighted normal modes properly or do I need to ensure they are in cartesian coordinates before performing the transformation?

1. Jacob, C.R & Reiher, M. J. Chem. Phys. 130, 084106 (2009); DOI: 10.1063/1.3077690

Answer 1

A partial, preliminary answer to collect my thoughts (or maybe just a test of Cunningham's Law).

I believe I have determined the reason for the discrepancy. In terms of calculating a derivative with respect to a given mode, my numerical derivative procedure seems to be incorrect. While this approach should work fine for modes that are all orthogonal to each other (which is the case for the normal modes), the localization procedure should generally lead to nonorthogonal modes.

If I take the directional derivative using the formula in my question, I won't just be getting the derivative along the desired local mode. Instead, I will get additional contributions from any of the other modes that overlap with it. For a simple example of this, consider the figure above. In a coordinate system where the $x$ axis is tilted by some amount $\phi$ towards the $z$ axis. The directional derivative of a function $f(x,y,z)$ along $x$ would no longer just be $\frac{\partial f}{\partial x}$, but would now have a component related to $\frac{\partial f}{\partial z}$.

I suspect (but have not proven) that the other procedure of transforming the derivative just gives me the component corresponding to the individual mode I'm interested in.