Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "Bias.h"
23 : #include "core/ActionRegister.h"
24 :
25 : namespace PLMD {
26 : namespace bias {
27 :
28 : //+PLUMEDOC BIAS UPPER_WALLS
29 : /*
30 : Defines a wall for the value of one or more collective variables,
31 : which limits the region of the phase space accessible during the simulation.
32 :
33 : If the inputs to the arguments are scalars the restraining potential starts acting on the system when the value of the CV is greater
34 : than a certain limit $a_i$ (AT) minus an offset $o_i$ (OFFSET). The expression for the bias due to the wall is given by:
35 :
36 : $$
37 : \sum_i {k_i}((x_i-a_i+o_i)/s_i)^e_i
38 : $$
39 :
40 : In this expression $k_i$ (KAPPA) is an energy constant in internal unit of the code, $s_i$ (EPS) a rescaling factor and
41 : $e_i$ (EXP) the exponent determining the power law. By default: EXP = 2, EPS = 1.0, OFFSET = 0.
42 :
43 : The following input tells plumed to add both a lower and an upper walls on the distance
44 : between atoms 3 and 5 and the distance between atoms 2 and 4. The lower and upper limits
45 : are defined at different values. The strength of the walls is the same for the four cases.
46 : It also tells plumed to print the energy of the walls.
47 :
48 : ```plumed
49 : d1: DISTANCE ATOMS=3,5
50 : d2: DISTANCE ATOMS=2,4
51 : uwall: UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0
52 : lwall: LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0
53 : PRINT ARG=uwall.bias,lwall.bias
54 : ```
55 :
56 : Alternatively, if the input to this action is a vector as in the example shown below:
57 :
58 : ```plumed
59 : d1: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
60 : uwall: UPPER_WALLS ARG=d1 AT=1.0 KAPPA=150.0 EXP=2 EPS=1 OFFSET=0
61 : PRINT ARG=uwall.bias
62 : ```
63 :
64 : The bias that acts is given by the following expression:
65 :
66 : $$
67 : \sum_i \sum_j {k_i}((x_{ij}-a_i+o_i)/s_i)^e_i
68 : $$
69 :
70 : The sum over $i$ here is runs over the two arguments, while the sum over $j$ runs over the three components of the input vectors.
71 : Notice, that regardless of whether the input is a scalar, vector or matrix the number of $k_i$, $a_i$, $o_i$, $e_i$ and $s_i$ values must be
72 : equal to the number of arguments to the action.
73 :
74 : ## Periodic arguments
75 :
76 : If an argument is periodic (e.g. a [TORSION](TORSION.md)), the wall acts on the value of the argument in its periodic domain
77 : and periodic boundary conditions are __not__ applied to the difference $x_i-a_i$. The value of AT must thus be within the
78 : periodic domain of the argument. The following input restrains a torsion within the interval $[-2,2]$:
79 :
80 : ```plumed
81 : t: TORSION ATOMS=1,2,3,4
82 : uwall: UPPER_WALLS ARG=t AT=2.0 KAPPA=150.0
83 : lwall: LOWER_WALLS ARG=t AT=-2.0 KAPPA=150.0
84 : PRINT ARG=t,uwall.bias,lwall.bias
85 : ```
86 :
87 : Notice that the bias is discontinuous at the boundary of the periodic domain (here at $\pm\pi$). You should thus make sure that
88 : the walls keep the argument far from this boundary, e.g. by using both an upper and a lower wall as in the input above.
89 :
90 : Versions of PLUMED before 2.9.7 and 2.10.3 applied periodic boundary conditions to the difference $x_i-a_i$. With that choice, each wall acts on half of the
91 : periodic domain, and the region left free by an upper and a lower wall is smaller than the one between AT values
92 : (see [#658](https://github.com/plumed/plumed2/issues/658)).
93 :
94 : */
95 : //+ENDPLUMEDOC
96 :
97 : //+PLUMEDOC BIAS LOWER_WALLS
98 : /*
99 : Defines a wall for the value of one or more collective variables,
100 : which limits the region of the phase space accessible during the simulation.
101 :
102 : If the inputs to the arguments are scalars the restraining potential starts acting on the system when the value of the CV is
103 : lower than a certain limit $a_i$ (AT) minus an offset $o_i$ (OFFSET). The expression for the bias due to the wall is given by:
104 :
105 : $$
106 : \sum_i {k_i}|(x_i-a_i-o_i)/s_i|^e_i
107 : $$
108 :
109 : In this expression $k_i$ (KAPPA) is an energy constant in internal unit of the code, $s_i$ (EPS) a rescaling factor and
110 : $e_i$ (EXP) the exponent determining the power law. By default: EXP = 2, EPS = 1.0, OFFSET = 0.
111 :
112 : The following input tells plumed to add both a lower and an upper walls on the distance
113 : between atoms 3 and 5 and the distance between atoms 2 and 4. The lower and upper limits
114 : are defined at different values. The strength of the walls is the same for the four cases.
115 : It also tells plumed to print the energy of the walls.
116 :
117 : ```plumed
118 : d1: DISTANCE ATOMS=3,5
119 : d2: DISTANCE ATOMS=2,4
120 : uwall: UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0
121 : lwall: LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0
122 : PRINT ARG=uwall.bias,lwall.bias
123 : ```
124 :
125 : Alternatively, if the input to this action is a vector as in the example shown below:
126 :
127 : ```plumed
128 : d1: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
129 : lwall: LOWER_WALLS ARG=d1 AT=1.0 KAPPA=150.0 EXP=2 EPS=1 OFFSET=0
130 : PRINT ARG=lwall.bias
131 : ```
132 :
133 : The bias that acts is given by the following expression:
134 :
135 : $$
136 : \sum_i \sum_j {k_i}|(x_{ij}-a_i-o_i)/s_i|^e_i
137 : $$
138 :
139 : The sum over $i$ here is runs over the two arguments, while the sum over $j$ runs over the three components of the input vectors.
140 : Notice, that regardless of whether the input is a scalar, vector or matrix the number of $k_i$, $a_i$, $o_i$, $e_i$ and $s_i$ values must be
141 : equal to the number of arguments to the action.
142 :
143 : ## Periodic arguments
144 :
145 : If an argument is periodic (e.g. a [TORSION](TORSION.md)), the wall acts on the value of the argument in its periodic domain
146 : and periodic boundary conditions are __not__ applied to the difference $x_i-a_i$. The value of AT must thus be within the
147 : periodic domain of the argument. The following input restrains a torsion within the interval $[-2,2]$:
148 :
149 : ```plumed
150 : t: TORSION ATOMS=1,2,3,4
151 : uwall: UPPER_WALLS ARG=t AT=2.0 KAPPA=150.0
152 : lwall: LOWER_WALLS ARG=t AT=-2.0 KAPPA=150.0
153 : PRINT ARG=t,uwall.bias,lwall.bias
154 : ```
155 :
156 : Notice that the bias is discontinuous at the boundary of the periodic domain (here at $\pm\pi$). You should thus make sure that
157 : the walls keep the argument far from this boundary, e.g. by using both an upper and a lower wall as in the input above.
158 :
159 : Versions of PLUMED before 2.9.7 and 2.10.3 applied periodic boundary conditions to the difference $x_i-a_i$. With that choice, each wall acts on half of the
160 : periodic domain, and the region left free by an upper and a lower wall is smaller than the one between AT values
161 : (see [#658](https://github.com/plumed/plumed2/issues/658)).
162 :
163 : */
164 : //+ENDPLUMEDOC
165 :
166 : enum class WallType {UPPER,LOWER};
167 :
168 : template<WallType W>
169 : class WallsScalar : public Bias {
170 : std::vector<double> at;
171 : std::vector<double> kappa;
172 : std::vector<double> exp;
173 : std::vector<double> eps;
174 : std::vector<double> offset;
175 : static constexpr auto force2="force2";
176 : inline double walloff(double cv, double off) {
177 : if constexpr(W==WallType::UPPER) {
178 14443 : return cv+off;
179 : } else {
180 297 : return cv-off;
181 : }
182 : }
183 : inline double wallpow(double scale, double exponent) {
184 : if constexpr(W==WallType::UPPER) {
185 145 : return std::pow(scale,exponent);
186 : } else {
187 109 : return std::pow(-scale,exponent);
188 : }
189 : }
190 : inline bool check(double scale) {
191 : if constexpr(W==WallType::UPPER) {
192 : return scale > 0.0;
193 : } else {
194 : return scale < 0.0;
195 : }
196 : }
197 : public:
198 : explicit WallsScalar(const ActionOptions&);
199 : void calculate() override;
200 : static void registerKeywords(Keywords& keys);
201 : };
202 :
203 : using UWalls = WallsScalar<WallType::UPPER>;
204 : PLUMED_REGISTER_ACTION(UWalls,"UPPER_WALLS_SCALAR")
205 :
206 : using LWalls = WallsScalar<WallType::LOWER>;
207 : PLUMED_REGISTER_ACTION(LWalls,"LOWER_WALLS_SCALAR")
208 :
209 : template<WallType W>
210 120 : void WallsScalar<W>::registerKeywords(Keywords& keys) {
211 120 : Bias::registerKeywords(keys);
212 : if constexpr(W==WallType::UPPER) {
213 202 : keys.setDisplayName("UPPER_WALLS");
214 : } else {
215 38 : keys.setDisplayName("LOWER_WALLS");
216 : }
217 :
218 120 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
219 120 : keys.add("compulsory","AT","the positions of the wall. The a_i in the expression for a wall.");
220 120 : keys.add("compulsory","KAPPA","the force constant for the wall. The k_i in the expression for a wall.");
221 120 : keys.add("compulsory","OFFSET","0.0","the offset for the start of the wall. The o_i in the expression for a wall.");
222 120 : keys.add("compulsory","EXP","2.0","the powers for the walls. The e_i in the expression for a wall.");
223 120 : keys.add("compulsory","EPS","1.0","the values for s_i in the expression for a wall");
224 240 : keys.addOutputComponent("force2","default","scalar","the instantaneous value of the squared force due to this bias potential");
225 120 : }
226 :
227 : template<WallType W>
228 57 : WallsScalar<W>::WallsScalar(const ActionOptions&ao):
229 : PLUMED_BIAS_INIT(ao),
230 114 : at(getNumberOfArguments(),0),
231 57 : kappa(getNumberOfArguments(),0.0),
232 57 : exp(getNumberOfArguments(),2.0),
233 57 : eps(getNumberOfArguments(),1.0),
234 114 : offset(getNumberOfArguments(),0.0) {
235 : // Note sizes of these vectors are automatically checked by parseVector :-)
236 57 : parseVector("OFFSET",offset);
237 57 : parseVector("EPS",eps);
238 57 : parseVector("EXP",exp);
239 57 : parseVector("KAPPA",kappa);
240 57 : parseVector("AT",at);
241 57 : checkRead();
242 :
243 57 : log.printf(" at");
244 120 : for(unsigned i=0; i<at.size(); i++) {
245 63 : log.printf(" %f",at[i]);
246 : }
247 57 : log.printf("\n");
248 57 : log.printf(" with an offset");
249 120 : for(unsigned i=0; i<offset.size(); i++) {
250 63 : log.printf(" %f",offset[i]);
251 : }
252 57 : log.printf("\n");
253 57 : log.printf(" with force constant");
254 120 : for(unsigned i=0; i<kappa.size(); i++) {
255 63 : log.printf(" %f",kappa[i]);
256 : }
257 57 : log.printf("\n");
258 57 : log.printf(" and exponent");
259 120 : for(unsigned i=0; i<exp.size(); i++) {
260 63 : log.printf(" %f",exp[i]);
261 : }
262 57 : log.printf("\n");
263 57 : log.printf(" rescaled");
264 120 : for(unsigned i=0; i<eps.size(); i++) {
265 63 : log.printf(" %f",eps[i]);
266 : }
267 57 : log.printf("\n");
268 :
269 : // The wall acts on the value of the argument in its domain, without periodic images
270 : // (see https://github.com/plumed/plumed2/issues/658).
271 120 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
272 : Value* arg=getPntrToArgument(i);
273 63 : if( arg->isPeriodic() ) {
274 : double min, max;
275 8 : arg->getDomain( min, max );
276 8 : if( at[i]<min || at[i]>max ) {
277 0 : error("the position of the wall (AT) for periodic argument " + arg->getName() + " should be within its periodic domain");
278 : }
279 8 : log.printf(" WARNING: argument %s is periodic. The wall acts on its value in the domain [%f,%f]\n"
280 : " without periodic images, so the bias is discontinuous at the boundary of the domain.\n"
281 : " Make sure that the walls keep the argument far from it (e.g. use both an upper and a lower wall).\n",
282 : arg->getName().c_str(), min, max );
283 : }
284 : }
285 :
286 114 : addComponent(force2);
287 57 : componentIsNotPeriodic(force2);
288 57 : }
289 :
290 :
291 :
292 : template<WallType W>
293 14492 : void WallsScalar<W>::calculate() {
294 : double ene = 0.0;
295 : double totf2 = 0.0;
296 29232 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
297 : double f = 0.0;
298 : // no periodic boundary conditions here, see the documentation and issue #658
299 14740 : const double cv=getArgument(i)-at[i];
300 14740 : const double off=offset[i];
301 14740 : const double epsilon=eps[i];
302 14740 : const double scale = walloff(cv,off)/epsilon;
303 14740 : if( check(scale) ) {
304 254 : const double k=kappa[i];
305 254 : const double exponent=exp[i];
306 : double power = wallpow( scale, exponent );
307 254 : f = -( k / epsilon ) * exponent * power / scale;
308 254 : ene += k * power;
309 254 : totf2 += f * f;
310 : }
311 14740 : setOutputForce(i,f);
312 : }
313 14492 : setBias(ene);
314 14492 : getPntrToComponent(force2)->set(totf2);
315 14492 : }
316 :
317 : }
318 : }
|
