Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "matrixtools/OuterProduct.h"
23 : #include "core/ActionRegister.h"
24 :
25 : //+PLUMEDOC MCOLVAR QUATERNION_PRODUCT_MATRIX
26 : /*
27 : Calculate the outer product matrix from two vectors of quaternions
28 :
29 : Calculate the outer product matrix from two vectors of quaternions. Quaternions are made of four numbers, one real number, and three imaginary numbers \textit{i}, \textit{j}, and \textit{k}. The imaginary numbers are not commutative:
30 :
31 : \[ij =k\] \[jk=i\] \[ki=j\] \[ji = -k\] \[ik = -j\] \[kj = -i\] \[ii = jj = kk = -1\]
32 :
33 : Thus multiplying two quaternions $\mathbf{q_1} = w_1 + x_1i + y_1j + z_1k $ and $\mathbf{q_2} = w_2 + x_2i + y_2j + z_2k$ yields the following four components:
34 :
35 : \[w_{12} = w_1w_2 - x_1x_2 - y_1y_2 -z_1z_2\]
36 : \[x_{12}i = (w_1x_2 + x_1w_2 + y_1z_2 - z_1y_2 )i\]
37 : \[y_{12}j = (w_1y_2 - x_1z_2 + y_1w_2 + z_1x_2)j\]
38 : \[z_{12}k = (w_1z_2 + x_1y_2 - y_1x_2 + z_1w_2)k\]
39 :
40 : Quaternions can also be multiplied by a real number, and the conjugate $\mathbf{\overline{q}} = w -xi - yj - zk$ is analogous to complex numbers.
41 :
42 :
43 : This action takes two equal sized vectors of quaternions of length $n$, and returns four $n\times n$ outer product matrices, corresponding to components w, x, y, and z in the above example. These matrices are real numbers, and will not behave with the usual nuances of quaternion algebra in any CUSTOMS, or other actions, that will need to be accounted for manually. The vectors of quaternions can be the same vectors, or different. Note, since the QUATERNION action returns unit quaternions, all quaternions returned here will also be unit quaternions.
44 :
45 : ```plumed
46 : #in a system of 4 molecules, each with 3 atoms
47 : #define quaternions 1-4
48 : quats12: QUATERNION ATOMS1=1,2,3 ATOMS2=4,5,6
49 : quats34: QUATERNION ATOMS1=7,8,9 ATOMS2=10,11,12
50 :
51 : #take the outer product of 1,2 and 3,4
52 : qp: QUATERNION_PRODUCT_MATRIX ARG=quats12.*,quats34.* #notice the asterisk
53 : #the components need to be passed in the order w1,x1,y1,z1,w2,x2,y2,z2
54 : #now use in an order parameter
55 : bpw: CUSTOM ARG=qp.* VAR=w,x,y,z FUNC=exp(w/2+x/2+y/2+z/2) PERIODIC=NO
56 : ```
57 :
58 : A common strategy when using this action is to multiply a function of elements of the quaternion product matrix by a contact matrix as shown below:
59 :
60 : ```plumed
61 : # Calculate some quaternions
62 : quats12: QUATERNION ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9 ATOMS4=10,11,12
63 : # Calculate a contact matrix
64 : cmap: CONTACT_MATRIX GROUP=1,4,7,10 SWITCH={RATIONAL R_0=0.2 D_MAX=0.7}
65 : # Calculate the quaternion product matrix
66 : qp: QUATERNION_PRODUCT_MATRIX ARG=quats12.*,quats12.* MASK=cmap
67 : # Now calculate a function of the quaternion elements (normally something more complicated than the following is performed
68 : func: CUSTOM ARG=cmap,qp.* VAR=x,w,i,j,k FUNC=x*(w+i+j+k) PERIODIC=NO
69 : # Now sum the rows of the above matrix
70 : ones: ONES SIZE=4
71 : op: MATRIX_VECTOR_PRODUCT ARG=func,ones
72 : # And output the values of the order parameter
73 : DUMPATOMS ATOMS=1,4,7,10 ARG=op FILE=func.yxz
74 : ```
75 :
76 : Notice how the MASK keyword is used in the CUSTOM in the input to QUATERNION_PRODUCT_MATRIX to prevent PLUMED from calculating the
77 : elements of the QUATERNION_PRODUCT_MATRIX that will be multiplied by the elements of that matrix `cmap` that is calculated in the
78 : CONTACT_MATRIX action when the calculations in the CUSTOM action with label `func` are performed. This is the same procedure that is performed
79 : in [TORSIONS_MATRIX](TORSIONS_MATRIX.md) and [MATRIX_PRODUCT](MATRIX_PRODUCT.md) to avoid computational expense. This procedure is automatically
80 : performed when you use the [ROPS](ROPS.md) shortcut.
81 :
82 : */
83 : //+ENDPLUMEDOC
84 :
85 : namespace PLMD {
86 : namespace crystdistrib {
87 :
88 : class QuaternionProductMatrix {
89 : public:
90 : static void registerKeywords( Keywords& keys );
91 : void setup( const std::vector<std::size_t>& shape, const std::string& func, matrixtools::OuterProductBase<QuaternionProductMatrix>* action );
92 : static void calculate( bool noderiv, const QuaternionProductMatrix& actdata, View<double> vals, MatrixElementOutput& output );
93 : };
94 :
95 : typedef matrixtools::OuterProductBase<QuaternionProductMatrix> qpop;
96 : PLUMED_REGISTER_ACTION(qpop,"QUATERNION_PRODUCT_MATRIX")
97 :
98 9 : void QuaternionProductMatrix::registerKeywords( Keywords& keys ) {
99 18 : keys.addInputKeyword("compulsory","ARG","vector","the labels of the quaternion vectors that you are outer product of");
100 18 : keys.addInputKeyword("optional","MASK","matrix","a matrix that is used to used to determine which elements of the output matrix to compute");
101 18 : keys.addOutputComponent("w","default","matrix","the real component of quaternion");
102 18 : keys.addOutputComponent("i","default","matrix","the i component of the quaternion");
103 18 : keys.addOutputComponent("j","default","matrix","the j component of the quaternion");
104 18 : keys.addOutputComponent("k","default","matrix","the k component of the quaternion");
105 9 : }
106 :
107 6 : void QuaternionProductMatrix::setup( const std::vector<std::size_t>& shape, const std::string& func, matrixtools::OuterProductBase<QuaternionProductMatrix>* action ) {
108 6 : unsigned nargs=action->getNumberOfArguments();
109 6 : if( action->getNumberOfMasks()>0 ) {
110 4 : nargs = nargs - action->getNumberOfMasks();
111 : }
112 6 : if( nargs!=8 ) {
113 0 : action->error("should be eight arguments to this action. Four quaternions for each set of atoms. You can repeat actions");
114 : }
115 54 : for(unsigned i=0; i<8; ++i) {
116 : Value* myarg=action->getPntrToArgument(i);
117 48 : if( myarg->getRank()!=1 ) {
118 0 : action->error("all arguments to this action should be vectors");
119 : }
120 48 : if( (myarg->getPntrToAction())->getName()!="QUATERNION_VECTOR" ) {
121 0 : action->error("all arguments to this action should be quaternions");
122 : }
123 48 : std::string mylab=(action->getPntrToArgument(i))->getName();
124 48 : std::size_t dot=mylab.find_first_of(".");
125 72 : if( (i==0 || i==4) && mylab.substr(dot+1)!="w" ) {
126 0 : action->error("quaternion arguments are in wrong order");
127 : }
128 72 : if( (i==1 || i==5) && mylab.substr(dot+1)!="i" ) {
129 0 : action->error("quaternion arguments are in wrong order");
130 : }
131 72 : if( (i==2 || i==6) && mylab.substr(dot+1)!="j" ) {
132 0 : action->error("quaternion arguments are in wrong order");
133 : }
134 72 : if( (i==3 || i==7) && mylab.substr(dot+1)!="k" ) {
135 0 : action->error("quaternion arguments are in wrong order");
136 : }
137 : }
138 6 : action->addComponent( "w", shape );
139 6 : action->componentIsNotPeriodic("w");
140 6 : action->addComponent( "i", shape );
141 6 : action->componentIsNotPeriodic("i");
142 6 : action->addComponent( "j", shape );
143 6 : action->componentIsNotPeriodic("j");
144 6 : action->addComponent( "k", shape );
145 6 : action->componentIsNotPeriodic("k");
146 6 : }
147 :
148 630393 : void QuaternionProductMatrix::calculate( bool noderiv, const QuaternionProductMatrix& actdata, View<double> vals, MatrixElementOutput& output ) {
149 630393 : std::vector<double> quat1(4), quat2(4);
150 :
151 3151965 : for(unsigned i=0; i<4; ++i) {
152 : // Retrieve the first quaternion
153 2521572 : quat1[i] = vals[i];
154 : // Retrieve the second quaternion
155 2521572 : quat2[i] = vals[4+i];
156 : }
157 :
158 : //make q1 the conjugate
159 630393 : quat1[1] *= -1;
160 630393 : quat1[2] *= -1;
161 630393 : quat1[3] *= -1;
162 :
163 : double pref=1;
164 : double pref2=1;
165 : double conj=1;
166 : //real part of q1*q2
167 3151965 : for(unsigned i=0; i<4; ++i) {
168 2521572 : if( i>0 ) {
169 : pref=-1;
170 : pref2=-1;
171 : }
172 2521572 : output.values[0] += pref*quat1[i]*quat2[i];
173 2521572 : if( noderiv ) {
174 1640388 : continue ;
175 : }
176 881184 : if (i>0) {
177 : conj=-1;
178 : }
179 881184 : output.derivs[0][i] = conj*pref*quat2[i];
180 881184 : output.derivs[0][4+i] = pref2*quat1[i];
181 : }
182 : //i component
183 : pref=1;
184 : conj=1;
185 : pref2=1;
186 3151965 : for(unsigned i=0; i<4; i++) {
187 2521572 : if(i==3) {
188 : pref=-1;
189 : } else {
190 : pref=1;
191 : }
192 2521572 : if(i==2) {
193 : pref2=-1;
194 : } else {
195 : pref2=1;
196 : }
197 2521572 : output.values[1] += pref*quat1[i]*quat2[(5-i)%4];
198 2521572 : if( noderiv ) {
199 1640388 : continue ;
200 : }
201 881184 : if (i>0) {
202 : conj=-1;
203 : }
204 881184 : output.derivs[1][i] = conj*pref*quat2[(5-i)%4];
205 881184 : output.derivs[1][4+i] = pref2*quat1[(5-i)%4];
206 : }
207 :
208 : //j component
209 : pref=1;
210 : conj=1;
211 : pref2=1;
212 3151965 : for (unsigned i=0; i<4; i++) {
213 2521572 : if(i==1) {
214 : pref=-1;
215 : } else {
216 : pref=1;
217 : }
218 2521572 : if (i==3) {
219 : pref2=-1;
220 : } else {
221 : pref2=1;
222 : }
223 2521572 : output.values[2] += pref*quat1[i]*quat2[(i+2)%4];
224 2521572 : if( noderiv ) {
225 1640388 : continue ;
226 : }
227 881184 : if (i>0) {
228 : conj=-1;
229 : }
230 881184 : output.derivs[2][i] = conj*pref*quat2[(i+2)%4];
231 881184 : output.derivs[2][4+i] = pref2*quat1[(i+2)%4];
232 : }
233 :
234 : //k component
235 : pref=1;
236 : conj=1;
237 : pref2=1;
238 3151965 : for (unsigned i=0; i<4; i++) {
239 2521572 : if(i==2) {
240 : pref=-1;
241 : } else {
242 : pref=1;
243 : }
244 2521572 : if(i==1) {
245 : pref2=-1;
246 : } else {
247 : pref2=1;
248 : }
249 2521572 : output.values[3] += pref*quat1[i]*quat2[(3-i)];
250 2521572 : if( noderiv ) {
251 1640388 : continue ;
252 : }
253 881184 : if (i>0) {
254 : conj=-1;
255 : }
256 881184 : output.derivs[3][i] = conj*pref*quat2[3-i];
257 881184 : output.derivs[3][4+i] = pref2*quat1[3-i];
258 : }
259 630393 : }
260 :
261 : }
262 : }
|
