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.use("ARG");
100 6 : keys.add("compulsory","TO","The array of target values");
101 6 : keys.add("compulsory","KAPPA","The array of force constants.");
102 6 : 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)");
103 6 : keys.add("optional","NOISE","Array of white noise intensities (add a temperature to the ABMD)");
104 6 : keys.add("optional","SEED","Array of seeds for the white noise (add a temperature to the ABMD)");
105 6 : keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
106 6 : keys.addOutputComponent("_min","default","one or multiple instances of this quantity can be referenced elsewhere in the input file. "
107 : " These quantities will be named with the arguments of the bias followed by "
108 : "the character string _min. These quantities tell the user the minimum value assumed by rho_m(t).");
109 3 : }
110 :
111 1 : ABMD::ABMD(const ActionOptions&ao):
112 : PLUMED_BIAS_INIT(ao),
113 2 : to(getNumberOfArguments(),0.0),
114 1 : min(getNumberOfArguments(),-1.0),
115 1 : kappa(getNumberOfArguments(),0.0),
116 1 : temp(getNumberOfArguments(),0.0),
117 1 : seed(getNumberOfArguments(),std::time(0)),
118 2 : random(getNumberOfArguments())
119 : {
120 : // Note : parseVector will check that number of arguments is correct
121 1 : parseVector("KAPPA",kappa);
122 2 : parseVector("MIN",min);
123 1 : if(min.size()==0) min.assign(getNumberOfArguments(),-1.0);
124 1 : if(min.size()!=getNumberOfArguments()) error("MIN array should have the same size as ARG array");
125 1 : parseVector("NOISE",temp);
126 1 : parseVector("SEED",seed);
127 1 : parseVector("TO",to);
128 1 : checkRead();
129 :
130 1 : log.printf(" min");
131 2 : for(unsigned i=0; i<min.size(); i++) log.printf(" %f",min[i]);
132 1 : log.printf("\n");
133 1 : log.printf(" to");
134 2 : for(unsigned i=0; i<to.size(); i++) log.printf(" %f",to[i]);
135 1 : log.printf("\n");
136 1 : log.printf(" with force constant");
137 2 : for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
138 1 : log.printf("\n");
139 :
140 2 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
141 1 : std::string str_min=getPntrToArgument(i)->getName()+"_min";
142 2 : addComponent(str_min); componentIsNotPeriodic(str_min);
143 1 : if(min[i]!=-1.0) getPntrToComponent(str_min)->set(min[i]);
144 : }
145 2 : for(unsigned i=0; i<getNumberOfArguments(); i++) {random[i].setSeed(-seed[i]);}
146 2 : addComponent("force2"); componentIsNotPeriodic("force2");
147 1 : }
148 :
149 :
150 5 : void ABMD::calculate() {
151 : double ene=0.0;
152 : double totf2=0.0;
153 10 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
154 5 : const double cv=difference(i,to[i],getArgument(i));
155 5 : const double cv2=cv*cv;
156 5 : const double k=kappa[i];
157 : double noise=0.;
158 5 : double diff=temp[i];
159 5 : if(diff>0) {
160 5 : noise = 2.*random[i].Gaussian()*diff;
161 5 : if(cv2<=diff) { diff=0.; temp[i]=0.; }
162 : }
163 :
164 : // min < 0 means that the variable has not been used in the input file, so the current position of the CV is used
165 : // cv2 < min means that the collective variable is nearer to the target value than at any other previous time so
166 : // min is set to the CV value
167 5 : if(min[i]<0.||cv2<min[i]) {
168 1 : min[i] = cv2;
169 : } else {
170 : // otherwise a noise is added to the minimum value
171 4 : min[i] += noise;
172 4 : const double f = -2.*k*(cv2-min[i])*cv;
173 4 : setOutputForce(i,f);
174 4 : ene += 0.5*k*(cv2-min[i])*(cv2-min[i]);
175 4 : totf2+=f*f;
176 : }
177 5 : getPntrToComponent(i+1)->set(min[i]);
178 : }
179 5 : setBias(ene);
180 5 : getPntrToComponent("force2")->set(totf2);
181 5 : }
182 :
183 : }
184 : }
185 :
186 :
|
