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