Bound angle between cuboids

[lastchance]

I presume the OP means "bond angle", not "bound angle", doesn't he? (Changed my replies to assume so.)

Also, it's presumably just polythene/polyethylene (saturated carbon chain) if all the angles are the same and equal to that for a tetrahedral cell.

[againtry]

bond and bound sort of mean the same thing when it comes to molecule chains connecting.

My guess is @OP is not concerned about what’s inside each cuboid, simply how the two connect knowing the allowable angle of 110 degrees which is determined by other chemistry means.

The random rotation at each ‘joint’ is a little ambiguous. It could be to simulate each joint finding a lock-in position for 2 cuboids or building a whole random string of 100 cuboids.

[alexas]

I greatly appreciate your help in this complicated task (for me) and advises. @lastchance becomes literally the last chance for me when similar obstacle is present. I can see that it should work, so will start to analyze it and adapt to my code. The intersection check is necessary so any chain that is physically impossible will be not generated.

I am coming from the field of Radiation Physics. And in the middle of my Phd which is about the simulation of the radiation effect on the polymer materials, I was frustrated that there is no any suitable chemistry model for my needs where I could replicate the polymer growth upon irradiation, i.e. I knew everything what I need for my radiation part, but what happens next, inside? And so I decided to create the simple coarse polymer growth model.

With zero skill in C++ and the knowledge of chemistry of 8 grade high school.

I have, of course, the support of actual chemist, but it is limited in time and sometimes I struggle with certain concepts, especially in 3d space ( I never had to use the similar knowledge previously). And to create really fast code (mandatory for MC simulation) that would continuously work with thousands of cuboids checking their interactions and intersections is a challenge itself for me too.

So no, I won't make any new vaccines (afterall it seems that Mother Nature will solve everything just by taking out the weakest - like it always did, no matter what we do)

[againtry]

@alexas

afterall it seems that Mother Nature will solve everything just by taking out the weakest - like it always did, no matter what we do

You've got to be joking with creationist crap like that!

[alexas]

:) i thought it's more about stoicism or at least fatalism :D just pls change Mother Nature = Laws of Physics, etc.

[alexas]

It seems that is really works!

However, I generated several chains, and it still sometimes gives very strange angles along it.
https://ibb.co/c1cFtwh

I also really want to learn here, so just few more question to be sure I understand what's happening:

        Pt perp = cross( old, ex );       // random perpendicular axis
        if ( len( perp ) < 0.01 ) perp = cross( old, ey );

We have previous axis, so we just look for the perpendicular axis for it either along x axis or along y axis? And if previous axis are almost perpendicular to any of the axis, we choose another one?

      Mat R = rotationMatrix( old, rnd( gen ) * 2 * PI );
      perp = R * perp;

So here we rotate that perpendicular axis randomly? Why ?

My function should return the new ExponentialMap only, so I return

	
	Pt axis = cos( delta ) * old + sin( delta ) * perp;  // new axis
	and random angle prng.get_float(0.0f, 2 * m_pi)

[alexas]

Ok I think I made mistake somewhere as your code with the last rotation = 0 gives nice angles of the cuboids while mine just go straight on the top of each other.

I think the edge maybe is the actual new rotation axis? What we actually return here?

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

[lastchance]

However, I generated several chains, and it still sometimes gives very strange angles along it.

They should definitely be cubes (all-equal-sides) - has your STL viewer got Covid-19?

     Pt perp = cross( old, ex );       // random perpendicular axis
        if ( len( perp ) < 0.01 ) perp = cross( old, ey );

We have previous axis, so we just look for the perpendicular axis for it either along x axis or along y axis? And if previous axis are almost perpendicular to any of the axis, we choose another one?

We need an axis perpendicular to the previous one (so that sin(70) worth of the original length can be directed that way). One easy way of generating a perpendicular vector is to use a cross product (with anything). The default just forms the cross product with the unit x-vector. This works unless the original was nearly lined up with the unit x-vector, in which case the cross product would give 0. In that circumstance it switches to a cross product with the unit y vector instead.

      Mat R = rotationMatrix( old, rnd( gen ) * 2 * PI );
      perp = R * perp;

So here we rotate that perpendicular axis randomly? Why ?

You asked for a new direction anywhere in a 70 degree cone. The "flat" of that cone is a circle (2.PI radians). So I rotated it by a random amount in [0,2.PI)



My function should return the new ExponentialMap only
I don't use an "ExponentialMap" (presumably a rotation-axis vector and angle), because it is more convenient for me to draw cubes if I know two of their side vectors and don't have to recover them later. Your method will use less storage for cubes, but will have to work harder to draw them later.

I think the edge maybe is the actual new rotation axis?

No. It is the side of a cube. You are storing cubes in a different way to me.
I suspect that if you create a rotation matrix with three COLUMNS [ Pt axis, Pt edge, Pt third ] (where third is the cross product of axis and edge) and then call the routine void axisAngle() on that matrix you would get the (different) axis and angle that you require to describe your cubes.

In general, a matrix effecting any linear transformation (which doesn't have to be a rotation) will have columns which are the transformed x, y and z unit vectors (the sides of an original, unrotated, cube).

[lastchance]

return ExponentialMap(axis.x, axis.y, axis.z, 0);

No! Please read my previous post. I repeat: we do not store cubes in the same way.

My "axis" is through the centre of a cube (in the direction of the polymer chain).
Your "axis" is an axis of rotation which defines the orientation of the cube (I think!).

I believe - but am not in a position to check - that if you first create a matrix with COLUMNS

    (  |     |       |   )
M = ( axis  edge   third )
    (  |     |       |   )

where Pt third = cross( axis, edge ); and then apply

axisAngle( const Mat &M, Pt &YOURaxis, double &YOURangle );
return ExponentialMap( YOURaxis.x, YOURaxis.y, YOURaxis.z, YOURangle);

you will probably get what you want. (But I can't check that.)

[alexas]

Thanks again. I'll need to mention this forum in my dissertation acknowledgements section.
I admit I still have very vague understanding how these rotations works but at the moment I would just need it to work.
Yes, my cube axis is an axis of rotation which defines the orientation of cube. The cube's starting position is based on the previous cube's top face.
So i implemented the proposed changes and set the edge rotation matrix angle to 0 with hope to see the the similar picture as in your code. However, I get smth like this: https://ibb.co/dLhFqVw

The function looks like this - maybe you can spot some probably obvious error? Maybe it is wrong that my previous axis is the orientation axis of the previous cuboid, and some conversion is needed ?

 auto prng = RandomGenerator(std::chrono::high_resolution_clock::now().time_since_epoch().count());

    const float delta = 1.22173f; //( 180 - 110 ) * PI / 180.0 = 1.22173
    const float pi = getPiValue();

    float px, py, pz, ptheta;
    previous_rot.get_rotation(px, py, pz, ptheta);
    Pt previous = {px, py, pz};   // previous axis  - (previous rotation axis!

    Pt ax = {1, 0, 0}, ay = {0, 1, 0}; // random perpendicular axis
    Pt perp = cross(previous, ax);
    if (len(perp) < 0.01) perp = cross(previous, ay);

    Mat R = rotationMatrix(previous, prng.get_float(0.0f, 1.0f));
    perp = R * perp;
    previous = previous / len(previous);//normalization
    perp = perp / len(perp);
    Pt axis = cos(delta) * previous + sin(delta) * perp;  // new axis

    Pt edge = sin(delta) * previous - cos(delta) * perp;
    R = rotationMatrix(axis, 0); // random angle
    edge = R * edge;

    Pt third = cross(axis, edge);

    /*(      |     |       |   )
      M = ( axis  edge   third )
      (      |     |       |   ) */

    Mat M = {{axis.x, edge.x, third.x},
             {axis.y, edge.y, third.y},
             {axis.z, edge.z, third.z}};

    Pt monomer_axis;
    float new_theta;
    axisAngle(M, monomer_axis, new_theta);

    return ExponentialMap(monomer_axis.x, monomer_axis.y, monomer_axis.z, new_theta);

[lastchance]

There are a few differences.

The most significant is that I don't really know exactly how your "Exponential map" defines the orientation of the cube. You can only try all possible permutations of axis, edge and third.


Next:
Mat R = rotationMatrix(previous, prng.get_float(0.0f, 1.0f));
Don't you want it between 0 and 2.PI radians, not 0 and 1?


It is unfortunate that you are using float rather than double.

    float new_theta;
    axisAngle(M, monomer_axis, new_theta);

In my code new_theta would have been a double. It depends how you are using my routines.

[lastchance]

@Alexas,
I have no way of testing your code, but does this work?

   auto prng = RandomGenerator(std::chrono::high_resolution_clock::now().time_since_epoch().count());

   const float delta = 1.22173f; //( 180 - 110 ) * PI / 180.0 = 1.22173
   const float pi = 3.14159;

   float px, py, pz, ptheta;
   previous_rot.get_rotation( px, py, pz, ptheta );
   Pt previous = {px, py, pz};   // previous axis  - (previous rotation axis! IS THIS UNIT LENGTH? ANY OTHER LENGTH?)

   Pt ax = {1, 0, 0}, ay = {0, 1, 0}; // unit vectors along cartesian axes
   Pt perp = cross(previous, ax);
   if (len(perp) < 0.01) perp = cross(previous, ay);

   Mat R = rotationMatrix(previous, prng.get_float(0.0f, 2 * pi));
   perp = R * perp;
   
   previous = previous / len(previous);  //normalise to length 1.0
   perp = perp / len(perp);              //normalise to length 1.0
   Pt axis = cos(delta) * previous + sin(delta) * perp;  // new axis: turned angle delta from previous (AND OF UNIT LENGTH)

   float new_theta = prng.get_float(0.0f, 2 * pi);
   return ExponentialMap( axis.x, axis.y, axis.z, new_theta );

We really could do with knowing exactly how your exponential map defines the orientation (and sewing-together) of cubes.

[alexas]

I am so sorry - those cuboids are non cubic only because I was using the different parameters regarding their size.
....and does it make any difference? Can I change the size of the cuboids parameters to see how it will look like?

If you remember my previous problem, we discussed, that cube's new position is based on the previous cube, but the rotation is absolute to the coordinate system. So exponential map of {1,0,0,1} would mean the rotated cuboid along its x axis for 1 rad, {0.5, 1, 0.5, 1.2} means 1.2 rad rotation about the axis(vector)0.5x,y,0.5z, and so on.

Yes, unfortunate typo, i fixed it now

Mat R = rotationMatrix(previous, prng.get_float(0.0f, 1.0f));

actually it kinda works, most of angles are what I would expect them to be, but some of them just getting very sharp, so don't kow why. The major problem is that I am not sure how the final result should look like.
The 2x4x8 cuboid should represent the approximation of Methacrylic acid (MAA) monomer. I was not able to find anywhere how the actual molecule chain should look like. I was just told about the bond angle "while the remaining rotations can be random" or smth like this.

[lastchance]

So exponential map of {1,0,0,1} would mean the rotated cuboid along its x axis for 1 rad, {0.5, 1, 0.5, 1.2} means 1.2 rad rotation about the axis(vector)0.5x,y,0.5z, and so on.

OK, that gives you the orientation of the cuboid. But how are successive cuboids stitched together? Can you refer me back to your previous thread from some months ago?
OK, previous thread was
http://www.cplusplus.com/forum/beginner/279658/2/#msg1208888
with a minor correction at
http://www.cplusplus.com/forum/beginner/279885/#msg1210032
Looks like they are stitched together along the rotated z axes.

Yes, unfortunate typo, i fixed it now
Mat R = rotationMatrix(previous, prng.get_float(0.0f, 1.0f));

Don't you mean
Mat R = rotationMatrix(previous, prng.get_float(0.0f, 2*pi));

[alexas]

Yes, I mean that i fixed that mistake to

Mat R = rotationMatrix(previous, prng.get_float(0.0f, 2 * pi));

The only question now remains if these lines are correct, i.e. if the previous orienational axis is suitable for it:

 float px, py, pz, ptheta;
    previous_rot.get_rotation(px, py, pz, ptheta);
    Pt previous = {px, py, pz};   // previous axis 

[lastchance]

@alexas,
This uses a ridiculous number of conversions, I'm afraid, because I only have experience of combining rotations using matrices. I've tried to get the vector along the monomer by converting your exponential map to an orientation matrix and multipling the ez unit vector by that, so giving a vector along the monomer. Then I have rotated this vector by 70 deg (180-110). Then created a new orientation matrix. Finally I've converted back to a new exponential map.

So, give or take some minor compile errors (which I have no way of checking) does this work?

   auto prng = RandomGenerator(std::chrono::high_resolution_clock::now().time_since_epoch().count());

   const float delta = 1.22173f; //( 180 - 110 ) * PI / 180.0 = 1.22173
   const float pi = 3.14159f;
   Pt ex = { 1, 0, 0 }, ey = { 0, 1, 0 }, ez = { 0, 0, 1 };

   float px, py, pz, ptheta;
   previous_rot.get_rotation( px, py, pz, ptheta );
   Pt previous = { px, py, pz };                 // previous ROTATION axis
   // Convert it to vector which is the rotated z axis
   Mat R = rotationMatrix( previous, ptheta );
   previous = R * ez;                            // Should be the alignment of the monomer


   Pt perp = cross(previous, ex);
      if (len(perp) < 0.01) perp = cross(previous, ey);

   Mat R = rotationMatrix( previous, prng.get_float( 0.0f, 2 * pi ) );
   perp = R * perp;
   previous = previous / len(previous);  // normalization to length 1
   perp = perp / len(perp);              // ditto
   Pt axis = cos(delta) * previous + sin(delta) * perp;  // new rotated z axis

   Pt edge = sin(delta) * previous - cos(delta) * perp;
   R = rotationMatrix( axis, prng.get_float( 0.0f, 2 * pi ) );
   edge = R * edge;

   Pt third = cross(axis, edge);

   /*        Rot(ex)  Rot(ey)  Rot(ez)
          (      |      |         |   )
      M = (    edge   third     axis  )
          (      |      |         |   )    */

   Mat M = { { edge.x, third.x, axis.x },
             { edge.y, third.y, axis.y },
             { edge.z, third.z, axis.z }};

   Pt monomer_axis;
   float new_theta;
   axisAngle( M, monomer_axis, new_theta );        // careful with float or double for new_theta

   return ExponentialMap( monomer_axis.x, monomer_axis.y, monomer_axis.z, new_theta );

[alexas]

It worked like a charm:
https://ibb.co/6rHM6D0
https://ibb.co/Yd4B684

Sir, you chose your nickname wisely, as you were actually the last chance for me to get over these obstacles. It is greatly appreciated and will not be forgotten.

[lastchance]

Ah, I'm glad that worked.

The method at the moment is long-winded, because it involved: convert the rotation to a matrix, create three new axes, create a new rotation matrix, convert it back to a rotation. Probably there's a much more efficient way going straight from rotation to rotation, but it's unfamiliar to me.

On the other hand, if it works, maybe best to leave it alone.

Good luck!

[alexas]

If I may again to ask hopefully the small question:

- Could I use your solution as well to connect two already existing chains, in such way:

1) supplying the rotation of chain's A end to the function above, just setting zero angles in R
matrix instead of random numbers;
2) the output is used as input to rotate the whole chain B, so it's first cuboid would be at the
desired angle with A's end cuboid

It seems to work in most cases, yet sometimes I still get sharp angles. However, it might be some other problems in the algorithm, so I am not sure.

[lastchance]

At the end of the first chain calculate the next rotation matrix in the usual way. Call this matrix A.

Find the matrix corresponding to the orientation of the first cuboid in the second chain. Call this matrix B.

Find the rotation matrix AB^-1 and apply it to the whole of the second chain. That can then be tagged to the end of the first chain.

Finding the inverse matrix B^-1 is easy, because B is a rotation matrix, and hence orthogonal. The inverse of a an orthogonal matrix is its transpose.

If you are doing rotations other than by using matrices then you will have to do the conversions first.