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 : \par Examples
30 :
31 : */
32 : //+ENDPLUMEDOC
33 :
34 : namespace PLMD {
35 : namespace crystdistrib {
36 :
37 : class QuaternionProductMatrix : public ActionWithMatrix {
38 : private:
39 : unsigned nderivatives;
40 : public:
41 : static void registerKeywords( Keywords& keys );
42 : explicit QuaternionProductMatrix(const ActionOptions&);
43 : unsigned getNumberOfDerivatives();
44 48 : unsigned getNumberOfColumns() const override { return getConstPntrToComponent(0)->getShape()[1]; }
45 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
46 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
47 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
48 : };
49 :
50 : PLUMED_REGISTER_ACTION(QuaternionProductMatrix,"QUATERNION_PRODUCT_MATRIX")
51 :
52 9 : void QuaternionProductMatrix::registerKeywords( Keywords& keys ) {
53 9 : ActionWithMatrix::registerKeywords(keys); keys.use("ARG");
54 18 : keys.addOutputComponent("w","default","the real component of quaternion");
55 18 : keys.addOutputComponent("i","default","the i component of the quaternion");
56 18 : keys.addOutputComponent("j","default","the j component of the quaternion");
57 18 : keys.addOutputComponent("k","default","the k component of the quaternion");
58 9 : }
59 :
60 6 : QuaternionProductMatrix::QuaternionProductMatrix(const ActionOptions&ao):
61 : Action(ao),
62 6 : ActionWithMatrix(ao)
63 : {
64 6 : if( getNumberOfArguments()!=8 ) error("should be eight arguments to this action. Four quaternions for each set of atoms. You can repeat actions");
65 6 : unsigned nquat = getPntrToArgument(0)->getNumberOfValues();
66 54 : for(unsigned i=0; i<8; ++i) {
67 48 : Value* myarg=getPntrToArgument(i); if( i==4 ) nquat = getPntrToArgument(i)->getNumberOfValues();
68 48 : if( myarg->getRank()!=1 ) error("all arguments to this action should be vectors");
69 48 : if( (myarg->getPntrToAction())->getName()!="QUATERNION_VECTOR" ) error("all arguments to this action should be quaternions");
70 48 : std::string mylab=getPntrToArgument(i)->getName(); std::size_t dot=mylab.find_first_of(".");
71 72 : if( (i==0 || i==4) && mylab.substr(dot+1)!="w" ) error("quaternion arguments are in wrong order");
72 72 : if( (i==1 || i==5) && mylab.substr(dot+1)!="i" ) error("quaternion arguments are in wrong order");
73 72 : if( (i==2 || i==6) && mylab.substr(dot+1)!="j" ) error("quaternion arguments are in wrong order");
74 72 : if( (i==3 || i==7) && mylab.substr(dot+1)!="k" ) error("quaternion arguments are in wrong order");
75 : }
76 6 : std::vector<unsigned> shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(4)->getShape()[0];
77 12 : addComponent( "w", shape ); componentIsNotPeriodic("w");
78 12 : addComponent( "i", shape ); componentIsNotPeriodic("i");
79 12 : addComponent( "j", shape ); componentIsNotPeriodic("j");
80 12 : addComponent( "k", shape ); componentIsNotPeriodic("k");
81 6 : nderivatives = buildArgumentStore(0);
82 6 : }
83 :
84 16 : unsigned QuaternionProductMatrix::getNumberOfDerivatives() {
85 16 : return nderivatives;
86 : }
87 :
88 18 : void QuaternionProductMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
89 18 : unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(4)->getShape()[0];
90 18 : if( indices.size()!=size_v+1 ) indices.resize( size_v+1 );
91 63 : for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + i;
92 : myvals.setSplitIndex( size_v + 1 );
93 18 : }
94 :
95 399598 : void QuaternionProductMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
96 399598 : unsigned ostrn, ind2=index2;
97 399598 : if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0];
98 :
99 399598 : std::vector<double> quat1(4), quat2(4);
100 :
101 : // Retrieve the first quaternion
102 1997990 : for(unsigned i=0; i<4; ++i) quat1[i] = getArgumentElement( i, index1, myvals );
103 : // Retrieve the second quaternion
104 1997990 : for(unsigned i=0; i<4; ++i) quat2[i] = getArgumentElement( 4+i, ind2, myvals );
105 :
106 : //make q1 the conjugate
107 399598 : quat1[1] *= -1;
108 399598 : quat1[2] *= -1;
109 399598 : quat1[3] *= -1;
110 :
111 :
112 : double pref=1;
113 : double pref2=1;
114 : double conj=1;
115 : //real part of q1*q2
116 1997990 : for(unsigned i=0; i<4; ++i) {
117 1598392 : if( i>0 ) {pref=-1; pref2=-1;}
118 1598392 : myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), pref*quat1[i]*quat2[i] );
119 1598392 : if( doNotCalculateDerivatives() ) continue ;
120 759240 : if (i>0) conj=-1;
121 759240 : addDerivativeOnVectorArgument( false, 0, i, index1, conj*pref*quat2[i], myvals );
122 759240 : addDerivativeOnVectorArgument( false, 0, 4+i, ind2, pref2*quat1[i], myvals );
123 : }
124 : //i component
125 : pref=1;
126 : conj=1;
127 : pref2=1;
128 1997990 : for (unsigned i=0; i<4; i++) {
129 1598392 : if(i==3) pref=-1;
130 : else pref=1;
131 1598392 : if(i==2) pref2=-1;
132 : else pref2=1;
133 1598392 : myvals.addValue( getConstPntrToComponent(1)->getPositionInStream(), pref*quat1[i]*quat2[(5-i)%4]);
134 1598392 : if( doNotCalculateDerivatives() ) continue ;
135 759240 : if (i>0) conj=-1;
136 759240 : addDerivativeOnVectorArgument( false, 1, i, index1, conj*pref*quat2[(5-i)%4], myvals );
137 759240 : addDerivativeOnVectorArgument( false, 1, 4+i, ind2, pref2*quat1[(5-i)%4], myvals );
138 : }
139 :
140 : //j component
141 : pref=1;
142 : conj=1;
143 : pref2=1;
144 1997990 : for (unsigned i=0; i<4; i++) {
145 1598392 : if(i==1) pref=-1;
146 : else pref=1;
147 1598392 : if (i==3) pref2=-1;
148 : else pref2=1;
149 1598392 : myvals.addValue( getConstPntrToComponent(2)->getPositionInStream(), pref*quat1[i]*quat2[(i+2)%4]);
150 1598392 : if( doNotCalculateDerivatives() ) continue ;
151 759240 : if (i>0) conj=-1;
152 759240 : addDerivativeOnVectorArgument( false, 2, i, index1, conj*pref*quat2[(i+2)%4], myvals );
153 759240 : addDerivativeOnVectorArgument( false, 2, 4+i, ind2, pref2*quat1[(i+2)%4], myvals );
154 : }
155 :
156 : //k component
157 : pref=1;
158 : conj=1;
159 : pref2=1;
160 1997990 : for (unsigned i=0; i<4; i++) {
161 1598392 : if(i==2) pref=-1;
162 : else pref=1;
163 1598392 : if(i==1) pref2=-1;
164 : else pref2=1;
165 1598392 : myvals.addValue( getConstPntrToComponent(3)->getPositionInStream(), pref*quat1[i]*quat2[(3-i)]);
166 1598392 : if( doNotCalculateDerivatives() ) continue ;
167 759240 : if (i>0) conj=-1;
168 759240 : addDerivativeOnVectorArgument( false, 3, i, index1, conj*pref*quat2[3-i], myvals );
169 759240 : addDerivativeOnVectorArgument( false, 3, 4+i, ind2, pref2*quat1[3-i], myvals );
170 :
171 : }
172 :
173 :
174 399598 : }
175 :
176 2028 : void QuaternionProductMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
177 2028 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) return ;
178 :
179 3915 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
180 3132 : unsigned nmat = getConstPntrToComponent(j)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
181 : std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntwo_atoms = myvals.getSplitIndex();
182 : // Quaternion for first molecule
183 15660 : unsigned base = 0; for(unsigned k=0; k<4; ++k) { matrix_indices[nmat_ind] = base + ival; base += getPntrToArgument(k)->getShape()[0]; nmat_ind++; }
184 : // Loop over row of matrix
185 762372 : for(unsigned i=1; i<ntwo_atoms; ++i) {
186 759240 : unsigned ind2 = indices[i]; if( ind2>=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
187 759240 : base = 4*getPntrToArgument(0)->getShape()[0];
188 : // Quaternion of second molecule
189 3796200 : for(unsigned k=0; k<4; ++k) { matrix_indices[nmat_ind] = base + ind2; base += getPntrToArgument(4+k)->getShape()[0]; nmat_ind++; }
190 : }
191 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
192 : }
193 :
194 : }
195 :
196 : }
197 : }
|