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/ActionWithValue.h"
      23             : #include "core/ActionWithArguments.h"
      24             : #include "core/ActionRegister.h"
      25             : 
      26             : //+PLUMEDOC MCOLVAR VSTACK
      27             : /*
      28             : Create a matrix by stacking vectors together
      29             : 
      30             : This action takes 2 or more vectors with the same number of elements as input and outputs a matrix.
      31             : This output matrix is constructed stacking the input vectors vertically.  The first row of the output
      32             : matrix thus contains the first elements of all the input vector, the second row contains the second elements
      33             : and so on.  In other words, the first input vector is in the first column of the output matrix, the second
      34             : input vector is in the second column and so.
      35             : 
      36             : The following shows an example of how this action operates in practise. The [DISTANCE](DISTANCE.md) command below calculates
      37             : the vectors containing four pairs of atoms.  The VSTACK command is then used to construct a $4 \times 3$ matrix
      38             : that contains all the vectors. The 1,1, 1,2 and 1,3 components in this matrix contain the
      39             : $x$, $y$ and $z$ components of the vector connecting atoms 1 and 2. The 2,1, 2,2 and 2,3 contain the
      40             : $x$, $y$ and $z$ components of the vector connecting atoms 3 and 4 and so on.
      41             : 
      42             : ```plumed
      43             : d1: DISTANCE ...
      44             :    COMPONENTS
      45             :    ATOMS1=1,2 ATOMS2=3,4
      46             :    ATOMS3=5,6 ATOMS4=7,8
      47             : ...
      48             : m: VSTACK ARG=d1.x,d1.y,d1.z
      49             : 
      50             : PRINT ARG=m FILE=matrix.dat
      51             : ```
      52             : 
      53             : */
      54             : //+ENDPLUMEDOC
      55             : 
      56             : namespace PLMD {
      57             : namespace valtools {
      58             : 
      59             : class VStack :
      60             :   public ActionWithValue,
      61             :   public ActionWithArguments {
      62             : public:
      63             :   static void registerKeywords( Keywords& keys );
      64             : /// Constructor
      65             :   explicit VStack(const ActionOptions&);
      66             : /// Get the number of derivatives
      67         317 :   unsigned getNumberOfDerivatives() override {
      68         317 :     return 0;
      69             :   }
      70             : ///
      71             :   void prepare() override ;
      72             : ///
      73             :   void calculate() override ;
      74             : ///
      75             :   void apply() override ;
      76             : ///
      77             :   void getMatrixColumnTitles( std::vector<std::string>& argnames ) const override ;
      78             : };
      79             : 
      80             : PLUMED_REGISTER_ACTION(VStack,"VSTACK")
      81             : 
      82         296 : void VStack::registerKeywords( Keywords& keys ) {
      83         296 :   Action::registerKeywords( keys );
      84         296 :   ActionWithValue::registerKeywords( keys );
      85         296 :   ActionWithArguments::registerKeywords( keys );
      86         296 :   keys.remove("NUMERICAL_DERIVATIVES");
      87         592 :   keys.addInputKeyword("compulsory","ARG","scalar/vector","the values that you would like to stack together to construct the output matrix");
      88         592 :   keys.setValueDescription("matrix","a matrix that contains the input vectors in its columns");
      89         296 : }
      90             : 
      91         161 : VStack::VStack(const ActionOptions& ao):
      92             :   Action(ao),
      93             :   ActionWithValue(ao),
      94         161 :   ActionWithArguments(ao) {
      95         161 :   if( getNumberOfArguments()==0 ) {
      96           0 :     error("no arguments were specificed");
      97             :   }
      98         161 :   if( getPntrToArgument(0)->getRank()>1 ) {
      99           0 :     error("all arguments should be vectors");
     100             :   }
     101             :   unsigned nvals=1;
     102             :   bool periodic=false;
     103             :   std::string smin, smax;
     104         161 :   if( getPntrToArgument(0)->getRank()==1 ) {
     105         141 :     nvals = getPntrToArgument(0)->getShape()[0];
     106             :   }
     107         161 :   if( getPntrToArgument(0)->isPeriodic() ) {
     108             :     periodic=true;
     109          24 :     getPntrToArgument(0)->getDomain( smin, smax );
     110             :   }
     111             : 
     112             :   bool derivbool=true;
     113        1354 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     114        1193 :     if( getPntrToArgument(i)->getRank()>1 || (getPntrToArgument(i)->getRank()==1 && getPntrToArgument(i)->hasDerivatives()) ) {
     115           0 :       error("all arguments should be vectors");
     116             :     }
     117        1193 :     if( getPntrToArgument(i)->getRank()==0 ) {
     118          41 :       if( nvals!=1 ) {
     119           0 :         error("all input vector should have same number of elements");
     120             :       }
     121        1152 :     } else if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
     122           0 :       error("all input vector should have same number of elements");
     123             :     }
     124        1193 :     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        1142 :     } else if( getPntrToArgument(i)->isPeriodic() ) {
     134           0 :       error("one argument is not periodic but " + getPntrToArgument(i)->getName() + " is periodic");
     135             :     }
     136        1193 :     if( !getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) {
     137             :       derivbool=false;
     138             :     }
     139             :   }
     140             :   // And create a value to hold the matrix
     141         161 :   std::vector<std::size_t> shape(2);
     142         161 :   shape[0]=nvals;
     143         161 :   shape[1]=getNumberOfArguments();
     144         161 :   addValue( shape );
     145         161 :   if( periodic ) {
     146          24 :     setPeriodic( smin, smax );
     147             :   } else {
     148         137 :     setNotPeriodic();
     149             :   }
     150             :   // And store this value
     151         161 :   getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
     152         161 :   if( derivbool ) {
     153          68 :     getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
     154             :   }
     155         322 : }
     156             : 
     157          25 : void VStack::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
     158          62 :   for(unsigned j=0; j<getNumberOfArguments(); ++j) {
     159          37 :     if( (getPntrToArgument(j)->getPntrToAction())->getName()=="COLLECT" ) {
     160          19 :       ActionWithArguments* aa = dynamic_cast<ActionWithArguments*>( getPntrToArgument(j)->getPntrToAction() );
     161          19 :       plumed_assert( aa && aa->getNumberOfArguments()==1 );
     162          19 :       argnames.push_back( (aa->getPntrToArgument(0))->getName() );
     163             :     } else {
     164          18 :       argnames.push_back( getPntrToArgument(j)->getName() );
     165             :     }
     166             :   }
     167          25 : }
     168             : 
     169        3759 : void VStack::prepare() {
     170        3759 :   if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->getShape()[0]==getPntrToComponent(0)->getShape()[0] ) {
     171        3737 :     return ;
     172             :   }
     173          22 :   std::vector<std::size_t> shape(2);
     174          22 :   shape[0] = getPntrToArgument(0)->getShape()[0];
     175          22 :   shape[1] = getNumberOfArguments();
     176          22 :   getPntrToComponent(0)->setShape(shape);
     177          22 :   getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
     178          22 : }
     179             : 
     180        3752 : void VStack::calculate() {
     181             :   unsigned nvals=1;
     182        3752 :   if( getPntrToArgument(0)->getRank()==1 ) {
     183        3187 :     nvals = getPntrToArgument(0)->getShape()[0];
     184             :   }
     185             : 
     186        3752 :   Value* valout = getPntrToComponent(0);
     187        3752 :   unsigned nargs = getNumberOfArguments();
     188       57664 :   for(unsigned i=0; i<nargs; ++i) {
     189             :     unsigned ipos = i;
     190             :     Value* myarg = getPntrToArgument(i);
     191    11358911 :     for(unsigned j=0; j<nvals; ++j) {
     192    11304999 :       valout->set( ipos, myarg->get(j) );
     193    11304999 :       ipos += nargs;
     194             :     }
     195             :   }
     196        3752 : }
     197             : 
     198        3700 : void VStack::apply() {
     199        3700 :   if( !getPntrToComponent(0)->forcesWereAdded() ) {
     200             :     return;
     201             :   }
     202             : 
     203             :   unsigned nvals=1;
     204        2970 :   if( getPntrToArgument(0)->getRank()==1 ) {
     205        2424 :     nvals = getPntrToArgument(0)->getShape()[0];
     206             :   }
     207        2970 :   Value* valout = getPntrToComponent(0);
     208        2970 :   unsigned nargs = getNumberOfArguments();
     209       52248 :   for(unsigned i=0; i<nargs; ++i) {
     210             :     unsigned ipos = i;
     211             :     Value* myarg = getPntrToArgument(i);
     212     2885268 :     for(unsigned j=0; j<nvals; ++j) {
     213     2835990 :       myarg->addForce( j, valout->getForce( ipos ) );
     214     2835990 :       ipos += nargs;
     215             :     }
     216             :   }
     217             : }
     218             : 
     219             : }
     220             : }