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