Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "Action.h"
      23             : #include "ActionAtomistic.h"
      24             : #include "ActionWithValue.h"
      25             : #include "ActionWithArguments.h"
      26             : #include "ActionWithVirtualAtom.h"
      27             : #include "ActionForInterface.h"
      28             : #include "DomainDecomposition.h"
      29             : #include "PbcAction.h"
      30             : #include "ActionToPutData.h"
      31             : #include "ActionToGetData.h"
      32             : #include "PlumedMain.h"
      33             : #include "tools/Log.h"
      34             : #include "tools/Exception.h"
      35             : #include "tools/Communicator.h"
      36             : #include "ActionSet.h"
      37             : #include <iomanip>
      38             : #include <iostream>
      39             : #include <algorithm>
      40             : namespace PLMD {
      41             : 
      42             : Keywords ActionOptions::emptyKeys;
      43             : 
      44       66035 : ActionOptions::ActionOptions(PlumedMain&p,const std::vector<std::string>&l):
      45       66035 :   plumed(p),
      46       66035 :   line(l),
      47       66035 :   keys(emptyKeys) {
      48       66035 : }
      49             : 
      50       66033 : ActionOptions::ActionOptions(const ActionOptions&ao,const Keywords&mykeys):
      51       66033 :   plumed(ao.plumed),
      52       66033 :   line(ao.line),
      53       66033 :   keys(mykeys) {
      54       66033 : }
      55             : 
      56       92056 : void Action::registerKeywords( Keywords& keys ) {
      57       92056 :   plumed_assert( keys.size()==0 );
      58       92056 :   keys.add( "hidden", "LABEL", "a label for the action so that its output can be referenced in the input to other actions.  Actions with scalar output are referenced using their label only.  Actions with vector output must have a separate label for every component.  Individual components are then referred to using label.component" );
      59       92056 :   keys.reserve("optional","UPDATE_FROM","Only update this action from this time");
      60       92056 :   keys.reserve("optional","UPDATE_UNTIL","Only update this action until this time");
      61       92056 :   keys.reserve("optional","RESTART","allows per-action setting of restart (YES/NO/AUTO)");
      62       92056 : }
      63             : 
      64       66033 : Action::Action(const ActionOptions&ao):
      65       66033 :   actionName(ao.line[0]),
      66             :   //skipping the name with the +1
      67       66033 :   linemap(ao.line.begin()+1,ao.line.end()),
      68       66033 :   update_from(std::numeric_limits<double>::max()),
      69       66033 :   update_until(std::numeric_limits<double>::max()),
      70       66033 :   timestep(0),
      71       66033 :   active(false),
      72       66033 :   restart(ao.plumed.getRestart()),
      73       66033 :   doCheckPoint(ao.plumed.getCPT()),
      74       66033 :   never_activate(actionName=="CONSTANT"),
      75       66033 :   plumed(ao.plumed),
      76       66033 :   log(plumed.getLog()),
      77       66033 :   comm(plumed.comm),
      78       66033 :   multi_sim_comm(plumed.multi_sim_comm),
      79       66033 :   keywords(ao.keys) {
      80             :   // Retrieve the timestep and save it
      81       66033 :   resetStoredTimestep();
      82             : 
      83             :   //line.erase(line.begin());
      84             :   //making the line map
      85       66033 :   if( !keywords.exists("NO_ACTION_LOG") ) {
      86       40773 :     log.printf("Action %s\n",actionName.c_str());
      87       40773 :     if(ao.fullPath.length()>0) {
      88          49 :       log<<"  from library: "<<ao.fullPath<<"\n";
      89             :     }
      90             :   }
      91             : 
      92       66033 :   if(comm.Get_rank()==0) {
      93       47826 :     replica_index=multi_sim_comm.Get_rank();
      94             :   }
      95       66033 :   comm.Bcast(replica_index,0);
      96             : 
      97       66033 :   if ( keywords.exists("LABEL") ) {
      98      131454 :     parse("LABEL",actionLabel);
      99             :   }
     100       66033 :   if(actionLabel.length()==0) {
     101             :     std::string s;
     102        4082 :     Tools::convert(plumed.getActionSet().size()-plumed.getActionSet().getNumberOfActionsForInterface(),s);
     103        8164 :     actionLabel="@"+s;
     104       61951 :   } else if ( actionLabel.find(".")!=std::string::npos ) {
     105          30 :     warning("using full stop in an action label should be avaoided as . has a special meaning in PLUMED action labels");
     106             :   }
     107       66033 :   if( plumed.getActionSet().selectWithLabel<Action*>(actionLabel) ) {
     108           0 :     error("label " + actionLabel + " has been already used");
     109             :   }
     110       66033 :   if( !keywords.exists("NO_ACTION_LOG") ) {
     111       40773 :     log.printf("  with label %s\n",actionLabel.c_str());
     112             :   }
     113             :   // Keywords from older versions of PLUMED are read as the keyword that replaces them,
     114             :   // so that they also satisfy compulsory keywords that are read by the base classes (e.g. GRID for ARG)
     115      722217 :   for(const auto & key : keywords.getKeys()) {
     116      656184 :     const std::string replacement=keywords.getDeprecatedReplacement(key);
     117      656184 :     if( replacement.length()>0 && linemap.renameKey(key,replacement) ) {
     118          55 :       warning("keyword " + key + " is deprecated and is read as " + replacement + ": use " + replacement + " instead");
     119             :     }
     120             :   }
     121       66033 :   if ( keywords.exists("UPDATE_FROM") ) {
     122        6140 :     parse("UPDATE_FROM",update_from);
     123             :   }
     124       66033 :   if( !keywords.exists("NO_ACTION_LOG") && update_from!=std::numeric_limits<double>::max()) {
     125          10 :     log.printf("  only update from time %f\n",update_from);
     126             :   }
     127       66033 :   if ( keywords.exists("UPDATE_UNTIL") ) {
     128        6140 :     parse("UPDATE_UNTIL",update_until);
     129             :   }
     130       66033 :   if( !keywords.exists("NO_ACTION_LOG") && update_until!=std::numeric_limits<double>::max()) {
     131          21 :     log.printf("  only update until time %f\n",update_until);
     132             :   }
     133       66033 :   if ( keywords.exists("RESTART") ) {
     134        3180 :     std::string srestart="AUTO";
     135        3179 :     parse("RESTART",srestart);
     136        3179 :     if( plumed.parseOnlyMode() ) {
     137           0 :       restart=false;
     138        3179 :     } else if(srestart=="YES") {
     139         239 :       restart=true;
     140        2940 :     } else if(srestart=="NO") {
     141          35 :       restart=false;
     142        2905 :     } else if(srestart=="AUTO") {
     143             :       // do nothing, this is the default
     144             :     } else {
     145           2 :       error("RESTART should be either YES, NO, or AUTO");
     146             :     }
     147             :   }
     148       66035 : }
     149             : 
     150      124973 : void Action::resetStoredTimestep() {
     151      124973 :   ActionWithValue* ts = plumed.getActionSet().selectWithLabel<ActionWithValue*>("timestep");
     152      124973 :   if( ts ) {
     153      113094 :     timestep = (ts->copyOutput(0))->get();
     154             :   }
     155      124973 : }
     156             : 
     157      132064 : Action::~Action() {
     158       66032 :   if(files.size()!=0) {
     159           0 :     std::cerr<<"WARNING: some files open in action "+getLabel()+" where not properly closed. This could lead to data loss!!\n";
     160             :   }
     161      198096 : }
     162             : 
     163         142 : FILE* Action::fopen(const char *path, const char *mode) {
     164             :   bool write(false);
     165         284 :   for(const char*p=mode; *p; p++)
     166         142 :     if(*p=='w' || *p=='a' || *p=='+') {
     167             :       write=true;
     168             :     }
     169             :   FILE* fp;
     170         142 :   if(write && comm.Get_rank()!=0) {
     171           0 :     fp=plumed.fopen("/dev/null",mode);
     172             :   } else {
     173         142 :     fp=plumed.fopen(path,mode);
     174             :   }
     175         141 :   files.insert(fp);
     176         141 :   return fp;
     177             : }
     178             : 
     179         187 : int Action::fclose(FILE*fp) {
     180             :   files.erase(fp);
     181         187 :   return plumed.fclose(fp);
     182             : }
     183             : 
     184      103835 : void Action::fflush() {
     185      103835 :   for(const auto & p : files) {
     186           0 :     std::fflush(p);
     187             :   }
     188      103835 : }
     189             : 
     190           3 : std::string Action::getKeyword(const std::string& key) {
     191             :   // Check keyword has been registered
     192           3 :   plumed_massert(keywords.exists(key), "keyword " + key + " has not been registered");
     193             : 
     194           3 :   std::string outkey=linemap.getKeyword(key);
     195           3 :   if(!outkey.empty()) {
     196           3 :     return outkey;
     197             :   }
     198           0 :   if( keywords.style(key,"compulsory") ) {
     199           0 :     if( keywords.getDefaultValue(key,outkey) ) {
     200           0 :       if( outkey.length()==0 ) {
     201           0 :         error("keyword " + key + " has weird default value");
     202             :       }
     203           0 :       return key + "=" +  outkey;
     204             :     } else {
     205           0 :       error("keyword " + key + " is compulsory for this action");
     206             :     }
     207             :   }
     208           0 :   return "";
     209             : }
     210             : 
     211      158904 : void Action::parseFlag(const std::string&key,bool & t) {
     212             :   // Check keyword has been registered
     213      158904 :   plumed_massert(keywords.exists(key), "keyword " + key + " has not been registered");
     214             :   // Check keyword is a flag
     215      325102 :   plumed_massert( keywords.style(key,"deprecated") || keywords.style(key,"flag") || keywords.style(key,"hidden"), "keyword " + key + " is not a flag");
     216             : 
     217             :   // Read in the flag otherwise get the default value from the keywords object
     218      158904 :   t = linemap.readAndRemoveFlag(key);
     219      158904 :   if(!t) {
     220      269884 :     if( keywords.style(key,"nohtml") ) {
     221           0 :       t=false;
     222      134942 :     } else if ( !keywords.getLogicalDefault(key,t) ) {
     223           0 :       log.printf("ERROR in action %s with label %s : flag %s has no default",actionName.c_str(),actionLabel.c_str(),key.c_str() );
     224           0 :       plumed_error();
     225             :     }
     226             :   }
     227      158904 : }
     228             : 
     229      218299 : void Action::addDependency(Action*action) {
     230             :   bool found=false;
     231    35635027 :   for(const auto & d : after ) {
     232    35495286 :     if( action==d ) {
     233             :       found=true;
     234             :       break;
     235             :     }
     236             :   }
     237      218299 :   if( !found ) {
     238      139741 :     after.push_back(action);
     239             :   }
     240      218299 : }
     241             : 
     242         242 : bool Action::checkForDependency( Action* action ) {
     243         467 :   for(const auto & d : after) {
     244         238 :     if( action==d ) {
     245             :       return true;
     246             :     }
     247         225 :     if( d->checkForDependency(action) ) {
     248             :       return true;
     249             :     }
     250             :   }
     251             :   return false;
     252             : }
     253             : 
     254     7389344 : void Action::activate() {
     255             : // This is set to true if actions are only need to be computed in setup (during checkRead)
     256     7389344 :   if( never_activate ) {
     257             :     return;
     258             :   }
     259             : // preparation step is called only the first time an Action is activated.
     260             : // since it could change its dependences (e.g. in an ActionAtomistic which is
     261             : // accessing to a virtual atom), this is done just before dependencies are
     262             : // activated
     263     7266706 :   if(!active) {
     264     3690340 :     this->unlockRequests();
     265     3690340 :     prepare();
     266     3690340 :     this->lockRequests();
     267             :   } else {
     268             :     return;
     269             :   }
     270     9553939 :   for(const auto & p : after) {
     271     5863599 :     p->activate();
     272             :   }
     273     3690340 :   active=true;
     274             : }
     275             : 
     276       26761 : void Action::setOption(const std::string &s) {
     277             : // This overloads the action and activate some options
     278       26761 :   actionOptions.insert(s);
     279       52507 :   for(const auto & p : after) {
     280       25746 :     p->setOption(s);
     281             :   }
     282       26761 : }
     283             : 
     284           0 : void Action::clearOptions() {
     285             : // This overloads the action and activate some options
     286             :   actionOptions.clear();
     287           0 : }
     288             : 
     289             : 
     290       28193 : void Action::clearDependencies() {
     291       28193 :   after.clear();
     292       28193 : }
     293             : 
     294       84447 : void Action::checkRead() {
     295       84447 :   if (!linemap.empty()) {
     296           2 :     std::string msg="cannot understand the following words from the input line : ";
     297           4 :     msg += linemap.keyList(", ");
     298           2 :     error(msg);
     299             :   }
     300             : 
     301       84445 :   setupConstantValues(false);
     302       84445 : }
     303             : 
     304      122869 : void Action::setupConstantValues( const bool& have_atoms ) {
     305      122869 :   if( have_atoms ) {
     306             :     // This ensures that we switch off actions that only depend on constant when passed from the
     307             :     // MD code on the first step
     308       38424 :     ActionAtomistic* at = castToActionAtomistic();
     309       38424 :     ActionWithValue* av = castToActionWithValue();
     310       38424 :     if( at && av ) {
     311       16656 :       never_activate=av->getNumberOfComponents()>0;
     312       16697 :       for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
     313       16656 :         if( !av->copyOutput(i)->isConstant() ) {
     314       16615 :           never_activate=false;
     315       16615 :           break;
     316             :         }
     317             :       }
     318             :     }
     319             :   }
     320      122869 :   ActionWithArguments* aa = castToActionWithArguments();
     321      122869 :   if( aa && aa->getNumberOfArguments()>0 && getName()!="BIASVALUE" ) {
     322       25085 :     never_activate = aa->calculateConstantValues( have_atoms );
     323             :   }
     324      122869 : }
     325             : 
     326        1749 : int Action::getNumberOfReplicas() const {
     327        1749 :   return plumed.getNumberOfReplicas();
     328             : }
     329             : 
     330        1739 : int Action::getReplicaIndex() const {
     331        1739 :   return plumed.getReplicaIndex();
     332             : }
     333             : 
     334     7372185 : long long int Action::getStep()const {
     335     7372185 :   return plumed.getStep();
     336             : }
     337             : 
     338     5035765 : double Action::getTime()const {
     339     5035765 :   return timestep*getStep();
     340             : }
     341             : 
     342      279706 : double Action::getTimeStep()const {
     343      279706 :   return timestep;
     344             : }
     345             : 
     346        1305 : double Action::getkBT() {
     347        1305 :   double temp=-1.0;
     348        1305 :   if( keywords.exists("TEMP") ) {
     349        2482 :     parse("TEMP",temp);
     350             :   }
     351        3147 :   if(temp>=0.0 && keywords.style("TEMP","optional") ) {
     352         717 :     return getKBoltzmann()*temp;
     353             :   }
     354         588 :   ActionForInterface* kb=plumed.getActionSet().selectWithLabel<ActionForInterface*>("kBT");
     355             :   double kbt=0;
     356         588 :   if(kb) {
     357         188 :     kbt=(kb->copyOutput(0))->get();
     358             :   }
     359         996 :   if( temp>=0 && keywords.style("TEMP","compulsory") ) {
     360         204 :     double kB=getKBoltzmann();
     361         204 :     if( kbt>0 && std::abs(kbt-kB*temp)>1e-4) {
     362             :       std::string strt1, strt2;
     363           0 :       Tools::convert( temp, strt1 );
     364           0 :       Tools::convert( kbt/kB, strt2 );
     365           0 :       warning("using TEMP=" + strt1 + " while MD engine uses " + strt2 + "\n");
     366             :     }
     367         204 :     kbt = kB*temp;
     368         204 :     plumed_massert(kbt>0,"your MD engine does not pass the temperature to plumed, you must specify it using TEMP");
     369             :     return kbt;
     370             :   }
     371             :   return kbt;
     372             : }
     373             : 
     374           0 : void Action::exit(int c) {
     375           0 :   plumed.exit(c);
     376           0 : }
     377             : 
     378           0 : void Action::calculateNumericalDerivatives( ActionWithValue* a ) {
     379           0 :   plumed_merror("if you get here it means that you are trying to use numerical derivatives for a class that does not implement them");
     380             : }
     381             : 
     382     2093856 : void Action::prepare() {
     383     2093856 :   return;
     384             : }
     385             : 
     386         162 : [[noreturn]] void Action::error( const std::string & msg ) const {
     387         162 :   log.printf("ERROR in input to action %s with label %s : %s \n \n", actionName.c_str(), actionLabel.c_str(), msg.c_str() );
     388         810 :   plumed_merror("ERROR in input to action " + actionName + " with label " + actionLabel + " : " + msg );
     389             : }
     390             : 
     391        4696 : void Action::warning( const std::string & msg ) {
     392        4696 :   log.printf("WARNING for action %s with label %s : %s \n", actionName.c_str(), actionLabel.c_str(), msg.c_str() );
     393        4696 : }
     394             : 
     395           0 : void Action::calculateFromPDB( const PDB& pdb ) {
     396           0 :   activate();
     397           0 :   for(const auto & p : after) {
     398           0 :     ActionWithValue*av=p->castToActionWithValue();
     399           0 :     if(av) {
     400           0 :       av->clearInputForces();
     401           0 :       av->clearDerivatives();
     402             :     }
     403           0 :     p->readAtomsFromPDB( pdb );
     404           0 :     p->calculate();
     405             :   }
     406           0 :   readAtomsFromPDB( pdb );
     407           0 :   calculate();
     408           0 : }
     409             : 
     410       56829 : bool Action::getExchangeStep()const {
     411       56829 :   return plumed.getExchangeStep();
     412             : }
     413             : 
     414          73 : std::string Action::cite(const std::string&s) {
     415          73 :   return plumed.cite(s);
     416             : }
     417             : 
     418             : /// Check if action should be updated.
     419     3527143 : bool Action::checkUpdate()const {
     420     3527143 :   double t=getTime();
     421     3527143 :   if(t<update_until && (update_from==std::numeric_limits<double>::max() || t>=update_from)) {
     422             :     return true;
     423             :   } else {
     424        1302 :     return false;
     425             :   }
     426             : }
     427             : 
     428        3929 : bool Action::getCPT() const {
     429        3929 :   return plumed.getCPT();
     430             : }
     431             : 
     432    15302850 : const Units& Action::getUnits() const {
     433    15302850 :   return plumed.getUnits();
     434             : }
     435             : 
     436        2977 : bool Action::usingNaturalUnits() const {
     437        2977 :   return plumed.usingNaturalUnits();
     438             : }
     439             : 
     440        1385 : double Action::getKBoltzmann() const {
     441        1385 :   if( usingNaturalUnits() ) {
     442             :     return 1.0;
     443             :   } else {
     444        1377 :     return kBoltzmann/getUnits().getEnergy();
     445             :   }
     446             : }
     447             : 
     448          56 : std::string Action::writeInGraph() const {
     449          56 :   std::string nam=getName();
     450          56 :   std::size_t u=nam.find_last_of("_");
     451          56 :   std::string sub=nam.substr(u+1);
     452          56 :   if( sub=="SCALAR" || sub=="VECTOR" || sub=="GRID" ) {
     453           8 :     return nam.substr(0,u);
     454             :   }
     455          48 :   return nam;
     456             : }
     457             : 
     458             : }
     459             :