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 3680 : keys.needsAction("SUM"); keys.needsAction("MEAN"); keys.needsAction("CUSTOM"); keys.needsAction("HIGHEST"); keys.needsAction("LOWEST");
69 2208 : keys.needsAction("LESS_THAN"); keys.needsAction("MORE_THAN"); keys.needsAction("BETWEEN");
70 736 : }
71 :
72 125 : void MultiColvarShortcuts::expandFunctions( const std::string& labout, const std::string& argin, const std::string& weights, ActionShortcut* action ) {
73 125 : std::map<std::string,std::string> keymap; readShortcutKeywords( keymap, action ); expandFunctions( labout, argin, weights, keymap, action );
74 125 : }
75 :
76 231 : void MultiColvarShortcuts::readShortcutKeywords( std::map<std::string,std::string>& keymap, ActionShortcut* action ) {
77 231 : Keywords keys; shortcutKeywords( keys ); action->readShortcutKeywords( keys, keymap );
78 231 : }
79 :
80 94 : void MultiColvarShortcuts::parseAtomList( const std::string& key, std::vector<std::string>& atoms, ActionShortcut* action ) {
81 94 : std::vector<std::string> astr; action->parseVector(key,astr); if( astr.size()==0 ) return ;
82 28 : Tools::interpretRanges( astr );
83 1265 : for(unsigned i=0; i<astr.size(); ++i) {
84 1237 : Group* mygr=action->plumed.getActionSet().selectWithLabel<Group*>(astr[i]);
85 1237 : if( mygr ) {
86 1 : std::vector<std::string> grstr( mygr->getGroupAtoms() );
87 101 : for(unsigned j=0; j<grstr.size(); ++j) atoms.push_back(grstr[j]);
88 1 : } else {
89 1236 : Group* mygr2=action->plumed.getActionSet().selectWithLabel<Group*>(astr[i] + "_grp");
90 1236 : if( mygr2 ) {
91 10 : std::vector<std::string> grstr( mygr2->getGroupAtoms() );
92 16082 : for(unsigned j=0; j<grstr.size(); ++j) atoms.push_back(grstr[j]);
93 1236 : } else atoms.push_back(astr[i]);
94 : }
95 : }
96 94 : }
97 :
98 240 : void MultiColvarShortcuts::expandFunctions( const std::string& labout, const std::string& argin, const std::string& weights,
99 : const std::map<std::string,std::string>& keymap, ActionShortcut* action ) {
100 240 : if( keymap.empty() ) return;
101 : // Parse LESS_THAN
102 294 : if( keymap.count("LESS_THAN") ) {
103 16 : std::string sum_arg = labout + "_lt", lt_string = keymap.find("LESS_THAN")->second;
104 16 : action->readInputLine( labout + "_lt: LESS_THAN ARG=" + argin + " SWITCH={" + lt_string + "}");
105 8 : if( weights.length()>0 ) {
106 1 : sum_arg = labout + "_wlt";
107 2 : action->readInputLine( labout + "_wlt: CUSTOM ARG=" + weights + "," + labout + "_lt FUNC=x*y PERIODIC=NO");
108 : }
109 16 : action->readInputLine( labout + "_lessthan: SUM ARG=" + sum_arg + " PERIODIC=NO");
110 : }
111 294 : if( keymap.count("LESS_THAN1") ) {
112 2 : for(unsigned i=1;; ++i) {
113 6 : std::string istr; Tools::convert( i, istr );
114 12 : if( !keymap.count("LESS_THAN" + istr ) ) { break; }
115 12 : std::string sum_arg = labout + "_lt" + istr, lt_string1 = keymap.find("LESS_THAN" + istr)->second;
116 8 : action->readInputLine( labout + "_lt" + istr + ": LESS_THAN ARG=" + argin + " SWITCH={" + lt_string1 + "}");
117 4 : if( weights.length()>0 ) {
118 0 : sum_arg = labout + "_wlt" + istr;
119 0 : action->readInputLine( labout + "_wlt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_lt" + istr + " FUNC=x*y PERIODIC=NO");
120 : }
121 8 : action->readInputLine( labout + "_lessthan-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
122 4 : }
123 : }
124 : // Parse MORE_THAN
125 294 : if( keymap.count("MORE_THAN") ) {
126 48 : std::string sum_arg=labout + "_mt", mt_string = keymap.find("MORE_THAN")->second;
127 48 : action->readInputLine( labout + "_mt: MORE_THAN ARG=" + argin + " SWITCH={" + mt_string + "}");
128 24 : if( weights.length()>0 ) {
129 1 : sum_arg = labout + "_wmt";
130 2 : action->readInputLine( labout + "_wmt: CUSTOM ARG=" + weights + "," + labout + "_mt FUNC=x*y PERIODIC=NO" );
131 : }
132 48 : action->readInputLine( labout + "_morethan: SUM ARG=" + sum_arg + " PERIODIC=NO");
133 : }
134 294 : if( keymap.count("MORE_THAN1") ) {
135 0 : for(unsigned i=1;; ++i) {
136 0 : std::string istr; Tools::convert( i, istr );
137 0 : if( !keymap.count("MORE_THAN" + istr ) ) { break; }
138 0 : std::string sum_arg = labout + "_mt" + istr, mt_string1 = keymap.find("MORE_THAN" + istr)->second;
139 0 : action->readInputLine( labout + "_mt" + istr + ": MORE_THAN ARG=" + argin + " SWITCH={" + mt_string1 + "}");
140 0 : if( weights.length()>0 ) {
141 0 : sum_arg = labout + "_wmt" + istr;
142 0 : action->readInputLine( labout + "_wmt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_lt" + istr + " FUNC=x*y PERIODIC=NO");
143 : }
144 0 : action->readInputLine( labout + "_morethan-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
145 0 : }
146 : }
147 : // Parse ALT_MIN
148 294 : if( keymap.count("ALT_MIN") ) {
149 1 : if( weights.length()>0 ) plumed_merror("cannot use ALT_MIN with this shortcut");
150 2 : std::string amin_string = keymap.find("ALT_MIN")->second;
151 1 : std::size_t dd = amin_string.find("BETA"); std::string beta_str = amin_string.substr(dd+5);
152 1 : beta_str.erase(std::remove_if(beta_str.begin(), beta_str.end(), ::isspace), beta_str.end());
153 2 : action->readInputLine( labout + "_me_altmin: CUSTOM ARG=" + argin + " FUNC=exp(-x*" + beta_str + ") PERIODIC=NO");
154 2 : action->readInputLine( labout + "_mec_altmin: SUM ARG=" + labout + "_me_altmin PERIODIC=NO");
155 2 : action->readInputLine( labout + "_altmin: CUSTOM ARG=" + labout + "_mec_altmin FUNC=-log(x)/" + beta_str + " PERIODIC=NO");
156 : }
157 : // Parse MIN
158 294 : if( keymap.count("MIN") ) {
159 1 : if( weights.length()>0 ) plumed_merror("cannot use MIN with this shortcut");
160 2 : std::string min_string = keymap.find("MIN")->second;
161 1 : std::size_t dd = min_string.find("BETA"); std::string beta_str = min_string.substr(dd+5);
162 1 : beta_str.erase(std::remove_if(beta_str.begin(), beta_str.end(), ::isspace), beta_str.end());
163 2 : action->readInputLine( labout + "_me_min: CUSTOM ARG=" + argin + " FUNC=exp(" + beta_str + "/x) PERIODIC=NO");
164 2 : action->readInputLine( labout + "_mec_min: SUM ARG=" + labout + "_me_min PERIODIC=NO");
165 2 : action->readInputLine( labout + "_min: CUSTOM ARG=" + labout + "_mec_min FUNC=" + beta_str + "/log(x) PERIODIC=NO");
166 : }
167 : // Parse MAX
168 294 : if( keymap.count("MAX") ) {
169 3 : if( weights.length()>0 ) plumed_merror("cannot use MAX with this shortcut");
170 6 : std::string max_string = keymap.find("MAX")->second;
171 3 : std::size_t dd = max_string.find("BETA"); std::string beta_str = max_string.substr(dd+5);
172 3 : beta_str.erase(std::remove_if(beta_str.begin(), beta_str.end(), ::isspace), beta_str.end());
173 6 : action->readInputLine( labout + "_me_max: CUSTOM ARG=" + argin + " FUNC=exp(x/" + beta_str + ") PERIODIC=NO");
174 6 : action->readInputLine( labout + "_mec_max: SUM ARG=" + labout + "_me_max PERIODIC=NO");
175 6 : action->readInputLine( labout + "_max: CUSTOM ARG=" + labout + "_mec_max FUNC=" + beta_str + "*log(x) PERIODIC=NO");
176 : }
177 : // Parse HIGHEST
178 294 : if( keymap.count("HIGHEST") ) {
179 8 : if( weights.length()>0 ) plumed_merror("cannot use HIGHEST with this shortcut");
180 16 : action->readInputLine( labout + "_highest: HIGHEST ARG=" + argin );
181 : }
182 : // Parse LOWEST
183 294 : if( keymap.count("LOWEST") ) {
184 3 : if( weights.length()>0 ) plumed_merror("cannot use LOWEST with this shortcut");
185 6 : action->readInputLine( labout + "_lowest: LOWEST ARG=" + argin );
186 : }
187 : // Parse SUM
188 294 : if( keymap.count("SUM") ) {
189 44 : std::string sum_arg=argin;
190 44 : if( weights.length()>0 ) {
191 22 : sum_arg = labout + "_wsum";
192 44 : action->readInputLine( labout + "_wsum: CUSTOM ARG=" + weights + "," + argin + " FUNC=x*y PERIODIC=NO");
193 : }
194 88 : action->readInputLine( labout + "_sum: SUM ARG=" + sum_arg + " PERIODIC=NO");
195 : }
196 : // Parse MEAN
197 294 : if( keymap.count("MEAN") ) {
198 64 : if( weights.length()>0 ) plumed_merror("cannot use MEAN with this shortcut");
199 128 : action->readInputLine( labout + "_mean: MEAN ARG=" + argin + " PERIODIC=NO");
200 : }
201 : // Parse BETWEEN
202 294 : if( keymap.count("BETWEEN") ) {
203 16 : std::string sum_arg=labout + "_bt", bt_string = keymap.find("BETWEEN")->second;
204 16 : action->readInputLine( labout + "_bt: BETWEEN ARG=" + argin + " SWITCH={" + bt_string + "}" );
205 8 : if( weights.length()>0 ) {
206 0 : sum_arg = labout + "_wbt";
207 0 : action->readInputLine( labout + "_wbt: CUSTOM ARG=" + weights + "," + labout + "_bt FUNC=x*y PERIODIC=NO");
208 : }
209 16 : action->readInputLine( labout + "_between: SUM ARG=" + sum_arg + " PERIODIC=NO");
210 : }
211 : std::string bt_string1;
212 294 : if( keymap.count("BETWEEN1") ) {
213 6 : for(unsigned i=1;; ++i) {
214 17 : std::string istr; Tools::convert( i, istr );
215 34 : if( !keymap.count("BETWEEN" + istr) ) break;
216 33 : std::string sum_arg=labout + "_bt" + istr, bt_string1 = keymap.find("BETWEEN" + istr)->second;
217 22 : action->readInputLine( labout + "_bt" + istr + ": BETWEEN ARG=" + argin + " SWITCH={" + bt_string1 + "}" );
218 11 : if( weights.length()>0 ) {
219 2 : sum_arg = labout + "_wbt" + istr;
220 2 : action->readInputLine( labout + "_wbt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_bt" + istr + " FUNC=x*y PERIODIC=NO");
221 : }
222 22 : action->readInputLine( labout + "_between-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
223 11 : }
224 : }
225 : // Parse HISTOGRAM
226 294 : if( keymap.count("HISTOGRAM") ) {
227 24 : std::vector<std::string> words=Tools::getWords( keymap.find("HISTOGRAM")->second );
228 12 : unsigned nbins; bool found=Tools::parse(words,"NBINS",nbins,0); // Need replica index
229 12 : if( !found ) plumed_merror("did not find NBINS in specification for HISTOGRAM");
230 12 : double lower; found=Tools::parse(words,"LOWER",lower,0);
231 12 : if( !found ) plumed_merror("did not find LOWER in specification for HISTOGRAM");
232 12 : double upper; found=Tools::parse(words,"UPPER",upper,0);
233 12 : if( !found ) plumed_merror("did not find UPPER in specification for HISTOGRAM");
234 12 : double delr = ( upper - lower ) / static_cast<double>( nbins );
235 12 : double smear=0.5; found=Tools::parse(words,"SMEAR",smear,0);
236 12 : if( !found ) smear = 0.5;
237 41 : for(unsigned i=0; i<nbins; ++i) {
238 58 : std::string smstr, istr; Tools::convert( i+1, istr ); Tools::convert( smear, smstr ); std::string sum_arg=labout + "_bt" + istr;
239 29 : std::string low_str, high_str; Tools::convert( lower + i*delr, low_str ); Tools::convert( lower + (i+1)*delr, high_str );
240 58 : action->readInputLine( labout + "_bt" + istr + ": BETWEEN ARG=" + argin + " SWITCH={" + words[0] + " LOWER=" + low_str + " UPPER=" + high_str + " SMEAR=" + smstr + "}");
241 29 : if( weights.length()>0 ) {
242 0 : sum_arg = labout + "_wbt" + istr;
243 0 : action->readInputLine( labout + "_wbt" + istr + ": CUSTOM ARG=" + weights + "," + labout + "_bt" + istr + " FUNC=x*y PERIODIC=NO");
244 : }
245 58 : action->readInputLine( labout + "_between-" + istr + ": SUM ARG=" + sum_arg + " PERIODIC=NO");
246 : }
247 12 : }
248 : }
249 :
250 : }
251 : }
|
