Conversión de la matriz Z en coordenadas cartesianas

Question

Para empezar, escribí un script que obtiene las coordenadas cartesianas de la molécula como entrada en el siguiente. Estas son $x,y,z$ coordenadas de la molécula de H2O2.

1 O          -1.7529   -0.5188    1.3324
2 O          -0.4737   -0.1091    0.7774
3 H          -2.2902    0.1678    0.8933
4 H           0.0636   -0.7957    1.2164

Luego, el script construye una matriz Z con ellos, así:

Z-mat :
O
O   1   1.45335189476
H   1   0.976176039452  2   96.5694760083
H   2   0.976131061897  1   96.5720363573   3   -179.995395182

Ahora necesito realizar la operación inversa y utilizar esta matriz Z como entrada y definir $x,y,z$ coordenadas para cada átomo. Esto es convertir una matriz Z en coordenadas cartesianas.

Mi pregunta es después de establecer el primer átomo como 0,0,0

1 O          0   0   0

y el segundo como 0,0,(distancia del primero) para ponerlo en el eje z

2 O          0   0   1.45335189476

¿Cómo debo tratar los átomos 3 y 4?

El tercer átomo debe tener coordenadas que son algo así si se lee correctamente:

3 H          0   distance*sin(angle)  z2+distance*cos(angle)

Tomando z2 como la coordenada z del átomo 2. No estoy seguro de si debo calcularlo como z2 + distance.cos(angle) o z2 - distance.cos(angle) y de qué depende, si ambas cosas son posibles.

Para el 4º átomo, utilizo coordenadas esféricas

r, theta, phi = (0.976, 96.572, -179.995)

para calcular $x,y,z$ valores de las fórmulas

x = r * sin(theta) * cos(phi)
y = r * sin(theta) * sin(phi)
z = r * cos(theta)

Si no hay ningún error hasta este punto, cómo voy a calcular las coordenadas cartesianas del 4º átomo utilizando estos $x,y,z$ ¿valores?

---

Mientras buscaba encontré El trabajo de TMPChem en GitHub y hace exactamente lo que quiero. Sin embargo, en su trabajo, hay una parte matemática que no entiendo:

# get local axis system from 3 coordinates
def get_local_axes(coords1, coords2, coords3):
    u21 = get_u12(coords1, coords2) #calculating vector between that points 1-2
    u23 = get_u12(coords2, coords3) #calculating vector between that points 2-3
    if (abs(get_udp(u21, u23)) >= 1.0):
        print('\nError: Co-linear atoms in an internal coordinate definition')
        sys.exit()
    u23c21 = get_ucp(u23, u21) # unit cross product
    u21c23c21 = get_ucp(u21, u23c21) # unit cross product
    z = u21
    y = u21c23c21
    x = get_ucp(y, z)
    local_axes = [x, y, z]
    return local_axes

¿Qué es "obtener un sistema de ejes local a partir de 3 coordenadas"?

A continuación se muestra el contexto de cómo se utiliza esta función:

bond_vector = get_bond_vector(atom.rval, atom.aval, atom.tval)
disp_vector = np.array(np.dot(bond_vector, self.atoms[i].local_axes))
for p in range(3):
    atom.coords[p] = self.atoms[atom.rnum].coords[p] + disp_vector[p]

La definición del vector de bonos está aquí:

def get_bond_vector(r, a, t):
    x = r * math.sin(a) * math.sin(t)
    y = r * math.sin(a) * math.cos(t)
    z = r * math.cos(a)
    bond_vector = [x, y, z]
    return bond_vector

De nuevo, la única parte que no entiendo es # get local axis system from 3 coordinates . ¿Qué hace esa función?

Answer 1

Intentaré responder a la pregunta en negrita desde una perspectiva matemática.

Básicamente, las dos primeras líneas de get_local_axes son para construir los vectores $\vec r_{12}$ y $\vec r_{23}$ . La siguiente línea construye el producto cruzado $$\frac{\vec r_{23}\times\vec r_{12}}{|\vec r_{23}\times\vec r_{12}|},$$

El siguiente vector es $$\frac{\vec r_{12}\times(\vec r_{23}\times\vec r_{12})}{|\vec r_{12}\times(\vec r_{23}\times\vec r_{12})|},$$

Se trata de un vector perpendicular a $\vec r_{12}$ y se encuentra en el plano donde $\vec r_{12}$ y $\vec r_{23}$ mentira.

Las siguientes líneas definen los tres ejes del sistema. $\vec e_z$ es sólo $\vec r_{12}$ normalizado, $\vec e_y$ se elige de manera que $yz$ es el plano donde se sitúan los tres puntos (átomos). y $\vec e_x$ se elige de manera que $(\vec e_x, \vec e_y, \vec e_z)$ forma un conjunto ortonormal y son los versores de un sistema de coordenadas de la derecha. Este sistema es el sistema de ejes local.

Este es el significado de "obtener un sistema de ejes local a partir de 3 coordenadas".