Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016,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 "PathProjectionCalculator.h"
      23             : #include "core/ActionWithValue.h"
      24             : #include "core/ActionWithArguments.h"
      25             : #include "core/ActionRegister.h"
      26             : #include "core/ActionSet.h"
      27             : #include "core/PlumedMain.h"
      28             : 
      29             : namespace PLMD {
      30             : namespace mapping {
      31             : 
      32          24 : void PathProjectionCalculator::registerKeywords(Keywords& keys) {
      33          48 :   keys.addInputKeyword("compulsory","ARG","matrix","the labels of the matrix that contains the vectors of displacements from each frame in the path");
      34          24 :   keys.add("compulsory","METRIC","the method to use for computing the displacement vectors between the reference frames");
      35          24 :   keys.add("compulsory","METRIC_COMPONENT","if the final action in your metric contains multiple components this keyword is used to specify the component that should be used");
      36          48 :   keys.addInputKeyword("compulsory","REFERENCE","vector","labels for actions that contain reference coordinates for each point on the path");
      37          24 : }
      38             : 
      39          10 : PathProjectionCalculator::PathProjectionCalculator( Action* act ):
      40          10 :   mypath_obj(NULL) {
      41          10 :   ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( act );
      42          10 :   if( aarg ) {
      43           6 :     mypath_obj = aarg->getPntrToArgument(0);
      44             :     // Check that we have only one argument as input
      45           6 :     if( aarg->getNumberOfArguments()!=1 ) {
      46           0 :       act->error("should only have one argument to this function");
      47             :     }
      48             :   }
      49             :   // Ensure that values are stored in base calculation and that PLUMED doesn't try to calculate this in the stream
      50          10 :   if( mypath_obj ) {
      51           6 :     mypath_obj->buildDataStore();
      52             :   }
      53             :   // Check that the input is a matrix
      54          10 :   if( mypath_obj )
      55           6 :     if( mypath_obj->getRank()!=2 ) {
      56           0 :       act->error("the input to this action should be a matrix");
      57             :     }
      58             :   // Get the labels for the reference points
      59             :   std::vector<std::string> reference_data;
      60          10 :   act->parseVector("REFERENCE", reference_data);
      61          10 :   std::vector<colvar::RMSDVector*> allrmsd = act->plumed.getActionSet().select<colvar::RMSDVector*>();
      62          10 :   ActionWithArguments::interpretArgumentList( reference_data, act->plumed.getActionSet(), act, refargs );
      63          26 :   for(unsigned i=0; i<refargs.size(); ++i ) {
      64          16 :     Action* thisact = refargs[i]->getPntrToAction();
      65          16 :     for(unsigned j=0; j<allrmsd.size(); ++j) {
      66           4 :       if( allrmsd[j]->checkForDependency(thisact) ) {
      67           4 :         rmsd_objects.push_back( allrmsd[j] );
      68             :         break;
      69             :       }
      70             :     }
      71          16 :     if( !refargs[i]->isConstant() ) {
      72           0 :       act->error("input" + refargs[i]->getName() + " is not constant");
      73             :     }
      74          16 :     if( refargs[i]->getRank()==2 && reference_data.size()!=1 ) {
      75           0 :       act->error("should only be one matrix in input to path projection object");
      76             :     }
      77          16 :     if( refargs[i]->getRank()>0 && refargs[i]->getShape()[0]!=refargs[0]->getShape()[0] ) {
      78           0 :       act->error("mismatch in number of reference frames in input to reference_data");
      79             :     }
      80             :   }
      81             :   // Create a plumed main object to compute distances between reference configurations
      82          10 :   int s=sizeof(double);
      83          10 :   metric.cmd("setRealPrecision",&s);
      84          10 :   metric.cmd("setMDEngine","plumed");
      85          10 :   int nat=0;
      86          10 :   metric.cmd("setNatoms",&nat);
      87          10 :   metric.cmd("setNoVirial");
      88          10 :   unsigned nargs=refargs.size();
      89          10 :   if( refargs[0]->getRank()==2 ) {
      90           5 :     nargs = refargs[0]->getShape()[1];
      91             :   }
      92             :   std::string str_nargs;
      93          10 :   Tools::convert( nargs, str_nargs );
      94          10 :   std::string period_str=" PERIODIC=NO";
      95          10 :   if( mypath_obj && mypath_obj->isPeriodic() ) {
      96             :     std::string min, max;
      97           1 :     mypath_obj->getDomain( min, max );
      98           2 :     period_str=" PERIODIC=" + min + "," + max;
      99             :   }
     100          20 :   metric.readInputLine("arg1: PUT UNIT=number SHAPE=" + str_nargs + period_str, true);
     101          20 :   metric.readInputLine("arg2: PUT UNIT=number SHAPE=" + str_nargs + period_str, true);
     102          10 :   double tstep=1.0;
     103          10 :   metric.cmd("setTimestep",&tstep);
     104             :   std::string inp;
     105          20 :   act->parse("METRIC",inp);
     106             :   inp += " ARG=arg2,arg1";
     107             :   const char* cinp=inp.c_str();
     108          10 :   std::vector<std::string> input=Tools::getWords(inp);
     109          10 :   if( input.size()==1 && !actionRegister().check(input[0]) ) {
     110           0 :     metric.cmd("setPlumedDat",cinp);
     111           0 :     metric.cmd("init");
     112             :   } else {
     113          10 :     metric.cmd("init");
     114          10 :     metric.cmd("readInputLine",cinp);
     115             :   }
     116             :   // Now setup stuff to retrieve the final displacement
     117          10 :   unsigned aind = metric.getActionSet().size()-1;
     118             :   while( true ) {
     119          15 :     const ActionShortcut* as=dynamic_cast<const ActionShortcut*>( metric.getActionSet()[aind].get() );
     120          15 :     if( !as ) {
     121             :       break ;
     122             :     }
     123           5 :     aind = aind - 1;
     124             :     plumed_assert( aind>=0 );
     125           5 :   }
     126          10 :   ActionWithValue* fav = dynamic_cast<ActionWithValue*>( metric.getActionSet()[aind].get() );
     127          10 :   if( !fav ) {
     128           0 :     act->error("final value should calculate relevant value that you want as reference");
     129             :   }
     130          10 :   std::string name = (fav->copyOutput(0))->getName();
     131          10 :   if( fav->getNumberOfComponents()>1 ) {
     132             :     std::string comp;
     133           5 :     act->parse("METRIC_COMPONENT",comp);
     134          10 :     name = fav->getLabel() + "." + comp;
     135             :   }
     136             :   long rank;
     137          20 :   metric.cmd("getDataRank " + name, &rank );
     138          10 :   if( rank==0 ) {
     139           0 :     rank=1;
     140             :   }
     141          10 :   std::vector<long> ishape( rank );
     142          20 :   metric.cmd("getDataShape " + name, &ishape[0] );
     143             :   unsigned nvals=1;
     144          20 :   for(unsigned i=0; i<ishape.size(); ++i) {
     145          10 :     nvals *= ishape[i];
     146             :   }
     147          10 :   data.resize( nvals );
     148          20 :   metric.cmd("setMemoryForData " + name, &data[0] );
     149          20 : }
     150             : 
     151        1114 : unsigned PathProjectionCalculator::getNumberOfFrames() const {
     152        1114 :   return refargs[0]->getShape()[0];
     153             : }
     154             : 
     155       11728 : void PathProjectionCalculator::computeVectorBetweenFrames( const unsigned& ifrom, const unsigned& ito ) {
     156       11728 :   int step = 1;
     157       11728 :   metric.cmd("setStep",&step);
     158       11728 :   std::vector<double> valdata1( data.size() ), valdata2( data.size() );
     159       11728 :   getReferenceConfiguration( ito, valdata2 );
     160       11728 :   getReferenceConfiguration( ifrom, valdata1 );
     161       11728 :   metric.cmd("setValue arg1", &valdata1[0] );
     162       11728 :   metric.cmd("setValue arg2", &valdata2[0] );
     163       11728 :   metric.cmd("calc");
     164       11728 : }
     165             : 
     166       11728 : void PathProjectionCalculator::getDisplaceVector( const unsigned& ifrom, const unsigned& ito, std::vector<double>& displace ) {
     167       11728 :   if( displace.size()!=data.size() ) {
     168        1669 :     displace.resize( data.size() );
     169             :   }
     170       11728 :   computeVectorBetweenFrames( ifrom, ito );
     171      404063 :   for(unsigned i=0; i<data.size(); ++i) {
     172      392335 :     displace[i] = data[i];
     173             :   }
     174       11728 : }
     175             : 
     176       27778 : void PathProjectionCalculator::getReferenceConfiguration( const unsigned& iframe, std::vector<double>& refpos ) const {
     177       27778 :   if( refpos.size()!=data.size() ) {
     178           2 :     refpos.resize( data.size() );
     179             :   }
     180       27778 :   if( refargs[0]->getRank()==2 ) {
     181      962880 :     for(unsigned i=0; i<refpos.size(); ++i) {
     182      938808 :       refpos[i] = refargs[0]->get( iframe*refpos.size() + i );
     183             :     }
     184             :   } else {
     185       11178 :     for(unsigned i=0; i<refpos.size(); ++i) {
     186        7472 :       refpos[i] = refargs[i]->get(iframe);
     187             :     }
     188             :   }
     189       27778 : }
     190             : 
     191        4322 : void PathProjectionCalculator::setReferenceConfiguration( const unsigned& iframe, std::vector<double>& refpos ) {
     192             :   plumed_dbg_assert( refpos.size()==data.size() );
     193        4322 :   if( refargs[0]->getRank()==2 ) {
     194      165360 :     for(unsigned i=0; i<refpos.size(); ++i) {
     195      161226 :       refargs[0]->set( iframe*refpos.size() + i, refpos[i] );
     196             :     }
     197             :   } else {
     198         572 :     for(unsigned i=0; i<refpos.size(); ++i) {
     199         384 :       refargs[i]->set( iframe, refpos[i] );
     200             :     }
     201             :   }
     202        4322 : }
     203             : 
     204          27 : void PathProjectionCalculator::updateDepedentRMSDObjects() {
     205          50 :   for(unsigned i=0; i<rmsd_objects.size(); ++i) {
     206          23 :     rmsd_objects[i]->setReferenceConfigurations();
     207             :   }
     208          27 : }
     209             : 
     210             : }
     211             : }
     212             :