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 ) error("No funnel file name was specified");
206 4 : parse("SCALE",scale_);
207 :
208 : //Reading optional arguments
209 4 : parse("KAPPA",KAPPA);
210 4 : parse("NBINS",NBINS);
211 4 : parse("NBINZ",NBINZ);
212 4 : parse("MINS",MINS);
213 4 : parse("RCYL",RCYL);
214 4 : parse("SLOPE",slope);
215 4 : parse("ALPHA",ALPHA);
216 : //Reading compulsory arguments
217 4 : parse("MAXS",MAXS);
218 4 : parse("ZCC",ZCC);
219 :
220 :
221 4 : checkRead();
222 :
223 4 : log.printf(" External potential from file %s\n",file.c_str());
224 4 : log.printf(" Multiplied by %lf\n",scale_);
225 4 : if(spline) {
226 4 : log.printf(" External potential uses spline interpolation\n");
227 : }
228 4 : if(sparsegrid) {
229 0 : log.printf(" External potential uses sparse grid\n");
230 : }
231 :
232 : // Non più necessario dalla 2.3
233 : // addComponent("bias"); componentIsNotPeriodic("bias");
234 :
235 4 : std::string funcl=getLabel() + ".bias";
236 :
237 : // int size = plumed.comm.Get_size();
238 : // int rank = plumed.comm.Get_rank();
239 4 : IFile I_hate_this;
240 4 : bool do_exist=I_hate_this.FileExist(file);
241 :
242 4 : if(walkers_mpi) {
243 2 : if(comm.Get_rank()==0 && multi_sim_comm.Get_rank()==0) {
244 1 : if(!do_exist) {
245 1 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
246 : }
247 : }
248 2 : multi_sim_comm.Barrier();
249 : } else {
250 2 : if(comm.Get_rank()==0) {
251 2 : if(!do_exist) {
252 2 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
253 : }
254 : }
255 : }
256 :
257 : /*
258 : if(comm.Get_rank()==0){
259 : if(multi_sim_comm.Get_rank()==0 && walkers_mpi){
260 : if(!do_exist){
261 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
262 : }
263 : } else {
264 : if(!do_exist){
265 : createBIAS(RCYL, ZCC, ALPHA, KAPPA, MINS, MAXS, NBINS, NBINZ, safety, sphere, slope, funcl, file);
266 : }
267 : }
268 : if(walkers_mpi) multi_sim_comm.Barrier();
269 : }
270 : */
271 4 : comm.Barrier();
272 :
273 : // read grid
274 4 : IFile gridfile;
275 4 : gridfile.open(file);
276 8 : BiasGrid_=Grid::create(funcl,getArguments(),gridfile,sparsegrid,spline,true);
277 : //not necessary anymore? gridfile.close();
278 4 : if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
279 12 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
280 16 : if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
281 : }
282 4 : comm.Barrier();
283 4 : if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier();
284 8 : log<<" Bibliography "<<plumed.cite("Limongelli, Bonomi, and Parrinello, PNAS 110, 6358 (2013)")<<"\n";
285 8 : }
286 :
287 :
288 3 : void Funnel::createBIAS(const double& R_cyl, const double& z_cc, const double& alpha,
289 : const double& KAPPA, const double& MIN_S, const double& MAX_S, const double& NBIN_S,
290 : const double& NBIN_Z, const double& safety, const bool& sphere, const double& slope,
291 : const string& funcl, const string& file) {
292 : //R_cyl and z_cc forms the parameters of the cylinder.
293 : //alpha defines the angle in degrees.
294 :
295 : //PARAMETERS OF THE CONE
296 3 : double tg_alpha= tan(alpha);
297 :
298 : //parameters for PROGRESSION
299 : //parameters for DISTANCE
300 : double MIN_Z=0;
301 : double MAX_Z;
302 3 : if (sphere==false) {
303 2 : MAX_Z=R_cyl + tg_alpha * (z_cc - MIN_S);
304 : }
305 : else {
306 1 : MAX_Z=z_cc+safety;
307 : }
308 :
309 : //bin size
310 3 : double DX_Z = (MAX_Z - MIN_Z) / NBIN_Z;
311 : double DX_S;
312 3 : if (sphere==false) {
313 2 : DX_S=(MAX_S - MIN_S) / NBIN_S;
314 : }
315 : else {
316 1 : DX_S=(MAX_S + z_cc + safety)/NBIN_S;
317 : }
318 :
319 : double SS, Zmax, ZZ, D, d;
320 : double POT, FZ, FS;
321 :
322 3 : PLMD::OFile pof;
323 3 : pof.open(file);
324 :
325 : //Write the header
326 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());
327 3 : if (sphere==false) pof.printf("#! SET min_%s %f\n", getPntrToArgument(0)->getName().c_str(), MIN_S);
328 1 : else pof.printf("#! SET min_%s %f\n", getPntrToArgument(0)->getName().c_str(), -z_cc-safety);
329 3 : pof.printf("#! SET max_%s %f\n", getPntrToArgument(0)->getName().c_str(), MAX_S);
330 3 : pof.printf("#! SET nbins_%s %f\n", getPntrToArgument(0)->getName().c_str(), NBIN_S);
331 3 : pof.printf("#! SET periodic_%s false\n", getPntrToArgument(0)->getName().c_str());
332 3 : pof.printf("#! SET min_%s %f\n", getPntrToArgument(1)->getName().c_str(), MIN_Z);
333 3 : pof.printf("#! SET max_%s %f\n", getPntrToArgument(1)->getName().c_str(), MAX_Z);
334 3 : pof.printf("#! SET nbins_%s %f\n", getPntrToArgument(1)->getName().c_str(), NBIN_Z);
335 3 : pof.printf("#! SET periodic_%s false\n", getPntrToArgument(1)->getName().c_str());
336 :
337 : //Calculate and write the GRID
338 : //Cone or cylinder?
339 :
340 1506 : for(int is=0; is <= NBIN_S; is++) {
341 1503 : if (sphere==false) {
342 1002 : SS = MIN_S + is * DX_S;
343 : }
344 : else {
345 501 : SS = - z_cc - safety + is * DX_S;
346 : }
347 : bool cone = false;
348 1503 : if (sphere==false) {
349 1002 : if(SS <= z_cc) cone = true;
350 : }
351 : else {
352 501 : if (SS <= sqrt(pow(z_cc,2)-pow(R_cyl,2))) cone = true;
353 : }
354 : //Set wall boundaries properly
355 : if(cone == true) {
356 902 : if(sphere==false) {
357 548 : Zmax = R_cyl + (z_cc - SS) * tg_alpha;
358 : }
359 : else {
360 354 : if (SS > -z_cc) {
361 277 : Zmax = sqrt(pow(z_cc,2) - pow(SS,2));
362 : }
363 : else {
364 : Zmax = 0;
365 : }
366 : }
367 : }
368 601 : else Zmax = R_cyl;
369 :
370 754506 : for(int iz=0; iz <= NBIN_Z; iz++) {
371 753003 : ZZ = MIN_Z + iz * DX_Z;
372 :
373 : //Inside or outside?
374 : bool inside;
375 753003 : if(ZZ < Zmax) inside = true;
376 : else inside = false;
377 :
378 : if(inside == true) {
379 : POT = 0;
380 : FS = 0;
381 : FZ = 0;
382 : }
383 : else {
384 518149 : if(cone == true) {
385 235130 : if(sphere==false) {
386 127824 : POT = 0.5 * KAPPA * (ZZ - Zmax) * (ZZ - Zmax);
387 127824 : FZ = - KAPPA * (ZZ - Zmax);
388 127824 : FS = - KAPPA * (ZZ - Zmax) * tg_alpha;
389 : }
390 : else {
391 107306 : D = sqrt(pow(ZZ,2)+pow(SS,2));
392 107306 : d = D - z_cc;
393 107306 : POT = 0.5 * KAPPA * pow(d,2);
394 107306 : FZ = - KAPPA * d * ZZ / D;
395 107306 : FS = - KAPPA * d * SS / D;
396 : }
397 : }
398 : else {
399 283019 : if(sphere==false) {
400 212018 : POT = 0.5 * KAPPA * (ZZ - Zmax) * (ZZ - Zmax);
401 212018 : FZ = - KAPPA * (ZZ - Zmax);
402 : FS = 0;
403 : }
404 : else {
405 71001 : D = sqrt(pow(ZZ,2)+pow(SS,2));
406 71001 : d = D - z_cc;
407 71001 : if(ZZ>=R_cyl+slope*(SS-z_cc)) {
408 45984 : POT = 0.5 * KAPPA * pow(d,2);
409 45984 : FZ = - KAPPA * d * ZZ / D;
410 45984 : FS = - KAPPA * d * SS / D;
411 : }
412 : else {
413 25017 : POT = 0.5 * KAPPA * pow(sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2))-
414 : z_cc,2);
415 25017 : FZ = - KAPPA*(sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2))-z_cc)*
416 25017 : ZZ/sqrt(pow((ZZ+slope*z_cc-R_cyl)/slope,2)+pow(ZZ,2));
417 : FS = 0;
418 : }
419 : }
420 : }
421 : }
422 753003 : pof.printf("%13.6lf %13.6lf %13.6lf %13.6lf %13.6lf\n", SS, ZZ, POT, FS, FZ);
423 : }
424 1503 : pof.printf("\n");
425 : }
426 3 : pof.close();
427 3 : }
428 :
429 120 : void Funnel::calculate()
430 : {
431 120 : unsigned ncv=getNumberOfArguments();
432 120 : vector<double> cv(ncv), der(ncv);
433 :
434 360 : for(unsigned i=0; i<ncv; ++i) {
435 240 : cv[i]=getArgument(i);
436 : }
437 :
438 : // log.printf(" In Funnel: %13.6lf %13.6lf\n", cv[0], cv[1]);
439 :
440 120 : double ene=scale_*BiasGrid_->getValueAndDerivatives(cv,der);
441 :
442 120 : setBias(ene);
443 :
444 : // set Forces
445 360 : for(unsigned i=0; i<ncv; ++i) {
446 240 : const double f=-scale_*der[i];
447 240 : setOutputForce(i,f);
448 : }
449 120 : }
450 :
451 : }
452 : }
|
