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 : #include "tools/LeptonCall.h"
25 :
26 : //+PLUMEDOC COLVAR OUTER_PRODUCT
27 : /*
28 : Calculate the outer product matrix of two vectors
29 :
30 : This action can be used to calculate the [outer product](https://en.wikipedia.org/wiki/Outer_product) of two
31 : vectors. As a (useless) example of what can be done with this action consider the following simple input:
32 :
33 : ```plumed
34 : d1: DISTANCE ATOMS1=1,2 ATOMS2=3,4
35 : d2: DISTANCE ATOMS1=5,6 ATOMS2=7,8 ATOMS3=9,10
36 : pp: OUTER_PRODUCT ARG=d1,d2
37 : PRINT ARG=pp FILE=colvar
38 : ```
39 :
40 : This input outputs a $2 \times 3$ matrix. If we call the 2 dimensional vector output by the first DISTANCE action
41 : $d$ and the 3 dimensional vector output by the second DISTANCE action $h$ then the $(i,j)$ element of the matrix
42 : output by the action with the label `pp` is given by:
43 :
44 : $$
45 : p_{ij} = d_i h_j
46 : $$
47 :
48 : These outer product matrices are useful if you are trying to calculate an adjacency matrix that says atoms are
49 : connected if they are within a certain distance of each other and if they satisfy a certain criterion. For example,
50 : consider the following input:
51 :
52 : ```plumed
53 : # Determine if atoms are within 5.3 nm of each other
54 : c1: CONTACT_MATRIX GROUP=1-100 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3}
55 : # Calculate the coordination numbers
56 : ones: ONES SIZE=100
57 : cc: MATRIX_VECTOR_PRODUCT ARG=c1,ones
58 : # Now use MORE_THAN to work out which atoms have a coordination number that is bigger than six
59 : cf: MORE_THAN ARG=cc SWITCH={RATIONAL D_0=5.5 R_0=0.5}
60 : # Now recalculate the contact matrix above as first step towards calculating adjacency matrix that measures if
61 : # atoms are close to each other and both have a coordination number that is bigger than six
62 : c2: CONTACT_MATRIX GROUP=1-100 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3}
63 : # 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
64 : cfm: OUTER_PRODUCT ARG=cf,cf
65 : # And multiply this by our contact matrix to determine the desired adjacency matrix
66 : m: CUSTOM ARG=c2,cfm FUNC=x*y PERIODIC=NO
67 : PRINT ARG=m FILE=colvar
68 : ```
69 :
70 : This input calculates a adjacency matrix which has element $(i,j)$ equal to one if atoms $i$ and $j$ have coordination numbers
71 : that are greater than 6 and if they are within 5.3 nm of each other.
72 :
73 : Notice that you can specify the function of the two input vectors that is to be calculated by using the `FUNC` keyword which accepts
74 : mathematical expressions of $x$ and $y$. In other words, the elements of the outer product are calculated using the lepton library
75 : 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
76 : the minimum of the two input variables or the maximum respectively.
77 :
78 : */
79 : //+ENDPLUMEDOC
80 :
81 : namespace PLMD {
82 : namespace matrixtools {
83 :
84 : class OuterProduct : public ActionWithMatrix {
85 : private:
86 : bool domin, domax, diagzero;
87 : LeptonCall function;
88 : unsigned nderivatives;
89 : bool stored_vector1, stored_vector2;
90 : public:
91 : static void registerKeywords( Keywords& keys );
92 : explicit OuterProduct(const ActionOptions&);
93 : unsigned getNumberOfDerivatives();
94 : void prepare() override ;
95 2298 : unsigned getNumberOfColumns() const override {
96 2298 : return getConstPntrToComponent(0)->getShape()[1];
97 : }
98 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
99 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
100 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
101 : };
102 :
103 : PLUMED_REGISTER_ACTION(OuterProduct,"OUTER_PRODUCT")
104 :
105 137 : void OuterProduct::registerKeywords( Keywords& keys ) {
106 137 : ActionWithMatrix::registerKeywords(keys);
107 274 : keys.addInputKeyword("compulsory","ARG","vector","the labels of the two vectors from which the outer product is being computed");
108 137 : keys.add("compulsory","FUNC","x*y","the function of the input vectors that should be put in the elements of the outer product");
109 137 : keys.addFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",false,"set all diagonal elements to zero");
110 274 : keys.setValueDescription("matrix","a matrix containing the outer product of the two input vectors that was obtained using the function that was input");
111 137 : }
112 :
113 77 : OuterProduct::OuterProduct(const ActionOptions&ao):
114 : Action(ao),
115 : ActionWithMatrix(ao),
116 77 : domin(false),
117 77 : domax(false) {
118 77 : if( getNumberOfArguments()!=2 ) {
119 0 : error("should be two arguments to this action, a matrix and a vector");
120 : }
121 77 : if( getPntrToArgument(0)->getRank()!=1 || getPntrToArgument(0)->hasDerivatives() ) {
122 0 : error("first argument to this action should be a vector");
123 : }
124 77 : if( getPntrToArgument(1)->getRank()!=1 || getPntrToArgument(1)->hasDerivatives() ) {
125 0 : error("first argument to this action should be a vector");
126 : }
127 :
128 : std::string func;
129 154 : parse("FUNC",func);
130 77 : if( func=="min") {
131 0 : domin=true;
132 0 : log.printf(" taking minimum of two input vectors \n");
133 77 : } else if( func=="max" ) {
134 2 : domax=true;
135 2 : log.printf(" taking maximum of two input vectors \n");
136 : } else {
137 75 : log.printf(" with function : %s \n", func.c_str() );
138 75 : std::vector<std::string> var(2);
139 : var[0]="x";
140 : var[1]="y";
141 75 : function.set( func, var, this );
142 75 : }
143 77 : parseFlag("ELEMENTS_ON_DIAGONAL_ARE_ZERO",diagzero);
144 77 : if( diagzero ) {
145 2 : log.printf(" setting diagonal elements equal to zero\n");
146 : }
147 :
148 77 : std::vector<unsigned> shape(2);
149 77 : shape[0]=getPntrToArgument(0)->getShape()[0];
150 77 : shape[1]=getPntrToArgument(1)->getShape()[0];
151 77 : addValue( shape );
152 77 : setNotPeriodic();
153 77 : nderivatives = buildArgumentStore(0);
154 77 : std::string headstr=getFirstActionInChain()->getLabel();
155 77 : stored_vector1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
156 77 : stored_vector2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
157 77 : if( getPntrToArgument(0)->isDerivativeZeroWhenValueIsZero() || getPntrToArgument(1)->isDerivativeZeroWhenValueIsZero() ) {
158 15 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
159 : }
160 77 : }
161 :
162 96 : unsigned OuterProduct::getNumberOfDerivatives() {
163 96 : return nderivatives;
164 : }
165 :
166 158 : void OuterProduct::prepare() {
167 158 : ActionWithVector::prepare();
168 158 : Value* myval=getPntrToComponent(0);
169 158 : if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0] && myval->getShape()[1]==getPntrToArgument(1)->getShape()[0] ) {
170 141 : return;
171 : }
172 17 : std::vector<unsigned> shape(2);
173 17 : shape[0] = getPntrToArgument(0)->getShape()[0];
174 17 : shape[1] = getPntrToArgument(1)->getShape()[0];
175 17 : myval->setShape( shape );
176 : }
177 :
178 27151 : void OuterProduct::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
179 27151 : unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[0];
180 27151 : if( diagzero ) {
181 990 : if( indices.size()!=size_v ) {
182 20 : indices.resize( size_v );
183 : }
184 : unsigned k=1;
185 99000 : for(unsigned i=0; i<size_v; ++i) {
186 98010 : if( task_index==i ) {
187 990 : continue ;
188 : }
189 97020 : indices[k] = size_v + i;
190 97020 : k++;
191 : }
192 : myvals.setSplitIndex( size_v );
193 : } else {
194 26161 : if( indices.size()!=size_v+1 ) {
195 249 : indices.resize( size_v+1 );
196 : }
197 1690193 : for(unsigned i=0; i<size_v; ++i) {
198 1664032 : indices[i+1] = start_n + i;
199 : }
200 : myvals.setSplitIndex( size_v + 1 );
201 : }
202 27151 : }
203 :
204 6874326 : void OuterProduct::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
205 6874326 : unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(), ind2=index2;
206 6874326 : if( index2>=getPntrToArgument(0)->getShape()[0] ) {
207 1976448 : ind2 = index2 - getPntrToArgument(0)->getShape()[0];
208 : }
209 6874326 : if( diagzero && index1==ind2 ) {
210 6511300 : return;
211 : }
212 :
213 : double fval;
214 6874326 : unsigned jarg = 0, kelem = index1;
215 6874326 : bool jstore=stored_vector1;
216 6874326 : std::vector<double> args(2);
217 6874326 : args[0] = getArgumentElement( 0, index1, myvals );
218 6874326 : args[1] = getArgumentElement( 1, ind2, myvals );
219 6874326 : if( domin ) {
220 0 : fval=args[0];
221 0 : if( args[1]<args[0] ) {
222 : fval=args[1];
223 0 : jarg=1;
224 0 : kelem=ind2;
225 0 : jstore=stored_vector2;
226 : }
227 6874326 : } else if( domax ) {
228 315192 : fval=args[0];
229 315192 : if( args[1]>args[0] ) {
230 : fval=args[1];
231 2055 : jarg=1;
232 2055 : kelem=ind2;
233 2055 : jstore=stored_vector2;
234 : }
235 : } else {
236 6559134 : fval=function.evaluate( args );
237 : }
238 :
239 6874326 : myvals.addValue( ostrn, fval );
240 6874326 : if( doNotCalculateDerivatives() ) {
241 : return ;
242 : }
243 :
244 366326 : if( domin || domax ) {
245 0 : addDerivativeOnVectorArgument( jstore, 0, jarg, kelem, 1.0, myvals );
246 : } else {
247 366326 : addDerivativeOnVectorArgument( stored_vector1, 0, 0, index1, function.evaluateDeriv( 0, args ), myvals );
248 366326 : addDerivativeOnVectorArgument( stored_vector2, 0, 1, ind2, function.evaluateDeriv( 1, args ), myvals );
249 : }
250 366326 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
251 : return ;
252 : }
253 363026 : unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
254 363026 : myvals.getMatrixRowDerivativeIndices( nmat )[nmat_ind] = arg_deriv_starts[1] + ind2;
255 363026 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+1 );
256 : }
257 :
258 39963 : void OuterProduct::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
259 39963 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
260 : return ;
261 : }
262 11402 : unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
263 11402 : myvals.getMatrixRowDerivativeIndices( nmat )[nmat_ind] = ival;
264 11402 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind+1 );
265 : }
266 :
267 : }
268 : }
|
