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