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 12 : void QuaternionProductMatrix::registerKeywords( Keywords& keys ) {
99 24 : keys.addInputKeyword("compulsory","ARG","vector","the labels of the quaternion vectors that you are outer product of");
100 24 : keys.addInputKeyword("optional","MASK","matrix","a matrix that is used to used to determine which elements of the output matrix to compute");
101 24 : keys.addOutputComponent("w","default","matrix","the real component of quaternion");
102 24 : keys.addOutputComponent("i","default","matrix","the i component of the quaternion");
103 24 : keys.addOutputComponent("j","default","matrix","the j component of the quaternion");
104 24 : keys.addOutputComponent("k","default","matrix","the k component of the quaternion");
105 12 : }
106 :
107 7 : void QuaternionProductMatrix::setup( const std::vector<std::size_t>& shape, const std::string& func, matrixtools::OuterProductBase<QuaternionProductMatrix>* action ) {
108 7 : unsigned nargs=action->getNumberOfArguments();
109 7 : if( action->getNumberOfMasks()>0 ) {
110 5 : nargs = nargs - action->getNumberOfMasks();
111 : }
112 7 : 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 63 : for(unsigned i=0; i<8; ++i) {
116 : Value* myarg=action->getPntrToArgument(i);
117 56 : if( myarg->getRank()!=1 ) {
118 0 : action->error("all arguments to this action should be vectors");
119 : }
120 56 : if( (myarg->getPntrToAction())->getName()!="QUATERNION_VECTOR" ) {
121 0 : action->error("all arguments to this action should be quaternions");
122 : }
123 56 : std::string mylab=(action->getPntrToArgument(i))->getName();
124 56 : std::size_t dot=mylab.find_first_of(".");
125 70 : if( (i==0 || i==4) && mylab.substr(dot+1)!="w" ) {
126 0 : action->error("quaternion arguments are in wrong order");
127 : }
128 70 : if( (i==1 || i==5) && mylab.substr(dot+1)!="i" ) {
129 0 : action->error("quaternion arguments are in wrong order");
130 : }
131 70 : if( (i==2 || i==6) && mylab.substr(dot+1)!="j" ) {
132 0 : action->error("quaternion arguments are in wrong order");
133 : }
134 70 : if( (i==3 || i==7) && mylab.substr(dot+1)!="k" ) {
135 0 : action->error("quaternion arguments are in wrong order");
136 : }
137 : }
138 7 : action->addComponent( "w", shape );
139 7 : action->componentIsNotPeriodic("w");
140 7 : action->addComponent( "i", shape );
141 7 : action->componentIsNotPeriodic("i");
142 7 : action->addComponent( "j", shape );
143 7 : action->componentIsNotPeriodic("j");
144 7 : action->addComponent( "k", shape );
145 7 : action->componentIsNotPeriodic("k");
146 7 : }
147 :
148 580943 : void QuaternionProductMatrix::calculate( bool noderiv, const QuaternionProductMatrix& actdata, View<double> vals, MatrixElementOutput& output ) {
149 580943 : std::vector<double> quat1(4), quat2(4);
150 :
151 2904715 : for(unsigned i=0; i<4; ++i) {
152 : // Retrieve the first quaternion
153 2323772 : quat1[i] = vals[i];
154 : // Retrieve the second quaternion
155 2323772 : quat2[i] = vals[4+i];
156 : }
157 :
158 : //make q1 the conjugate
159 580943 : quat1[1] *= -1;
160 580943 : quat1[2] *= -1;
161 580943 : quat1[3] *= -1;
162 :
163 : double pref=1;
164 : double pref2=1;
165 : double conj=1;
166 : //real part of q1*q2
167 2904715 : for(unsigned i=0; i<4; ++i) {
168 2323772 : if( i>0 ) {
169 : pref=-1;
170 : pref2=-1;
171 : }
172 2323772 : output.values[0] += pref*quat1[i]*quat2[i];
173 2323772 : if( noderiv ) {
174 1541380 : continue ;
175 : }
176 782392 : if (i>0) {
177 : conj=-1;
178 : }
179 782392 : output.derivs[0][i] = conj*pref*quat2[i];
180 782392 : output.derivs[0][4+i] = pref2*quat1[i];
181 : }
182 : //i component
183 : pref=1;
184 : conj=1;
185 : pref2=1;
186 2904715 : for(unsigned i=0; i<4; i++) {
187 2323772 : if(i==3) {
188 : pref=-1;
189 : } else {
190 : pref=1;
191 : }
192 1742829 : if(i==2) {
193 : pref2=-1;
194 : } else {
195 : pref2=1;
196 : }
197 2323772 : output.values[1] += pref*quat1[i]*quat2[(5-i)%4];
198 2323772 : if( noderiv ) {
199 1541380 : continue ;
200 : }
201 782392 : if (i>0) {
202 : conj=-1;
203 : }
204 782392 : output.derivs[1][i] = conj*pref*quat2[(5-i)%4];
205 782392 : 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 2904715 : for (unsigned i=0; i<4; i++) {
213 2323772 : if(i==1) {
214 : pref=-1;
215 : } else {
216 : pref=1;
217 : }
218 1742829 : if (i==3) {
219 : pref2=-1;
220 : } else {
221 : pref2=1;
222 : }
223 2323772 : output.values[2] += pref*quat1[i]*quat2[(i+2)%4];
224 2323772 : if( noderiv ) {
225 1541380 : continue ;
226 : }
227 782392 : if (i>0) {
228 : conj=-1;
229 : }
230 782392 : output.derivs[2][i] = conj*pref*quat2[(i+2)%4];
231 782392 : 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 2904715 : for (unsigned i=0; i<4; i++) {
239 2323772 : if(i==2) {
240 : pref=-1;
241 : } else {
242 : pref=1;
243 : }
244 1742829 : if(i==1) {
245 : pref2=-1;
246 : } else {
247 : pref2=1;
248 : }
249 2323772 : output.values[3] += pref*quat1[i]*quat2[(3-i)];
250 2323772 : if( noderiv ) {
251 1541380 : continue ;
252 : }
253 782392 : if (i>0) {
254 : conj=-1;
255 : }
256 782392 : output.derivs[3][i] = conj*pref*quat2[3-i];
257 782392 : output.derivs[3][4+i] = pref2*quat1[3-i];
258 : }
259 580943 : }
260 :
261 : }
262 : }
|
