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