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.addFlag("NOSPLINE",false,"specifies that no spline interpolation is to be used when calculating the energy and forces due to the external potential");
148 6 : keys.addFlag("SPARSE",false,"specifies that the external potential uses a sparse grid");
149 6 : keys.addFlag("SPHERE",false, "The Funnel potential including the binding site can be spherical instead of a cone");
150 6 : keys.add("compulsory","SCALE","1.0","a factor that multiplies the external potential, useful to invert free energies");
151 : // old stuff?
152 : // componentsAreNotOptional(keys);
153 : // keys.addOutputComponent("bias","default","the instantaneous value of the bias potential");
154 :
155 : //Defining optional arguments
156 6 : keys.add("optional","NBINS","number of bins along fps.lp");
157 6 : keys.add("optional","NBINZ","number of bins along fps.ld");
158 6 : keys.add("optional","MINS","minimum value assumed by fps.lp, if the ligand is able to go beyond this value the simulation will crash");
159 6 : keys.add("optional","KAPPA","constant to be used for the funnel-shape restraint potential");
160 6 : keys.add("optional","RCYL","radius of the cylindrical section");
161 6 : keys.add("optional","SAFETY","To be used in case the SPHERE flag is chosen, it regulates how much the potential extends (in nm)");
162 6 : 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");
163 6 : keys.add("optional","ALPHA","angle to change the width of the cone section");
164 : //Defining compulsory arguments
165 6 : keys.add("compulsory","MAXS","MAXS","maximum value assumed by fps.lp");
166 6 : keys.add("compulsory","ZCC","ZCC","switching point between cylinder and cone");
167 6 : keys.add("compulsory","FILE","name of the Funnel potential file");
168 6 : keys.addFlag("WALKERS_MPI",false,"To be used when gromacs + multiple walkers are used");
169 6 : }
170 :
171 : // Old version 2.3
172 : //Funnel::~Funnel(){
173 : // delete BiasGrid_;
174 : //}
175 :
176 4 : Funnel::Funnel(const ActionOptions& ao):
177 : PLUMED_BIAS_INIT(ao),
178 : // Old version 2.3
179 : // BiasGrid_(NULL),
180 4 : NBINS(500.0),
181 4 : NBINZ(500.0),
182 4 : MINS(0.0),
183 4 : KAPPA(84.0),
184 4 : RCYL(0.1),
185 4 : safety(1.0),
186 4 : slope(1.0),
187 4 : ALPHA(1.413)
188 :
189 : {
190 4 : bool sparsegrid=false;
191 4 : parseFlag("SPARSE",sparsegrid);
192 4 : bool nospline=false;
193 4 : parseFlag("NOSPLINE",nospline);
194 4 : bool spline=!nospline;
195 4 : bool walkers_mpi=false;
196 4 : parseFlag("WALKERS_MPI",walkers_mpi);
197 : // bool components=false;
198 : // parseFlag("POINTS",components);
199 4 : bool sphere=false;
200 4 : parseFlag("SPHERE",sphere);
201 8 : parse("SAFETY",safety);
202 : string file;
203 8 : parse("FILE",file);
204 4 : if( file.length()==0 ) {
205 0 : error("No funnel file name was specified");
206 : }
207 4 : parse("SCALE",scale_);
208 :
209 : //Reading optional arguments
210 4 : parse("KAPPA",KAPPA);
211 4 : parse("NBINS",NBINS);
212 4 : parse("NBINZ",NBINZ);
213 4 : parse("MINS",MINS);
214 4 : parse("RCYL",RCYL);
215 4 : parse("SLOPE",slope);
216 4 : parse("ALPHA",ALPHA);
217 : //Reading compulsory arguments
218 4 : parse("MAXS",MAXS);
219 4 : parse("ZCC",ZCC);
220 :
221 :
222 4 : checkRead();
223 :
224 4 : log.printf(" External potential from file %s\n",file.c_str());
225 4 : log.printf(" Multiplied by %lf\n",scale_);
226 4 : if(spline) {
227 4 : log.printf(" External potential uses spline interpolation\n");
228 : }
229 4 : if(sparsegrid) {
230 0 : log.printf(" External potential uses sparse grid\n");
231 : }
232 :
233 : // Non più necessario dalla 2.3
234 : // addComponent("bias"); componentIsNotPeriodic("bias");
235 :
236 4 : std::string funcl=getLabel() + ".bias";
237 :
238 : // int size = plumed.comm.Get_size();
239 : // int rank = plumed.comm.Get_rank();
240 4 : IFile I_hate_this;
241 4 : bool do_exist=I_hate_this.FileExist(file);
242 :
243 4 : if(walkers_mpi) {
244 2 : if(comm.Get_rank()==0 && multi_sim_comm.Get_rank()==0) {
245 1 : if(!do_exist) {
246 1 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
247 : }
248 : }
249 2 : multi_sim_comm.Barrier();
250 : } else {
251 2 : if(comm.Get_rank()==0) {
252 2 : if(!do_exist) {
253 2 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
254 : }
255 : }
256 : }
257 :
258 : /*
259 : if(comm.Get_rank()==0){
260 : if(multi_sim_comm.Get_rank()==0 && walkers_mpi){
261 : if(!do_exist){
262 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
263 : }
264 : } else {
265 : if(!do_exist){
266 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
267 : }
268 : }
269 : if(walkers_mpi) multi_sim_comm.Barrier();
270 : }
271 : */
272 4 : comm.Barrier();
273 :
274 : // read grid
275 4 : IFile gridfile;
276 4 : gridfile.open(file);
277 8 : BiasGrid_=Grid::create(funcl,getArguments(),gridfile,sparsegrid,spline,true);
278 : //not necessary anymore? gridfile.close();
279 4 : if(BiasGrid_->getDimension()!=getNumberOfArguments()) {
280 0 : error("mismatch between dimensionality of input grid and number of arguments");
281 : }
282 12 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
283 16 : if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) {
284 0 : error("periodicity mismatch between arguments and input bias");
285 : }
286 : }
287 4 : comm.Barrier();
288 4 : if(comm.Get_rank()==0 && walkers_mpi) {
289 2 : multi_sim_comm.Barrier();
290 : }
291 8 : log<<" Bibliography "<<plumed.cite("Limongelli, Bonomi, and Parrinello, PNAS 110, 6358 (2013)")<<"\n";
292 8 : }
293 :
294 :
295 3 : void Funnel::createBIAS(const double& R_cyl, const double& z_cc, const double& alpha,
296 : const double& KAPPA, const double& MIN_S, const double& MAX_S, const double& NBIN_S,
297 : const double& NBIN_Z, const double& safety, const bool& sphere, const double& slope,
298 : const string& funcl, const string& file) {
299 : //R_cyl and z_cc forms the parameters of the cylinder.
300 : //alpha defines the angle in degrees.
301 :
302 : //PARAMETERS OF THE CONE
303 3 : double tg_alpha= tan(alpha);
304 :
305 : //parameters for PROGRESSION
306 : //parameters for DISTANCE
307 : double MIN_Z=0;
308 : double MAX_Z;
309 3 : if (sphere==false) {
310 2 : MAX_Z=R_cyl + tg_alpha * (z_cc - MIN_S);
311 : } else {
312 1 : MAX_Z=z_cc+safety;
313 : }
314 :
315 : //bin size
316 3 : double DX_Z = (MAX_Z - MIN_Z) / NBIN_Z;
317 : double DX_S;
318 3 : if (sphere==false) {
319 2 : DX_S=(MAX_S - MIN_S) / NBIN_S;
320 : } else {
321 1 : DX_S=(MAX_S + z_cc + safety)/NBIN_S;
322 : }
323 :
324 : double SS, Zmax, ZZ, D, d;
325 : double POT, FZ, FS;
326 :
327 3 : PLMD::OFile pof;
328 3 : pof.open(file);
329 :
330 : //Write the header
331 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());
332 3 : if (sphere==false) {
333 2 : pof.printf("#! SET min_%s %f\n", getPntrToArgument(0)->getName().c_str(), MIN_S);
334 : } else {
335 1 : pof.printf("#! SET min_%s %f\n", getPntrToArgument(0)->getName().c_str(), -z_cc-safety);
336 : }
337 3 : pof.printf("#! SET max_%s %f\n", getPntrToArgument(0)->getName().c_str(), MAX_S);
338 3 : pof.printf("#! SET nbins_%s %f\n", getPntrToArgument(0)->getName().c_str(), NBIN_S);
339 3 : pof.printf("#! SET periodic_%s false\n", getPntrToArgument(0)->getName().c_str());
340 3 : pof.printf("#! SET min_%s %f\n", getPntrToArgument(1)->getName().c_str(), MIN_Z);
341 3 : pof.printf("#! SET max_%s %f\n", getPntrToArgument(1)->getName().c_str(), MAX_Z);
342 3 : pof.printf("#! SET nbins_%s %f\n", getPntrToArgument(1)->getName().c_str(), NBIN_Z);
343 3 : pof.printf("#! SET periodic_%s false\n", getPntrToArgument(1)->getName().c_str());
344 :
345 : //Calculate and write the GRID
346 : //Cone or cylinder?
347 :
348 1506 : for(int is=0; is <= NBIN_S; is++) {
349 1503 : if (sphere==false) {
350 1002 : SS = MIN_S + is * DX_S;
351 : } else {
352 501 : SS = - z_cc - safety + is * DX_S;
353 : }
354 : bool cone = false;
355 1503 : if (sphere==false) {
356 1002 : if(SS <= z_cc) {
357 : cone = true;
358 : }
359 : } else {
360 501 : if (SS <= sqrt(pow(z_cc,2)-pow(R_cyl,2))) {
361 : cone = true;
362 : }
363 : }
364 : //Set wall boundaries properly
365 : if(cone == true) {
366 902 : if(sphere==false) {
367 548 : Zmax = R_cyl + (z_cc - SS) * tg_alpha;
368 : } else {
369 354 : if (SS > -z_cc) {
370 277 : Zmax = sqrt(pow(z_cc,2) - pow(SS,2));
371 : } else {
372 : Zmax = 0;
373 : }
374 : }
375 : } else {
376 601 : Zmax = R_cyl;
377 : }
378 :
379 754506 : for(int iz=0; iz <= NBIN_Z; iz++) {
380 753003 : ZZ = MIN_Z + iz * DX_Z;
381 :
382 : //Inside or outside?
383 : bool inside;
384 753003 : if(ZZ < Zmax) {
385 : inside = true;
386 : } else {
387 : inside = false;
388 : }
389 :
390 : if(inside == true) {
391 : POT = 0;
392 : FS = 0;
393 : FZ = 0;
394 : } else {
395 518149 : if(cone == true) {
396 235130 : if(sphere==false) {
397 127824 : POT = 0.5 * KAPPA * (ZZ - Zmax) * (ZZ - Zmax);
398 127824 : FZ = - KAPPA * (ZZ - Zmax);
399 127824 : FS = - KAPPA * (ZZ - Zmax) * tg_alpha;
400 : } else {
401 107306 : D = sqrt(pow(ZZ,2)+pow(SS,2));
402 107306 : d = D - z_cc;
403 107306 : POT = 0.5 * KAPPA * pow(d,2);
404 107306 : FZ = - KAPPA * d * ZZ / D;
405 107306 : FS = - KAPPA * d * SS / D;
406 : }
407 : } else {
408 283019 : if(sphere==false) {
409 212018 : POT = 0.5 * KAPPA * (ZZ - Zmax) * (ZZ - Zmax);
410 212018 : FZ = - KAPPA * (ZZ - Zmax);
411 : FS = 0;
412 : } else {
413 71001 : D = sqrt(pow(ZZ,2)+pow(SS,2));
414 71001 : d = D - z_cc;
415 71001 : if(ZZ>=R_cyl+slope*(SS-z_cc)) {
416 45984 : POT = 0.5 * KAPPA * pow(d,2);
417 45984 : FZ = - KAPPA * d * ZZ / D;
418 45984 : FS = - KAPPA * d * SS / D;
419 : } else {
420 25017 : POT = 0.5 * KAPPA * pow(sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2))-
421 : z_cc,2);
422 25017 : FZ = - KAPPA*(sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2))-z_cc)*
423 25017 : ZZ/sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2));
424 : FS = 0;
425 : }
426 : }
427 : }
428 : }
429 753003 : pof.printf("%13.6lf %13.6lf %13.6lf %13.6lf %13.6lf\n", SS, ZZ, POT, FS, FZ);
430 : }
431 1503 : pof.printf("\n");
432 : }
433 3 : pof.close();
434 3 : }
435 :
436 120 : void Funnel::calculate() {
437 120 : unsigned ncv=getNumberOfArguments();
438 120 : vector<double> cv(ncv), der(ncv);
439 :
440 360 : for(unsigned i=0; i<ncv; ++i) {
441 240 : cv[i]=getArgument(i);
442 : }
443 :
444 : // log.printf(" In Funnel: %13.6lf %13.6lf\n", cv[0], cv[1]);
445 :
446 120 : double ene=scale_*BiasGrid_->getValueAndDerivatives(cv,der);
447 :
448 120 : setBias(ene);
449 :
450 : // set Forces
451 360 : for(unsigned i=0; i<ncv; ++i) {
452 240 : const double f=-scale_*der[i];
453 240 : setOutputForce(i,f);
454 : }
455 120 : }
456 :
457 : }
458 : }
|
