Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2017 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 VSTACK
26 : /*
27 : Create a matrix by stacking vectors together
28 :
29 : This action takes 2 or more vectors with the same number of elements as input and outputs a matrix.
30 : This output matrix is constructed stacking the input vectors vertically. The first row of the output
31 : matrix thus contains the first elements of all the input vector, the second row contains the second elements
32 : and so on. In other words, the first input vector is in the first column of the output matrix, the second
33 : input vector is in the second column and so.
34 :
35 : The following shows an example of how this action operates in practise. The [DISTANCE](DISTANCE.md) command below calculates
36 : the vectors containing four pairs of atoms. The VSTACK command is then used to construct a $4 \times 3$ matrix
37 : that contains all the vectors. The 1,1, 1,2 and 1,3 components in this matrix contain the
38 : $x$, $y$ and $z$ components of the vector connecting atoms 1 and 2. The 2,1, 2,2 and 2,3 contain the
39 : $x$, $y$ and $z$ components of the vector connecting atoms 3 and 4 and so on.
40 :
41 : ```plumed
42 : d1: DISTANCE ...
43 : COMPONENTS
44 : ATOMS1=1,2 ATOMS2=3,4
45 : ATOMS3=5,6 ATOMS4=7,8
46 : ...
47 : m: VSTACK ARG=d1.x,d1.y,d1.z
48 :
49 : PRINT ARG=m FILE=matrix.dat
50 : ```
51 :
52 : */
53 : //+ENDPLUMEDOC
54 :
55 : namespace PLMD {
56 : namespace valtools {
57 :
58 : class VStack : public ActionWithMatrix {
59 : private:
60 : std::vector<bool> stored;
61 : public:
62 : static void registerKeywords( Keywords& keys );
63 : /// Constructor
64 : explicit VStack(const ActionOptions&);
65 : /// Get the number of derivatives
66 294 : unsigned getNumberOfDerivatives() override {
67 294 : return 0;
68 : }
69 : ///
70 : void prepare() override ;
71 : ///
72 1270505 : unsigned getNumberOfColumns() const override {
73 1270505 : return getNumberOfArguments();
74 : }
75 : ///
76 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const override ;
77 : ///
78 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override ;
79 : ///
80 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
81 : ///
82 : void getMatrixColumnTitles( std::vector<std::string>& argnames ) const override ;
83 : };
84 :
85 : PLUMED_REGISTER_ACTION(VStack,"VSTACK")
86 :
87 254 : void VStack::registerKeywords( Keywords& keys ) {
88 254 : ActionWithMatrix::registerKeywords( keys );
89 508 : keys.addInputKeyword("compulsory","ARG","scalar/vector","the values that you would like to stack together to construct the output matrix");
90 508 : keys.setValueDescription("matrix","a matrix that contains the input vectors in its columns");
91 254 : }
92 :
93 139 : VStack::VStack(const ActionOptions& ao):
94 : Action(ao),
95 139 : ActionWithMatrix(ao) {
96 139 : if( getNumberOfArguments()==0 ) {
97 0 : error("no arguments were specificed");
98 : }
99 139 : if( getPntrToArgument(0)->getRank()>1 ) {
100 0 : error("all arguments should be vectors");
101 : }
102 : unsigned nvals=1;
103 : bool periodic=false;
104 : std::string smin, smax;
105 139 : if( getPntrToArgument(0)->getRank()==1 ) {
106 119 : nvals = getPntrToArgument(0)->getShape()[0];
107 : }
108 139 : if( getPntrToArgument(0)->isPeriodic() ) {
109 : periodic=true;
110 24 : getPntrToArgument(0)->getDomain( smin, smax );
111 : }
112 :
113 1228 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
114 1089 : if( getPntrToArgument(i)->getRank()>1 || (getPntrToArgument(i)->getRank()==1 && getPntrToArgument(i)->hasDerivatives()) ) {
115 0 : error("all arguments should be vectors");
116 : }
117 1089 : if( getPntrToArgument(i)->getRank()==0 ) {
118 41 : if( nvals!=1 ) {
119 0 : error("all input vector should have same number of elements");
120 : }
121 1048 : } else if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
122 0 : error("all input vector should have same number of elements");
123 : }
124 1089 : if( periodic ) {
125 51 : if( !getPntrToArgument(i)->isPeriodic() ) {
126 0 : error("one argument is periodic but " + getPntrToArgument(i)->getName() + " is not periodic");
127 : }
128 : std::string tmin, tmax;
129 51 : getPntrToArgument(i)->getDomain( tmin, tmax );
130 51 : if( tmin!=smin || tmax!=smax ) {
131 0 : error("domain of argument " + getPntrToArgument(i)->getName() + " is different from domain for all other arguments");
132 : }
133 1038 : } else if( getPntrToArgument(i)->isPeriodic() ) {
134 0 : error("one argument is not periodic but " + getPntrToArgument(i)->getName() + " is periodic");
135 : }
136 : }
137 : // And create a value to hold the matrix
138 139 : std::vector<unsigned> shape(2);
139 139 : shape[0]=nvals;
140 139 : shape[1]=getNumberOfArguments();
141 139 : addValue( shape );
142 139 : if( periodic ) {
143 24 : setPeriodic( smin, smax );
144 : } else {
145 115 : setNotPeriodic();
146 : }
147 : // And store this value
148 139 : getPntrToComponent(0)->buildDataStore();
149 139 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
150 : // Setup everything so we can build the store
151 139 : done_in_chain=true;
152 139 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(0)->getPntrToAction() );
153 139 : if( av ) {
154 99 : const ActionWithVector* head0 = av->getFirstActionInChain();
155 996 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
156 928 : ActionWithVector* avv=dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
157 928 : if( !avv ) {
158 4 : continue;
159 : }
160 924 : if( head0!=avv->getFirstActionInChain() ) {
161 31 : done_in_chain=false;
162 31 : break;
163 : }
164 : }
165 : } else {
166 40 : done_in_chain=false;
167 : }
168 139 : unsigned nder = buildArgumentStore(0);
169 : // This checks which values have been stored
170 139 : stored.resize( getNumberOfArguments() );
171 139 : std::string headstr=getFirstActionInChain()->getLabel();
172 1228 : for(unsigned i=0; i<stored.size(); ++i) {
173 1089 : stored[i] = getPntrToArgument(i)->ignoreStoredValue( headstr );
174 : }
175 139 : }
176 :
177 23 : void VStack::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
178 60 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
179 37 : if( (getPntrToArgument(j)->getPntrToAction())->getName()=="COLLECT" ) {
180 17 : ActionWithArguments* aa = dynamic_cast<ActionWithArguments*>( getPntrToArgument(j)->getPntrToAction() );
181 17 : plumed_assert( aa && aa->getNumberOfArguments()==1 );
182 17 : argnames.push_back( (aa->getPntrToArgument(0))->getName() );
183 : } else {
184 20 : argnames.push_back( getPntrToArgument(j)->getName() );
185 : }
186 : }
187 23 : }
188 :
189 3100 : void VStack::prepare() {
190 3100 : ActionWithVector::prepare();
191 3100 : if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->getShape()[0]==getPntrToComponent(0)->getShape()[0] ) {
192 3082 : return ;
193 : }
194 18 : std::vector<unsigned> shape(2);
195 18 : shape[0] = getPntrToArgument(0)->getShape()[0];
196 18 : shape[1] = getNumberOfArguments();
197 18 : getPntrToComponent(0)->setShape(shape);
198 18 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
199 : }
200 :
201 214418 : void VStack::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
202 214418 : unsigned nargs = getNumberOfArguments();
203 214418 : unsigned nvals = getConstPntrToComponent(0)->getShape()[0];
204 214418 : if( indices.size()!=nargs+1 ) {
205 51622 : indices.resize( nargs+1 );
206 : }
207 4287565 : for(unsigned i=0; i<nargs; ++i) {
208 4073147 : indices[i+1] = nvals + i;
209 : }
210 : myvals.setSplitIndex( nargs + 1 );
211 214418 : }
212 :
213 4073147 : void VStack::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
214 4073147 : unsigned ind2 = index2;
215 4073147 : if( index2>=getConstPntrToComponent(0)->getShape()[0] ) {
216 4073147 : ind2 = index2 - getConstPntrToComponent(0)->getShape()[0];
217 : }
218 4073147 : myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), getArgumentElement( ind2, index1, myvals ) );
219 :
220 4073147 : if( doNotCalculateDerivatives() ) {
221 380548 : return;
222 : }
223 3692599 : addDerivativeOnVectorArgument( stored[ind2], 0, ind2, index1, 1.0, myvals );
224 : }
225 :
226 214418 : void VStack::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
227 214418 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
228 : return ;
229 : }
230 :
231 3 : unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
232 : std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
233 3 : plumed_assert( nmat_ind<matrix_indices.size() );
234 12 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
235 : bool found=false;
236 9 : ActionWithValue* iav = getPntrToArgument(i)->getPntrToAction();
237 9 : for(unsigned j=0; j<i; ++j) {
238 6 : if( iav==getPntrToArgument(j)->getPntrToAction() ) {
239 : found=true;
240 : break;
241 : }
242 : }
243 9 : if( found ) {
244 6 : continue ;
245 : }
246 :
247 : unsigned istrn = getPntrToArgument(i)->getPositionInStream();
248 48 : for(unsigned k=0; k<myvals.getNumberActive(istrn); ++k) {
249 45 : matrix_indices[nmat_ind] = myvals.getActiveIndex(istrn,k);
250 45 : nmat_ind++;
251 : }
252 : }
253 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
254 : }
255 :
256 : }
257 : }
|
