1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-2021 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 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() override;
94 static void registerKeywords(Keywords& keys);
95 };
96
97 PLUMED_REGISTER_ACTION(ABMD,"ABMD")
98
registerKeywords(Keywords & keys)99 void ABMD::registerKeywords(Keywords& keys) {
100 Bias::registerKeywords(keys);
101 keys.use("ARG");
102 keys.add("compulsory","TO","The array of target values");
103 keys.add("compulsory","KAPPA","The array of force constants.");
104 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 keys.add("optional","NOISE","Array of white noise intensities (add a temperature to the ABMD)");
106 keys.add("optional","SEED","Array of seeds for the white noise (add a temperature to the ABMD)");
107 keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
108 keys.addOutputComponent("_min","default","one or multiple instances of this quantity can be referenced 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 }
112
ABMD(const ActionOptions & ao)113 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 seed(getNumberOfArguments(),time(0)),
120 random(getNumberOfArguments())
121 {
122 // Note : parseVector will check that number of arguments is correct
123 parseVector("KAPPA",kappa);
124 parseVector("MIN",min);
125 if(min.size()==0) min.assign(getNumberOfArguments(),-1.0);
126 if(min.size()!=getNumberOfArguments()) error("MIN array should have the same size as ARG array");
127 parseVector("NOISE",temp);
128 parseVector("SEED",seed);
129 parseVector("TO",to);
130 checkRead();
131
132 log.printf(" min");
133 for(unsigned i=0; i<min.size(); i++) log.printf(" %f",min[i]);
134 log.printf("\n");
135 log.printf(" to");
136 for(unsigned i=0; i<to.size(); i++) log.printf(" %f",to[i]);
137 log.printf("\n");
138 log.printf(" with force constant");
139 for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
140 log.printf("\n");
141
142 for(unsigned i=0; i<getNumberOfArguments(); i++) {
143 std::string str_min=getPntrToArgument(i)->getName()+"_min";
144 addComponent(str_min); componentIsNotPeriodic(str_min);
145 if(min[i]!=-1.0) getPntrToComponent(str_min)->set(min[i]);
146 }
147 for(unsigned i=0; i<getNumberOfArguments(); i++) {random[i].setSeed(-seed[i]);}
148 addComponent("force2"); componentIsNotPeriodic("force2");
149 }
150
151
calculate()152 void ABMD::calculate() {
153 double ene=0.0;
154 double totf2=0.0;
155 for(unsigned i=0; i<getNumberOfArguments(); ++i) {
156 const double cv=difference(i,to[i],getArgument(i));
157 const double cv2=cv*cv;
158 const double k=kappa[i];
159 double noise=0.;
160 double diff=temp[i];
161 if(diff>0) {
162 noise = 2.*random[i].Gaussian()*diff;
163 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 if(min[i]<0.||cv2<min[i]) {
170 min[i] = cv2;
171 } else {
172 // otherwise a noise is added to the minimum value
173 min[i] += noise;
174 const double f = -2.*k*(cv2-min[i])*cv;
175 setOutputForce(i,f);
176 ene += 0.5*k*(cv2-min[i])*(cv2-min[i]);
177 totf2+=f*f;
178 }
179 getPntrToComponent(i+1)->set(min[i]);
180 }
181 setBias(ene);
182 getPntrToComponent("force2")->set(totf2);
183 }
184
185 }
186 }
187
188
189