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