Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 "ActionRegister.h"
24 :
25 :
26 : using namespace std;
27 :
28 :
29 : namespace PLMD {
30 : namespace bias {
31 :
32 : //+PLUMEDOC BIAS MOVINGRESTRAINT
33 : /*
34 : Add a time-dependent, harmonic restraint on one or more variables.
35 :
36 : This form of bias can be used to performed steered MD \cite Grubmuller3
37 : and Jarzynski sampling \cite jarzynski.
38 :
39 : The harmonic restraint on your system is given by:
40 :
41 : \f[
42 : V(\vec{s},t) = \frac{1}{2} \kappa(t) ( \vec{s} - \vec{s}_0(t) )^2
43 : \f]
44 :
45 : The time dependence of \f$\kappa\f$ and \f$\vec{s}_0\f$ are specified by a list of
46 : STEP, KAPPA and AT keywords. These keywords tell plumed what values \f$\kappa\f$ and \f$\vec{s}_0\f$
47 : should have at the time specified by the corresponding STEP keyword. Inbetween these times
48 : the values of \f$\kappa\f$ and \f$\vec{s}_0\f$ are linearly interpolated.
49 :
50 : Additional material and examples can be also found in the tutorial \ref belfast-5
51 :
52 : \par Examples
53 :
54 : The following input is dragging the distance between atoms 2 and 4
55 : from 1 to 2 in the first 1000 steps, then back in the next 1000 steps.
56 : In the following 500 steps the restraint is progressively switched off.
57 : \plumedfile
58 : DISTANCE ATOMS=2,4 LABEL=d
59 : MOVINGRESTRAINT ...
60 : ARG=d
61 : STEP0=0 AT0=1.0 KAPPA0=100.0
62 : STEP1=1000 AT1=2.0
63 : STEP2=2000 AT2=1.0
64 : STEP3=2500 KAPPA3=0.0
65 : ... MOVINGRESTRAINT
66 : \endplumedfile
67 : The following input is progressively building restraints
68 : distances between atoms 1 and 5 and between atoms 2 and 4
69 : in the first 1000 steps. Afterwards, the restraint is kept
70 : static.
71 : \plumedfile
72 : DISTANCE ATOMS=1,5 LABEL=d1
73 : DISTANCE ATOMS=2,4 LABEL=d2
74 : MOVINGRESTRAINT ...
75 : ARG=d1,d2
76 : STEP0=0 AT0=1.0,1.5 KAPPA0=0.0,0.0
77 : STEP1=1000 AT1=1.0,1.5 KAPPA1=1.0,1.0
78 : ... MOVINGRESTRAINT
79 : \endplumedfile
80 : The following input is progressively bringing atoms 1 and 2
81 : close to each other with an upper wall
82 : \plumedfile
83 : DISTANCE ATOMS=1,2 LABEL=d1
84 : MOVINGRESTRAINT ...
85 : ARG=d1
86 : VERSE=U
87 : STEP0=0 AT0=1.0 KAPPA0=10.0
88 : STEP1=1000 AT1=0.0
89 : ... MOVINGRESTRAINT
90 : \endplumedfile
91 :
92 : By default the Action is issuing some values which are
93 : the work on each degree of freedom, the center of the harmonic potential,
94 : the total bias deposited
95 :
96 : (See also \ref DISTANCE).
97 :
98 : \attention Work is not computed properly when KAPPA is time dependent.
99 :
100 : */
101 : //+ENDPLUMEDOC
102 :
103 :
104 12 : class MovingRestraint : public Bias {
105 : std::vector<std::vector<double> > at;
106 : std::vector<std::vector<double> > kappa;
107 : std::vector<long int> step;
108 : std::vector<double> oldaa;
109 : std::vector<double> oldk;
110 : std::vector<double> olddpotdk;
111 : std::vector<double> oldf;
112 : std::vector<string> verse;
113 : std::vector<double> work;
114 : double tot_work;
115 : public:
116 : explicit MovingRestraint(const ActionOptions&);
117 : void calculate();
118 : static void registerKeywords( Keywords& keys );
119 : };
120 :
121 6455 : PLUMED_REGISTER_ACTION(MovingRestraint,"MOVINGRESTRAINT")
122 :
123 4 : void MovingRestraint::registerKeywords( Keywords& keys ) {
124 4 : Bias::registerKeywords(keys);
125 8 : keys.use("ARG");
126 20 : keys.add("compulsory","VERSE","B","Tells plumed whether the restraint is only acting for CV larger (U) or smaller (L) than "
127 : "the restraint or whether it is acting on both sides (B)");
128 16 : keys.add("numbered","STEP","This keyword appears multiple times as STEPx with x=0,1,2,...,n. Each value given represents "
129 : "the MD step at which the restraint parameters take the values KAPPAx and ATx.");
130 12 : keys.reset_style("STEP","compulsory");
131 16 : keys.add("numbered","AT","ATx is equal to the position of the restraint at time STEPx. For intermediate times this parameter "
132 : "is linearly interpolated. If no ATx is specified for STEPx then the values of AT are kept constant "
133 : "during the interval of time between STEPx-1 and STEPx.");
134 12 : keys.reset_style("AT","compulsory");
135 16 : keys.add("numbered","KAPPA","KAPPAx is equal to the value of the force constants at time STEPx. For intermediate times this "
136 : "parameter is linearly interpolated. If no KAPPAx is specified for STEPx then the values of KAPPAx "
137 : "are kept constant during the interval of time between STEPx-1 and STEPx.");
138 12 : keys.reset_style("KAPPA","compulsory");
139 16 : keys.addOutputComponent("work","default","the total work performed changing this restraint");
140 16 : keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
141 16 : keys.addOutputComponent("_cntr","default","one or multiple instances of this quantity will be refereceable elsewhere in the input file. "
142 : "these quantities will named with the arguments of the bias followed by "
143 : "the character string _cntr. These quantities give the instantaneous position "
144 : "of the center of the harmonic potential.");
145 16 : keys.addOutputComponent("_work","default","one or multiple instances of this quantity will be refereceable elsewhere in the input file. "
146 : "These quantities will named with the arguments of the bias followed by "
147 : "the character string _work. These quantities tell the user how much work has "
148 : "been done by the potential in dragging the system along the various colvar axis.");
149 16 : keys.addOutputComponent("_kappa","default","one or multiple instances of this quantity will be refereceable elsewhere in the input file. "
150 : "These quantities will named with the arguments of the bias followed by "
151 : "the character string _kappa. These quantities tell the user the time dependent value of kappa.");
152 4 : }
153 :
154 3 : MovingRestraint::MovingRestraint(const ActionOptions&ao):
155 : PLUMED_BIAS_INIT(ao),
156 15 : verse(getNumberOfArguments())
157 : {
158 6 : parseVector("VERSE",verse);
159 6 : vector<long int> ss(1); ss[0]=-1;
160 6 : std::vector<double> kk( getNumberOfArguments() ), aa( getNumberOfArguments() );
161 8 : for(int i=0;; i++) {
162 : // Read in step
163 22 : if( !parseNumberedVector("STEP",i,ss) ) break;
164 40 : for(unsigned j=0; j<step.size(); j++) {
165 16 : if(ss[0]<step[j]) error("in moving restraint step number must always increase");
166 : }
167 8 : step.push_back(ss[0]);
168 :
169 : // Try to read kappa
170 18 : if( !parseNumberedVector("KAPPA",i,kk) ) kk=kappa[i-1];
171 8 : kappa.push_back(kk);
172 :
173 : // Now read AT
174 17 : if( !parseNumberedVector("AT",i,aa) ) aa=at[i-1];
175 8 : at.push_back(aa);
176 8 : }
177 3 : checkRead();
178 :
179 30 : for(unsigned i=0; i<step.size(); i++) {
180 16 : log.printf(" step%u %ld\n",i,step[i]);
181 8 : log.printf(" at");
182 46 : for(unsigned j=0; j<at[i].size(); j++) log.printf(" %f",at[i][j]);
183 8 : log.printf("\n");
184 8 : log.printf(" with force constant");
185 46 : for(unsigned j=0; j<kappa[i].size(); j++) log.printf(" %f",kappa[i][j]);
186 8 : log.printf("\n");
187 : };
188 :
189 9 : addComponent("force2"); componentIsNotPeriodic("force2");
190 :
191 : // add the centers of the restraint as additional components that can be retrieved (useful for debug)
192 :
193 : std::string comp;
194 11 : for(unsigned i=0; i< getNumberOfArguments() ; i++) {
195 8 : comp=getPntrToArgument(i)->getName()+"_cntr"; // each spring has its own center
196 4 : addComponent(comp); componentIsNotPeriodic(comp);
197 8 : comp=getPntrToArgument(i)->getName()+"_work"; // each spring has its own work
198 4 : addComponent(comp); componentIsNotPeriodic(comp);
199 8 : comp=getPntrToArgument(i)->getName()+"_kappa"; // each spring has its own kappa
200 4 : addComponent(comp); componentIsNotPeriodic(comp);
201 8 : work.push_back(0.); // initialize the work value
202 : }
203 9 : addComponent("work"); componentIsNotPeriodic("work");
204 3 : tot_work=0.0;
205 :
206 3 : log<<" Bibliography ";
207 9 : log<<cite("Grubmuller, Heymann, and Tavan, Science 271, 997 (1996)")<<"\n";
208 :
209 3 : }
210 :
211 :
212 561 : void MovingRestraint::calculate() {
213 : double ene=0.0;
214 : double totf2=0.0;
215 561 : unsigned narg=getNumberOfArguments();
216 561 : long int now=getStep();
217 561 : std::vector<double> kk(narg),aa(narg),f(narg),dpotdk(narg);
218 561 : if(now<=step[0]) {
219 2 : kk=kappa[0];
220 2 : aa=at[0];
221 2 : oldaa=at[0];
222 2 : oldk=kappa[0];
223 2 : olddpotdk.resize(narg);
224 2 : oldf.resize(narg);
225 1118 : } else if(now>=step[step.size()-1]) {
226 47 : kk=kappa[step.size()-1];
227 94 : aa=at[step.size()-1];
228 : } else {
229 : unsigned i=0;
230 1034 : for(i=1; i<step.size(); i++) if(now<step[i]) break;
231 1536 : double c2=(now-step[i-1])/double(step[i]-step[i-1]);
232 512 : double c1=1.0-c2;
233 2592 : for(unsigned j=0; j<narg; j++) kk[j]=(c1*kappa[i-1][j]+c2*kappa[i][j]);
234 1552 : for(unsigned j=0; j<narg; j++) {
235 1040 : const double a1=at[i-1][j];
236 520 : const double a2=at[i][j];
237 1560 : aa[j]=(c1*a1+c2*(a1+difference(j,a1,a2)));
238 : }
239 : }
240 561 : tot_work=0.0;
241 1703 : for(unsigned i=0; i<narg; ++i) {
242 571 : const double cv=difference(i,aa[i],getArgument(i)); // this gives: getArgument(i) - aa[i]
243 2284 : getPntrToComponent(getPntrToArgument(i)->getName()+"_cntr")->set(aa[i]);
244 571 : const double k=kk[i];
245 1142 : f[i]=-k*cv;
246 571 : if(verse[i]=="U" && cv<0) continue;
247 571 : if(verse[i]=="L" && cv>0) continue;
248 1713 : plumed_assert(verse[i]=="U" || verse[i]=="L" || verse[i]=="B");
249 1142 : dpotdk[i]=0.5*cv*cv;
250 6271 : if(oldaa.size()==aa.size() && oldf.size()==f.size()) work[i]+=0.5*(oldf[i]+f[i])*(aa[i]-oldaa[i]) + 0.5*( dpotdk[i]+olddpotdk[i] )*(kk[i]-oldk[i]);
251 1713 : getPntrToComponent(getPntrToArgument(i)->getName()+"_work")->set(work[i]);
252 1713 : getPntrToComponent(getPntrToArgument(i)->getName()+"_kappa")->set(kk[i]);
253 571 : tot_work+=work[i];
254 571 : ene+=0.5*k*cv*cv;
255 571 : setOutputForce(i,f[i]);
256 571 : totf2+=f[i]*f[i];
257 : };
258 1122 : getPntrToComponent("work")->set(tot_work);
259 561 : oldf=f;
260 561 : oldaa=aa;
261 561 : oldk=kk;
262 561 : olddpotdk=dpotdk;
263 : setBias(ene);
264 1122 : getPntrToComponent("force2")->set(totf2);
265 561 : }
266 :
267 : }
268 4839 : }
269 :
270 :
|
