Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "core/ActionRegister.h"
24 : #include "tools/Random.h"
25 : #include "core/PlumedMain.h"
26 :
27 : namespace PLMD {
28 : namespace bias {
29 :
30 : //+PLUMEDOC BIAS EXTENDED_LAGRANGIAN
31 : /*
32 : Add extended Lagrangian.
33 :
34 : This action can be used to create fictitious collective variables that are coupled to the real ones.
35 : In this method potential and kinetic contributions are added to the energy of the system using
36 :
37 : $$
38 : V=\sum_i \frac{k_i}{2} (x_i-s_i)^2 + \sum_i \frac{\dot{s}_i^2}{2m_i}
39 : $$.
40 :
41 : In this expression, $x_i$ is the $i$th scalar-argument that is input to the bias potential so
42 : the resulting potential is similar to a [RESTRAINT](RESTRAINT.md) that is forced to take scalars in input.
43 : However, in this method the restraint center moves with time following Hamiltonian
44 : dynamics with mass $m_i$.
45 :
46 : This bias potential accepts thus vectorial keywords (one element per argument)
47 : to define the coupling constant (KAPPA) and a relaxation time $tau$ (TAU).
48 : The mass is them computed as $m=k(\frac{\tau}{2\pi})^2$.
49 :
50 : Notice that this action creates several components.
51 : The ones named XX_fict are the fictitious coordinates. It is possible
52 : to add further forces on them by means of other bias potential,
53 : e.g. to obtain an indirect [METAD](METAD.md) as in the first paper cited below.
54 : Also notice that the velocities of the fictitious coordinates
55 : are reported (XX_vfict). However, printed velocities are the ones
56 : at the previous step.
57 :
58 : It is also possible to provide a non-zero friction (one value per component).
59 : This is then used to implement a Langevin thermostat, so as to implement the
60 : TAMD/dAFED method that is discussed in the second two papers cited below. Notice that
61 : here a massive Langevin thermostat is used, whereas usually
62 : TAMD employs an overamped Langevin dynamics while dAFED employs
63 : a Gaussian thermostat.
64 :
65 : !!!! warning "bias not campatible with replica exchange"
66 :
67 : The bias potential is reported in the component bias.
68 : Notice that this bias potential, although formally compatible with
69 : replica exchange framework, probably does not work as expected in that case.
70 : Indeed, since fictitious coordinates are not swapped upon exchange,
71 : acceptace can be expected to be extremely low unless (by chance) two neighboring
72 : replicas have the fictitious variables located properly in space.
73 :
74 : !!! warning "restart not supported"
75 :
76 : [RESTART](RESTART.md) is not properly supported by this action. Indeed,
77 : at every start the position of the fictitious variable is reset to the value
78 : of the real variable, and its velocity is set to zero.
79 : This is not expected to introduce big errors, but certainly is
80 : introducing a small inconsistency between a single long run
81 : and many shorter runs.
82 :
83 : ## Examples
84 :
85 : The following input tells plumed to perform a metadynamics
86 : with an extended Lagrangian on two torsional angles.
87 :
88 : ```plumed
89 : phi: TORSION ATOMS=5,7,9,15
90 : psi: TORSION ATOMS=7,9,15,17
91 : ex: EXTENDED_LAGRANGIAN ARG=phi,psi KAPPA=20,20.0 TAU=0.1,0.1
92 : METAD ARG=ex.phi_fict,ex.psi_fict PACE=100 SIGMA=0.35,0.35 HEIGHT=0.1
93 : # monitor the two variables
94 : PRINT STRIDE=10 ARG=phi,psi,ex.phi_fict,ex.psi_fict FILE=COLVAR
95 : ```
96 :
97 : The following input tells plumed to perform a TAMD (or dAFED)
98 : calculation on two torsional angles, keeping the two variables
99 : at a fictitious temperature of 3000K with a Langevin thermostat
100 : with friction 10
101 :
102 : ```plumed
103 : phi: TORSION ATOMS=5,7,9,15
104 : psi: TORSION ATOMS=7,9,15,17
105 : ex: EXTENDED_LAGRANGIAN ARG=phi,psi KAPPA=20,20.0 TAU=0.1,0.1 FRICTION=10,10 TEMP=3000
106 : # monitor the two variables
107 : PRINT STRIDE=10 ARG=phi,psi,ex.phi_fict,ex.psi_fict FILE=COLVAR
108 : ```
109 :
110 : */
111 : //+ENDPLUMEDOC
112 :
113 : class ExtendedLagrangian : public Bias {
114 : bool firsttime;
115 : std::vector<double> fict;
116 : std::vector<double> vfict;
117 : std::vector<double> vfict_laststep;
118 : std::vector<double> ffict;
119 : std::vector<double> kappa;
120 : std::vector<double> tau;
121 : std::vector<double> friction;
122 : std::vector<Value*> fictValue;
123 : std::vector<Value*> vfictValue;
124 : double kbt;
125 : Random rand;
126 : public:
127 : explicit ExtendedLagrangian(const ActionOptions&);
128 : void calculate() override;
129 : void update() override;
130 : static void registerKeywords(Keywords& keys);
131 : };
132 :
133 : PLUMED_REGISTER_ACTION(ExtendedLagrangian,"EXTENDED_LAGRANGIAN")
134 :
135 6 : void ExtendedLagrangian::registerKeywords(Keywords& keys) {
136 6 : Bias::registerKeywords(keys);
137 6 : keys.add("compulsory","KAPPA","specifies that the restraint is harmonic and what the values of the force constants on each of the variables are");
138 6 : keys.add("compulsory","TAU","specifies that the restraint is harmonic and what the values of the force constants on each of the variables are");
139 6 : keys.add("compulsory","FRICTION","0.0","add a friction to the variable");
140 6 : keys.add("optional","TEMP","the system temperature - needed when FRICTION is present. If not provided will be taken from MD code (if available)");
141 12 : keys.addOutputComponent("_fict","default","scalar","one or multiple instances of this quantity can be referenced elsewhere in the input file. "
142 : "These quantities will named with the arguments of the bias followed by "
143 : "the character string _tilde. It is possible to add forces on these variable.");
144 12 : keys.addOutputComponent("_vfict","default","scalar","one or multiple instances of this quantity can be referenced elsewhere in the input file. "
145 : "These quantities will named with the arguments of the bias followed by "
146 : "the character string _tilde. It is NOT possible to add forces on these variable.");
147 6 : keys.addDOI("10.1103/PhysRevLett.90.238302");
148 6 : keys.addDOI("10.1016/j.cplett.2006.05.062");
149 6 : keys.addDOI("10.1021/jp805039u");
150 6 : }
151 :
152 2 : ExtendedLagrangian::ExtendedLagrangian(const ActionOptions&ao):
153 : PLUMED_BIAS_INIT(ao),
154 2 : firsttime(true),
155 4 : fict(getNumberOfArguments(),0.0),
156 2 : vfict(getNumberOfArguments(),0.0),
157 2 : vfict_laststep(getNumberOfArguments(),0.0),
158 2 : ffict(getNumberOfArguments(),0.0),
159 2 : kappa(getNumberOfArguments(),0.0),
160 2 : tau(getNumberOfArguments(),0.0),
161 2 : friction(getNumberOfArguments(),0.0),
162 2 : fictValue(getNumberOfArguments(),NULL),
163 2 : vfictValue(getNumberOfArguments(),NULL),
164 4 : kbt(0.0) {
165 2 : parseVector("TAU",tau);
166 2 : parseVector("FRICTION",friction);
167 2 : parseVector("KAPPA",kappa);
168 2 : kbt=getkBT();
169 2 : checkRead();
170 :
171 : bool hasFriction=false;
172 6 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
173 : // the mass of the fictitious particle is KAPPA*TAU^2/(4 pi^2)
174 4 : if(!(kappa[i]>0.0)) {
175 0 : error("KAPPA should be positive");
176 : }
177 4 : if(!(tau[i]>0.0)) {
178 0 : error("TAU should be positive");
179 : }
180 4 : if(friction[i]<0.0) {
181 0 : error("FRICTION should not be negative");
182 : }
183 4 : if(friction[i]>0.0) {
184 : hasFriction=true;
185 : }
186 : }
187 2 : if(hasFriction && !(kbt>0.0)) {
188 2 : log.printf(" WARNING: FRICTION without a temperature (TEMP, or passed by the MD code): the fictitious particles are cooled to zero temperature\n");
189 : }
190 :
191 2 : log.printf(" with harmonic force constant");
192 6 : for(unsigned i=0; i<kappa.size(); i++) {
193 4 : log.printf(" %f",kappa[i]);
194 : }
195 2 : log.printf("\n");
196 :
197 2 : log.printf(" with relaxation time");
198 6 : for(unsigned i=0; i<tau.size(); i++) {
199 4 : log.printf(" %f",tau[i]);
200 : }
201 2 : log.printf("\n");
202 :
203 2 : if(hasFriction) {
204 2 : log.printf(" with friction");
205 6 : for(unsigned i=0; i<friction.size(); i++) {
206 4 : log.printf(" %f",friction[i]);
207 : }
208 2 : log.printf("\n");
209 : }
210 :
211 2 : log.printf(" and kbt");
212 2 : log.printf(" %f",kbt);
213 2 : log.printf("\n");
214 :
215 6 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
216 4 : std::string comp=getPntrToArgument(i)->getName()+"_fict";
217 4 : addComponentWithDerivatives(comp);
218 4 : if(getPntrToArgument(i)->isPeriodic()) {
219 : std::string a,b;
220 4 : getPntrToArgument(i)->getDomain(a,b);
221 4 : componentIsPeriodic(comp,a,b);
222 : } else {
223 0 : componentIsNotPeriodic(comp);
224 : }
225 4 : fictValue[i]=getPntrToComponent(comp);
226 8 : comp=getPntrToArgument(i)->getName()+"_vfict";
227 4 : addComponent(comp);
228 4 : componentIsNotPeriodic(comp);
229 4 : vfictValue[i]=getPntrToComponent(comp);
230 : }
231 :
232 2 : log<<" Bibliography "<<plumed.cite("Iannuzzi, Laio, and Parrinello, Phys. Rev. Lett. 90, 238302 (2003)");
233 2 : if(hasFriction) {
234 2 : log<<plumed.cite("Maragliano and Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006)");
235 4 : log<<plumed.cite("Abrams and Tuckerman, J. Phys. Chem. B 112, 15742 (2008)");
236 : }
237 2 : log<<"\n";
238 2 : }
239 :
240 :
241 24 : void ExtendedLagrangian::calculate() {
242 :
243 24 : if(firsttime) {
244 6 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
245 4 : fict[i]=getArgument(i);
246 : }
247 2 : firsttime=false;
248 : }
249 : double ene=0.0;
250 72 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
251 48 : const double cv=difference(i,fict[i],getArgument(i));
252 48 : const double k=kappa[i];
253 48 : const double f=-k*cv;
254 48 : ene+=0.5*k*cv*cv;
255 48 : setOutputForce(i,f);
256 48 : ffict[i]=-f;
257 : };
258 24 : setBias(ene);
259 72 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
260 48 : fict[i]=fictValue[i]->bringBackInPbc(fict[i]);
261 48 : fictValue[i]->set(fict[i]);
262 48 : vfictValue[i]->set(vfict_laststep[i]);
263 : }
264 24 : }
265 :
266 24 : void ExtendedLagrangian::update() {
267 24 : double dt=getTimeStep()*getStride();
268 72 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
269 48 : double mass=kappa[i]*tau[i]*tau[i]/(4*pi*pi); // should be k/omega**2
270 48 : double c1=std::exp(-0.5*friction[i]*dt);
271 48 : double c2=std::sqrt(kbt*(1.0-c1*c1)/mass);
272 : // consider additional forces on the fictitious particle
273 : // (e.g. MetaD stuff)
274 48 : ffict[i]+=fictValue[i]->getForce();
275 :
276 : // update velocity (half step)
277 48 : vfict[i]+=ffict[i]*0.5*dt/mass;
278 : // thermostat (half step)
279 48 : vfict[i]=c1*vfict[i]+c2*rand.Gaussian();
280 : // save full step velocity to be dumped at next step
281 48 : vfict_laststep[i]=vfict[i];
282 : // thermostat (half step)
283 48 : vfict[i]=c1*vfict[i]+c2*rand.Gaussian();
284 : // update velocity (half step)
285 48 : vfict[i]+=ffict[i]*0.5*dt/mass;
286 : // update position (full step)
287 48 : fict[i]+=vfict[i]*dt;
288 : }
289 24 : }
290 :
291 : }
292 :
293 : }
|
