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 "OuterProduct.h"
23 : #include "core/ActionShortcut.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/LeptonCall.h"
26 :
27 : //+PLUMEDOC COLVAR OUTER_PRODUCT
28 : /*
29 : Calculate the outer product matrix of two vectors
30 :
31 : This action can be used to calculate the [outer product](https://en.wikipedia.org/wiki/Outer_product) of two
32 : vectors. As a (useless) example of what can be done with this action consider the following simple input:
33 :
34 : ```plumed
35 : d1: DISTANCE ATOMS1=1,2 ATOMS2=3,4
36 : d2: DISTANCE ATOMS1=5,6 ATOMS2=7,8 ATOMS3=9,10
37 : pp: OUTER_PRODUCT ARG=d1,d2
38 : PRINT ARG=pp FILE=colvar
39 : ```
40 :
41 : This input outputs a $2 \times 3$ matrix. If we call the 2 dimensional vector output by the first DISTANCE action
42 : $d$ and the 3 dimensional vector output by the second DISTANCE action $h$ then the $(i,j)$ element of the matrix
43 : output by the action with the label `pp` is given by:
44 :
45 : $$
46 : p_{ij} = d_i h_j
47 : $$
48 :
49 : These outer product matrices are useful if you are trying to calculate an adjacency matrix that says atoms are
50 : connected if they are within a certain distance of each other and if they satisfy a certain criterion. For example,
51 : consider the following input:
52 :
53 : ```plumed
54 : # Determine if atoms are within 5.3 nm of each other
55 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3}
56 : # Calculate the coordination numbers
57 : ones: ONES SIZE=100
58 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
59 : # Now use MORE_THAN to work out which atoms have a coordination number that is bigger than six
60 : cf: MORE_THAN ARG=cc SWITCH={RATIONAL D_0=5.5 R_0=0.5}
61 : # Now recalculate the contact matrix above as first step towards calculating adjacency matrix that measures if
62 : # atoms are close to each other and both have a coordination number that is bigger than six
63 : c2: CONTACT_MATRIX GROUP=1-100 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3}
64 : # Now make a matrix in which element i,j is one if atom i and atom j both have a coordination number that is greater than 6
65 : cfm: OUTER_PRODUCT ARG=cf,cf MASK=c2
66 : # And multiply this by our contact matrix to determine the desired adjacency matrix
67 : m: CUSTOM ARG=c2,cfm FUNC=x*y PERIODIC=NO
68 : f: SUM ARG=m PERIODIC=NO
69 : PRINT ARG=f FILE=colvar
70 : ```
71 :
72 : This input calculates a adjacency matrix which has element $(i,j)$ equal to one if atoms $i$ and $j$ have coordination numbers
73 : that are greater than 6 and if they are within 5.3 nm of each other. Notice how the `MASK` keyword is used in the input to the
74 : OUTER_PRODUCT action here to ensure that do not calculate elements of the `cfm` matrix that will be mulitplied by elements of the
75 : matrix `c2` that are zero. The final quantity output is equal to two times the number of pairs of atoms that are within 5.3 nm of each
76 : and which both have coordination numbers of six.
77 :
78 : Notice that you can specify the function of the two input vectors that is to be calculated by using the `FUNC` keyword which accepts
79 : mathematical expressions of $x$ and $y$. In other words, the elements of the outer product are calculated using the lepton library
80 : that is used in the [CUSTOM](CUSTOM.md) action. In addition, you can set `FUNC=min` or `FUNC=max` to set the elements of the outer product equal to
81 : the minimum of the two input variables or the maximum respectively.
82 :
83 : ## Calculating angles in the first coordination sphere
84 :
85 : We can use OUTER_PRODUCT to calculate a matrix of angles between bonds as shown below:
86 :
87 : ```plumed
88 : # Calculate the directors for a set of vectors
89 : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ATOMS5=1,6
90 : dm: DISTANCE ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ATOMS5=1,6
91 : dx: CUSTOM ARG=d.x,dm FUNC=x/y PERIODIC=NO
92 : dy: CUSTOM ARG=d.y,dm FUNC=x/y PERIODIC=NO
93 : dz: CUSTOM ARG=d.z,dm FUNC=x/y PERIODIC=NO
94 : # Construct a matrix that contains all the directors of the vectors calculated
95 : v: VSTACK ARG=dx,dy,dz
96 : # Transpose v
97 : vT: TRANSPOSE ARG=v
98 : # Transform the distances by a switching functions to determine pairs of atoms that are bonded
99 : sw: LESS_THAN ARG=dm SWITCH={RATIONAL R_0=0.2}
100 : # Calculate the matrix of dot products between the input directors
101 : dpmat: MATRIX_PRODUCT ELEMENTS_ON_DIAGONAL_ARE_ZERO ARG=v,vT
102 : # Use the transformed distances to determine which triples of atoms are bonded
103 : swmat: OUTER_PRODUCT ELEMENTS_ON_DIAGONAL_ARE_ZERO ARG=sw,sw
104 : # And calculate the angles
105 : angles: CUSTOM ARG=swmat,dpmat FUNC=x*acos(y) PERIODIC=NO
106 : # Print the matrix of angles
107 : PRINT ARG=angles FILE=colvar
108 : ```
109 :
110 : Notice that we have to use the `ELEMENTS_ON_DIAGONAL_ARE_ZERO` flag here to avoid numerical issues in the calculation.
111 :
112 : */
113 : //+ENDPLUMEDOC
114 :
115 : namespace PLMD {
116 : namespace matrixtools {
117 :
118 : class OuterProduct : public ActionShortcut {
119 : public:
120 : static void getKeywords( Keywords& keys );
121 : static void registerKeywords( Keywords& keys );
122 : explicit OuterProduct(const ActionOptions&);
123 : };
124 :
125 : PLUMED_REGISTER_ACTION(OuterProduct,"OUTER_PRODUCT")
126 :
127 311 : void OuterProduct::getKeywords( Keywords& keys ) {
128 311 : keys.setDisplayName("OUTER_PRODUCT");
129 622 : keys.addInputKeyword("compulsory","ARG","vector","the labels of the two vectors from which the outer product is being computed");
130 622 : keys.addInputKeyword("optional","MASK","matrix","a matrix that is used to used to determine which elements of the output matrix to compute");
131 311 : keys.add("compulsory","FUNC","x*y","the function of the input vectors that should be put in the elements of the outer product");
132 311 : keys.addFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",false,"set all diagonal elements to zero");
133 622 : keys.setValueDescription("matrix","a matrix containing the outer product of the two input vectors that was obtained using the function that was input");
134 311 : }
135 :
136 143 : void OuterProduct::registerKeywords( Keywords& keys ) {
137 143 : ActionShortcut::registerKeywords( keys );
138 143 : getKeywords( keys );
139 143 : keys.addActionNameSuffix("_MIN");
140 143 : keys.addActionNameSuffix("_MAX");
141 143 : keys.addActionNameSuffix("_FUNC");
142 143 : }
143 :
144 81 : OuterProduct::OuterProduct(const ActionOptions&ao):
145 : Action(ao),
146 81 : ActionShortcut(ao) {
147 :
148 : std::string func;
149 162 : parse("FUNC",func);
150 81 : if( func=="min") {
151 0 : readInputLine( getShortcutLabel() + ": OUTER_PRODUCT_MIN FUNC=" + func + " " + convertInputLineToString() );
152 81 : } else if( func=="max" ) {
153 4 : readInputLine( getShortcutLabel() + ": OUTER_PRODUCT_MAX FUNC=" + func + " " + convertInputLineToString() );
154 : } else {
155 158 : readInputLine( getShortcutLabel() + ": OUTER_PRODUCT_FUNC FUNC=" + func + " " + convertInputLineToString() );
156 : }
157 81 : }
158 :
159 : class OutputProductMin {
160 : public:
161 : static void registerKeywords( Keywords& keys );
162 : void setup( const std::vector<std::size_t>& shape,
163 : const std::string& func,
164 : OuterProductBase<OutputProductMin>* action );
165 : static void calculate( bool noderiv,
166 : const OutputProductMin& actdata,
167 : View<double> vals,
168 : MatrixElementOutput& output );
169 : };
170 :
171 : typedef OuterProductBase<OutputProductMin> opmin;
172 : PLUMED_REGISTER_ACTION(opmin,"OUTER_PRODUCT_MIN")
173 :
174 2 : void OutputProductMin::registerKeywords( Keywords& keys ) {
175 2 : OuterProduct::getKeywords( keys );
176 2 : }
177 :
178 0 : void OutputProductMin::setup( const std::vector<std::size_t>& shape,
179 : const std::string& func,
180 : OuterProductBase<OutputProductMin>* action ) {
181 0 : plumed_assert( func=="min" );
182 0 : action->log.printf(" taking minimum of two input vectors \n");
183 0 : }
184 :
185 0 : void OutputProductMin::calculate( bool noderiv,
186 : const OutputProductMin& actdata,
187 : View<double> vals,
188 : MatrixElementOutput& output ) {
189 0 : if( vals[0]<vals[1] ) {
190 0 : output.derivs[0][0] = 1;
191 0 : output.derivs[0][1] = 0;
192 0 : output.values[0] = vals[0];
193 0 : return;
194 : }
195 0 : output.derivs[0][0] = 0;
196 0 : output.derivs[0][1] = 1;
197 0 : output.values[0] = vals[1];
198 : }
199 :
200 : class OutputProductMax {
201 : public:
202 : static void registerKeywords( Keywords& keys );
203 : void setup( const std::vector<std::size_t>& shape,
204 : const std::string& func,
205 : OuterProductBase<OutputProductMax>* action );
206 : static void calculate( bool noderiv,
207 : const OutputProductMax& actdata,
208 : View<const double> vals,
209 : MatrixElementOutput& output );
210 : };
211 :
212 : typedef OuterProductBase<OutputProductMax> opmax;
213 : PLUMED_REGISTER_ACTION(opmax,"OUTER_PRODUCT_MAX")
214 :
215 6 : void OutputProductMax::registerKeywords( Keywords& keys ) {
216 6 : OuterProduct::getKeywords( keys );
217 6 : }
218 :
219 2 : void OutputProductMax::setup( const std::vector<std::size_t>& shape,
220 : const std::string& func,
221 : OuterProductBase<OutputProductMax>* action ) {
222 2 : plumed_assert( func=="max" );
223 2 : action->log.printf(" taking maximum of two input vectors \n");
224 2 : }
225 :
226 315192 : void OutputProductMax::calculate( bool noderiv,
227 : const OutputProductMax& actdata,
228 : View<const double> vals,
229 : MatrixElementOutput& output ) {
230 315192 : if( vals[0]>vals[1] ) {
231 2055 : output.derivs[0][0] = 1;
232 2055 : output.derivs[0][1] = 0;
233 2055 : output.values[0] = vals[0];
234 2055 : return;
235 : }
236 313137 : output.derivs[0][0] = 0;
237 313137 : output.derivs[0][1] = 1;
238 313137 : output.values[0] = vals[1];
239 : }
240 :
241 : class OutputProductFunc {
242 : public:
243 : std::string inputf;
244 : LeptonCall function;
245 : static void registerKeywords( Keywords& keys );
246 : void setup( const std::vector<std::size_t>& shape,
247 : const std::string& func,
248 : OuterProductBase<OutputProductFunc>* action );
249 : static void calculate( bool noderiv,
250 : const OutputProductFunc& actdata,
251 : View<const double> vals,
252 : MatrixElementOutput& output );
253 79 : OutputProductFunc& operator=( const OutputProductFunc& m ) {
254 79 : inputf = m.inputf;
255 79 : std::vector<std::string> var(2);
256 : var[0]="x";
257 : var[1]="y";
258 79 : function.set( inputf, var );
259 79 : return *this;
260 79 : }
261 : };
262 :
263 : typedef OuterProductBase<OutputProductFunc> opfunc;
264 : PLUMED_REGISTER_ACTION(opfunc,"OUTER_PRODUCT_FUNC")
265 :
266 160 : void OutputProductFunc::registerKeywords( Keywords& keys ) {
267 160 : OuterProduct::getKeywords( keys );
268 160 : }
269 :
270 79 : void OutputProductFunc::setup( const std::vector<std::size_t>& shape,
271 : const std::string& func,
272 : OuterProductBase<OutputProductFunc>* action ) {
273 79 : action->log.printf(" with function : %s \n", func.c_str() );
274 79 : inputf = func;
275 79 : std::vector<std::string> var(2);
276 : var[0]="x";
277 : var[1]="y";
278 79 : function.set( func, var, action );
279 79 : }
280 :
281 10497626 : void OutputProductFunc::calculate( bool noderiv,
282 : const OutputProductFunc& actdata,
283 : View<const double> vals,
284 : MatrixElementOutput& output ) {
285 10497626 : output.values[0] = actdata.function.evaluate( vals );
286 10497626 : output.derivs[0][0] = actdata.function.evaluateDeriv( 0, vals );
287 10497626 : output.derivs[0][1] = actdata.function.evaluateDeriv( 1, vals );
288 10497626 : }
289 :
290 : }
291 : }
|
