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 : #ifdef __PLUMED_HAS_OPENACC
23 : #define __PLUMED_USE_OPENACC 1
24 : #endif //__PLUMED_HAS_OPENACC
25 : #include "MatrixTimesMatrix.h"
26 : #include "core/ActionRegister.h"
27 :
28 : //+PLUMEDOC MCOLVAR MATRIX_PRODUCT
29 : /*
30 : Calculate the product of two matrices
31 :
32 : This action allows you to [multiply](https://en.wikipedia.org/wiki/Matrix_multiplication) two matrices. The following input shows an
33 : example where two contact matrices are multiplied together.
34 :
35 : ```plumed
36 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.1}
37 : c2: CONTACT_MATRIX GROUP=1-100 SWITCH={RATIONAL R_0=0.2}
38 : m: MATRIX_PRODUCT ARG=c1,c2
39 : PRINT ARG=m FILE=colvar
40 : ```
41 :
42 : The functionality in this action is useful for claculating the relative orientations of large numbers of molecules. For example in
43 : his input the [DISTANCE](DISTANCE.md) command is used to calculate the orientation of a collection of molecules. We then can then use the [VSTACK](VSTACK.md), [TRANSPOSE](TRANSPOSE.md) and the
44 : MATRIX_PRODUCT commands to calculate the dot products between all these vectors
45 :
46 : ```plumed
47 : # Calculate the vectors connecting these three pairs of atoms
48 : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
49 : # Construct a matrix that contains all the components of the vectors calculated
50 : v: VSTACK ARG=d.x,d.y,d.z
51 : # Transpose v
52 : vT: TRANSPOSE ARG=v
53 : # And now calculate the 3x3 matrix of dot products
54 : m: MATRIX_PRODUCT ARG=v,vT
55 : # And output the matrix product to a file
56 : PRINT ARG=m FILE=colvar
57 : ```
58 :
59 : ## The MASK keyword
60 :
61 : We can use MATRIX_PRODUCT to calculate a measure of local order as shown below:
62 :
63 : ```plumed
64 : # Calculate the vectors connecting these three pairs of atoms
65 : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
66 : # Normalize the distance vectors
67 : dm: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
68 : dx: CUSTOM ARG=d.x,dm FUNC=x/y PERIODIC=NO
69 : dy: CUSTOM ARG=d.y,dm FUNC=x/y PERIODIC=NO
70 : dz: CUSTOM ARG=d.z,dm FUNC=x/y PERIODIC=NO
71 : # Construct a matrix that contains all the directors of the vectors calculated
72 : v: VSTACK ARG=dx,dy,dz
73 : # Transpose v
74 : vT: TRANSPOSE ARG=v
75 : # Calculate a contact matrix between pairs of atoms
76 : cmap: CONTACT_MATRIX GROUP=1,3,5 SWITCH={RATIONAL R_0=0.2 D_MAX=0.5}
77 : # Calculate the matrix of dot products
78 : prod: MATRIX_PRODUCT ARG=v,vT MASK=cmap
79 : # Take the product of the two matrices we computed above
80 : p: CUSTOM ARG=cmap,prod FUNC=x*y PERIODIC=NO
81 : # And compute the average of the dot product for the molecules in the
82 : # first coordination sphere around each of the atoms
83 : ones: ONES SIZE=3
84 : numer: MATRIX_VECTOR_PRODUCT ARG=p,ones
85 : denom: MATRIX_VECTOR_PRODUCT ARG=cmap,ones
86 : order: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
87 : # Print the order parameter values
88 : PRINT ARG=order FILE=colvar
89 : ```
90 :
91 : In the above input the [DISTANCE](DISTANCE.md) action is used to calculate the orientations of some molecules (you would normally use more distances than the three
92 : that are used in this input when doing this sort of calculation). We then construct a matrix called `v` that contains the directors of the vectors that define the
93 : orientation of these molecules. The final vector of values `order` that is output here measures the average for the dot product between the orientations of the molecules
94 : and the molecules in their first coordination sphere.
95 :
96 : Ideass similar to this are used in the calculation of [LOCAL_Q3](LOCAL_Q3.md), [LOCAL_Q4](LOCAL_Q4.md) and [LOCAL_Q6](LOCAL_Q6.md). Notice, that it is very important to
97 : use `MASK` keyword in the MATRIX_PRODUCT action when you are doing this type calculation. Using this keyword ensures that PLUMED does not calculate the $(i,j)$ matrix of
98 : the matrix `prod` if the corresponding element in `cmap` is zero. Using this keyword thus ensures that we can avoid a large number of unecessary calculations and improves
99 : the performance of the calculation considerably.
100 :
101 : ## Calculating angles
102 :
103 : MATRIX_PRODUCT can also be used to calculate a matrix of angles between bonds as shown below:
104 :
105 : ```plumed
106 : # Calculate the directors for a set of vectors
107 : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ATOMS5=1,6
108 : dm: DISTANCE ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ATOMS5=1,6
109 : dx: CUSTOM ARG=d.x,dm FUNC=x/y PERIODIC=NO
110 : dy: CUSTOM ARG=d.y,dm FUNC=x/y PERIODIC=NO
111 : dz: CUSTOM ARG=d.z,dm FUNC=x/y PERIODIC=NO
112 : # Construct a matrix that contains all the directors of the vectors calculated
113 : v: VSTACK ARG=dx,dy,dz
114 : # Transpose v
115 : vT: TRANSPOSE ARG=v
116 : # Calculate the matrix of dot products between the input directors
117 : dpmat: MATRIX_PRODUCT ELEMENTS_ON_DIAGONAL_ARE_ZERO ARG=v,vT
118 : # And calculate the angles
119 : angles: CUSTOM ARG=dpmat FUNC=acos(x) PERIODIC=NO
120 : # Print the matrix of angles
121 : PRINT ARG=angles FILE=colvar
122 : ```
123 :
124 : Notice that we have to use the `ELEMENTS_ON_DIAGONAL_ARE_ZERO` flag here to avoid numerical issues in the calculation.
125 :
126 : */
127 : //+ENDPLUMEDOC
128 :
129 : //+PLUMEDOC ANALYSIS DISSIMILARITIES
130 : /*
131 : Calculate the matrix of dissimilarities between a trajectory of atomic configurations.
132 :
133 : This action allows you to calculate a dissimilarity matrix, which is a square matrix tht describes the
134 : pairwise distinction between a set of objects. This action is routinely used in dimensionality reduction
135 : calculations. The example shown below shows how we might get a matrix of dissimilarities upon which we can run a
136 : dimensionality reduction calculation.
137 :
138 : ```plumed
139 : d1: DISTANCE ATOMS=1,2
140 : d2: DISTANCE ATOMS=3,4
141 : # Collect the values of the two distances and store them for later analysis
142 : ff: COLLECT_FRAMES ARG=d1,d2 STRIDE=1
143 : # Transpose the matrix that is collected above
144 : ff_dataT: TRANSPOSE ARG=ff_data
145 : # Now compute all the dissimilarities between the collected frames
146 : ss1: DISSIMILARITIES ARG=ff_data,ff_dataT
147 : DUMPVECTOR ARG=ss1 FILE=mymatrix.dat
148 : ```
149 :
150 : Some dimensionality reduction algorithms take the squares of the dissimilarities rather than the dissimilarities.
151 : To calculate these quantities using PLUMED you use the SQUARED keyword as shown below:
152 :
153 : ```plumed
154 : d1: DISTANCE ATOMS=1,2
155 : d2: DISTANCE ATOMS=3,4
156 : # Collect the values of the two distances and store them for later analysis
157 : ff: COLLECT_FRAMES ARG=d1,d2 STRIDE=1
158 : # Transpose the matrix that is collected above
159 : ff_dataT: TRANSPOSE ARG=ff_data
160 : # Now compute all the dissimilarities between the collected frames
161 : ss1: DISSIMILARITIES ARG=ff_data,ff_dataT SQUARED
162 : DUMPVECTOR ARG=ss1 FILE=mymatrix.dat
163 : ```
164 :
165 : ## The MASK keyword
166 :
167 : The MASK keyword for DISSIMILARITIES works in the same way that it does for [MATRIX_PRODUCT](MATRIX_PRODUCT.md). This keyword thus
168 : expects a matrix in input and will only evaluate elements the elements of the dissimilarity matrix whose corresponding elements in
169 : the matrix that was input to MASK are non-zero. The following input shows an example of how this might be used in an example input
170 : that has not been used in any production calculation.
171 :
172 : ```plumed
173 : # Calculate and store three vectors connecting three pairs of atoms
174 : d1: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
175 : ff1: VSTACK ARG=d1.x,d1.y,d1.z
176 : # Calculate the magnitude of the vector and transform this magnitude using a
177 : # switching function
178 : d1m: CUSTOM ARG=d1.x,d1.y,d1.z FUNC=sqrt(x*x+y*y+z*z) PERIODIC=NO
179 : sw1: LESS_THAN ARG=d1m SWITCH={RATIONAL R_0=0.2 D_MAX=0.5}
180 : # Now perform the same operations for a different set of three atoms
181 : d2: DISTANCE COMPONENTS ATOMS1=7,8 ATOMS2=9,10 ATOMS3=11,12
182 : ff2: VSTACK ARG=d2.x,d2.y,d2.z
183 : ff2T: TRANSPOSE ARG=ff2
184 : d2m: CUSTOM ARG=d2.x,d2.y,d2.z FUNC=sqrt(x*x+y*y+z*z) PERIODIC=NO
185 : sw2: LESS_THAN ARG=d2m SWITCH={RATIONAL R_0=0.2 D_MAX=0.5}
186 : # Use OUTER_PRODUCT to work out which pairs of atoms are less than a certain cutoff
187 : swmat: OUTER_PRODUCT ARG=sw1,sw2
188 : # Now calculate the dissimilarities between the two sets of vectors for those pairs of atoms that are within a certain cutoff
189 : d: DISSIMILARITIES ARG=ff1,ff2T MASK=swmat
190 : mat: CUSTOM ARG=swmat,d FUNC=x*y PERIODIC=NO
191 : # And compute the average dissimilarity
192 : numer: SUM ARG=mat PERIODIC=NO
193 : denom: SUM ARG=swmat PERIODIC=NO
194 : v: CUSTOM ARG=numer,denom FUNC=x/y PERIODIC=NO
195 : PRINT ARG=v FILE=colvar
196 : ```
197 :
198 : We have not provided a more relevant example because we haven't really used this functionality for anything. If you have a better example
199 : for a way of using this functionality please get in touch so we can update this part of the manual with your example.
200 :
201 : */
202 : //+ENDPLUMEDOC
203 :
204 : namespace PLMD {
205 : namespace matrixtools {
206 :
207 : struct MatrixProduct {
208 : static void registerKeywords( Keywords& keys );
209 : void setup( MatrixTimesMatrix<MatrixProduct>* action, const Value* myval ) {}
210 : static void calculate( bool noderiv,
211 : const MatrixProduct& actdata,
212 : const InputVectors& vectors,
213 : MatrixElementOutput& output );
214 : #ifdef __PLUMED_USE_OPENACC
215 : void toACCDevice() const {
216 : #pragma acc enter data copyin(this[0:1])
217 : }
218 : void removeFromACCDevice() const {
219 : #pragma acc exit data delete(this[0:1])
220 : }
221 : #endif //__PLUMED_USE_OPENACC
222 : };
223 :
224 : typedef MatrixTimesMatrix<MatrixProduct> mtimes;
225 : PLUMED_REGISTER_ACTION(mtimes,"MATRIX_PRODUCT")
226 :
227 77 : void MatrixProduct::registerKeywords( Keywords& keys ) {
228 154 : keys.setValueDescription("matrix","the product of the two input matrices");
229 77 : }
230 :
231 5196607 : void MatrixProduct::calculate( bool noderiv,
232 : const MatrixProduct& actdata,
233 : const InputVectors& vectors,
234 : MatrixElementOutput& output ) {
235 : std::size_t pp = vectors.arg1.size();
236 137736753 : for(unsigned i=0; i<vectors.nelem; ++i) {
237 132540146 : output.values[0] += vectors.arg1[i]*vectors.arg2[i];
238 132540146 : output.derivs[0][i] = vectors.arg2[i];
239 132540146 : output.derivs[0][pp+i] = vectors.arg1[i];
240 : }
241 5196607 : }
242 :
243 24 : struct Dissimilarities {
244 : bool squared{false};
245 : bool periodic{false};
246 : double min{0};
247 : double max{0};
248 : double max_minus_min{0};
249 : double inv_max_minus_min{0};
250 : static void registerKeywords( Keywords& keys );
251 : void setup( MatrixTimesMatrix<Dissimilarities>* action,
252 : const Value* myval );
253 : static void calculate( bool noderiv,
254 : const Dissimilarities& actdata,
255 : const InputVectors& vectors,
256 : MatrixElementOutput& output );
257 : #ifdef __PLUMED_USE_OPENACC
258 : void toACCDevice() const {
259 : #pragma acc enter data copyin(this[0:1], squared, periodic, min, max, \
260 : max_minus_min, inv_max_minus_min)
261 :
262 : }
263 : void removeFromACCDevice() const {
264 : #pragma acc exit data delete(inv_max_minus_min, max_minus_min, max, \
265 : min, periodic, squared, this[0:1])
266 : }
267 : #endif //__PLUMED_USE_OPENACC
268 : };
269 :
270 : typedef MatrixTimesMatrix<Dissimilarities> dissims;
271 : PLUMED_REGISTER_ACTION(dissims,"DISSIMILARITIES")
272 :
273 21 : void Dissimilarities::registerKeywords( Keywords& keys ) {
274 21 : keys.addFlag("SQUARED",false,"calculate the squares of the dissimilarities (this option cannot be used with MATRIX_PRODUCT)");
275 42 : keys.setValueDescription("matrix","the dissimilarity matrix");
276 21 : }
277 :
278 12 : void Dissimilarities::setup( MatrixTimesMatrix<Dissimilarities>* action, const Value* myval ) {
279 12 : action->parseFlag("SQUARED",squared);
280 12 : if( squared ) {
281 7 : action->log.printf(" calculating the squares of the dissimilarities \n");
282 : }
283 12 : periodic = myval->isPeriodic();
284 12 : if( periodic ) {
285 : std::string strmin, strmax;
286 3 : myval->getDomain( strmin, strmax );
287 3 : Tools::convert( strmin, min );
288 3 : Tools::convert( strmax, max );
289 3 : max_minus_min = max - min;
290 3 : plumed_assert( max_minus_min>0 );
291 3 : inv_max_minus_min=1.0/max_minus_min;
292 : }
293 12 : }
294 :
295 545753 : void Dissimilarities::calculate( bool noderiv, const Dissimilarities& actdata, const InputVectors& vectors, MatrixElementOutput& output ) {
296 : std::size_t pp = vectors.arg1.size();
297 2175346 : for(unsigned i=0; i<vectors.nelem; ++i) {
298 1629593 : double tmp1 = vectors.arg1[i] - vectors.arg2[i];
299 1629593 : if( actdata.periodic ) {
300 1530000 : tmp1 = tmp1*actdata.inv_max_minus_min;
301 1530000 : tmp1 = Tools::pbc(tmp1);
302 1530000 : tmp1 = tmp1*actdata.max_minus_min;
303 : }
304 1629593 : output.values[0] += tmp1*tmp1;
305 1629593 : output.derivs[0][i] = 2*tmp1;
306 1629593 : output.derivs[0][pp+i] = -2*tmp1;
307 : }
308 545753 : if( actdata.squared ) {
309 : return ;
310 : }
311 388 : output.values[0] = sqrt( output.values[0] );
312 :
313 1151 : for(unsigned i=0; i<vectors.nelem; ++i) {
314 763 : output.derivs[0][i] = output.derivs[0][i] / (2*output.values[0]);
315 763 : output.derivs[0][pp+i] = output.derivs[0][vectors.nelem+i] / (2*output.values[0]);
316 : }
317 : }
318 :
319 : }
320 : }
|
