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 [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
• 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? – taciteloquence May 21 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 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. – Nike Dattani May 21 at 3:37

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