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