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 "core/ActionWithMatrix.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 : */
59 : //+ENDPLUMEDOC
60 :
61 : namespace PLMD {
62 : namespace crystdistrib {
63 :
64 : class QuaternionProductMatrix : public ActionWithMatrix {
65 : private:
66 : unsigned nderivatives;
67 : public:
68 : static void registerKeywords( Keywords& keys );
69 : explicit QuaternionProductMatrix(const ActionOptions&);
70 : unsigned getNumberOfDerivatives();
71 48 : unsigned getNumberOfColumns() const override {
72 48 : return getConstPntrToComponent(0)->getShape()[1];
73 : }
74 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
75 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
76 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
77 : };
78 :
79 : PLUMED_REGISTER_ACTION(QuaternionProductMatrix,"QUATERNION_PRODUCT_MATRIX")
80 :
81 9 : void QuaternionProductMatrix::registerKeywords( Keywords& keys ) {
82 9 : ActionWithMatrix::registerKeywords(keys);
83 18 : keys.addInputKeyword("compulsory","ARG","vector","the labels of the quaternion vectors that you are outer product of");
84 18 : keys.addOutputComponent("w","default","matrix","the real component of quaternion");
85 18 : keys.addOutputComponent("i","default","matrix","the i component of the quaternion");
86 18 : keys.addOutputComponent("j","default","matrix","the j component of the quaternion");
87 18 : keys.addOutputComponent("k","default","matrix","the k component of the quaternion");
88 9 : }
89 :
90 6 : QuaternionProductMatrix::QuaternionProductMatrix(const ActionOptions&ao):
91 : Action(ao),
92 6 : ActionWithMatrix(ao) {
93 6 : if( getNumberOfArguments()!=8 ) {
94 0 : error("should be eight arguments to this action. Four quaternions for each set of atoms. You can repeat actions");
95 : }
96 6 : unsigned nquat = getPntrToArgument(0)->getNumberOfValues();
97 54 : for(unsigned i=0; i<8; ++i) {
98 : Value* myarg=getPntrToArgument(i);
99 48 : if( i==4 ) {
100 6 : nquat = getPntrToArgument(i)->getNumberOfValues();
101 : }
102 48 : if( myarg->getRank()!=1 ) {
103 0 : error("all arguments to this action should be vectors");
104 : }
105 48 : if( (myarg->getPntrToAction())->getName()!="QUATERNION_VECTOR" ) {
106 0 : error("all arguments to this action should be quaternions");
107 : }
108 48 : std::string mylab=getPntrToArgument(i)->getName();
109 48 : std::size_t dot=mylab.find_first_of(".");
110 72 : if( (i==0 || i==4) && mylab.substr(dot+1)!="w" ) {
111 0 : error("quaternion arguments are in wrong order");
112 : }
113 72 : if( (i==1 || i==5) && mylab.substr(dot+1)!="i" ) {
114 0 : error("quaternion arguments are in wrong order");
115 : }
116 72 : if( (i==2 || i==6) && mylab.substr(dot+1)!="j" ) {
117 0 : error("quaternion arguments are in wrong order");
118 : }
119 72 : if( (i==3 || i==7) && mylab.substr(dot+1)!="k" ) {
120 0 : error("quaternion arguments are in wrong order");
121 : }
122 : }
123 6 : std::vector<unsigned> shape(2);
124 6 : shape[0]=getPntrToArgument(0)->getShape()[0];
125 6 : shape[1]=getPntrToArgument(4)->getShape()[0];
126 6 : addComponent( "w", shape );
127 6 : componentIsNotPeriodic("w");
128 6 : addComponent( "i", shape );
129 6 : componentIsNotPeriodic("i");
130 6 : addComponent( "j", shape );
131 6 : componentIsNotPeriodic("j");
132 6 : addComponent( "k", shape );
133 6 : componentIsNotPeriodic("k");
134 6 : nderivatives = buildArgumentStore(0);
135 6 : }
136 :
137 16 : unsigned QuaternionProductMatrix::getNumberOfDerivatives() {
138 16 : return nderivatives;
139 : }
140 :
141 18 : void QuaternionProductMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
142 18 : unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(4)->getShape()[0];
143 18 : if( indices.size()!=size_v+1 ) {
144 6 : indices.resize( size_v+1 );
145 : }
146 63 : for(unsigned i=0; i<size_v; ++i) {
147 45 : indices[i+1] = start_n + i;
148 : }
149 : myvals.setSplitIndex( size_v + 1 );
150 18 : }
151 :
152 399598 : void QuaternionProductMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
153 399598 : unsigned ostrn, ind2=index2;
154 399598 : if( index2>=getPntrToArgument(0)->getShape()[0] ) {
155 45 : ind2 = index2 - getPntrToArgument(0)->getShape()[0];
156 : }
157 :
158 399598 : std::vector<double> quat1(4), quat2(4);
159 :
160 : // Retrieve the first quaternion
161 1997990 : for(unsigned i=0; i<4; ++i) {
162 1598392 : quat1[i] = getArgumentElement( i, index1, myvals );
163 : }
164 : // Retrieve the second quaternion
165 1997990 : for(unsigned i=0; i<4; ++i) {
166 1598392 : quat2[i] = getArgumentElement( 4+i, ind2, myvals );
167 : }
168 :
169 : //make q1 the conjugate
170 399598 : quat1[1] *= -1;
171 399598 : quat1[2] *= -1;
172 399598 : quat1[3] *= -1;
173 :
174 :
175 : double pref=1;
176 : double pref2=1;
177 : double conj=1;
178 : //real part of q1*q2
179 1997990 : for(unsigned i=0; i<4; ++i) {
180 1598392 : if( i>0 ) {
181 : pref=-1;
182 : pref2=-1;
183 : }
184 1598392 : myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), pref*quat1[i]*quat2[i] );
185 1598392 : if( doNotCalculateDerivatives() ) {
186 839152 : continue ;
187 : }
188 759240 : if (i>0) {
189 : conj=-1;
190 : }
191 759240 : addDerivativeOnVectorArgument( false, 0, i, index1, conj*pref*quat2[i], myvals );
192 759240 : addDerivativeOnVectorArgument( false, 0, 4+i, ind2, pref2*quat1[i], myvals );
193 : }
194 : //i component
195 : pref=1;
196 : conj=1;
197 : pref2=1;
198 1997990 : for (unsigned i=0; i<4; i++) {
199 1598392 : if(i==3) {
200 : pref=-1;
201 : } else {
202 : pref=1;
203 : }
204 1598392 : if(i==2) {
205 : pref2=-1;
206 : } else {
207 : pref2=1;
208 : }
209 1598392 : myvals.addValue( getConstPntrToComponent(1)->getPositionInStream(), pref*quat1[i]*quat2[(5-i)%4]);
210 1598392 : if( doNotCalculateDerivatives() ) {
211 839152 : continue ;
212 : }
213 759240 : if (i>0) {
214 : conj=-1;
215 : }
216 759240 : addDerivativeOnVectorArgument( false, 1, i, index1, conj*pref*quat2[(5-i)%4], myvals );
217 759240 : addDerivativeOnVectorArgument( false, 1, 4+i, ind2, pref2*quat1[(5-i)%4], myvals );
218 : }
219 :
220 : //j component
221 : pref=1;
222 : conj=1;
223 : pref2=1;
224 1997990 : for (unsigned i=0; i<4; i++) {
225 1598392 : if(i==1) {
226 : pref=-1;
227 : } else {
228 : pref=1;
229 : }
230 1598392 : if (i==3) {
231 : pref2=-1;
232 : } else {
233 : pref2=1;
234 : }
235 1598392 : myvals.addValue( getConstPntrToComponent(2)->getPositionInStream(), pref*quat1[i]*quat2[(i+2)%4]);
236 1598392 : if( doNotCalculateDerivatives() ) {
237 839152 : continue ;
238 : }
239 759240 : if (i>0) {
240 : conj=-1;
241 : }
242 759240 : addDerivativeOnVectorArgument( false, 2, i, index1, conj*pref*quat2[(i+2)%4], myvals );
243 759240 : addDerivativeOnVectorArgument( false, 2, 4+i, ind2, pref2*quat1[(i+2)%4], myvals );
244 : }
245 :
246 : //k component
247 : pref=1;
248 : conj=1;
249 : pref2=1;
250 1997990 : for (unsigned i=0; i<4; i++) {
251 1598392 : if(i==2) {
252 : pref=-1;
253 : } else {
254 : pref=1;
255 : }
256 1598392 : if(i==1) {
257 : pref2=-1;
258 : } else {
259 : pref2=1;
260 : }
261 1598392 : myvals.addValue( getConstPntrToComponent(3)->getPositionInStream(), pref*quat1[i]*quat2[(3-i)]);
262 1598392 : if( doNotCalculateDerivatives() ) {
263 839152 : continue ;
264 : }
265 759240 : if (i>0) {
266 : conj=-1;
267 : }
268 759240 : addDerivativeOnVectorArgument( false, 3, i, index1, conj*pref*quat2[3-i], myvals );
269 759240 : addDerivativeOnVectorArgument( false, 3, 4+i, ind2, pref2*quat1[3-i], myvals );
270 :
271 : }
272 :
273 :
274 399598 : }
275 :
276 2028 : void QuaternionProductMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
277 2028 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
278 : return ;
279 : }
280 :
281 3915 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
282 3132 : unsigned nmat = getConstPntrToComponent(j)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
283 : std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
284 : unsigned ntwo_atoms = myvals.getSplitIndex();
285 : // Quaternion for first molecule
286 : unsigned base = 0;
287 15660 : for(unsigned k=0; k<4; ++k) {
288 12528 : matrix_indices[nmat_ind] = base + ival;
289 12528 : base += getPntrToArgument(k)->getShape()[0];
290 12528 : nmat_ind++;
291 : }
292 : // Loop over row of matrix
293 762372 : for(unsigned i=1; i<ntwo_atoms; ++i) {
294 759240 : unsigned ind2 = indices[i];
295 759240 : if( ind2>=getPntrToArgument(0)->getShape()[0] ) {
296 0 : ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
297 : }
298 759240 : base = 4*getPntrToArgument(0)->getShape()[0];
299 : // Quaternion of second molecule
300 3796200 : for(unsigned k=0; k<4; ++k) {
301 3036960 : matrix_indices[nmat_ind] = base + ind2;
302 3036960 : base += getPntrToArgument(4+k)->getShape()[0];
303 3036960 : nmat_ind++;
304 : }
305 : }
306 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
307 : }
308 :
309 : }
310 :
311 : }
312 : }
|
