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 LOWER_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 : The restraining potential starts acting on the system when the value of the CV is greater
34 : (in the case of UPPER_WALLS) or lower (in the case of LOWER_WALLS) than a certain limit \f$a_i\f$ (AT)
35 : minus an offset \f$o_i\f$ (OFFSET).
36 : The expression for the bias due to the wall is given by:
37 :
38 : for UPPER_WALLS:
39 : \f$
40 : \sum_i {k_i}((x_i-a_i+o_i)/s_i)^e_i
41 : \f$
42 :
43 : for LOWER_WALLS:
44 : \f$
45 : \sum_i {k_i}|(x_i-a_i-o_i)/s_i|^e_i
46 : \f$
47 :
48 : \f$k_i\f$ (KAPPA) is an energy constant in internal unit of the code, \f$s_i\f$ (EPS) a rescaling factor and
49 : \f$e_i\f$ (EXP) the exponent determining the power law. By default: EXP = 2, EPS = 1.0, OFFSET = 0.
50 :
51 :
52 : \par Examples
53 :
54 : The following input tells plumed to add both a lower and an upper walls on the distance
55 : between atoms 3 and 5 and the distance between atoms 2 and 4. The lower and upper limits
56 : are defined at different values. The strength of the walls is the same for the four cases.
57 : It also tells plumed to print the energy of the walls.
58 : \plumedfile
59 : DISTANCE ATOMS=3,5 LABEL=d1
60 : DISTANCE ATOMS=2,4 LABEL=d2
61 : UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=uwall
62 : LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=lwall
63 : PRINT ARG=uwall.bias,lwall.bias
64 : \endplumedfile
65 : (See also \ref DISTANCE and \ref PRINT).
66 :
67 :
68 : \par Periodic arguments
69 :
70 : If an argument is periodic (e.g. a \ref TORSION), the wall acts on the value of the argument in its periodic domain
71 : and periodic boundary conditions are __not__ applied to the difference \f$x_i-a_i\f$. The value of AT must thus be within the
72 : periodic domain of the argument. This is the case both when the arguments are scalars and when they are vectors.
73 : The following input restrains a torsion within the interval \f$[-2,2]\f$:
74 : \plumedfile
75 : t: TORSION ATOMS=1,2,3,4
76 : uwall: UPPER_WALLS ARG=t AT=2.0 KAPPA=150.0
77 : lwall: LOWER_WALLS ARG=t AT=-2.0 KAPPA=150.0
78 : PRINT ARG=t,uwall.bias,lwall.bias
79 : \endplumedfile
80 :
81 : Notice that the bias is discontinuous at the boundary of the periodic domain (here at \f$\pm\pi\f$). You should thus make sure that
82 : 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.
83 :
84 : Versions of PLUMED before 2.9.7 applied periodic boundary conditions to the difference \f$x_i-a_i\f$. With that choice, each wall acts on half of the
85 : periodic domain, and the region left free by an upper and a lower wall is smaller than the one between AT values
86 : (see \issue{658}).
87 :
88 : */
89 : //+ENDPLUMEDOC
90 :
91 : //+PLUMEDOC BIAS LOWER_WALLS_SCALAR
92 : /*
93 : Defines a wall for the value of one or more collective variables,
94 : which limits the region of the phase space accessible during the simulation.
95 :
96 : \par Examples
97 :
98 : */
99 : //+ENDPLUMEDOC
100 :
101 : class LWalls : public Bias {
102 : std::vector<double> at;
103 : std::vector<double> kappa;
104 : std::vector<double> exp;
105 : std::vector<double> eps;
106 : std::vector<double> offset;
107 : public:
108 : explicit LWalls(const ActionOptions&);
109 : void calculate() override;
110 : static void registerKeywords(Keywords& keys);
111 : };
112 :
113 : PLUMED_REGISTER_ACTION(LWalls,"LOWER_WALLS_SCALAR")
114 :
115 20 : void LWalls::registerKeywords(Keywords& keys) {
116 20 : Bias::registerKeywords(keys);
117 20 : keys.setDisplayName("LOWER_WALLS");
118 20 : keys.use("ARG");
119 40 : keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
120 40 : keys.add("compulsory","AT","the positions of the wall. The a_i in the expression for a wall.");
121 40 : keys.add("compulsory","KAPPA","the force constant for the wall. The k_i in the expression for a wall.");
122 40 : keys.add("compulsory","OFFSET","0.0","the offset for the start of the wall. The o_i in the expression for a wall.");
123 40 : keys.add("compulsory","EXP","2.0","the powers for the walls. The e_i in the expression for a wall.");
124 40 : keys.add("compulsory","EPS","1.0","the values for s_i in the expression for a wall");
125 40 : keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
126 20 : }
127 :
128 8 : LWalls::LWalls(const ActionOptions&ao):
129 : PLUMED_BIAS_INIT(ao),
130 16 : at(getNumberOfArguments(),0),
131 8 : kappa(getNumberOfArguments(),0.0),
132 8 : exp(getNumberOfArguments(),2.0),
133 8 : eps(getNumberOfArguments(),1.0),
134 16 : offset(getNumberOfArguments(),0.0) {
135 : // Note sizes of these vectors are automatically checked by parseVector :-)
136 8 : parseVector("OFFSET",offset);
137 8 : parseVector("EPS",eps);
138 8 : parseVector("EXP",exp);
139 8 : parseVector("KAPPA",kappa);
140 8 : parseVector("AT",at);
141 8 : checkRead();
142 :
143 8 : log.printf(" at");
144 18 : for(unsigned i=0; i<at.size(); i++) {
145 10 : log.printf(" %f",at[i]);
146 : }
147 8 : log.printf("\n");
148 8 : log.printf(" with an offset");
149 18 : for(unsigned i=0; i<offset.size(); i++) {
150 10 : log.printf(" %f",offset[i]);
151 : }
152 8 : log.printf("\n");
153 8 : log.printf(" with force constant");
154 18 : for(unsigned i=0; i<kappa.size(); i++) {
155 10 : log.printf(" %f",kappa[i]);
156 : }
157 8 : log.printf("\n");
158 8 : log.printf(" and exponent");
159 18 : for(unsigned i=0; i<exp.size(); i++) {
160 10 : log.printf(" %f",exp[i]);
161 : }
162 8 : log.printf("\n");
163 8 : log.printf(" rescaled");
164 18 : for(unsigned i=0; i<eps.size(); i++) {
165 10 : log.printf(" %f",eps[i]);
166 : }
167 8 : log.printf("\n");
168 :
169 : // The wall acts on the value of the argument in its domain, without periodic images
170 : // (see https://github.com/plumed/plumed2/issues/658).
171 18 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
172 : Value* arg=getPntrToArgument(i);
173 10 : if( arg->isPeriodic() ) {
174 : double min, max;
175 4 : arg->getDomain( min, max );
176 4 : if( at[i]<min || at[i]>max ) {
177 0 : error("the position of the wall (AT) for periodic argument " + arg->getName() + " should be within its periodic domain");
178 : }
179 4 : log.printf(" WARNING: argument %s is periodic. The wall acts on its value in the domain [%f,%f]\n"
180 : " without periodic images, so the bias is discontinuous at the boundary of the domain.\n"
181 : " Make sure that the walls keep the argument far from it (e.g. use both an upper and a lower wall).\n",
182 : arg->getName().c_str(), min, max );
183 : }
184 : }
185 :
186 16 : addComponent("force2");
187 8 : componentIsNotPeriodic("force2");
188 8 : }
189 :
190 117 : void LWalls::calculate() {
191 : double ene = 0.0;
192 : double totf2 = 0.0;
193 298 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
194 : double f = 0.0;
195 : // no periodic boundary conditions here, see the documentation and issue #658
196 181 : const double cv=getArgument(i)-at[i];
197 181 : const double off=offset[i];
198 181 : const double epsilon=eps[i];
199 181 : const double lscale = (cv-off)/epsilon;
200 181 : if( lscale < 0.) {
201 76 : const double k=kappa[i];
202 76 : const double exponent=exp[i];
203 76 : double power = std::pow( -lscale, exponent );
204 76 : f = -( k / epsilon ) * exponent * power / lscale;
205 76 : ene += k * power;
206 76 : totf2 += f * f;
207 : }
208 181 : setOutputForce(i,f);
209 : }
210 117 : setBias(ene);
211 117 : getPntrToComponent("force2")->set(totf2);
212 117 : }
213 :
214 : }
215 : }
|
