Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-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 "ActionWithAveraging.h"
      23             : #include "analysis/DataCollectionObject.h"
      24             : #include "analysis/ReadAnalysisFrames.h"
      25             : #include "core/PlumedMain.h"
      26             : #include "core/ActionSet.h"
      27             : #include "bias/ReweightBase.h"
      28             : 
      29             : namespace PLMD {
      30             : namespace vesselbase {
      31             : 
      32         113 : void ActionWithAveraging::registerKeywords( Keywords& keys ) {
      33         113 :   Action::registerKeywords( keys );
      34         113 :   ActionPilot::registerKeywords( keys );
      35         113 :   ActionAtomistic::registerKeywords( keys );
      36         113 :   ActionWithArguments::registerKeywords( keys );
      37         113 :   ActionWithValue::registerKeywords( keys );
      38         113 :   ActionWithVessel::registerKeywords( keys );
      39         226 :   keys.add("compulsory","STRIDE","1","the frequency with which the data should be collected and added to the quantity being averaged");
      40         226 :   keys.add("compulsory","CLEAR","0","the frequency with which to clear all the accumulated data.  The default value "
      41             :            "of 0 implies that all the data will be used and that the grid will never be cleared");
      42         226 :   keys.add("optional","LOGWEIGHTS","list of actions that calculates log weights that should be used to weight configurations when calculating averages");
      43         226 :   keys.add("compulsory","NORMALIZATION","true","This controls how the data is normalized it can be set equal to true, false or ndata.  The differences between these options are explained in the manual page for \\ref HISTOGRAM");
      44         113 :   keys.remove("NUMERICAL_DERIVATIVES");
      45         113 : }
      46             : 
      47          73 : ActionWithAveraging::ActionWithAveraging( const ActionOptions& ao ):
      48             :   Action(ao),
      49             :   ActionPilot(ao),
      50             :   ActionAtomistic(ao),
      51             :   ActionWithArguments(ao),
      52             :   ActionWithValue(ao),
      53             :   ActionWithVessel(ao),
      54          73 :   myaverage(NULL),
      55          73 :   activated(false),
      56          73 :   my_analysis_object(NULL),
      57          73 :   normalization(t),
      58          73 :   useRunAllTasks(false),
      59          73 :   clearstride(0),
      60          73 :   lweight(0),cweight(0) {
      61         146 :   if( keywords.exists("CLEAR") ) {
      62          58 :     parse("CLEAR",clearstride);
      63          58 :     if( clearstride>0 ) {
      64          12 :       if( clearstride%getStride()!=0 ) {
      65           0 :         error("CLEAR parameter must be a multiple of STRIDE");
      66             :       }
      67          12 :       log.printf("  clearing grid every %u steps \n",clearstride);
      68             :     }
      69             :   }
      70          73 :   if( ActionWithAveraging::getNumberOfArguments()>0 ) {
      71          29 :     my_analysis_object=dynamic_cast<analysis::AnalysisBase*>( getPntrToArgument(0)->getPntrToAction() );
      72          38 :     for(unsigned i=1; i<ActionWithAveraging::getNumberOfArguments(); i++) {
      73           9 :       if( my_analysis_object && my_analysis_object->getLabel()!=(getPntrToArgument(i)->getPntrToAction())->getLabel() ) {
      74           0 :         error("all arguments should be from one single analysis object");
      75             :       }
      76             :     }
      77          29 :     if( my_analysis_object ) {
      78           7 :       if( getStride()!=1 ) {
      79           0 :         error("stride should not have been set when calculating average from analysis data");
      80             :       }
      81           7 :       setStride(0);
      82           7 :       addDependency( my_analysis_object );
      83             :     }
      84             :   }
      85         146 :   if( keywords.exists("LOGWEIGHTS") ) {
      86             :     std::vector<std::string> wwstr;
      87         116 :     parseVector("LOGWEIGHTS",wwstr);
      88          58 :     if( wwstr.size()>0 ) {
      89           7 :       log.printf("  reweighting using weights from ");
      90             :     }
      91          58 :     std::vector<Value*> arg( getArguments() );
      92          66 :     for(unsigned i=0; i<wwstr.size(); ++i) {
      93           8 :       ActionWithValue* val = plumed.getActionSet().selectWithLabel<ActionWithValue*>(wwstr[i]);
      94           8 :       if( !val ) {
      95           0 :         error("could not find value named");
      96             :       }
      97           8 :       bias::ReweightBase* iswham=dynamic_cast<bias::ReweightBase*>( val );
      98           8 :       if( iswham && iswham->buildsWeightStore() ) {
      99           0 :         error("to use wham you must gather data using COLLECT_FRAMES");
     100             :       }
     101           8 :       weights.push_back( val->copyOutput(val->getLabel()) );
     102           8 :       arg.push_back( val->copyOutput(val->getLabel()) );
     103           8 :       log.printf("%s ",wwstr[i].c_str() );
     104             :     }
     105          58 :     if( wwstr.size()>0 ) {
     106           7 :       log.printf("\n");
     107             :     } else {
     108          51 :       log.printf("  weights are all equal to one\n");
     109             :     }
     110          58 :     requestArguments( arg );
     111          58 :   }
     112         146 :   if( keywords.exists("NORMALIZATION") ) {
     113             :     std::string normstr;
     114         116 :     parse("NORMALIZATION",normstr);
     115          58 :     if( normstr=="true" ) {
     116          23 :       normalization=t;
     117          35 :     } else if( normstr=="false" ) {
     118           9 :       normalization=f;
     119          26 :     } else if( normstr=="ndata" ) {
     120          26 :       normalization=ndata;
     121             :     } else {
     122           0 :       error("invalid instruction for NORMALIZATION flag should be true, false, or ndata");
     123             :     }
     124             :   }
     125          73 : }
     126             : 
     127          53 : bool ActionWithAveraging::ignoreNormalization() const {
     128          53 :   if( normalization==f ) {
     129           9 :     return true;
     130             :   }
     131             :   return false;
     132             : }
     133             : 
     134          68 : void ActionWithAveraging::setAveragingAction( std::unique_ptr<AveragingVessel> av_vessel, const bool& usetasks ) {
     135             :   // cppcheck-suppress danglingLifetime
     136          68 :   myaverage=av_vessel.get();
     137          68 :   addVessel( std::move(av_vessel) );
     138          68 :   useRunAllTasks=usetasks;
     139          68 :   resizeFunctions();
     140          68 : }
     141             : 
     142        1240 : void ActionWithAveraging::lockRequests() {
     143             :   ActionAtomistic::lockRequests();
     144             :   ActionWithArguments::lockRequests();
     145        1240 : }
     146             : 
     147        1240 : void ActionWithAveraging::unlockRequests() {
     148             :   ActionAtomistic::unlockRequests();
     149             :   ActionWithArguments::unlockRequests();
     150        1240 : }
     151             : 
     152         374 : unsigned ActionWithAveraging::getNumberOfQuantities() const {
     153         374 :   if( my_analysis_object ) {
     154          26 :     return getNumberOfArguments()+2;
     155             :   }
     156             :   return 2;
     157             : }
     158             : 
     159           0 : void ActionWithAveraging::calculateNumericalDerivatives(PLMD::ActionWithValue*) {
     160           0 :   error("not possible to compute numerical derivatives for this action");
     161             : }
     162             : 
     163        1258 : void ActionWithAveraging::update() {
     164        1258 :   if( (clearstride!=1 && getStep()==0) || (!onStep() && !my_analysis_object) ) {
     165          81 :     return;
     166             :   }
     167        1177 :   if( my_analysis_object ) {
     168           7 :     analysis::ReadAnalysisFrames* myfram = dynamic_cast<analysis::ReadAnalysisFrames*>( my_analysis_object );
     169           7 :     if( !activated && !myfram && !onStep() ) {
     170             :       return ;
     171           7 :     } else if( !activated && !my_analysis_object->onStep() ) {
     172             :       return ;
     173             :     }
     174             :   }
     175             : 
     176             :   // Clear if it is time to reset
     177        1177 :   if( myaverage ) {
     178        1176 :     if( myaverage->wasreset() ) {
     179          82 :       clearAverage();
     180             :     }
     181             :   }
     182             :   // Calculate the weight for all reweighting
     183        1177 :   if ( weights.size()>0 && !my_analysis_object ) {
     184             :     double sum=0;
     185        3018 :     for(unsigned i=0; i<weights.size(); ++i) {
     186        2009 :       sum+=weights[i]->get();
     187             :     }
     188        1009 :     lweight=sum;
     189        1009 :     cweight = exp( sum );
     190             :   } else {
     191         168 :     lweight=0;
     192         168 :     cweight=1.0;
     193             :   }
     194             :   // Prepare the task list for averaging
     195        1177 :   if( my_analysis_object ) {
     196        2115 :     for(unsigned i=getFullNumberOfTasks(); i<my_analysis_object->getNumberOfDataPoints(); ++i) {
     197        2108 :       addTaskToList(i);
     198             :     }
     199           7 :     deactivateAllTasks();
     200           7 :     cweight=0;
     201        2115 :     for(unsigned i=0; i<my_analysis_object->getNumberOfDataPoints(); ++i) {
     202        2108 :       taskFlags[i]=1;
     203        2108 :       cweight += my_analysis_object->getWeight(i);
     204             :     }
     205           7 :     lockContributors();
     206             :   }
     207             :   // Prepare to do the averaging
     208        1177 :   prepareForAveraging();
     209             :   // Run all the tasks (if required
     210        1177 :   if( my_analysis_object || useRunAllTasks ) {
     211         145 :     runAllTasks();
     212             :   }
     213             :   // This the averaging if it is not done using task list
     214             :   else {
     215        1032 :     performOperations( true );
     216             :   }
     217             :   // Update the norm
     218        1177 :   double normt = cweight;
     219        1177 :   if( !my_analysis_object && normalization==ndata ) {
     220             :     normt = 1;
     221             :   }
     222        1177 :   if( myaverage && my_analysis_object ) {
     223             :     myaverage->setNorm( normt );
     224        1170 :   } else if( myaverage ) {
     225        1169 :     myaverage->setNorm( normt + myaverage->getNorm() );
     226             :   }
     227             :   // Finish the averaging
     228        1177 :   finishAveraging();
     229             :   // By resetting here we are ensuring that the grid will be cleared at the start of the next step
     230        1177 :   if( myaverage ) {
     231        1176 :     if( getStride()==0 || (clearstride>0 && getStep()%clearstride==0) ) {
     232          49 :       myaverage->reset();
     233             :     }
     234             :   }
     235             : }
     236             : 
     237       60010 : void ActionWithAveraging::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
     238       60010 :   if( my_analysis_object ) {
     239        2108 :     const analysis::DataCollectionObject& mystore=my_analysis_object->getStoredData( current, false );
     240        4216 :     for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     241        2108 :       myvals.setValue( 1+i, mystore.getArgumentValue( ActionWithArguments::getArguments()[i]->getName() ) );
     242             :     }
     243        2108 :     myvals.setValue( 0, my_analysis_object->getWeight(current) );
     244        2108 :     if( normalization==f ) {
     245           0 :       myvals.setValue( 1+getNumberOfArguments(), 1.0 );
     246             :     } else {
     247        2108 :       myvals.setValue( 1+getNumberOfArguments(), 1.0 / cweight );
     248             :     }
     249        2108 :     accumulateAverage( myvals );
     250             :   } else {
     251       57902 :     runTask( current, myvals );
     252             :   }
     253       60010 : }
     254             : 
     255         102 : void ActionWithAveraging::clearAverage() {
     256         102 :   plumed_assert( myaverage->wasreset() );
     257         102 :   myaverage->clear();
     258         102 : }
     259             : 
     260           0 : void ActionWithAveraging::performOperations( const bool& from_update ) {
     261           0 :   plumed_error();
     262             : }
     263             : 
     264          58 : void ActionWithAveraging::runFinalJobs() {
     265          58 :   if( my_analysis_object && getStride()==0 ) {
     266           7 :     activated=true;
     267           7 :     update();
     268             :   }
     269          58 : }
     270             : 
     271             : }
     272             : }