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 "tools/Random.h"
24 : #include "ActionRegister.h"
25 : #include <ctime>
26 :
27 : using namespace std;
28 :
29 : namespace PLMD {
30 : namespace bias {
31 :
32 : //+PLUMEDOC BIAS ABMD
33 : /*
34 : Adds a ratchet-and-pawl like restraint on one or more variables.
35 :
36 : This action can be used to evolve a system towards a target value in
37 : CV space using an harmonic potential moving with the thermal fluctuations of the CV
38 : \cite ballone \cite provasi10abmd \cite camilloni11abmd. The biasing potential in this
39 : method is as follows:
40 :
41 : \f$
42 : V(\rho(t)) = \left \{ \begin{array}{ll} \frac{K}{2}\left(\rho(t)-\rho_m(t)\right)^2, &\rho(t)>\rho_m(t)\\
43 : 0, & \rho(t)\le\rho_m(t), \end{array} \right .
44 : \f$
45 :
46 :
47 : where
48 :
49 :
50 : \f$
51 : \rho(t)=\left(CV(t)-TO\right)^2
52 : \f$
53 :
54 :
55 : and
56 :
57 :
58 : \f$
59 : \rho_m(t)=\min_{0\le\tau\le t}\rho(\tau)+\eta(t)
60 : \f$.
61 :
62 : The method is based on the introduction of a biasing potential which is zero when
63 : the system is moving towards the desired arrival point and which damps the
64 : fluctuations when the system attempts to move in the opposite direction. As in the
65 : case of the ratchet and pawl system, propelled by thermal motion of the solvent
66 : molecules, the biasing potential does not exert work on the system. \f$\eta(t)\f$ is
67 : an additional white noise acting on the minimum position of the bias.
68 :
69 : \par Examples
70 :
71 : The following input sets up two biases, one on the distance between atoms 3 and 5
72 : and another on the distance between atoms 2 and 4. The two target values are defined
73 : using TO and the two strength using KAPPA. The total energy of the bias is printed.
74 : \plumedfile
75 : DISTANCE ATOMS=3,5 LABEL=d1
76 : DISTANCE ATOMS=2,4 LABEL=d2
77 : ABMD ARG=d1,d2 TO=1.0,1.5 KAPPA=5.0,5.0 LABEL=abmd
78 : PRINT ARG=abmd.bias,abmd.d1_min,abmd.d2_min
79 : \endplumedfile
80 :
81 : */
82 : //+ENDPLUMEDOC
83 :
84 3 : class ABMD : public Bias {
85 : std::vector<double> to;
86 : std::vector<double> min;
87 : std::vector<double> kappa;
88 : std::vector<double> temp;
89 : std::vector<int> seed;
90 : vector<Random> random;
91 : public:
92 : explicit ABMD(const ActionOptions&);
93 : void calculate();
94 : static void registerKeywords(Keywords& keys);
95 : };
96 :
97 6453 : PLUMED_REGISTER_ACTION(ABMD,"ABMD")
98 :
99 2 : void ABMD::registerKeywords(Keywords& keys) {
100 2 : Bias::registerKeywords(keys);
101 4 : keys.use("ARG");
102 8 : keys.add("compulsory","TO","The array of target values");
103 8 : keys.add("compulsory","KAPPA","The array of force constants.");
104 8 : keys.add("optional","MIN","Array of starting values for the bias (set rho_m(t), otherwise it is set using the current value of ARG)");
105 8 : keys.add("optional","NOISE","Array of white noise intensities (add a temperature to the ABMD)");
106 8 : keys.add("optional","SEED","Array of seeds for the white noise (add a temperature to the ABMD)");
107 8 : keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
108 8 : keys.addOutputComponent("_min","default","one or multiple instances of this quantity will be refereceable elsewhere in the input file. "
109 : " These quantities will be named with the arguments of the bias followed by "
110 : "the character string _min. These quantities tell the user the minimum value assumed by rho_m(t).");
111 2 : }
112 :
113 1 : ABMD::ABMD(const ActionOptions&ao):
114 : PLUMED_BIAS_INIT(ao),
115 : to(getNumberOfArguments(),0.0),
116 : min(getNumberOfArguments(),-1.0),
117 : kappa(getNumberOfArguments(),0.0),
118 : temp(getNumberOfArguments(),0.0),
119 1 : seed(getNumberOfArguments(),time(0)),
120 8 : random(getNumberOfArguments())
121 : {
122 : // Note : parseVector will check that number of arguments is correct
123 2 : parseVector("KAPPA",kappa);
124 2 : parseVector("MIN",min);
125 3 : if(min.size()==0) min.assign(getNumberOfArguments(),-1.0);
126 1 : if(min.size()!=getNumberOfArguments()) error("MIN array should have the same size as ARG array");
127 2 : parseVector("NOISE",temp);
128 2 : parseVector("SEED",seed);
129 2 : parseVector("TO",to);
130 1 : checkRead();
131 :
132 1 : log.printf(" min");
133 5 : for(unsigned i=0; i<min.size(); i++) log.printf(" %f",min[i]);
134 1 : log.printf("\n");
135 1 : log.printf(" to");
136 5 : for(unsigned i=0; i<to.size(); i++) log.printf(" %f",to[i]);
137 1 : log.printf("\n");
138 1 : log.printf(" with force constant");
139 5 : for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
140 1 : log.printf("\n");
141 :
142 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
143 1 : std::string str_min=getPntrToArgument(i)->getName()+"_min";
144 1 : addComponent(str_min); componentIsNotPeriodic(str_min);
145 1 : if(min[i]!=-1.0) getPntrToComponent(str_min)->set(min[i]);
146 : }
147 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {random[i].setSeed(-seed[i]);}
148 3 : addComponent("force2"); componentIsNotPeriodic("force2");
149 1 : }
150 :
151 :
152 5 : void ABMD::calculate() {
153 : double ene=0.0;
154 : double totf2=0.0;
155 15 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
156 5 : const double cv=difference(i,to[i],getArgument(i));
157 5 : const double cv2=cv*cv;
158 5 : const double k=kappa[i];
159 : double noise=0.;
160 5 : double diff=temp[i];
161 5 : if(diff>0) {
162 5 : noise = 2.*random[i].Gaussian()*diff;
163 5 : if(cv2<=diff) { diff=0.; temp[i]=0.; }
164 : }
165 :
166 : // min < 0 means that the variable has not been used in the input file, so the current position of the CV is used
167 : // cv2 < min means that the collective variable is nearer to the target value than at any other previous time so
168 : // min is set to the CV value
169 5 : if(min[i]<0.||cv2<min[i]) {
170 1 : min[i] = cv2;
171 : } else {
172 : // otherwise a noise is added to the minimum value
173 4 : min[i] += noise;
174 8 : const double f = -2.*k*(cv2-min[i])*cv;
175 : setOutputForce(i,f);
176 8 : ene += 0.5*k*(cv2-min[i])*(cv2-min[i]);
177 4 : totf2+=f*f;
178 : }
179 10 : getPntrToComponent(i+1)->set(min[i]);
180 : }
181 : setBias(ene);
182 10 : getPntrToComponent("force2")->set(totf2);
183 5 : }
184 :
185 : }
186 4839 : }
187 :
188 :
|
