# Discrepancy between numerical and transformed derivatives

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 : $$\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
• I just want to be clear: your problem is that taking the derivative and then doing the transformation is not giving you the same result at doing the transformation and then taking the derivative? May 21, 2020 at 2:43
• @taciteloquence basically, yes. I have the derivatives from outputted from a Gaussian calculation, as well as the normal modes themselves. If I transform the derivatives, I get a different result then if I transform the modes and recompute the derivatives.
– Tyberius
May 21, 2020 at 3:27
• I wonder if "numerical-convergence" might be appropriate here? It's not exactly talking about numerical-convergence like in the sense of SCF, but there's a numerical discrepancy between derivatives calculated numerically, so people "watching" the convergence tag might be interested. The people interested in "scientific computation" or "numerical techniques" might be interested. May 21, 2020 at 3:37
• I commented on your other question, but the issue is that you have U instead of the inverse of U. The gradient is in the dual (reciprocal) space, since it is d/dQ rather than Q (ie dimensions of 1/Q) Mar 21 at 2:36

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}$$.