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