Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 MOVINGRESTRAINT
29 : /*
30 : Add a time-dependent, harmonic restraint on one or more variables.
31 :
32 : This form of bias can be used to perform the steered MD
33 : and Jarzynski sampling calculations that are discussed in the papers cited below.
34 :
35 : The harmonic restraint on your system is given by:
36 :
37 : $$
38 : V(\vec{s},t) = \frac{1}{2} \kappa(t) ( \vec{s} - \vec{s}_0(t) )^2
39 : $$
40 :
41 : The time dependence of $\kappa$ and $\vec{s}_0$ are specified by a list of
42 : STEP, KAPPA and AT keywords. These keywords tell plumed what values $\kappa$ and $\vec{s}_0$
43 : should have at the time specified by the corresponding STEP keyword. In between these times
44 : the values of $\kappa$ and $\vec{s}_0$ are linearly interpolated.
45 :
46 : ## Examples
47 :
48 : The following input is dragging the distance between atoms 2 and 4
49 : from 1 to 2 in the first 1000 steps, then back in the next 1000 steps.
50 : In the following 500 steps the restraint is progressively switched off.
51 :
52 : ```plumed
53 : DISTANCE ATOMS=2,4 LABEL=d
54 : MOVINGRESTRAINT ...
55 : ARG=d
56 : STEP0=0 AT0=1.0 KAPPA0=100.0
57 : STEP1=1000 AT1=2.0
58 : STEP2=2000 AT2=1.0
59 : STEP3=2500 KAPPA3=0.0
60 : ... MOVINGRESTRAINT
61 : ```
62 :
63 : The following input is progressively building restraints
64 : distances between atoms 1 and 5 and between atoms 2 and 4
65 : in the first 1000 steps. Afterwards, the restraint is kept
66 : static.
67 :
68 : ```plumed
69 : DISTANCE ATOMS=1,5 LABEL=d1
70 : DISTANCE ATOMS=2,4 LABEL=d2
71 : MOVINGRESTRAINT ...
72 : ARG=d1,d2
73 : STEP0=0 AT0=1.0,1.5 KAPPA0=0.0,0.0
74 : STEP1=1000 AT1=1.0,1.5 KAPPA1=1.0,1.0
75 : ... MOVINGRESTRAINT
76 : ```
77 :
78 : The following input is progressively bringing atoms 1 and 2
79 : close to each other with an upper wall
80 :
81 : ```plumed
82 : DISTANCE ATOMS=1,2 LABEL=d1
83 : MOVINGRESTRAINT ...
84 : ARG=d1
85 : VERSE=U
86 : STEP0=0 AT0=1.0 KAPPA0=10.0
87 : STEP1=1000 AT1=0.0
88 : ... MOVINGRESTRAINT
89 : ```
90 :
91 : By default the Action is issuing some values which are
92 : the work on each degree of freedom, the center of the harmonic potential,
93 : the total bias deposited
94 :
95 : (See also [DISTANCE](DISTANCE.md)).
96 :
97 : !!! attention ""
98 :
99 : Work is not computed properly when KAPPA is time dependent.
100 :
101 : */
102 : //+ENDPLUMEDOC
103 :
104 :
105 : class MovingRestraint : public Bias {
106 : std::vector<std::vector<double> > at;
107 : std::vector<std::vector<double> > kappa;
108 : std::vector<long long int> step;
109 : std::vector<double> oldaa;
110 : std::vector<double> oldk;
111 : std::vector<double> olddpotdk;
112 : std::vector<double> oldf;
113 : std::vector<std::string> verse;
114 : std::vector<double> work;
115 : //memeber used as temporary working objects
116 : std::vector<double> kk,aa,f,dpotdk;
117 : std::vector<Value*> valueCntr;
118 : std::vector<Value*> valueWork;
119 : std::vector<Value*> valueKappa;
120 : Value* valueTotWork=nullptr;
121 : Value* valueForce2=nullptr;
122 : public:
123 : explicit MovingRestraint(const ActionOptions&);
124 : void calculate() override;
125 : static void registerKeywords( Keywords& keys );
126 : };
127 :
128 : PLUMED_REGISTER_ACTION(MovingRestraint,"MOVINGRESTRAINT")
129 :
130 12 : void MovingRestraint::registerKeywords( Keywords& keys ) {
131 12 : Bias::registerKeywords(keys);
132 12 : keys.add("compulsory","VERSE","B","Tells plumed whether the restraint is only acting for CV larger (U) or smaller (L) than "
133 : "the restraint or whether it is acting on both sides (B)");
134 12 : keys.add("numbered","STEP","This keyword appears multiple times as STEPx with x=0,1,2,...,n. Each value given represents "
135 : "the MD step at which the restraint parameters take the values KAPPAx and ATx.");
136 24 : keys.reset_style("STEP","compulsory");
137 12 : keys.add("numbered","AT","ATx is equal to the position of the restraint at time STEPx. For intermediate times this parameter "
138 : "is linearly interpolated. If no ATx is specified for STEPx then the values of AT are kept constant "
139 : "during the interval of time between STEP(x-1) and STEPx.");
140 24 : keys.reset_style("AT","compulsory");
141 12 : keys.add("numbered","KAPPA","KAPPAx is equal to the value of the force constants at time STEPx. For intermediate times this "
142 : "parameter is linearly interpolated. If no KAPPAx is specified for STEPx then the values of KAPPAx "
143 : "are kept constant during the interval of time between STEP(x-1) and STEPx.");
144 24 : keys.reset_style("KAPPA","compulsory");
145 24 : keys.addOutputComponent("work","default","scalar","the total work performed changing this restraint");
146 24 : keys.addOutputComponent("force2","default","scalar","the instantaneous value of the squared force due to this bias potential");
147 24 : keys.addOutputComponent("_cntr","default","scalar","one or multiple instances of this quantity can be referenced elsewhere in the input file. "
148 : "these quantities will named with the arguments of the bias followed by "
149 : "the character string _cntr. These quantities give the instantaneous position "
150 : "of the center of the harmonic potential.");
151 24 : keys.addOutputComponent("_work","default","scalar","one or multiple instances of this quantity can be referenced elsewhere in the input file. "
152 : "These quantities will named with the arguments of the bias followed by "
153 : "the character string _work. These quantities tell the user how much work has "
154 : "been done by the potential in dragging the system along the various colvar axis.");
155 24 : keys.addOutputComponent("_kappa","default","scalar","one or multiple instances of this quantity can be referenced elsewhere in the input file. "
156 : "These quantities will named with the arguments of the bias followed by "
157 : "the character string _kappa. These quantities tell the user the time dependent value of kappa.");
158 12 : keys.addDOI("10.1126/science.271.5251.997");
159 12 : keys.addDOI("10.1103/PhysRevLett.78.2690");
160 12 : }
161 :
162 8 : MovingRestraint::MovingRestraint(const ActionOptions&ao):
163 : PLUMED_BIAS_INIT(ao),
164 10 : verse(getNumberOfArguments()),
165 10 : kk(getNumberOfArguments()),
166 10 : aa(getNumberOfArguments()),
167 10 : f(getNumberOfArguments()),
168 18 : dpotdk(getNumberOfArguments()) {
169 16 : parseVector("VERSE",verse);
170 16 : for(const auto & v : verse) {
171 9 : if(v!="U" && v!="L" && v!="B") {
172 2 : error("VERSE should be U, L or B, found " + v);
173 : }
174 : }
175 9 : std::vector<long long int> ss(1);
176 7 : ss[0]=-1;
177 14 : for(int i=0;; i++) {
178 : // Read in step
179 42 : if( !parseNumberedVector("STEP",i,ss) ) {
180 : break;
181 : }
182 26 : for(unsigned j=0; j<step.size(); j++) {
183 11 : if(ss[0]<step[j]) {
184 0 : error("in moving restraint step number must always increase");
185 : }
186 : }
187 15 : step.push_back(ss[0]);
188 :
189 : // Try to read kappa
190 30 : if( !parseNumberedVector("KAPPA",i,kk) ) {
191 4 : if(i==0) {
192 1 : error("KAPPA0 is required");
193 : }
194 3 : kk=kappa[i-1];
195 : }
196 14 : kappa.push_back(kk);
197 :
198 : // Now read AT
199 28 : if( !parseNumberedVector("AT",i,aa) ) {
200 2 : if(i==0) {
201 0 : error("AT0 is required");
202 : }
203 2 : aa=at[i-1];
204 : }
205 14 : at.push_back(aa);
206 14 : }
207 6 : if(step.empty()) {
208 0 : error("STEP0 is required");
209 : }
210 6 : checkRead();
211 :
212 20 : for(unsigned i=0; i<step.size(); i++) {
213 14 : log.printf(" step%u %lld\n",i,step[i]);
214 14 : log.printf(" at");
215 30 : for(unsigned j=0; j<at[i].size(); j++) {
216 16 : log.printf(" %f",at[i][j]);
217 : }
218 14 : log.printf("\n");
219 14 : log.printf(" with force constant");
220 30 : for(unsigned j=0; j<kappa[i].size(); j++) {
221 16 : log.printf(" %f",kappa[i][j]);
222 : }
223 14 : log.printf("\n");
224 : };
225 :
226 12 : addComponent("force2");
227 6 : componentIsNotPeriodic("force2");
228 6 : valueForce2=getPntrToComponent("force2");
229 :
230 : // add the centers of the restraint as additional components that can be retrieved (useful for debug)
231 :
232 : std::string comp;
233 13 : for(unsigned i=0; i< getNumberOfArguments() ; i++) {
234 14 : comp=getPntrToArgument(i)->getName()+"_cntr"; // each spring has its own center
235 7 : addComponent(comp);
236 7 : componentIsNotPeriodic(comp);
237 7 : valueCntr.push_back(getPntrToComponent(comp));
238 14 : comp=getPntrToArgument(i)->getName()+"_work"; // each spring has its own work
239 7 : addComponent(comp);
240 7 : componentIsNotPeriodic(comp);
241 7 : valueWork.push_back(getPntrToComponent(comp));
242 14 : comp=getPntrToArgument(i)->getName()+"_kappa"; // each spring has its own kappa
243 7 : addComponent(comp);
244 7 : componentIsNotPeriodic(comp);
245 7 : valueKappa.push_back(getPntrToComponent(comp));
246 7 : work.push_back(0.); // initialize the work value
247 : }
248 12 : addComponent("work");
249 6 : componentIsNotPeriodic("work");
250 6 : valueTotWork=getPntrToComponent("work");
251 :
252 6 : log<<" Bibliography ";
253 12 : log<<cite("Grubmuller, Heymann, and Tavan, Science 271, 997 (1996)")<<"\n";
254 :
255 19 : }
256 :
257 :
258 686 : void MovingRestraint::calculate() {
259 : double ene=0.0;
260 : double totf2=0.0;
261 686 : unsigned narg=getNumberOfArguments();
262 686 : long long int now=getStep();
263 686 : if(now<=step[0]) {
264 5 : kk=kappa[0];
265 5 : aa=at[0];
266 5 : oldaa=at[0];
267 5 : oldk=kappa[0];
268 5 : olddpotdk.resize(narg);
269 5 : oldf.resize(narg);
270 681 : } else if(now>=step[step.size()-1]) {
271 47 : kk=kappa[step.size()-1];
272 47 : aa=at[step.size()-1];
273 : } else {
274 : unsigned i=0;
275 639 : for(i=1; i<step.size()-1; i++)
276 10 : if(now<step[i]) {
277 : break;
278 : }
279 634 : double c2=(now-step[i-1])/double(step[i]-step[i-1]);
280 634 : double c1=1.0-c2;
281 1276 : for(unsigned j=0; j<narg; j++) {
282 642 : kk[j]=(c1*kappa[i-1][j]+c2*kappa[i][j]);
283 : }
284 1276 : for(unsigned j=0; j<narg; j++) {
285 642 : const double a1=at[i-1][j];
286 642 : const double a2=at[i][j];
287 642 : aa[j]=(c1*a1+c2*(a1+difference(j,a1,a2)));
288 : }
289 : }
290 : double tot_work=0.0;
291 1382 : for(unsigned i=0; i<narg; ++i) {
292 696 : const double cv=difference(i,aa[i],getArgument(i)); // this gives: getArgument(i) - aa[i]
293 696 : valueCntr[i]->set(aa[i]);
294 696 : const double k=kk[i];
295 : // with VERSE=U (L) the restraint acts only when the CV is larger (smaller) than the center:
296 : // otherwise the potential, the force and the derivative with respect to kappa are zero
297 696 : const bool acting=!( (verse[i]=="U" && cv<0) || (verse[i]=="L" && cv>0) );
298 : if(acting) {
299 636 : f[i]=-k*cv;
300 636 : dpotdk[i]=0.5*cv*cv;
301 636 : ene+=0.5*k*cv*cv;
302 : } else {
303 60 : f[i]=0.0;
304 60 : dpotdk[i]=0.0;
305 : }
306 696 : if(oldaa.size()==aa.size() && oldf.size()==f.size()) {
307 695 : work[i]+=0.5*(oldf[i]+f[i])*(aa[i]-oldaa[i]) + 0.5*( dpotdk[i]+olddpotdk[i] )*(kk[i]-oldk[i]);
308 : }
309 696 : valueWork[i]->set(work[i]);
310 696 : valueKappa[i]->set(kk[i]);
311 696 : tot_work+=work[i];
312 696 : setOutputForce(i,f[i]);
313 696 : totf2+=f[i]*f[i];
314 : }
315 686 : valueTotWork->set(tot_work);
316 686 : oldf=f;
317 686 : oldaa=aa;
318 686 : oldk=kk;
319 686 : olddpotdk=dpotdk;
320 686 : setBias(ene);
321 686 : valueForce2->set(totf2);
322 686 : }
323 :
324 : }
325 : }
326 :
327 :
|
