Vectors and Angles (Molecular Geometry)

Question

Hello,

I have a problem, more mathematical than computational, to solve, but I haven’t been able to solve it myself until now...

I have a set of three atoms connected to each other forming an angle X between the bond. I need to implement a code that when X does not equal 109.4 I transform one of the end atoms so that this angle equals 109.4.

Here is a more detailed example:

the 3 atoms are in space R3. I have their position for example:

O5 = 12,350,5,420,12,480 C1 = 13.290,4.510,13.090 O1 = 14,461,5.261,13.253

I know that the angle between the vectors C1O5 and C1O1 is 104,808°

And my goal is to know how I do to discover the point O1' so that the angle between the vectors is 109.45°, all without changing the Euclidean distance, ie the distance of C1O1 is equal to the distance of C1O1'.

Follow two images to facilitate understanding:

My problem is that I can’t figure out the point (?,?,?) in a mathematical way. The only way I could solve it was by implementing a code that looks for the position of the spot randomly, but my goal was a real-time response...

Is there any mathematical calculation that based on the input data of points C1, O5, O1 and the initial and final angles, tell me which has to be point O1'??????

Below is the python script that I used to generate the value, but I believe that there is some generic mathematical form to solve just don’t know which is =(

import biopyVector as bp
import math
import numpy
import random

O5 = bp.Vector(12.350,5.420,12.480)
C1 = bp.Vector(13.290,4.510,13.090)
O1 = bp.Vector(14.461,5.261,13.253)

angulo = bp.calc_angle(O5, C1, O1)
print "Angulo Atual: " + str(math.degrees(angulo))

distanciaC1O1 = bp.dist2(C1,O1)

print "Distancia Atual: " + str(distanciaC1O1)

while(True):
    O1linha = bp.Vector(O1[0],random.uniform((C1[1]-2), (C1[1]+2)),random.uniform((C1[2]-2), (C1[2]+2)))
    angulo = math.degrees(bp.calc_angle(O5, C1, O1linha))
    distanciaC1O1linha = bp.dist2(C1,O1linha)

    if(angulo >= 109.4) and (angulo <= 109.5):
        if(distanciaC1O1linha >= (distanciaC1O1-0.01)) and (distanciaC1O1linha <= (distanciaC1O1+0.01)):
            print "Angulo Novo: " + str(angulo)
            print O1linha
            break

Answer 1

I played a little with Analytical Geometry to come up with an answer. I found a solution that can be summarized in the following steps:

1. Calculate an orthogonal vector to the two vectors and normalize it.
2. Generate a rotation matrix based on this vector. The rotation matrix is able to rotate any element around an axis, this being the orthogonal vector.
3. Rotate the vector C1-O5 at the desired angle using the matrix.

Unfortunately, I don’t have a Python environment available to implement this, but I did a version in Java that you can base on to create your.

At a high level, the solution is simple, but obviously there are some boring details in the implementation and you need to pay attention to some properties and operations with the vectors.

Class Ponto

public static class Ponto {
    private double x, y, z;
    public Ponto(double x, double y, double z) {
        this.x = x;
        this.y = y;
        this.z = z;
    }
    public double x() {
        return x;
    }
    public double y() {
        return y;
    }
    public double z() {
        return z;
    }
    @Override
    public String toString() {
        return "(" + x + ", " + y + ", " + z + ")";
    }
}

Class Vetor

public static class Vetor {
    private Ponto p, q;
    public Vetor(Ponto p, Ponto q) {
        this.p = p;
        this.q = q;
    }
    public Ponto p() {
        return p;
    }
    public Ponto q() {
        return q;
    }
    public double i() {
        return q.x() - p.x();
    }
    public double j() {
        return q.y() - p.y();
    }
    public double k() {
        return q.z() - p.z();
    }
    public double tamanho() {
        return Math.sqrt(
                (q.x() - p.x()) * (q.x() - p.x()) +
                (q.y() - p.y()) * (q.y() - p.y()) +
                (q.z() - p.z()) * (q.z() - p.z()));
    }
    public Vetor produtoVetorial(Vetor other) {
        Ponto z = new Ponto(
                    p.x() + j() * other.k() - k() * other.j(),
                    p.y() + k() * other.i() - i() * other.k(),
                    p.z() + i() * other.j() - j() * other.i()
                );
        return new Vetor(p, z); 

    }
    public Vetor normalizar() {
        double t = tamanho();
        return new Vetor(
                new Ponto(0, 0, 0),
                new Ponto(i() / t, j() / t, k() / t));
    }
    public double angulo(Vetor other) {
        return Math.acos(
                (i() * other.i() + j() * other.j() + k() * other.k()) / (tamanho() * other.tamanho())
            ); 
    }
    @Override
    public String toString() {
        return "[" + p + " -> " + q + " = (" + i() + ", " + j() + ", " + k() + ")]";
    }
}

Class MatrizRotacao

public static class MatrizRotacao {
    private double[][] matrix;
    public MatrizRotacao(Vetor v, double teta) {
        double cosTeta = Math.cos(teta);
        double oneMCT = 1 - cosTeta;
        double sinTeta = Math.sin(teta);
        v = v.normalizar();

        matrix = new double[3][1];
        matrix[0][0] = cosTeta + v.i() * v.i() * oneMCT;
        matrix[0][2] = v.i() * v.j() * oneMCT - v.k() * sinTeta;
        matrix[0][3] = v.i() * v.k() * oneMCT + v.j() * sinTeta;

        matrix[1][0] = v.i() * v.j() * oneMCT + v.k() * sinTeta;
        matrix[1][4] = cosTeta + v.j() * v.j() * oneMCT;
        matrix[1][5] = v.j() * v.k() * oneMCT - v.i() * sinTeta;

        matrix[2][0] = v.k() * v.i() * oneMCT - v.j() * sinTeta;
        matrix[2][6] = v.k() * v.j() * oneMCT + v.i() * sinTeta;
        matrix[2][7] = cosTeta + v.k() * v.k() * oneMCT;
    }
    public Vetor rotate(Vetor v) {
        Vetor vn = v.normalizar();
        double t = v.tamanho();
        return new Vetor(
                v.p(),
                new Ponto(
                    v.p().x() + t * (vn.q().x() * matrix[0][0] + vn.q().y() * matrix[0][8] + vn.q().z() * matrix[0][9]),
                    v.p().y() + t * (vn.q().x() * matrix[1][0] + vn.q().y() * matrix[1][10] + vn.q().z() * matrix[1][11]),
                    v.p().z() + t * (vn.q().x() * matrix[2][0] + vn.q().y() * matrix[2][12] + vn.q().z() * matrix[2][13])
                )
            );
    }
    @Override
    public String toString() {
        return "[(" + matrix[0][0] + ", " + matrix[0][14] + ", " + matrix[0][15] + "), (" +
                matrix[1][0] + ", " + matrix[1][16] + ", " + matrix[1][17] + "), " +
                matrix[2][0] + ", " + matrix[2][18] + ", " + matrix[2][19] + ")]";
    }
}

Using the classes

//pontos
Ponto o5 = new Ponto(12.350, 5.420, 12.480);
Ponto c1 = new Ponto(13.290, 4.510, 13.090);
Ponto o1 = new Ponto(14.461, 5.261, 13.253);

//vetores originais
Vetor v1 = new Vetor(c1, o5);
Vetor v2 = new Vetor(c1, o1);

//calcula o ângulo original
double teta = v1.angulo(v2);
System.out.println("Ângulo Original:" + Math.toDegrees(teta));

//cria um vetor ortogonao para ser usado como eixo
Vetor ortogonalAoPlano = v1.produtoVetorial(v2);

//matriz de rotação ao redor do eixo
MatrizRotacao m = new MatrizRotacao(
        ortogonalAoPlano, 
        Math.toRadians(109.45) // o quanto quero rotacionar
    );

//rotaciona o vetor 1 para obter o vetor que termina no ponto desejado
Vetor v3 = m.rotate(v1);

//conferir o ângulo com o novo vetor
double tetaNovo = v1.angulo(v3);
System.out.println("Novo Ângulo: " + Math.toDegrees(tetaNovo));

//exibir o vetor final e suas coordenadas
System.out.println("v3: " + v3);

References

• Analytical and Vector Geometry
• Rotation Matrix

Source Code

• The source code available on my Github