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