```python
import psi4
from psi4 import *
from psi4.core import *
import numpy as np
import os
sys.path.append('os.getcwd()')
from opt_helper import stre, bend, intcosMisc, linearAlgebra
```

## The Step Back-transformation
An optimization algorithm carried out in internal coordinates (see, e.g., the RFO tutorial) will generate a displacement step to be taken in internal coordinates.  The conversion of the step into Cartesian coordinates is here called the "back-transformation" ("back", since the original gradient was computed in Cartesians).

As shown in the tutorial on coordinates and the B-matrix,

$$\textbf {B} \Delta x = \Delta q $$

and the $\textbf A^T$ matrix defined by

$$ \textbf A^T \equiv (\textbf{B} \textbf{u} \textbf {B}^T)^{-1} \textbf {B} \textbf{u}$$

was shown to be the left inverse of $\textbf B^T$ where __u__ is an arbitrary symmetric matrix.  Attention must be paid to the non-square nature of __A__ and __B__.  Here, we have

\begin{align}
\textbf B \Delta x &= \Delta q \\
\\
\textbf B \Delta x &= \big( \textbf{BuB}^T \big) \big( \textbf{BuB}^T\big)^{-1} \Delta q \\
\\
\textbf B \Delta x &= \textbf B \big[ \textbf{uB}^T \big( \textbf{BuB}^T\big)^{-1}\big] \Delta q \\
\\
\Delta x &= \textbf{uB}^T \big( \textbf{BuB}^T\big)^{-1} \Delta q = \textbf A \Delta q \\
\end{align}

The __u__ matrix may be chosen to be the unit matrix which gives

$$\Delta x = \textbf B^T (\textbf B \textbf B^T)^{-1} \Delta q$$

where redundant coordinates can be accommodated simply by using the generalized inverse.  It is common to introduce $ \textbf{G} = \textbf B \textbf B^T $ and write the expression as

$$ \Delta x = \textbf{B}^T \textbf{G}^{-1} \Delta q$$

Note the __G__ matrix is a square matrix of dimension (number of internals) by (number of internals).  This equation is exact only for infinitesimal displacements, because the B-matrix elements depend upon the molecular geometry (i.e., the Cartesian coordinates).  Thus, the back-transformation is carried out iteratively.

To converge on a Cartesian geometry with the desired internal coordinate values, we repeatedly compute the difference between the current internal coordinate values and the desired ones (generating repeated $\Delta q$'s) and using the equation above to compute a new Cartesian geometry.

### Illustration of back-transformation
The back-transformation will now be demonstrated by taking a 0.2 au step increase in the bond lengths and a 5 degree increase in the bond angle of a water molecule.


```python
# Setup the water molecule and coordinates.
mol = psi4.geometry("""
O
H 1 1.7
H 1 1.7 2 104
unit au
""")
# We'll use cc-pVDZ RHF.
psi4.set_options({"basis": "cc-pvdz"})
mol.update_geometry()
xyz_0 = np.array( mol.geometry() )

# Generate the internal coordinates manually.  Show their values.
intcos = [stre.STRE(0,1), stre.STRE(0,2), bend.BEND(1,0,2)]
print("%15s%15s" % ('Coordinate', 'Value'))
for I in intcos:
  print("%15s = %15.8f %15.8f" % (I, I.q(xyz_0), I.qShow(xyz_0)))

# Handy variables for later.
Natom = mol.natom()
Nintco = len(intcos)
Ncart = 3*Natom
```


```python
# Create an internal coordinate displacement of +0.2au in bond lengths,
# and +5 degrees in the bond angle.
dq = np.array( [0.2, 0.2, 5.0/180*np.pi], float)

B = intcosMisc.Bmat(intcos, xyz_0)
G = np.dot(B, B.T)
G_inv = linearAlgebra.symmMatInv(G, redundant=True)

# Dx = B^T G^(-1) Dq
dx = np.dot(B.T, np.dot(G_inv, dq))
print("Displacement in Cartesians")
print(dx)

# Add Dx to original geometry.
xyz_1 = np.add(np.reshape(dx, (3, 3)), xyz_0)
print("New geometry in cartesians")
print(xyz_1)

# Compute internal coordinate values of new geometry.
print("\n%15s%15s" % ('Coordinate', 'Value'))
for I in intcos:
  print("%15s = %15.8f %15.8f" % (I, I.q(xyz_1), I.qShow(xyz_1)))
```

You see that the desired internal coordinate value is not _exactly_ achieved.  You can play with the desired displacement and observe more diverse behavior.  For water, if you displace only the bond lengths, then the result will be exact, because if the bond angle is fixed then the direction of the displacements (_s_-vectors on each atom) are constant wrt to the bond lengths.  On the other hand, the displacement directions for the bend depend upon the value of the angle.  So if you displace only along a bend, the result will not be exact.  In general, the result is reasonable but only approximate for small displacements.

### Illustration of iterative back-transformation
Finally, we demonstrate how convergence to the desired internal coordinate displacement can be achieved by an interative process.


```python
# Create array of target internal coordinate values.
dq_target = np.array( [0.2, 0.2, 5.0/180*np.pi], float)
q_target = np.zeros( (len(intcos)), float)
for i, intco in enumerate(intcos):
    q_target[i] = intco.q(xyz_0) + dq_target[i]
    
xyz = xyz_0.copy()
rms_dq = 1
niter = 1

while rms_dq > 1e-10:
    print("Iteration %d" % niter)
    dq = dq_target.copy()

    # Compute distance from target in internal coordinates.
    for i, intco in enumerate(intcos):
        dq[i] = q_target[i] - intco.q(xyz)
    rms_dq = np.sqrt(np.mean(dq**2))
    print("\tRMS(dq) = %10.5e" % rms_dq)
    

    # Dx = B^T G^(-1) Dq
    B = intcosMisc.Bmat(intcos, xyz)
    G = np.dot(B, B.T)
    G_inv = linearAlgebra.symmMatInv(G, redundant=True)
    dx = np.dot(B.T, np.dot(G_inv, dq))
    print("\tRMS(dx) = %10.5e" % np.sqrt(np.mean(dx**2)))

    # Compute new Cartesian geometry.
    xyz[:] += np.reshape(dx, (3,3))
    niter += 1

print("\nFinal converged geometry.")
print(xyz)

# Compute internal coordinate values of new geometry.
print("\n%15s%15s" % ('Coordinate', 'Value'))
for I in intcos:
  print("%15s = %15.8f %15.8f" % (I, I.q(xyz), I.qShow(xyz)))
```

The exact desired displacement is achieved.

Due to the non-orthogonal nature of the coordinates, the iterations may not always converge.  In this case, common tactics include using the Cartesian geometry generated by the first back-transformation step, or using the Cartesian geometry that was closest to the desired internal coordinates.  Hopefully, as a geometry optimization proceeds, the forces and displacements get smaller and convergence occurs.

A serious complication in procedures such as this one are discontinuities in the values of the internal coordinates.  In some way, the internal coordinate values must be canonicalized so that, e.g., an increase in a torsion from 179 degrees to -178 degrees is interpreted as an increase of 3 degrees.  Similar problems present for bond angles near 180 degrees.  (The consistent computation of these changes in values and forces is also critical in Hessian update schemes.)

2ffa88155d34e02cd25aa73d1c3f8a031c9991e2
ipynb
Jupyter Notebook
Example/Psi4Numpy/13-GeometryOptimization/13e_Step-Backtransformation.ipynb
yychuang/109-2-compchem-lite
cbf17e542f9447e89fb48de1b28759419ffff956
BSD-3-Clause
2017-03-01T08:04:48.000Z
2022-03-23T08:52:04.000Z
Example/Psi4Numpy/13-GeometryOptimization/13e_Step-Backtransformation.ipynb
yychuang/109-2-compchem-lite
cbf17e542f9447e89fb48de1b28759419ffff956
BSD-3-Clause
2017-03-03T13:20:20.000Z
2022-03-05T18:20:27.000Z
Example/Psi4Numpy/13-GeometryOptimization/13e_Step-Backtransformation.ipynb
yychuang/109-2-compchem-lite
cbf17e542f9447e89fb48de1b28759419ffff956
BSD-3-Clause
2017-02-17T19:44:47.000Z
2022-03-22T05:52:43.000Z
Qwen/Qwen-72B
1. YES
2. YES