Bound angle between cuboids

[alexas]

Dear Experts,
Angles have worn me down yet again.
Now I am trying to set up the required connection between my molecules. So Molecules are represented by cuboids which are described by three points (bottom center, center, top center) and rotation (in AxisAngle notion, x,y,z, theta).

The molecules(cuboids) are randomly rotated but should be connected to each other at ~110 deg angle if we draw the plane between the n-th molecule beginning, the connection point with n+1 molecule, and n+1 -th moelcule endpoint.

I tried two approach, both does not bring expected results. If any talented man with good visual thinking could help, I would be grateful.

//function takes the rotation input of the previous molecule
ExponentialMap
StructureGenerator:: _random_rotation_bound(ExponentialMap& prot) const
{
    const auto m_pi = getPiValue();
    auto prng = RandomGenerator(std::chrono::high_resolution_clock::now().time_since_epoch().count());

    float px, py, pz, theta;
    prot.get_rotation(px, py, pz, theta);
// here we create the random rotation axis, with intended 70 degrees rotation

    auto rand = prng.get_float(0.0f, 1.0f);
    ExponentialMap nrot = {prng.get_float(0.0f, 1.0f),prng.get_float(0.0f, 1.0f), prng.get_float(0.0f, 1.0f), 1.22173f};
    if (rand > 0.5) nrot = {prng.get_float(0.0f, 1.0f),prng.get_float(0.0f, 1.0f), prng.get_float(0.0f, 1.0f), -1.22173f};

// combining rotations (rotation matrix multiplication)

    prot+=nrot;

// now making random rotation with random angle
    ExponentialMap rot = ExponentialMap(prng.get_float(0.0f, 1.0f), prng.get_float(0.0f, 1.0f), prng.get_float(0.0f, 1.0f), prng.get_float(0.0f, 2 * m_pi));

// and combining them again
    rot+=prot;

    return rot;
}

second try:

//function takes n-th cuboid's beginning point and intersection point (first vector)
ExponentialMap
StructureGenerator:: _random_rotation_bound(float fx, float fy, float fz,
                                            float sx, float sy, float sz) const
{
    // calculation of the first vector and its normalization
    float dfx = sx - fx;
    float dfy = sy - fy;
    float dfz = sz - fz;
    auto fnorm = 1.0f / std::sqrt(dfx * dfx + dfy * dfy + dfz * dfz);
    dfx *= fnorm;
    dfy *= fnorm;
    dfz *= fnorm;

    // Second vector
    float dx = std::sin(1.9111355f);
    float dy = std::cos(1.9111355f);
    float dz = -(fz);
    if (dz > 0)
    {
        dx = std::cos(1.9111355f);
        dy = std::sin(1.9111355f);
    }

    // creating expmap with first vector and random rotation
    const auto m_pi = getPiValue();
    auto prng = RandomGenerator(std::chrono::high_resolution_clock::now().time_since_epoch().count());
    ExponentialMap rot = ExponentialMap(dfx, dfy, dfz, prng.get_float(0.0f, 2 * m_pi));

    //rotating the second vector
    float refx, refy, refz;
    rot.rotate(refx, refy, refz, dx, dy, dz);

    return ExponentialMap(refx, refy, refz, prng.get_float(0.0f, 2 * m_pi));

[keskiverto]

Cuboid "molecule"? I've heard of molecules that consist of atoms and bonds.

Where is the "beginning" in a cuboid?
Where is the "connection point" in a cuboid?
Where is the "end point" in a cuboid?

Nevertheless, there are three points (beginning, shared connection, end).
Therefore, two vectors: from connection to beginning and from connection to end.
Angle between two vectors ... you do know how to compute?


On the other hand, if you have one cuboid -- two points -- you have axis defined by them. There are infinite number of vectors that start from (connection) point on the axis and all have same angle from the axis.

[alexas]

Thank you for your attention,
I attached the schematics, it will explain better: https://ibb.co/3SHydmH
Cuboid is the approximation of the monomer.

I don't know how to have the constant bound angle between vectors while they should be able to have the random rotation.

It's not even easy to explain, and I myself don't have SUCH good understanding that could it explain easily. Just imagine bunch of randomly rotated cuboids, but all they must be connected by the certain angle in the plane which is made by that points (like in the picture)

[keskiverto]

A is beginning of first cuboid. B is end of first cuboid. C is end of second cuboid.
It should be trivial to move (translate) the second cuboid to have its beginning on B.

Vector vA = B-A. Vector vB = C-B. You want vB at certain angle compared to vA.

Pick vector vN that is orthogonal to vA. There are infinite number to choose from. Perhaps pick one and then rotate it by random amount around vA. The result should still be orthogonal to vA.

Put vA and vN on B. Rotate vA around vN by (Pi-110 degrees). Make vA as long as vB. The new location for C is B+vA.

[alexas]

yes , that's exactly what I was doing. There is no problem to connect cuboids just in that angle.

The problem is that I need to rotate additionally each cuboid in the remaining directions randomly yet preserving that bound angle that should remain as constant.

I admit again, I don't get the whole concept very well ( I am not a chemist). However I know what result should look like. You can see in the picture below and you see the chains goes straight only, while actually they should bend at some point, i.e. should not go so straight all the time.
https://ibb.co/0qgXLWS

I admit that my explanation and purpose may sound vague.

[againtry]

The (local) plane formed by the 2 end points and common connection point remains unchanged as long as the intersection angle remains unchanged. That's effectively what @kbw describes and you agree with.

The only other possible rotations are:
a) rotations of that plane about any number of 3D axes
b) rotations of each cuboid about it's own axis

[alexas]

Thank you very much. Your a) statement just opened my eyes.

The whole algorithm/my thinking is wrong. The third point should be random! I am not good in geometry obviously, i have lack of visual imagination for that thing.

As far as I understand now:
- I randomly choose the third point c,
- we create equation of a plane and extract vector c-a
- an we rotate it by
- and we possibly can extract a-b-c vector from that equation of plane
-and then this vector should be translated to AxisAngle notion where angle could be calculated again from the equation of plane to perpendicular y-z coordinate plane.

However, no actual experience in doing that. Maybe someone with more abilities could guide me through (in case someone easily understand this), would be most appreciated.

[alexas]

Hello again,

Here is my code that would address the problem. It gives the actual angles between the cuboids, but they just go round and round, i.e. randomness disappears. If anyone has any advice, kindly please.
Here is the result image: https://ibb.co/KmC53BD

//function gets the coordinates of the first cuboid beginning and the end 
//selecting random 3rd point
auto prng = RandomGenerator(std::chrono::high_resolution_clock::now().time_since_epoch().count());
    float x3 = prng.get_float(-1.0f, 1.0f);
    float y3 = prng.get_float(-1.0f, 1.0f);
    float z3 = prng.get_float(-1.0f, 1.0f);

//creating the new rotation axis with the intended angle of ~110 degree to vector 1->2
   auto cp = get_cross_product(x1, y1, z1,
                                x2, y2, z2,
                                x3, y3, z3);

   
 ExponentialMap r (cp[0], cp[1], cp[2], 1.91099476f);

    float vec23x = x3 - x2;
    float vec23y = y3 - y2;
    float vec23z = z3 - z2;

//generating the new AxisAngle for the intended rotation of 3cuboid
    float newvec23x, newvec23y, newvec23z;
    r.rotate(newvec23x, newvec23y, newvec23z, vec23x, vec23y, vec23z);

    float vec12x = x3 - x2;
    float vec12y = y3 - y2;
    float vec12z = z3 - z2;

    auto angle = acos((vec12x * newvec23x + vec12y * newvec23y + vec12z * newvec23z) /
                              ((std::sqrt(vec12x * vec12x + vec12y * vec12y + vec12z * vec12z)) *
                      (std::sqrt(newvec23x * newvec23x + newvec23y * newvec23y + newvec23z * newvec23z))));

    return ExponentialMap(newvec23x, newvec23y, newvec23z, (3.14136125f- angle));

[againtry]

I suggest you mark your imgur picture with the 3 x,y,z points.

[lastchance]

@alexas,
Unless your clock is really high speed I suspect you are seeding your random-number generator with the same number every time, and consequently you end up with the same three "random" numbers every time.

Seeding operations based on the clock should be done only once.

[keskiverto]

    float vec23x = x3 - x2;
    float vec23y = y3 - y2;
    float vec23z = z3 - z2;
    float newvec23x, newvec23y, newvec23z;
    r.rotate(newvec23x, newvec23y, newvec23z, vec23x, vec23y, vec23z);

C++ has classes and operator overloading, which allow encapsulation.
Those take some effort to set up, but then the "program" has less repetition, which could lead to typos and does obscure the logic.

    MyVec vec23 = p3 - p2;
    MyVec newvec23;
    r.rotate( newvec23, vec23 );

Or, perhaps even:

vector<Point> residue = ...;
Transform t = ...; // rotates around pivot and translates
for ( auto& point : residue ) {
  point = t.rotation * ( point - t.pivot ) + t.pivot + t.translation;
}

[alexas]

thank you very much for your attention and input.

@againtry, https://ibb.co/NVv0CzS if I got you correctly. We generate the first cuboid. Its bottom center is point 1, its top center is point 2. The next cuboid will start at its top center, and we are looking randomly for point 3, at which the next cuboid should position its top center, so creating the needed angle between the cuboids.

@lastchance, good remark, but I checked it, the numbers are being generated randomly with no correlation.

@keskiverto, thank you, I will certainly try to learn more advanced c++ in the near future.

There is just something wrong here. Bound angle should remain the same but the cuboid may be rotated in any other direction in the remaining axis. I tried to introduce additional randomness with the codeline below, and it changes the picture, but not as it should be.

    float rand = prng.get_float(0, 1.0f);
    ExponentialMap r (cp[0], cp[1], cp[2], 1.91099476f);
    if (rand > 0.5f)
        r = {cp[0], cp[1], cp[2], -1.91099476f};

[againtry]

@alexas
More or less what I imagined you are doing.

I'd concentrate on a single 'strand' of 2 cuboids and get the geometry of that sorted before going to long and even multiple strands.

Rotating pt 3 about the 1-2 (local) axis results in a cone - the surface of revolution line 2-3 describes. The height of the cone is L*sin(70) and radius L*sin(70), L is the length of line 2-3.

These are all constants.

The circular path of pt 3 and hence the coordinates of pt 3 can be described by a (t) parametric function where t is the angle of rotation around the circle.

From that you get the local coordinates of 3 relative to the line 1-2 axis, which in turn can be transformed to global coordinates.

Of course there are other molecular constraints that stop pt 3 positions mapping a complete circle. See foldit.it if you haven't already.

[keskiverto]

Perhaps,
* Create a "pre-template" cuboid. Beginning is at (0,0,0), end is at (x,0,0), and other points as they should.

* Create a "template" cuboid by rotating the pre-template by 70 degrees around Y-axis. The beginning is on axis of rotation and thus stays on (0,0,0) and the end moves to (x*cos(70deg), x*sin(70deg), 0). The other points accordingly. The angle between -X and template is now 110 degrees

* For a new cuboid, make copy of the template.
- Then rotate the copy by random amount around X-axis
- Determine rotation R that would make X-axis point to same direction as the latest cuboid in the chain
- Rotate the copy by R
- Translate the copy by endpoint of latest cuboid. Now the beginning of the copy should be on end of the latest

[againtry]

https://imgur.com/ivSeE3P

[lastchance]

100-cube polymer, 110-degree bond angle.

This code creates an image file (stl.stl). Open directly in Windows 10, or use https://www.viewstl.com/

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <random>
#include <cmath>
#include <ctime>
using namespace std;

const double PI = 4.0 * atan( 1.0 );
using Mat = vector< vector<double> >;
mt19937 gen( time( 0 ) );
uniform_real_distribution<double> rnd( 0.0, 1.0 );

//====

struct Pt
{
   double x, y, z;
   Pt( double x = 0, double y = 0, double z = 0 ) : x( x ), y( y ), z( z ) {}
};

struct Triangle
{
   Pt v1, v2, v3;
   Triangle( const Pt &v1, const Pt &v2, const Pt &v3 ) : v1( v1 ), v2( v2 ), v3( v3 ) {}
};

struct Polymer
{
   double side;
   Pt start;
   vector< pair<Pt,Pt> > cubes;

   Polymer( double side, const Pt &start, const vector< pair<Pt,Pt> > &cubes ) : side( side ), start( start ), cubes( cubes ) {}
};

//====

Pt operator + ( const Pt &a, const Pt &b ){ return { a.x + b.x, a.y + b.y, a.z + b.z }; }
Pt operator - ( const Pt &a, const Pt &b ){ return { a.x - b.x, a.y - b.y, a.z - b.z }; }
Pt operator / ( const Pt &a, double d    ){ return { a.x / d, a.y / d, a.z / d }; }
Pt operator * ( double d   , const Pt &a ){ return { d * a.x, d * a.y, d * a.z }; }
Pt cross      ( const Pt &a, const Pt &b ){ return { a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x }; }
double dot    ( const Pt &a, const Pt &b ){ return a.x * b.x + a.y * b.y + a.z * b.z; }
double normsq ( const Pt &a              ){ return dot( a, a ); }
double len    ( const Pt &a              ){ return sqrt( normsq( a ) ); }
Pt unit       ( const Pt &a              ){ return a / len( a ); }
ostream & operator << ( ostream &out, const Pt a ){ return out << a.x << " " << a.y << " " << a.z << " "; }

Pt operator * ( const Mat &M, const Pt &v )
{
   return { M[0][0] * v.x + M[0][1] * v.y + M[0][2] * v.z,
            M[1][0] * v.x + M[1][1] * v.y + M[1][2] * v.z,
            M[2][0] * v.x + M[2][1] * v.y + M[2][2] * v.z };
}

Mat operator * ( const Mat &M, const Mat &N )
{
   Mat R( M.size(), vector<double>( N[0].size(), 0.0 ) );
   for ( int i = 0; i < M.size(); i++ )
   {
      for ( int j = 0; j < N[0].size(); j++ )
      {
         for ( int k = 0; k < N.size(); k++ ) R[i][j] += M[i][k] * N[k][j];
      }
   }
   return R;
}

//====

Mat rotationMatrix( const Pt &axis, double radians )
{
   Mat result( 3, vector<double>( 3, 0.0 ) );
   Pt n = axis / len( axis );
   double cosA = cos( radians ), sinA = sin( radians ), cosA1 = 1.0 - cosA;
   result[0][0] = cosA + n.x * n.x * cosA1;
   result[0][1] =      + n.x * n.y * cosA1 - n.z * sinA;
   result[0][2] =      + n.x * n.z * cosA1 + n.y * sinA;
   result[1][0] =      + n.y * n.x * cosA1 + n.z * sinA;
   result[1][1] = cosA + n.y * n.y * cosA1;
   result[1][2] =      + n.y * n.z * cosA1 - n.x * sinA;
   result[2][0] =      + n.z * n.x * cosA1 - n.y * sinA;
   result[2][1] =      + n.z * n.y * cosA1 + n.x * sinA;
   result[2][2] = cosA + n.z * n.z * cosA1;
   return result;
}

//---

void axisAngle( const Mat &R, Pt &axis, double &angle )
{
   angle = acos( 0.5 * ( R[0][0] + R[1][1] + R[2][2] - 1 ) );
   double s = sin( angle );
   if ( abs( s ) < 1.0e-6 )       // Edge case sin(theta)=0: painful!
   {
      // Find largest diagonal element and use that column
      int i = 0;
      if ( R[1][1] > R[i][i] ) i = 1;
      if ( R[2][2] > R[i][i] ) i = 2;
      axis = Pt( R[0][i], R[1][i], R[2][i] );
      if ( i == 0 ) axis.x += 1;
      if ( i == 1 ) axis.y += 1;
      if ( i == 2 ) axis.z += 1;
      // Note: axis must be normalised below; it currently has length 2n_i
   }
   else
   {
      // Antisymmetric part should determine axis without problems
      axis = ( 0.5 / s ) * Pt( R[2][1] - R[1][2], R[0][2] - R[2][0], R[1][0] - R[0][1] );
   }
   axis = axis / len( axis );     // Normalise to length 1
}

//====




//====

class Stl;

class Shape
{
public:
   virtual void addToPlot( Stl &stl ) = 0;
};

//====

class Stl
{
   vector<Triangle> triangles;
public:
   void add( Shape *S ) { S->addToPlot( *this ); }
   void addTriangle( const Pt &v1, const Pt &v2, const Pt &v3 ) { triangles.push_back( Triangle( v1, v2, v3 ) ); }
   void addRectangle( const Pt &v1, const Pt &v2, const Pt &v3, const Pt &v4 ) { addTriangle( v1, v2, v3 );   addTriangle( v1, v3, v4 ); }
   void insert( const Stl &stl ) { triangles.insert( triangles.end(), stl.triangles.begin(), stl.triangles.end() ); }
   void translate( const Pt &v ) { for ( Triangle &t : triangles ) t = Triangle( t.v1 + v, t.v2 + v, t.v3 + v ); }
   void scale( double s ) { for ( Triangle &t : triangles ) t = Triangle( s * t.v1, s * t.v2, s * t.v3 ); }
   void scale( double s, const Pt &centre ) { translate( -1 * centre );   scale( s );   translate( centre ); }
   void rotate( const Mat &M ) { for ( Triangle &t : triangles ) t = Triangle( M * t.v1, M * t.v2, M * t.v3 ); }
   void rotate( const Mat &M, const Pt &centre ) { translate( -1 * centre );   rotate( M );   translate( centre ); }
   void draw( const string &filename );
};

//---

void Stl::draw( const string &filename )
{
   ofstream out( filename );
   out << "solid\n";

   for ( Triangle t : triangles )
   {
      Pt n = cross( t.v2 - t.v1, t.v3 - t.v1 );
      n = n / len( n );

      out << "   facet normal " << n << '\n';
      out << "      outer loop\n";
      out << "         vertex " << t.v1 << '\n';
      out << "         vertex " << t.v2 << '\n';
      out << "         vertex " << t.v3 << '\n';
      out << "      endloop\n";
      out << "   endfacet\n";
   }

   out << "endsolid\n";
}

//====

class Cuboid : public Shape
{
   Pt centre;
   Pt side1, side2, side3;

public:
   Cuboid( const Pt &c, double L ) : centre( c ), side1( Pt( L, 0, 0 ) ), side2( Pt( 0, L, 0 ) ), side3( Pt( 0, 0, L ) ) {}
   Cuboid( const Pt &c, const Pt &s1, const Pt &s2, const Pt &s3 ) : centre( c ), side1( s1 ), side2( s2 ), side3( s3 ) {}

   void addToPlot( Stl &stl );
};

//---

void Cuboid::addToPlot( Stl &stl )
{
   Pt v1 = centre - 0.5 * ( side1 + side2 + side3 );
   Pt v2 = v1 + side1, v3 = v2 + side2, v4 = v3 - side1;
   Pt v5 = v1 + side3, v6 = v2 + side3, v7 = v3 + side3, v8 = v4 + side3;
   stl.addRectangle( v1, v2, v6, v5 );
   stl.addRectangle( v2, v3, v7, v6 );
   stl.addRectangle( v3, v4, v8, v7 );
   stl.addRectangle( v4, v1, v5, v8 );
   stl.addRectangle( v1, v4, v3, v2 );
   stl.addRectangle( v5, v6, v7, v8 );
}

//====

void plot( const string &filename, const Polymer &P )
{
   Stl stl;
   Pt start = P.start;
   for ( auto &C : P.cubes )
   {
      Pt side1 = C.first;               // axis
      Pt side2 = C.second;              // perpendicular edge
      Pt side3 = cross( side1, side2 ); // other perpendicular edge
      side1 = P.side * unit( side1 );
      side2 = P.side * unit( side2 );
      side3 = P.side * unit( side3 );

      Pt centre = start + 0.5 * side1;
      Cuboid cub( centre, side1, side2, side3 );
      stl.add( &cub );
      start = start + side1;
   }
   stl.draw( filename );
}

//====

int main()
{
   const int Ncubes = 100;
   const double delta = ( 180 - 110 ) * PI / 180.0;
   Pt ex = { 1, 0, 0 }, ey = { 0, 1, 0 };

   Polymer P( 1.0, { 0, 0, 0 }, vector< pair<Pt,Pt> >{ { ex, ey } } );
   for ( int i = 1; i < Ncubes; i++ ) 
   {
      Pt old  = P.cubes.back().first;   // previous axis

      Pt perp = cross( old, ex );       // random perpendicular axis
         if ( len( perp ) < 0.01 ) perp = cross( old, ey );
      Mat R = rotationMatrix( old, rnd( gen ) * 2 * PI );
      perp = R * perp;

      old = old / len( old );   perp = perp / len( perp ); // normalise
      Pt axis = cos( delta ) * old + sin( delta ) * perp;  // new axis

      Pt edge = sin( delta ) * old - cos( delta ) * perp;
      R = rotationMatrix( axis, rnd( gen ) * 2 * PI );
      edge = R * edge;                                     // new edge

      P.cubes.push_back( { axis, edge } );
   }

   plot( "stl.stl", P );
}

[againtry]

Small point. Add #include <fstream> and it works.

[lastchance]

Small point. Add #include <fstream> and it works.

Thanks, @Againtry. I was struggling to get the code under the character limit and took it out, forgetting that I needed it.

Code edited so that it should work.

The polymer gets very contorted unless you increase the bond angle (say, to 150 deg).


Wonder if they'll use this to develop new vaccines to Covid-19 - workplace has just slapped another facemask mandate in place.

[againtry]

Who knows, we are now at variant Omicron so when we get to the end of the road with variant Omega the Lastchance Vaccine will save the world.

[keskiverto]

Wonder if they'll use this to develop new vaccines to Covid-19

I hope not. There are plenty of both open source and proprietary molecular modeling applications and suites already, which do a more detailed job. Force fields, molecular docking, dynamic simulations, etc.

Besides, if these cuboids represent chain of residues of a protein, they are only an abstraction of the main chain. One could go more plain by tracing just the alpha carbons. However, as already stated, there are steric hindrances -- all random rotations are not feasible, and sidechains of amino acid residues do play important role in function; most molecular interactions are via atoms of sidechains.