Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "core/ActionPilot.h"
      23             : #include "core/ActionWithArguments.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "core/PlumedMain.h"
      26             : #include "core/ActionSet.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace generic {
      30             : 
      31             : //+PLUMEDOC PRINTANALYSIS UPDATE_IF
      32             : /*
      33             : Conditional update of other actions.
      34             : 
      35             : This action can be used to enable and disable the update step for the following actions
      36             : depending on the value of its arguments. This allows one to extract snapshots
      37             : with value of some CVs in a given range.
      38             : 
      39             : When called with MORE_THAN and/or LESS_THAN keywords, this action starts an if block.
      40             : The block is executed if all the arguments are less than all the respective values
      41             : in the LESS_THAN keyword (if present) and all the arguments are more than all the
      42             : respective values in the MORE_THAN keyword (if present).
      43             : 
      44             : When called with the END flag, this action ends the corresponding IF block.
      45             : Notice that in this case one should also provide the ARG keyword. It is recommended to
      46             : use the same ARG keyword that was used to begin the block, so as to make the input more readable.
      47             : 
      48             : Of course, blocks can be nested at will.
      49             : 
      50             : There are many potential usages for this keyword. One might e.g. decide to analyze some variable
      51             : only when another variable is within a given range.
      52             : 
      53             : !!! caution ""
      54             : 
      55             :     Notice that not all the possible usage make
      56             :     particular sense. For example, conditionally updating a [METAD](METAD.md) keyword
      57             :     (that is: adding hills only if a variable is within a given range)
      58             :     can lead to unexpected results.
      59             : 
      60             : ## Examples
      61             : 
      62             : The following input instructs plumed to dump all the snapshots in which atoms 1 and 2 are within
      63             : 0.3 nm of each other:
      64             : 
      65             : ```plumed
      66             : d: DISTANCE ATOMS=1,2
      67             : UPDATE_IF ARG=d LESS_THAN=0.3 STRIDE=1
      68             : DUMPATOMS ATOMS=@mdatoms FILE=output.xyz
      69             : UPDATE_IF ARG=d END
      70             : ```
      71             : 
      72             : Notice, that the STRIDE is set equal to one by default.  You will likely always want to set it to this value.
      73             : 
      74             : The following input instructs plumed to dump all the snapshots where an atom is in touch with
      75             : the solute.
      76             : 
      77             : ```plumed
      78             : solute: GROUP ATOMS=1-124
      79             : coord: COORDINATION GROUPA=solute GROUPB=500 R_0=0.5
      80             : 
      81             : # A coordination number higher than 0.5 indicate that there is at least one
      82             : # atom of group `solute` at less than 5 A from atom number 500
      83             : 
      84             : UPDATE_IF ARG=coord MORE_THAN=0.5
      85             : DUMPATOMS ATOMS=solute,500 FILE=output.xyz
      86             : UPDATE_IF ARG=coord END
      87             : ```
      88             : 
      89             : */
      90             : //+ENDPLUMEDOC
      91             : 
      92             : class UpdateIf:
      93             :   public ActionPilot,
      94             :   public ActionWithArguments {
      95             :   std::vector<double> lower;
      96             :   std::vector<double> upper;
      97             :   bool on;
      98             :   bool end;
      99             : public:
     100             :   void prepare() override;
     101             :   void calculate() override;
     102             :   void beforeUpdate() override;
     103             :   explicit UpdateIf(const ActionOptions&);
     104             :   static void registerKeywords(Keywords& keys);
     105          28 :   void apply() override {}
     106             :   ~UpdateIf();
     107             : };
     108             : 
     109             : PLUMED_REGISTER_ACTION(UpdateIf,"UPDATE_IF")
     110             : 
     111          11 : void UpdateIf::registerKeywords(Keywords& keys) {
     112          11 :   Action::registerKeywords(keys);
     113          11 :   ActionPilot::registerKeywords(keys);
     114          11 :   ActionWithArguments::registerKeywords(keys);
     115          22 :   keys.addInputKeyword("compulsory","ARG","scalar","the labels of values that should be used to make the decision on whether to update or not");
     116          11 :   keys.add("compulsory","STRIDE","1","the frequency with which the quantities of interest should be output");
     117          11 :   keys.addFlag("END",false,"end");
     118          11 :   keys.add("optional","LESS_THAN","upper bound");
     119          11 :   keys.add("optional","MORE_THAN","lower bound");
     120          11 : }
     121             : 
     122           7 : UpdateIf::UpdateIf(const ActionOptions&ao):
     123             :   Action(ao),
     124             :   ActionPilot(ao),
     125             :   ActionWithArguments(ao),
     126           7 :   on(false),
     127           7 :   end(false) {
     128           7 :   parseFlag("END",end);
     129           7 :   parseVector("LESS_THAN",upper);
     130           7 :   parseVector("MORE_THAN",lower);
     131           7 :   if(end) {
     132             :     // an END without a matching opening line would remove the flag of the
     133             :     // enclosing level at every step (undefined behavior in PlumedMain::update)
     134             :     int open=0;
     135          13 :     for(const auto & u : plumed.getActionSet().select<UpdateIf*>()) {
     136          17 :       open += ( u->end ? -1 : 1 );
     137             :     }
     138           3 :     if(open<=0) {
     139           0 :       error("END without a matching UPDATE_IF");
     140             :     }
     141             :   }
     142           7 :   if(end && upper.size()!=0) {
     143           0 :     error("END and LESS_THAN are not compatible");
     144             :   }
     145           7 :   if(end && lower.size()!=0) {
     146           0 :     error("END and MORE_THAN are not compatible");
     147             :   }
     148           7 :   if(upper.size()==0) {
     149           4 :     upper.assign(getNumberOfArguments(),+std::numeric_limits<double>::max());
     150             :   }
     151           7 :   if(lower.size()==0) {
     152           5 :     lower.assign(getNumberOfArguments(),-std::numeric_limits<double>::max());
     153             :   }
     154           7 :   if(upper.size()!=getNumberOfArguments()) {
     155           0 :     error("LESS_THAN should have the same size as ARG");
     156             :   }
     157           7 :   if(lower.size()!=getNumberOfArguments()) {
     158           0 :     error("MORE_THAN should have the same size as ARG");
     159             :   }
     160          15 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     161           8 :     log<<"  boundaries for argument "<<i<<"    "<<lower[i]<<" "<<upper[i]<<"\n";
     162             :   }
     163           7 :   checkRead();
     164           7 : }
     165             : 
     166          28 : void UpdateIf::prepare() {
     167          28 :   on=false;
     168          28 : }
     169             : 
     170          28 : void UpdateIf::calculate() {
     171          28 :   on=true;
     172          60 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     173          32 :     if(getArgument(i)>=upper[i] || getArgument(i)<=lower[i]) {
     174          10 :       on=false;
     175             :     }
     176             :   }
     177          28 : }
     178             : 
     179          28 : void UpdateIf::beforeUpdate() {
     180          28 :   if(end) {
     181          12 :     plumed.updateFlagsPop();
     182             :   } else {
     183          16 :     if(on) {
     184           7 :       plumed.updateFlagsPush(plumed.updateFlagsTop());
     185             :     } else {
     186           9 :       plumed.updateFlagsPush(false);
     187             :     }
     188             :   }
     189          28 : }
     190             : 
     191             : 
     192          14 : UpdateIf::~UpdateIf() {
     193          14 : }
     194             : 
     195             : }
     196             : 
     197             : 
     198             : }