Skip to content

Commit

Permalink
Added derivatives of quaternion CVs
Browse files Browse the repository at this point in the history
  • Loading branch information
Gareth Aneurin Tribello committed Sep 28, 2023
1 parent 4bd2575 commit 0fe398e
Show file tree
Hide file tree
Showing 6 changed files with 206 additions and 45 deletions.
1 change: 1 addition & 0 deletions regtest/crystdistrib/rt-quaternion-derivs/Makefile
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
include ../../scripts/test.make
2 changes: 2 additions & 0 deletions regtest/crystdistrib/rt-quaternion-derivs/config
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
type=driver
arg="--plumed plumed.dat --ixyz small.xyz --timestep=0.001 --trajectory-stride=25 --dump-forces forces --dump-forces-fmt=%8.4f --debug-forces=forces.num"
55 changes: 55 additions & 0 deletions regtest/crystdistrib/rt-quaternion-derivs/deriv.reference
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
#! FIELDS time parameter q1.w q1.i q1.j q1.k
0.000000 0 -0.8575 -0.5001 -0.2945 -0.6478
0.000000 1 -0.5855 0.5974 -0.0322 -0.2183
0.000000 2 -0.3754 1.4125 0.1645 0.1055
0.000000 3 0.8913 0.4167 0.1661 0.6655
0.000000 4 0.5254 -0.4494 0.2601 0.1870
0.000000 5 0.2458 -1.0930 0.3277 -0.1733
0.000000 6 -0.0338 0.0834 0.1284 -0.0177
0.000000 7 0.0600 -0.1480 -0.2279 0.0314
0.000000 8 0.1296 -0.3195 -0.4922 0.0678
0.000000 9 0.1023 0.1058 0.0979 0.0808
0.000000 10 0.1071 -0.1481 -0.0977 0.0440
0.000000 11 0.1095 -0.3367 -0.2430 0.0159
0.000000 12 -0.2998 -0.0921 0.0094 -0.2202
0.000000 13 -0.1380 0.0712 -0.1933 -0.0442
0.000000 14 -0.0153 0.1926 -0.3429 0.0880
0.000000 15 0.1655 0.0703 0.0212 0.1231
0.000000 16 0.0918 -0.0716 0.0640 0.0320
0.000000 17 0.0357 -0.1770 0.0954 -0.0366
0.025000 0 -0.8871 -0.5108 -0.3143 -0.6879
0.025000 1 -0.5781 0.6589 0.0028 -0.2164
0.025000 2 -0.4085 1.5087 0.2204 0.0909
0.025000 3 0.9195 0.4281 0.1918 0.7056
0.025000 4 0.5152 -0.4982 0.2353 0.1819
0.025000 5 0.2757 -1.1695 0.2823 -0.1638
0.025000 6 -0.0324 0.0827 0.1225 -0.0178
0.025000 7 0.0629 -0.1606 -0.2381 0.0345
0.025000 8 0.1327 -0.3391 -0.5026 0.0729
0.025000 9 0.0939 0.0989 0.0925 0.0761
0.025000 10 0.0983 -0.1514 -0.1056 0.0417
0.025000 11 0.1085 -0.3339 -0.2491 0.0211
0.025000 12 -0.3038 -0.0931 0.0006 -0.2296
0.025000 13 -0.1301 0.0764 -0.1915 -0.0398
0.025000 14 -0.0206 0.1982 -0.3372 0.0873
0.025000 15 0.1668 0.0682 0.0223 0.1273
0.025000 16 0.0856 -0.0731 0.0650 0.0290
0.025000 17 0.0362 -0.1753 0.0990 -0.0362
0.050000 0 -0.9332 -0.4279 -0.2709 -0.7220
0.050000 1 -0.4990 0.6339 0.0280 -0.2148
0.050000 2 -0.3755 1.6251 0.2511 0.0661
0.050000 3 0.9544 0.3615 0.1775 0.7350
0.050000 4 0.4489 -0.4769 0.1930 0.1839
0.050000 5 0.2561 -1.2507 0.2759 -0.1398
0.050000 6 -0.0212 0.0664 0.0934 -0.0131
0.050000 7 0.0501 -0.1570 -0.2210 0.0309
0.050000 8 0.1193 -0.3744 -0.5270 0.0737
0.050000 9 0.0941 0.0774 0.0721 0.0743
0.050000 10 0.0779 -0.1410 -0.1026 0.0378
0.050000 11 0.0954 -0.3488 -0.2650 0.0258
0.050000 12 -0.3204 -0.0853 -0.0124 -0.2452
0.050000 13 -0.1216 0.0788 -0.1699 -0.0448
0.050000 14 -0.0254 0.2249 -0.3452 0.0812
0.050000 15 0.1511 0.0495 0.0180 0.1160
0.050000 16 0.0648 -0.0581 0.0531 0.0255
0.050000 17 0.0275 -0.1562 0.0978 -0.0295
2 changes: 2 additions & 0 deletions regtest/crystdistrib/rt-quaternion-derivs/plumed.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
q1: QUATERNION ATOMS=1,2,3
DUMPDERIVATIVES ARG=q1.* FILE=deriv FMT=%8.4f
33 changes: 33 additions & 0 deletions regtest/crystdistrib/rt-quaternion-derivs/small.xyz
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
9
8.268800 8.271100 8.270000
X 3.640000 1.382000 1.364000
X 3.503000 1.700000 1.181000
X 3.055000 0.897000 1.436000
X 3.108000 1.381000 1.357000
X 2.940000 1.310000 1.522000
X 3.821000 1.482000 1.373000
X 3.752000 1.694000 1.289000
X 3.239000 0.984000 1.287000
X 2.916000 0.937000 1.614000
9
8.268800 8.271100 8.270000
X 3.637000 1.385000 1.358000
X 3.515000 1.698000 1.180000
X 3.073000 0.890000 1.455000
X 3.126000 1.386000 1.389000
X 2.952000 1.297000 1.538000
X 3.821000 1.480000 1.384000
X 3.757000 1.681000 1.267000
X 3.246000 0.990000 1.294000
X 2.889000 0.946000 1.581000
9
8.268800 8.271100 8.270000
X 3.637000 1.379000 1.353000
X 3.526000 1.704000 1.197000
X 3.076000 0.894000 1.457000
X 3.121000 1.374000 1.375000
X 2.935000 1.297000 1.518000
X 3.813000 1.480000 1.403000
X 3.752000 1.688000 1.286000
X 3.242000 0.976000 1.290000
X 2.861000 0.960000 1.575000
158 changes: 113 additions & 45 deletions src/crystdistrib/Quaternion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ Quaternion::Quaternion(const ActionOptions&ao):
derivs(4),
virial(4)
{
for(unsigned i=0; i<4; ++i) derivs[i].resize(2);
for(unsigned i=0; i<4; ++i) derivs[i].resize(3);
std::vector<AtomNumber> atoms;
parseAtomList(-1,atoms,this);
if(atoms.size()!=3) error("Number of specified atoms should be 3");
Expand Down Expand Up @@ -163,77 +163,145 @@ void Quaternion::calculateCV( const unsigned& mode, const std::vector<double>& m
//declarations
Vector vec1_comp = delta( pos[0], pos[1] ); //components between atom 1 and 2
Vector vec2_comp = delta( pos[0], pos[2] ); //components between atom 1 and 3

double x1, x2, x3; //vector 1
double y1, y2, y3; //vector 2
double z1, z2, z3; //vector 3

////////////////////////////////////
////////x-vector calculations///////
double vec1_magnitude = sqrt((vec1_comp[0]*vec1_comp[0])+(vec1_comp[1]*vec1_comp[1])+(vec1_comp[2]*vec1_comp[2])); //normalization factor for atom 2 - atom 1

x1 = vec1_comp[0]/vec1_magnitude;
x2 = vec1_comp[1]/vec1_magnitude;
x3 = vec1_comp[2]/vec1_magnitude;
double magx = vec1_comp.modulo(); Vector x = vec1_comp / magx;
std::vector<Tensor> dx(3); double magx3= magx*magx*magx;
dx[0](0,0) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 ); // dx/dx
dx[0](0,1) = ( vec1_comp[0]*vec1_comp[1]/magx3 ); // dx/dy
dx[0](0,2) = ( vec1_comp[0]*vec1_comp[2]/magx3 ); // dx/dz
dx[0](1,0) = ( vec1_comp[1]*vec1_comp[0]/magx3 ); // dy/dx
dx[0](1,1) = ( -(vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 ); // dy/dy
dx[0](1,2) = ( vec1_comp[1]*vec1_comp[2]/magx3 );
dx[0](2,0) = ( vec1_comp[2]*vec1_comp[0]/magx3 );
dx[0](2,1) = ( vec1_comp[2]*vec1_comp[1]/magx3 );
dx[0](2,2) = ( -(vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );

dx[1](0,0) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[2]*vec1_comp[2])/magx3 );
dx[1](0,1) = ( -vec1_comp[0]*vec1_comp[1]/magx3 );
dx[1](0,2) = ( -vec1_comp[0]*vec1_comp[2]/magx3 );
dx[1](1,0) = ( -vec1_comp[1]*vec1_comp[0]/magx3 );
dx[1](1,1) = ( (vec1_comp[0]*vec1_comp[0]+vec1_comp[2]*vec1_comp[2])/magx3 );
dx[1](1,2) = ( -vec1_comp[1]*vec1_comp[2]/magx3 );
dx[1](2,0) = ( -vec1_comp[2]*vec1_comp[0]/magx3 );
dx[1](2,1) = ( -vec1_comp[2]*vec1_comp[1]/magx3 );
dx[1](2,2) = ( (vec1_comp[1]*vec1_comp[1]+vec1_comp[0]*vec1_comp[0])/magx3 );
dx[2].zero();

/////////////////////////////////////
////////y-vector calculations////////
//project vec2_comp on to vec1_comp
//first do dot product of unormalized x and unormed y, divided by magnitude of x^2
double fac = ((vec1_comp[0]*vec2_comp[0])+(vec1_comp[1]*vec2_comp[1])+(vec1_comp[2]*vec2_comp[2]))/(vec1_magnitude*vec1_magnitude); //fac meaning factor on front
double dp = dotProduct( vec1_comp, vec2_comp ); double magx2=magx*magx;
std::vector<Vector> fac_derivs(3); double magx4=magx2*magx2, fac = dp/magx2; //fac meaning factor on front
fac_derivs[0] = (-vec2_comp - vec1_comp)/magx2 + 2*dp*vec1_comp / magx4;
fac_derivs[1] = (vec2_comp)/(magx2) - 2*dp*vec1_comp / magx4;
fac_derivs[2] = (vec1_comp)/(magx2);
//now multiply fac by unormed x, and subtract it from unormed y, then normalize
y1 = vec2_comp[0] - fac*vec1_comp[0];
y2 = vec2_comp[1] - fac*vec1_comp[1];
y3 = vec2_comp[2] - fac*vec1_comp[2];
Vector y = vec2_comp - fac*vec1_comp; std::vector<Tensor> dy(3);
dy[0](0,0) = -1 - fac_derivs[0][0]*vec1_comp[0] + fac; // dx/dx
dy[0](0,1) = -fac_derivs[0][0]*vec1_comp[1]; // dx/dy
dy[0](0,2) = -fac_derivs[0][0]*vec1_comp[2];
dy[0](1,0) = -fac_derivs[0][1]*vec1_comp[0];
dy[0](1,1) = -1 - fac_derivs[0][1]*vec1_comp[1] + fac;
dy[0](1,2) = -fac_derivs[0][1]*vec1_comp[2];
dy[0](2,0) = -fac_derivs[0][2]*vec1_comp[0];
dy[0](2,1) = -fac_derivs[0][2]*vec1_comp[1];
dy[0](2,2) = -1 - fac_derivs[0][2]*vec1_comp[2] + fac;

dy[1](0,0) = -fac_derivs[1][0]*vec1_comp[0] - fac;
dy[1](0,1) = -fac_derivs[1][0]*vec1_comp[1];
dy[1](0,2) = -fac_derivs[1][0]*vec1_comp[2];
dy[1](1,0) = -fac_derivs[1][1]*vec1_comp[0];
dy[1](1,1) = -fac_derivs[1][1]*vec1_comp[1] - fac;
dy[1](1,2) = -fac_derivs[1][1]*vec1_comp[2];
dy[1](2,0) = -fac_derivs[1][2]*vec1_comp[0];
dy[1](2,1) = -fac_derivs[1][2]*vec1_comp[1];
dy[1](2,2) = -fac_derivs[1][2]*vec1_comp[2] - fac;

dy[2](0,0) = 1 - fac_derivs[2][0]*vec1_comp[0];
dy[2](0,1) = -fac_derivs[2][0]*vec1_comp[1];
dy[2](0,2) = -fac_derivs[2][0]*vec1_comp[2];
dy[2](1,0) = -fac_derivs[2][1]*vec1_comp[0];
dy[2](1,1) = 1 - fac_derivs[2][1]*vec1_comp[1];
dy[2](1,2) = -fac_derivs[2][1]*vec1_comp[2];
dy[2](2,0) = -fac_derivs[2][2]*vec1_comp[0];
dy[2](2,1) = -fac_derivs[2][2]*vec1_comp[1];
dy[2](2,2) = 1 - fac_derivs[2][2]*vec1_comp[2];
//now normalize, and we have our y vector
double y_magnitude = sqrt((y1*y1)+(y2*y2)+(y3*y3));
y1 = y1/y_magnitude;
y2 = y2/y_magnitude;
y3 = y3/y_magnitude;
double magy = y.modulo(); double imagy = 1/magy, magy3 = magy*magy*magy;
Tensor multi_y; for(unsigned i=0; i<3; ++i) multi_y.setRow(i, y[i]*y);
y = y / magy; for(unsigned i=0; i<3; ++i) dy[i] = dy[i] / magy - matmul( dy[i], multi_y ) / magy3;
/////////////////////////////////////
///////z-vector calculations/////////

//simply cross them to make z
Vector cross_xy = crossProduct(Vector(x1,x2,x3),Vector(y1,y2,y3));
Vector z = crossProduct(x,y); std::vector<Tensor> dz(3);
dz[0].setRow( 0, crossProduct( dx[0].getRow(0), y ) + crossProduct( x, dy[0].getRow(0) ) );
dz[0].setRow( 1, crossProduct( dx[0].getRow(1), y ) + crossProduct( x, dy[0].getRow(1) ) );
dz[0].setRow( 2, crossProduct( dx[0].getRow(2), y ) + crossProduct( x, dy[0].getRow(2) ) );

dz[1].setRow( 0, crossProduct( dx[1].getRow(0), y ) + crossProduct( x, dy[1].getRow(0) ) );
dz[1].setRow( 1, crossProduct( dx[1].getRow(1), y ) + crossProduct( x, dy[1].getRow(1) ) );
dz[1].setRow( 2, crossProduct( dx[1].getRow(2), y ) + crossProduct( x, dy[1].getRow(2) ) );

z1 = cross_xy[0];
z2 = cross_xy[1];
z3 = cross_xy[2];
dz[2].setRow( 0, crossProduct( x, dy[2].getRow(0) ) );
dz[2].setRow( 1, crossProduct( x, dy[2].getRow(1) ) );
dz[2].setRow( 2, crossProduct( x, dy[2].getRow(2) ) );

//the above 9 components form an orthonormal basis, centered on the molecule in question
//the rotation matrix is generally the inverse of this matrix, and in this case since it is 1) orthogonal and 2) its determinant is 1
//the inverse is simply the transpose
//we can now convert it to a quaternion, which will then be used to rotate the vector separating 2 molecules - > q1_bar * r_21 * q1, where r21 is the vector separating the atoms
//r21 is a quaternion with a 0 real part ^^^
//essentially this rotates the vector into the frame of molecule 1
double tr = x1 + y2 + z3; //trace of the rotation matrix
double tr = x[0] + y[1] + z[2]; //trace of the rotation matrix

//this converts to quaternion
std::vector<Vector> dS(3);
if (tr > 0) {
double S = sqrt(tr+1.0) * 2; // S=4*qw
vals[0] = 0.25 * S;
vals[1] = (z2 - y3) / S;
vals[2] = (x3 - z1) / S;
vals[3] = (y1 - x2) / S;
} else if ((x1 > y2)&(x1 > z3)) {
float S = sqrt(1.0 + x1 - y2 - z3) * 2; // S=4*qx
vals[0] = (z2 - y3) / S;
vals[1] = 0.25 * S;
vals[2] = (x2 + y1) / S;
vals[3] = (x3 + z1) / S;
} else if (y2 > z3) {
float S = sqrt(1.0 + y2 - x1 - z3) * 2; // S=4*qy
vals[0] = (x3 - z1) / S;
vals[1] = (x2 + y1) / S;
vals[2] = 0.25 * S;
vals[3] = (y3 + z2) / S;
double S = sqrt(tr+1.0) * 2; // S=4*qw
for(unsigned i=0; i<3; ++i) dS[i] = (2/S)*(dx[i].getCol(0) + dy[i].getCol(1) + dz[i].getCol(2));
vals[0] = 0.25 * S; for(unsigned i=0; i<3; ++i) derivs[0][i] =0.25*dS[i];
vals[1] = (z[1] - y[2]) / S;
for(unsigned i=0; i<3; ++i) derivs[1][i] = (1/S)*(dz[i].getCol(1) - dy[i].getCol(2)) - (vals[1]/S)*dS[i];
vals[2] = (x[2] - z[0]) / S;
for(unsigned i=0; i<3; ++i) derivs[2][i] = (1/S)*(dx[i].getCol(2) - dz[i].getCol(0)) - (vals[2]/S)*dS[i];
vals[3] = (y[0] - x[1]) / S;
for(unsigned i=0; i<3; ++i) derivs[3][i] = (1/S)*(dy[i].getCol(0) - dx[i].getCol(1)) - (vals[3]/S)*dS[i];
} else if ((x[0] > y[1])&(x[0] > z[2])) {
float S = sqrt(1.0 + x[0] - y[1] - z[2]) * 2; // S=4*qx
for(unsigned i=0; i<3; ++i) dS[i] = (2/S)*(dx[i].getCol(0) - dy[i].getCol(1) - dz[i].getCol(2));
vals[0] = (z[1] - y[2]) / S;
for(unsigned i=0; i<3; ++i) derivs[0][i] = (1/S)*(dz[i].getCol(1) - dy[i].getCol(2)) - (vals[0]/S)*dS[i];
vals[1] = 0.25 * S; for(unsigned i=0; i<3; ++i) derivs[1][i] =0.25*dS[i];
vals[2] = (x[1] + y[0]) / S;
for(unsigned i=0; i<3; ++i) derivs[2][i] = (1/S)*(dx[i].getCol(1) + dy[i].getCol(0)) - (vals[2]/S)*dS[i];
vals[3] = (x[2] + z[0]) / S;
for(unsigned i=0; i<3; ++i) derivs[3][i] = (1/S)*(dx[i].getCol(2) + dz[i].getCol(0)) - (vals[3]/S)*dS[i];
} else if (y[1] > z[2]) {
float S = sqrt(1.0 + y[1] - x[0] - z[2]) * 2; // S=4*qy
for(unsigned i=0; i<3; ++i) dS[i] = (2/S)*( -dx[i].getCol(0) + dy[i].getCol(1) - dz[i].getCol(2));
vals[0] = (x[2] - z[0]) / S;
for(unsigned i=0; i<3; ++i) derivs[0][i] = (1/S)*(dx[i].getCol(2) - dz[i].getCol(0)) - (vals[0]/S)*dS[i];
vals[1] = (x[1] + y[0]) / S;
for(unsigned i=0; i<3; ++i) derivs[1][i] = (1/S)*(dx[i].getCol(1) + dy[i].getCol(0)) - (vals[1]/S)*dS[i];
vals[2] = 0.25 * S; for(unsigned i=0; i<3; ++i) derivs[2][i] =0.25*dS[i];
vals[3] = (y[2] + z[1]) / S;
for(unsigned i=0; i<3; ++i) derivs[3][i] = (1/S)*(dy[i].getCol(2) + dz[i].getCol(1)) - (vals[3]/S)*dS[i];
} else {
float S = sqrt(1.0 + z3 - x1 - y2) * 2; // S=4*qz
vals[0] = (y1 - x2) / S;
vals[1] = (x3 + z1) / S;
vals[2] = (y3 + z2) / S;
vals[3] = 0.25 * S;
float S = sqrt(1.0 + z[2] - x[0] - y[1]) * 2; // S=4*qz
for(unsigned i=0; i<3; ++i) dS[i] = (2/S)*(-dx[i].getCol(0) - dy[i].getCol(1) + dz[i].getCol(2));
vals[0] = (y[0] - x[1]) / S;
for(unsigned i=0; i<3; ++i) derivs[0][i] = (1/S)*(dy[i].getCol(0) - dx[i].getCol(1)) - (vals[0]/S)*dS[i];
vals[1] = (x[2] + z[0]) / S;
for(unsigned i=0; i<3; ++i) derivs[1][i] = (1/S)*(dx[i].getCol(2) + dz[i].getCol(0)) - (vals[1]/S)*dS[i];
vals[2] = (y[2] + z[1]) / S;
for(unsigned i=0; i<3; ++i) derivs[2][i] = (1/S)*(dy[i].getCol(2) + dz[i].getCol(1)) - (vals[2]/S)*dS[i];
vals[3] = 0.25 * S; for(unsigned i=0; i<3; ++i) derivs[3][i] =0.25*dS[i];
}
setBoxDerivativesNoPbc( pos, derivs, virial );
//debugging vector calcs
//log.printf("%f ", x1);
//log.printf("%f ", x2);
Expand Down

0 comments on commit 0fe398e

Please sign in to comment.