Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2019-2020
3 :
4 : This file is part of funnel code module.
5 :
6 : The FM code respects the CC BY-NC license.
7 : Users are free to download, adapt and use the code as long as it is not for commercial purposes.
8 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
9 : #include "bias/Bias.h"
10 : #include "core/ActionRegister.h"
11 : #include "tools/Grid.h"
12 : #include "tools/Exception.h"
13 : #include "tools/File.h"
14 : #include <cstring>
15 : #include "tools/Communicator.h"
16 : #include "core/ActionSet.h"
17 : #include "tools/FileBase.h"
18 : #include <memory>
19 : #include "core/PlumedMain.h"
20 :
21 : using namespace std;
22 : using namespace PLMD::bias;
23 :
24 :
25 : namespace PLMD {
26 : namespace funnel {
27 :
28 : //+PLUMEDOC BIAS FUNNEL
29 : /*
30 : Calculate a funnel-shape restraint potential that is defined on a grid that is read during the setup.
31 :
32 : If the input file is not already present, it will create one with the name specified in the FILE flag.
33 : The potential has a two-dimensional resolution since it has been devised to be used with the two
34 : components of [FUNNEL_PS](FUNNEL_PS.md) (i.e., fps.lp and fps.ld) and it is divided in two sections, a cone shape
35 : attached to a cylindrical one. The user can customize the shape of both the sections by modifying a
36 : number of flags. In particular the cone section of the funnel is calculated with the following formula:
37 :
38 : $$
39 : MAX_Z=R_{cyl} + tg_{alpha} * (z_{cc} - MIN_S)
40 : $$
41 :
42 : where $ MAX_Z $ is the radius of the cone base, $ R_{cyl} $ is the radius of the cylinder part,
43 : $ tg_{alpha} $ is the angle regulating how steep the cone is, $ z_{cc} $ is the switching point
44 : between cone and cylinder, and $ MIN_S $ is the lowest possible value assumed by fps.lp of [FUNNEL_PS](FUNNEL_PS.md).
45 : As for the cylinder, it starts from the value of $ z_{cc} $ and stops at the value of $ MAX_S $
46 : with a section of $ pi*r_{cyl}^2 $.
47 :
48 : There is the option of transforming the cone region into a sphere with the use of the SPHERE flag. In this
49 : case, the new shape will have a radius of $ z_{cc} $. It might be necessary tuning the SAFETY option
50 : to select how much the potential extends from the sphere.
51 :
52 : ## Examples
53 :
54 : The following is an input for a calculation with a funnel potential that is defined in the file BIAS
55 : and that acts on the collective variables defined by FUNNEL_PS.
56 :
57 : ```plumed
58 : #SETTINGS INPUTFILES=regtest/funnel/rt-funnel/BIAS.reference INPUTFILELINES=1-15
59 :
60 : lig: COM ATOMS=3221,3224,3225,3228,3229,3231,3233,3235,3237
61 : fps: FUNNEL_PS LIGAND=lig REFERENCE=start.pdb ANCHOR=2472 POINTS=4.724,5.369,4.069,4.597,5.721,4.343
62 :
63 : funnel: FUNNEL ...
64 : ARG=fps.lp,fps.ld ZCC=1.8 ALPHA=0.55
65 : RCYL=0.1 MINS=-0.5 MAXS=3.7 KAPPA=35100
66 : NBINS=500 NBINZ=500
67 : FILE=regtest/funnel/rt-funnel/BIAS.reference
68 : ...
69 : ```
70 :
71 : The Funnel potential should always be used in combination with the collective variable [FUNNEL_PS](FUNNEL_PS.md), since it
72 : is constructed to take as inputs fps.lp and fps.ld (the former linepos and linedist of Funnel-Metadynamics FM).
73 : In the first block of data the value of fps.lp (the value in the first column) is kept fixed
74 : and the value of the function is given at 500 equally spaced values for fps.ld between 0 and 1.51. In
75 : the second block of data fps.lp is fixed at $-0.5 + \frac{4.2}{500}$ and the value of the function
76 : is given at 500 equally spaced values for fps.ld between 0 and 1.51. In the third block of data the same
77 : is done but fps.lp is fixed at $-0.5 + \frac{8.4}{100}$ and so on until you get to the five hundredth
78 : block of data where fps.lp is fixed at $3.7$.
79 :
80 : It is possible to switch the shape of the cone region, transforming it in a sphere, with the flag SPHERE.
81 :
82 : ```plumed
83 : #SETTINGS INPUTFILES=regtest/funnel/rt-funnel/BIAS.reference INPUTFILELINES=1-15
84 :
85 : lig: COM ATOMS=545,546,547,548,549,550,551,552,553
86 : fps: FUNNEL_PS LIGAND=lig REFERENCE=ref.pdb ANCHOR=52 POINTS=2.793,3.696,3.942,3.607,4.298,3.452
87 :
88 : funnel: FUNNEL ...
89 : ARG=fps.lp,fps.ld ZCC=4.0 RCYL=0.1 MINS=0.2
90 : MAXS=4.9 KAPPA=100000 NBINS=500
91 : NBINZ=500 SPHERE SAFETY=1.0
92 : FILE=regtest/funnel/rt-funnel/BIAS.reference
93 : ...
94 : ```
95 :
96 : */
97 : //+ENDPLUMEDOC
98 : class Funnel : public Bias {
99 :
100 : private:
101 : std::unique_ptr<GridBase> BiasGrid_;
102 :
103 : /////////////////////
104 : // old 2.3
105 : //Grid* BiasGrid_;
106 : /////////////////////
107 : //Optional parameters
108 : double NBINS;
109 : double NBINZ;
110 : double MINS;
111 : double KAPPA;
112 : double RCYL;
113 : double safety;
114 : double slope;
115 : double ALPHA;
116 : //Compulsory parameters
117 : double MAXS;
118 : double ZCC;
119 : double scale_;
120 :
121 :
122 : public:
123 : explicit Funnel(const ActionOptions&);
124 :
125 : // old Funnel-2.3
126 : // ~Funnel();
127 :
128 : void calculate();
129 : static void registerKeywords(Keywords& keys);
130 : void createBIAS(const double& R_cyl, const double& z_cc, const double& alpha, const double& KAPPA,
131 : const double& MIN_S, const double& MAX_S, const double& NBIN_S, const double& NBIN_Z,
132 : const double& safety, const bool& sphere, const double& slope, const string& funcl,
133 : const string& file);
134 : // void createBIAS3D(const double& R_cyl, const double& z_cc, const double& alpha,
135 : // const double& KAPPA, const double& MIN_S, const double& MAX_S, const double& NBIN_S,
136 : // const double& NBIN_Z, const double& safety, const bool& sphere, const double& slope,
137 : // const string& funcl);
138 : };
139 :
140 : PLUMED_REGISTER_ACTION(Funnel,"FUNNEL")
141 :
142 8 : void Funnel::registerKeywords(Keywords& keys) {
143 8 : Bias::registerKeywords(keys);
144 8 : keys.addFlag("NOSPLINE",false,"specifies that no spline interpolation is to be used when calculating the energy and forces due to the external potential");
145 8 : keys.addFlag("SPARSE",false,"specifies that the external potential uses a sparse grid");
146 8 : keys.addFlag("SPHERE",false, "The Funnel potential including the binding site can be spherical instead of a cone");
147 8 : keys.add("compulsory","SCALE","1.0","a factor that multiplies the external potential, useful to invert free energies");
148 : // old stuff?
149 : // componentsAreNotOptional(keys);
150 : // keys.addOutputComponent("bias","default","the instantaneous value of the bias potential");
151 :
152 : //Defining optional arguments
153 8 : keys.add("optional","NBINS","number of bins along fps.lp");
154 8 : keys.add("optional","NBINZ","number of bins along fps.ld");
155 8 : keys.add("optional","MINS","minimum value assumed by fps.lp, if the ligand is able to go beyond this value the simulation will crash");
156 8 : keys.add("optional","KAPPA","constant to be used for the funnel-shape restraint potential");
157 8 : keys.add("optional","RCYL","radius of the cylindrical section");
158 8 : keys.add("optional","SAFETY","To be used in case the SPHERE flag is chosen, it regulates how much the potential extends (in nm)");
159 8 : keys.add("optional","SLOPE","Adjust the behavior of the potential outside the funnel, greater values than 1.0 will tend to push the ligand more towards the cylinder and vice versa");
160 8 : keys.add("optional","ALPHA","angle to change the width of the cone section");
161 : //Defining compulsory arguments
162 8 : keys.add("compulsory","MAXS","MAXS","maximum value assumed by fps.lp");
163 8 : keys.add("compulsory","ZCC","ZCC","switching point between cylinder and cone");
164 8 : keys.add("compulsory","FILE","name of the Funnel potential file");
165 8 : keys.addFlag("WALKERS_MPI",false,"To be used when gromacs + multiple walkers are used");
166 8 : }
167 :
168 : // Old version 2.3
169 : //Funnel::~Funnel(){
170 : // delete BiasGrid_;
171 : //}
172 :
173 4 : Funnel::Funnel(const ActionOptions& ao):
174 : PLUMED_BIAS_INIT(ao),
175 : // Old version 2.3
176 : // BiasGrid_(NULL),
177 4 : NBINS(500.0),
178 4 : NBINZ(500.0),
179 4 : MINS(0.0),
180 4 : KAPPA(84.0),
181 4 : RCYL(0.1),
182 4 : safety(1.0),
183 4 : slope(1.0),
184 4 : ALPHA(1.413)
185 :
186 : {
187 4 : bool sparsegrid=false;
188 4 : parseFlag("SPARSE",sparsegrid);
189 4 : bool nospline=false;
190 4 : parseFlag("NOSPLINE",nospline);
191 4 : bool spline=!nospline;
192 4 : bool walkers_mpi=false;
193 4 : parseFlag("WALKERS_MPI",walkers_mpi);
194 4 : bool sphere=false;
195 4 : parseFlag("SPHERE",sphere);
196 8 : parse("SAFETY",safety);
197 : string file;
198 8 : parse("FILE",file);
199 4 : if( file.length()==0 ) {
200 0 : error("No funnel file name was specified");
201 : }
202 4 : parse("SCALE",scale_);
203 :
204 : //Reading optional arguments
205 4 : parse("KAPPA",KAPPA);
206 4 : parse("NBINS",NBINS);
207 4 : parse("NBINZ",NBINZ);
208 4 : parse("MINS",MINS);
209 4 : parse("RCYL",RCYL);
210 4 : parse("SLOPE",slope);
211 4 : parse("ALPHA",ALPHA);
212 : //Reading compulsory arguments
213 4 : parse("MAXS",MAXS);
214 4 : parse("ZCC",ZCC);
215 :
216 :
217 4 : checkRead();
218 :
219 4 : log.printf(" External potential from file %s\n",file.c_str());
220 4 : log.printf(" Multiplied by %lf\n",scale_);
221 4 : if(spline) {
222 4 : log.printf(" External potential uses spline interpolation\n");
223 : }
224 4 : if(sparsegrid) {
225 0 : log.printf(" External potential uses sparse grid\n");
226 : }
227 :
228 : // Non più necessario dalla 2.3
229 : // addComponent("bias"); componentIsNotPeriodic("bias");
230 :
231 4 : std::string funcl=getLabel() + ".bias";
232 :
233 4 : IFile I_hate_this;
234 4 : bool do_exist=I_hate_this.FileExist(file);
235 :
236 4 : if(walkers_mpi) {
237 2 : if(comm.Get_rank()==0 && multi_sim_comm.Get_rank()==0) {
238 1 : if(!do_exist) {
239 1 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
240 : }
241 : }
242 2 : multi_sim_comm.Barrier();
243 : } else {
244 2 : if(comm.Get_rank()==0) {
245 2 : if(!do_exist) {
246 2 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
247 : }
248 : }
249 : }
250 :
251 4 : comm.Barrier();
252 :
253 : // read grid
254 4 : IFile gridfile;
255 4 : gridfile.open(file);
256 8 : BiasGrid_=Grid::create(funcl,getArguments(),gridfile,sparsegrid,spline,true);
257 4 : if(BiasGrid_->getDimension()!=getNumberOfArguments()) {
258 0 : error("mismatch between dimensionality of input grid and number of arguments");
259 : }
260 12 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
261 16 : if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) {
262 0 : error("periodicity mismatch between arguments and input bias");
263 : }
264 : }
265 4 : comm.Barrier();
266 4 : if(comm.Get_rank()==0 && walkers_mpi) {
267 2 : multi_sim_comm.Barrier();
268 : }
269 4 : log<<" Bibliography "<<plumed.cite("Limongelli, Bonomi, and Parrinello, PNAS 110, 6358 (2013)")<<"\n";
270 8 : }
271 :
272 :
273 3 : void Funnel::createBIAS(const double& R_cyl, const double& z_cc, const double& alpha,
274 : const double& KAPPA, const double& MIN_S, const double& MAX_S, const double& NBIN_S,
275 : const double& NBIN_Z, const double& safety, const bool& sphere, const double& slope,
276 : const string& funcl, const string& file) {
277 : //R_cyl and z_cc forms the parameters of the cylinder.
278 : //alpha defines the angle in degrees.
279 :
280 : //PARAMETERS OF THE CONE
281 3 : double tg_alpha= tan(alpha);
282 :
283 : //parameters for PROGRESSION
284 : //parameters for DISTANCE
285 : double MIN_Z=0;
286 : double MAX_Z;
287 3 : if (sphere==false) {
288 2 : MAX_Z=R_cyl + tg_alpha * (z_cc - MIN_S);
289 : } else {
290 1 : MAX_Z=z_cc+safety;
291 : }
292 :
293 : //bin size
294 3 : double DX_Z = (MAX_Z - MIN_Z) / NBIN_Z;
295 : double DX_S;
296 3 : if (sphere==false) {
297 2 : DX_S=(MAX_S - MIN_S) / NBIN_S;
298 : } else {
299 1 : DX_S=(MAX_S + z_cc + safety)/NBIN_S;
300 : }
301 :
302 : double SS, Zmax, ZZ, D, d;
303 : double POT, FZ, FS;
304 :
305 3 : PLMD::OFile pof;
306 3 : pof.open(file);
307 :
308 : //Write the header
309 3 : pof.printf("#! FIELDS %s %s %s der_%s der_%s \n", getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str(), funcl.c_str(), getPntrToArgument(0)->getName().c_str(), getPntrToArgument(1)->getName().c_str());
310 3 : if (sphere==false) {
311 2 : pof.printf("#! SET min_%s %f\n", getPntrToArgument(0)->getName().c_str(), MIN_S);
312 : } else {
313 1 : pof.printf("#! SET min_%s %f\n", getPntrToArgument(0)->getName().c_str(), -z_cc-safety);
314 : }
315 3 : pof.printf("#! SET max_%s %f\n", getPntrToArgument(0)->getName().c_str(), MAX_S);
316 3 : pof.printf("#! SET nbins_%s %f\n", getPntrToArgument(0)->getName().c_str(), NBIN_S);
317 3 : pof.printf("#! SET periodic_%s false\n", getPntrToArgument(0)->getName().c_str());
318 3 : pof.printf("#! SET min_%s %f\n", getPntrToArgument(1)->getName().c_str(), MIN_Z);
319 3 : pof.printf("#! SET max_%s %f\n", getPntrToArgument(1)->getName().c_str(), MAX_Z);
320 3 : pof.printf("#! SET nbins_%s %f\n", getPntrToArgument(1)->getName().c_str(), NBIN_Z);
321 3 : pof.printf("#! SET periodic_%s false\n", getPntrToArgument(1)->getName().c_str());
322 :
323 : //Calculate and write the GRID
324 : //Cone or cylinder?
325 :
326 1506 : for(int is=0; is <= NBIN_S; is++) {
327 1503 : if (sphere==false) {
328 1002 : SS = MIN_S + is * DX_S;
329 : } else {
330 501 : SS = - z_cc - safety + is * DX_S;
331 : }
332 : bool cone = false;
333 1503 : if (sphere==false) {
334 1002 : if(SS <= z_cc) {
335 : cone = true;
336 : }
337 : } else {
338 501 : if (SS <= sqrt(pow(z_cc,2)-pow(R_cyl,2))) {
339 : cone = true;
340 : }
341 : }
342 : //Set wall boundaries properly
343 : if(cone == true) {
344 902 : if(sphere==false) {
345 548 : Zmax = R_cyl + (z_cc - SS) * tg_alpha;
346 : } else {
347 354 : if (SS > -z_cc) {
348 277 : Zmax = sqrt(pow(z_cc,2) - pow(SS,2));
349 : } else {
350 : Zmax = 0;
351 : }
352 : }
353 : } else {
354 601 : Zmax = R_cyl;
355 : }
356 :
357 754506 : for(int iz=0; iz <= NBIN_Z; iz++) {
358 753003 : ZZ = MIN_Z + iz * DX_Z;
359 :
360 : //Inside or outside?
361 : bool inside;
362 753003 : if(ZZ < Zmax) {
363 : inside = true;
364 : } else {
365 : inside = false;
366 : }
367 :
368 : if(inside == true) {
369 : POT = 0;
370 : FS = 0;
371 : FZ = 0;
372 : } else {
373 518149 : if(cone == true) {
374 235130 : if(sphere==false) {
375 127824 : POT = 0.5 * KAPPA * (ZZ - Zmax) * (ZZ - Zmax);
376 127824 : FZ = - KAPPA * (ZZ - Zmax);
377 127824 : FS = - KAPPA * (ZZ - Zmax) * tg_alpha;
378 : } else {
379 107306 : D = sqrt(pow(ZZ,2)+pow(SS,2));
380 107306 : d = D - z_cc;
381 107306 : POT = 0.5 * KAPPA * pow(d,2);
382 107306 : FZ = - KAPPA * d * ZZ / D;
383 107306 : FS = - KAPPA * d * SS / D;
384 : }
385 : } else {
386 283019 : if(sphere==false) {
387 212018 : POT = 0.5 * KAPPA * (ZZ - Zmax) * (ZZ - Zmax);
388 212018 : FZ = - KAPPA * (ZZ - Zmax);
389 : FS = 0;
390 : } else {
391 71001 : D = sqrt(pow(ZZ,2)+pow(SS,2));
392 71001 : d = D - z_cc;
393 71001 : if(ZZ>=R_cyl+slope*(SS-z_cc)) {
394 45984 : POT = 0.5 * KAPPA * pow(d,2);
395 45984 : FZ = - KAPPA * d * ZZ / D;
396 45984 : FS = - KAPPA * d * SS / D;
397 : } else {
398 25017 : POT = 0.5 * KAPPA * pow(sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2))-
399 : z_cc,2);
400 25017 : FZ = - KAPPA*(sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2))-z_cc)*
401 25017 : ZZ/sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2));
402 : FS = 0;
403 : }
404 : }
405 : }
406 : }
407 753003 : pof.printf("%13.6lf %13.6lf %13.6lf %13.6lf %13.6lf\n", SS, ZZ, POT, FS, FZ);
408 : }
409 1503 : pof.printf("\n");
410 : }
411 3 : pof.close();
412 3 : }
413 :
414 120 : void Funnel::calculate() {
415 120 : unsigned ncv=getNumberOfArguments();
416 120 : vector<double> cv(ncv), der(ncv);
417 :
418 360 : for(unsigned i=0; i<ncv; ++i) {
419 240 : cv[i]=getArgument(i);
420 : }
421 :
422 :
423 120 : double ene=scale_*BiasGrid_->getValueAndDerivatives(cv,der);
424 :
425 120 : setBias(ene);
426 :
427 : // set Forces
428 360 : for(unsigned i=0; i<ncv; ++i) {
429 240 : const double f=-scale_*der[i];
430 240 : setOutputForce(i,f);
431 : }
432 120 : }
433 :
434 : }
435 : }
|
