Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-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 "MultiColvarShortcuts.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "core/ActionSet.h"
      25             : #include "core/Group.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace multicolvar {
      29             : 
      30             : namespace {
      31             : // Read the value of BETA from the specification of MIN, MAX or ALT_MIN.
      32             : // The value is returned as it was written, so that it can be used in a function.
      33           7 : std::string readBeta( const std::string& what, const std::string& spec, ActionShortcut* action ) {
      34           7 :   std::vector<std::string> words=Tools::getWords( spec );
      35             :   std::string beta;
      36          14 :   if( !Tools::parse(words,"BETA",beta,0) || beta.length()==0 ) {
      37           3 :     action->error("did not find BETA in the specification of " + what + ".  Use " + what + "={BETA=...}");
      38             :   }
      39           6 :   if( words.size()>0 ) {
      40           5 :     action->error("found " + words[0] + " in the specification of " + what + " but the only option is BETA");
      41             :   }
      42           5 :   return beta;
      43           7 : }
      44             : }
      45             : 
      46         888 : void MultiColvarShortcuts::shortcutKeywords( Keywords& keys ) {
      47         888 :   keys.add("numbered","LESS_THAN","calculate the number of variables that are less than a certain target value. "
      48             :            "This quantity is calculated using \\f$\\sum_i \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ "
      49             :            "is a switching function.");
      50        1776 :   keys.linkActionInDocs("LESS_THAN","LESS_THAN");
      51        1776 :   keys.reset_style("LESS_THAN","deprecated");
      52        1776 :   keys.addOutputComponent("lessthan","LESS_THAN","scalar","the number of colvars that have a value less than a threshold");
      53         888 :   keys.add("numbered","MORE_THAN","calculate the number of variables that are more than a certain target value. "
      54             :            "This quantity is calculated using \\f$\\sum_i 1 - \\sigma(s_i)\\f$, where \\f$\\sigma(s)\\f$ "
      55             :            "is a switching function.");
      56        1776 :   keys.linkActionInDocs("MORE_THAN","MORE_THAN");
      57        1776 :   keys.reset_style("MORE_THAN","deprecated");
      58        1776 :   keys.addOutputComponent("morethan","MORE_THAN","scalar","the number of colvars that have a value more than a threshold");
      59         888 :   keys.add("optional","ALT_MIN","calculate the minimum value. "
      60             :            "To make this quantity continuous the minimum is calculated using "
      61             :            "\\f$ \\textrm{min} = -\\frac{1}{\\beta} \\log \\sum_i \\exp\\left( -\\beta s_i \\right)  \\f$ "
      62             :            "The value of \\f$\\beta\\f$ in this function is specified using (BETA=\\f$\\beta\\f$).");
      63        1776 :   keys.reset_style("ALT_MIN","deprecated");
      64        1776 :   keys.addOutputComponent("altmin","ALT_MIN","scalar","the minimum value of the cv");
      65         888 :   keys.add("optional","MIN","calculate the minimum value. "
      66             :            "To make this quantity continuous the minimum is calculated using "
      67             :            "\\f$ \\textrm{min} = \\frac{\\beta}{ \\log \\sum_i \\exp\\left( \\frac{\\beta}{s_i} \\right) } \\f$ "
      68             :            "The value of \\f$\\beta\\f$ in this function is specified using (BETA=\\f$\\beta\\f$)");
      69        1776 :   keys.reset_style("MIN","deprecated");
      70        1776 :   keys.addOutputComponent("min","MIN","scalar","the minimum colvar");
      71         888 :   keys.add("optional","MAX","calculate the maximum value. "
      72             :            "To make this quantity continuous the maximum is calculated using "
      73             :            "\\f$ \\textrm{max} = \\beta \\log \\sum_i \\exp\\left( \\frac{s_i}{\\beta}\\right) \\f$ "
      74             :            "The value of \\f$\\beta\\f$ in this function is specified using (BETA=\\f$\\beta\\f$)");
      75        1776 :   keys.reset_style("MAX","deprecated");
      76        1776 :   keys.addOutputComponent("max","MAX","scalar","the maximum colvar");
      77         888 :   keys.add("numbered","BETWEEN","calculate the number of values that are within a certain range. "
      78             :            "These quantities are calculated using kernel density estimation as described in "
      79             :            "the documentation for KDE.");
      80        1776 :   keys.linkActionInDocs("BETWEEN","BETWEEN");
      81        1776 :   keys.reset_style("BETWEEN","deprecated");
      82        1776 :   keys.addOutputComponent("between","BETWEEN","scalar","the number of colvars that have a value that lies in a particular interval");
      83         888 :   keys.addFlag("HIGHEST",false,"this flag allows you to recover the highest of these variables.");
      84        1776 :   keys.reset_style("HIGHEST","deprecated");
      85        1776 :   keys.addOutputComponent("highest","HIGHEST","scalar","the largest of the colvars");
      86         888 :   keys.add("optional","HISTOGRAM","calculate a discretized histogram of the distribution of values. "
      87             :            "This shortcut allows you to calculates NBIN quantites like BETWEEN.");
      88        1776 :   keys.reset_style("HISTOGRAM","deprecated");
      89         888 :   keys.addFlag("LOWEST",false,"this flag allows you to recover the lowest of these variables.");
      90        1776 :   keys.reset_style("LOWEST","deprecated");
      91        1776 :   keys.addOutputComponent("lowest","LOWEST","scalar","the smallest of the colvars");
      92         888 :   keys.addFlag("SUM",false,"calculate the sum of all the quantities.");
      93        1776 :   keys.reset_style("SUM","deprecated");
      94        1776 :   keys.addOutputComponent("sum","SUM","scalar","the sum of the colvars");
      95         888 :   keys.addFlag("MEAN",false,"calculate the mean of all the quantities.");
      96        1776 :   keys.reset_style("MEAN","deprecated");
      97        1776 :   keys.addOutputComponent("mean","MEAN","scalar","the mean of the colvars");
      98         888 :   keys.needsAction("SUM");
      99         888 :   keys.needsAction("MEAN");
     100         888 :   keys.needsAction("CUSTOM");
     101         888 :   keys.needsAction("HIGHEST");
     102         888 :   keys.needsAction("LOWEST");
     103         888 :   keys.needsAction("LESS_THAN");
     104         888 :   keys.needsAction("MORE_THAN");
     105         888 :   keys.needsAction("BETWEEN");
     106         888 : }
     107             : 
     108         140 : void MultiColvarShortcuts::expandFunctions( const std::string& labout, const std::string& argin, const std::string& weights, ActionShortcut* action ) {
     109             :   std::map<std::string,std::string> keymap;
     110         140 :   readShortcutKeywords( keymap, action );
     111         140 :   expandFunctions( labout, argin, weights, keymap, action );
     112         135 : }
     113             : 
     114         260 : void MultiColvarShortcuts::readShortcutKeywords( std::map<std::string,std::string>& keymap, ActionShortcut* action ) {
     115         260 :   Keywords keys;
     116         260 :   shortcutKeywords( keys );
     117         260 :   action->readShortcutKeywords( keys, keymap );
     118         260 : }
     119             : 
     120         127 : void MultiColvarShortcuts::parseAtomList( const std::string& key, std::vector<std::string>& atoms, ActionShortcut* action ) {
     121         127 :   std::vector<std::string> astr;
     122         127 :   action->parseVector(key,astr);
     123         127 :   if( astr.size()==0 ) {
     124             :     return ;
     125             :   }
     126          37 :   Tools::interpretRanges( astr );
     127        1287 :   for(unsigned i=0; i<astr.size(); ++i) {
     128        1250 :     Group* mygr=action->plumed.getActionSet().selectWithLabel<Group*>(astr[i]);
     129        1250 :     if( mygr ) {
     130           5 :       std::vector<std::string> grstr( mygr->getGroupAtoms() );
     131         121 :       for(unsigned j=0; j<grstr.size(); ++j) {
     132         116 :         atoms.push_back(grstr[j]);
     133             :       }
     134           5 :     } else {
     135        1245 :       Group* mygr2=action->plumed.getActionSet().selectWithLabel<Group*>(astr[i] + "_grp");
     136        1245 :       if( mygr2 ) {
     137          10 :         std::vector<std::string> grstr( mygr2->getGroupAtoms() );
     138       16082 :         for(unsigned j=0; j<grstr.size(); ++j) {
     139       16072 :           atoms.push_back(grstr[j]);
     140             :         }
     141          10 :       } else {
     142        1235 :         atoms.push_back(astr[i]);
     143             :       }
     144             :     }
     145             :   }
     146         127 : }
     147             : 
     148         281 : void MultiColvarShortcuts::expandFunctions( const std::string& labout, const std::string& argin, const std::string& weights,
     149             :     const std::map<std::string,std::string>& keymap, ActionShortcut* action, std::vector<std::string>* outputs ) {
     150         281 :   if( keymap.empty() ) {
     151             :     return;
     152             :   }
     153             :   // Store the end of the label (labout + suffix) of each scalar that is the sum of the elements of a vector or matrix
     154             :   auto record = [outputs]( const std::string& suffix ) {
     155         168 :     if( outputs ) {
     156          15 :       outputs->push_back( suffix );
     157             :     }
     158             :   };
     159             :   // Parse LESS_THAN
     160         350 :   if( keymap.count("LESS_THAN") ) {
     161          26 :     std::string sum_arg = labout + "_lt", lt_string = keymap.find("LESS_THAN")->second;
     162          52 :     action->readInputLine( labout + "_lt: LESS_THAN ARG=" + argin + " SWITCH={" + lt_string + "}");
     163          13 :     if( weights.length()>0 ) {
     164           5 :       sum_arg = labout + "_wlt";
     165          20 :       action->readInputLine( labout + "_wlt: CUSTOM ARG=" + weights + "," + labout + "_lt FUNC=x*y PERIODIC=NO");
     166             :     }
     167          26 :     action->readInputLine( labout + "_lessthan: SUM ARG=" + sum_arg + " PERIODIC=NO");
     168          26 :     record("_lessthan");
     169             :   }
     170         350 :   if( keymap.count("LESS_THAN1") ) {
     171           3 :     for(unsigned i=1;; ++i) {
     172             :       std::string istr;
     173           8 :       Tools::convert( i, istr );
     174          16 :       if( !keymap.count("LESS_THAN" + istr ) ) {
     175             :         break;
     176             :       }
     177          15 :       std::string sum_arg = labout + "_lt" + istr, lt_string1 = keymap.find("LESS_THAN" + istr)->second;
     178          25 :       action->readInputLine( labout + "_lt" + istr + ": LESS_THAN ARG=" + argin + " SWITCH={" + lt_string1 + "}");
     179           5 :       if( weights.length()>0 ) {
     180           1 :         sum_arg = labout + "_wlt" + istr;
     181           6 :         action->readInputLine( labout + "_wlt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_lt" + istr + " FUNC=x*y PERIODIC=NO");
     182             :       }
     183          15 :       action->readInputLine( labout + "_lessthan-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
     184          10 :       record("_lessthan-" + istr);
     185           5 :     }
     186             :   }
     187             :   // Parse MORE_THAN
     188         350 :   if( keymap.count("MORE_THAN") ) {
     189          52 :     std::string sum_arg=labout + "_mt", mt_string = keymap.find("MORE_THAN")->second;
     190         104 :     action->readInputLine( labout + "_mt: MORE_THAN ARG=" + argin + " SWITCH={" + mt_string + "}");
     191          26 :     if( weights.length()>0 ) {
     192           3 :       sum_arg = labout + "_wmt";
     193          12 :       action->readInputLine( labout + "_wmt: CUSTOM ARG=" + weights + "," + labout + "_mt FUNC=x*y PERIODIC=NO" );
     194             :     }
     195          52 :     action->readInputLine( labout + "_morethan: SUM ARG=" + sum_arg + " PERIODIC=NO");
     196          52 :     record("_morethan");
     197             :   }
     198         350 :   if(  keymap.count("MORE_THAN1") ) {
     199           5 :     for(unsigned i=1;; ++i) {
     200             :       std::string istr;
     201          15 :       Tools::convert( i, istr );
     202          30 :       if( !keymap.count("MORE_THAN" + istr ) ) {
     203             :         break;
     204             :       }
     205          30 :       std::string sum_arg = labout + "_mt" + istr, mt_string1 = keymap.find("MORE_THAN" + istr)->second;
     206          50 :       action->readInputLine( labout + "_mt" + istr + ": MORE_THAN ARG=" + argin + " SWITCH={" + mt_string1 + "}");
     207          10 :       if( weights.length()>0 ) {
     208          10 :         sum_arg = labout + "_wmt" + istr;
     209          60 :         action->readInputLine( labout + "_wmt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_mt" + istr + " FUNC=x*y PERIODIC=NO");
     210             :       }
     211          30 :       action->readInputLine( labout + "_morethan-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
     212          20 :       record("_morethan-" + istr);
     213          10 :     }
     214             :   }
     215             :   // Parse ALT_MIN
     216         350 :   if( keymap.count("ALT_MIN") ) {
     217           1 :     if( weights.length()>0 ) {
     218           0 :       plumed_merror("cannot use ALT_MIN with this shortcut");
     219             :     }
     220           3 :     std::string beta_str = readBeta( "ALT_MIN", keymap.find("ALT_MIN")->second, action );
     221           3 :     action->readInputLine( labout + "_me_altmin: CUSTOM ARG=" + argin + " FUNC=exp(-x*" + beta_str + ") PERIODIC=NO");
     222           2 :     action->readInputLine( labout + "_mec_altmin: SUM ARG=" + labout + "_me_altmin PERIODIC=NO");
     223           4 :     action->readInputLine( labout + "_altmin: CUSTOM ARG=" + labout + "_mec_altmin FUNC=-log(x)/" + beta_str + " PERIODIC=NO");
     224             :   }
     225             :   // Parse MIN
     226         350 :   if( keymap.count("MIN") ) {
     227           2 :     if( weights.length()>0 ) {
     228           0 :       plumed_merror("cannot use MIN with this shortcut");
     229             :     }
     230           6 :     std::string beta_str = readBeta( "MIN", keymap.find("MIN")->second, action );
     231           3 :     action->readInputLine( labout + "_me_min: CUSTOM ARG=" + argin + " FUNC=exp(" + beta_str + "/x) PERIODIC=NO");
     232           2 :     action->readInputLine( labout + "_mec_min: SUM ARG=" + labout + "_me_min PERIODIC=NO");
     233           4 :     action->readInputLine( labout + "_min: CUSTOM ARG=" + labout + "_mec_min FUNC=" + beta_str + "/log(x) PERIODIC=NO");
     234             :   }
     235             :   // Parse MAX
     236         348 :   if( keymap.count("MAX") ) {
     237           4 :     if( weights.length()>0 ) {
     238           0 :       plumed_merror("cannot use MAX with this shortcut");
     239             :     }
     240          12 :     std::string beta_str = readBeta( "MAX", keymap.find("MAX")->second, action );
     241           9 :     action->readInputLine( labout + "_me_max: CUSTOM ARG=" + argin + " FUNC=exp(x/" + beta_str + ") PERIODIC=NO");
     242           6 :     action->readInputLine( labout + "_mec_max: SUM ARG=" + labout + "_me_max PERIODIC=NO");
     243          12 :     action->readInputLine( labout + "_max: CUSTOM ARG=" + labout + "_mec_max FUNC=" + beta_str  + "*log(x) PERIODIC=NO");
     244             :   }
     245             :   // Parse HIGHEST
     246         346 :   if( keymap.count("HIGHEST") ) {
     247           8 :     if( weights.length()>0 ) {
     248           0 :       plumed_merror("cannot use HIGHEST with this shortcut");
     249             :     }
     250          16 :     action->readInputLine( labout + "_highest: HIGHEST ARG=" + argin );
     251             :   }
     252             :   // Parse LOWEST
     253         346 :   if( keymap.count("LOWEST") ) {
     254           3 :     if( weights.length()>0 ) {
     255           0 :       plumed_merror("cannot use LOWEST with this shortcut");
     256             :     }
     257           6 :     action->readInputLine( labout + "_lowest: LOWEST ARG=" + argin );
     258             :   }
     259             :   // Parse SUM
     260         346 :   if( keymap.count("SUM") ) {
     261          62 :     std::string sum_arg=argin;
     262          62 :     if( weights.length()>0 ) {
     263          34 :       sum_arg = labout + "_wsum";
     264         136 :       action->readInputLine( labout + "_wsum: CUSTOM ARG=" + weights + "," + argin + " FUNC=x*y PERIODIC=NO");
     265             :     }
     266         124 :     action->readInputLine( labout + "_sum: SUM ARG=" + sum_arg + " PERIODIC=NO");
     267         124 :     record("_sum");
     268             :   }
     269             :   // Parse MEAN
     270         346 :   if( keymap.count("MEAN") ) {
     271          65 :     if( weights.length()>0 ) {
     272           0 :       plumed_merror("cannot use MEAN with this shortcut");
     273             :     }
     274         195 :     action->readInputLine( labout + "_mean: MEAN ARG=" + argin + " PERIODIC=NO");
     275             :   }
     276             :   // Parse BETWEEN
     277         346 :   if( keymap.count("BETWEEN") ) {
     278          12 :     std::string sum_arg=labout + "_bt";
     279          24 :     std::string bt_string = keymap.find("BETWEEN")->second;
     280          48 :     action->readInputLine( labout + "_bt: BETWEEN ARG=" + argin + " SWITCH={" + bt_string + "}" );
     281          12 :     if( weights.length()>0 ) {
     282           4 :       sum_arg = labout + "_wbt";
     283          16 :       action->readInputLine( labout + "_wbt: CUSTOM ARG=" + weights + "," + labout + "_bt FUNC=x*y PERIODIC=NO");
     284             :     }
     285          24 :     action->readInputLine( labout + "_between: SUM ARG=" + sum_arg + " PERIODIC=NO");
     286          24 :     record("_between");
     287             :   }
     288             : 
     289         346 :   if( keymap.count("BETWEEN1") ) {
     290           6 :     for(unsigned i=1;; ++i) {
     291             :       std::string istr;
     292          17 :       Tools::convert( i, istr );
     293          34 :       if( !keymap.count("BETWEEN" + istr) ) {
     294             :         break;
     295             :       }
     296          11 :       std::string sum_arg=labout + "_bt" + istr;
     297          22 :       std::string bt_string1 = keymap.find("BETWEEN" + istr)->second;
     298          55 :       action->readInputLine( labout + "_bt" + istr + ": BETWEEN ARG=" + argin + " SWITCH={" + bt_string1 + "}" );
     299          11 :       if( weights.length()>0 ) {
     300           1 :         sum_arg = labout + "_wbt" + istr;
     301           6 :         action->readInputLine( labout + "_wbt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_bt" + istr + " FUNC=x*y PERIODIC=NO");
     302             :       }
     303          33 :       action->readInputLine( labout + "_between-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
     304          22 :       record("_between-" + istr);
     305          11 :     }
     306             :   }
     307             :   // Parse HISTOGRAM
     308         346 :   if( keymap.count("HISTOGRAM") ) {
     309          33 :     std::vector<std::string> words=Tools::getWords( keymap.find("HISTOGRAM")->second );
     310             :     unsigned nbins;
     311          15 :     bool found=Tools::parse(words,"NBINS",nbins,0); // Need replica index
     312          15 :     if( !found ) {
     313           0 :       action->error("did not find NBINS in specification for HISTOGRAM");
     314             :     }
     315          15 :     if( nbins==0 ) {
     316           1 :       action->error("NBINS in specification for HISTOGRAM must be larger than zero");
     317             :     }
     318             :     double lower;
     319          14 :     found=Tools::parse(words,"LOWER",lower,0);
     320          14 :     if( !found ) {
     321           0 :       action->error("did not find LOWER in specification for HISTOGRAM");
     322             :     }
     323             :     double upper;
     324          14 :     found=Tools::parse(words,"UPPER",upper,0);
     325          14 :     if( !found ) {
     326           0 :       action->error("did not find UPPER in specification for HISTOGRAM");
     327             :     }
     328          14 :     if( upper<=lower ) {
     329           1 :       action->error("UPPER must be larger than LOWER in specification for HISTOGRAM");
     330             :     }
     331             : 
     332          13 :     double delr = ( upper - lower ) / static_cast<double>( nbins );
     333          13 :     double smear=0.5;
     334          13 :     found=Tools::parse(words,"SMEAR",smear,0);
     335          13 :     if( !found ) {
     336           6 :       smear = 0.5;
     337             :     }
     338          13 :     if( words.size()==0 ) {
     339           1 :       action->error("did not find the kernel type (e.g. GAUSSIAN) in specification for HISTOGRAM");
     340             :     }
     341          41 :     for(unsigned i=0; i<nbins; ++i) {
     342             :       std::string smstr, istr;
     343          29 :       Tools::convert( i+1, istr );
     344          29 :       Tools::convert( smear, smstr );
     345          58 :       std::string sum_arg=labout + "_bt" + istr;
     346             :       std::string low_str, high_str;
     347          29 :       Tools::convert( lower + i*delr, low_str );
     348          29 :       Tools::convert( lower + (i+1)*delr, high_str );
     349         232 :       action->readInputLine( labout + "_bt" + istr + ": BETWEEN ARG=" + argin + " SWITCH={" + words[0] + " LOWER=" + low_str + " UPPER=" + high_str + " SMEAR=" + smstr + "}");
     350          29 :       if( weights.length()>0 ) {
     351           0 :         sum_arg = labout + "_wbt" + istr;
     352           0 :         action->readInputLine( labout + "_wbt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_bt" + istr + " FUNC=x*y PERIODIC=NO");
     353             :       }
     354          87 :       action->readInputLine( labout + "_between-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
     355          58 :       record("_between-" + istr);
     356             :     }
     357          15 :   }
     358             : }
     359             : 
     360             : }
     361             : }