Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "core/ActionPilot.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/PlumedMain.h"
27 :
28 : namespace PLMD {
29 : namespace analysis {
30 :
31 : //+PLUMEDOC PRINTANALYSIS COMMITTOR
32 : /*
33 : Does a committor analysis.
34 :
35 : \par Examples
36 :
37 : The following input monitors two torsional angles during a simulation,
38 : defines two basins (A and B) as a function of the two torsions and
39 : stops the simulation when it falls in one of the two. In the log
40 : file will be shown the latest values for the CVs and the basin reached.
41 : \plumedfile
42 : TORSION ATOMS=1,2,3,4 LABEL=r1
43 : TORSION ATOMS=2,3,4,5 LABEL=r2
44 : COMMITTOR ...
45 : ARG=r1,r2
46 : STRIDE=10
47 : BASIN_LL1=0.15,0.20
48 : BASIN_UL1=0.25,0.40
49 : BASIN_LL2=-0.25,-0.40
50 : BASIN_UL2=-0.15,-0.20
51 : ... COMMITTOR
52 : \endplumedfile
53 :
54 : */
55 : //+ENDPLUMEDOC
56 :
57 8 : class Committor :
58 : public ActionPilot,
59 : public ActionWithArguments
60 : {
61 : private:
62 : std::string file;
63 : OFile ofile;
64 : std::string fmt;
65 : std::vector< std::vector<double> > lowerlimits;
66 : std::vector< std::vector<double> > upperlimits;
67 : unsigned nbasins;
68 : unsigned basin;
69 : bool doNotStop;
70 : public:
71 : static void registerKeywords( Keywords& keys );
72 : explicit Committor(const ActionOptions&ao);
73 : void calculate();
74 24 : void apply() {}
75 : };
76 :
77 6454 : PLUMED_REGISTER_ACTION(Committor,"COMMITTOR")
78 :
79 3 : void Committor::registerKeywords( Keywords& keys ) {
80 3 : Action::registerKeywords(keys);
81 3 : ActionPilot::registerKeywords(keys);
82 3 : ActionWithArguments::registerKeywords(keys);
83 6 : keys.use("ARG");
84 12 : keys.add("numbered", "BASIN_LL","List of lower limits for basin #");
85 12 : keys.add("numbered", "BASIN_UL","List of upper limits for basin #");
86 15 : keys.reset_style("BASIN_LL","compulsory"); keys.reset_style("BASIN_UL","compulsory");
87 15 : keys.add("compulsory","STRIDE","1","the frequency with which the CVs are analysed");
88 12 : keys.add("optional","FILE","the name of the file on which to output the reached basin");
89 12 : keys.add("optional","FMT","the format that should be used to output real numbers");
90 9 : keys.addFlag("NOSTOP",false,"if true do not stop the simulation when reaching a basin but just keep track of it");
91 3 : }
92 :
93 2 : Committor::Committor(const ActionOptions&ao):
94 : Action(ao),
95 : ActionPilot(ao),
96 : ActionWithArguments(ao),
97 : fmt("%f"),
98 : basin(0),
99 8 : doNotStop(false)
100 : {
101 2 : ofile.link(*this);
102 4 : parse("FILE",file);
103 2 : if(file.length()>0) {
104 2 : ofile.open(file);
105 4 : log.printf(" on file %s\n",file.c_str());
106 : } else {
107 0 : log.printf(" on plumed log file\n");
108 0 : ofile.link(log);
109 : }
110 4 : parse("FMT",fmt);
111 4 : fmt=" "+fmt;
112 4 : log.printf(" with format %s\n",fmt.c_str());
113 :
114 6 : for(unsigned b=1;; ++b ) {
115 : std::vector<double> tmpl, tmpu;
116 16 : parseNumberedVector("BASIN_LL", b, tmpl );
117 16 : parseNumberedVector("BASIN_UL", b, tmpu );
118 10 : if( tmpl.empty() && tmpu.empty() ) break;
119 6 : if( tmpl.size()!=getNumberOfArguments()) error("Wrong number of values for BASIN_LL: they should be equal to the number of arguments");
120 6 : if( tmpu.size()!=getNumberOfArguments()) error("Wrong number of values for BASIN_UL: they should be equal to the number of arguments");
121 6 : lowerlimits.push_back(tmpl);
122 6 : upperlimits.push_back(tmpu);
123 6 : nbasins=b;
124 6 : }
125 :
126 4 : parseFlag("NOSTOP", doNotStop);
127 :
128 2 : checkRead();
129 :
130 :
131 8 : for(unsigned b=0; b<nbasins; b++) {
132 6 : log.printf(" BASIN %u definition:\n", b+1);
133 30 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
134 36 : if(lowerlimits[b][i]>upperlimits[b][i]) error("COMMITTOR: UPPER bounds must always be greater than LOWER bounds");
135 36 : log.printf(" %f - %f\n", lowerlimits[b][i], upperlimits[b][i]);
136 : }
137 6 : if(doNotStop) log.printf(" COMMITOR will keep track of the visited basins without stopping the simulations\n");
138 : }
139 :
140 10 : for(unsigned i=0; i<getNumberOfArguments(); ++i) ofile.setupPrintValue( getPntrToArgument(i) );
141 2 : }
142 :
143 24 : void Committor::calculate() {
144 : std::vector<unsigned> inbasin;
145 48 : inbasin.assign (nbasins,1);
146 :
147 : // check if current configuration belongs to a basin
148 168 : for(unsigned b=0; b<nbasins; ++b) {
149 360 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
150 372 : if(getArgument(i)>lowerlimits[b][i]&&getArgument(i)<upperlimits[b][i]) {
151 : inbasin[b]*=1;
152 : } else {
153 86 : inbasin[b]*=0;
154 : }
155 : }
156 : }
157 :
158 : // check in which basin we are if any and if this is the same or a new one
159 : bool inonebasin = false;
160 112 : for(unsigned b=0; b<nbasins; ++b) {
161 114 : if(inbasin[b]==1) {
162 13 : if(basin!=(b+1)) {
163 8 : basin = b+1;
164 16 : ofile.fmtField(" %f");
165 16 : ofile.printField("time",getTime());
166 40 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
167 16 : ofile.fmtField(fmt);
168 32 : ofile.printField( getPntrToArgument(i), getArgument(i) );
169 : }
170 16 : ofile.printField("basin", static_cast<int> (b+1));
171 8 : ofile.printField();
172 : }
173 : inonebasin = true;
174 : break;
175 : }
176 : }
177 11 : if(!inonebasin) basin = 0;
178 :
179 : // then check if the simulation should be stopped
180 24 : if(inonebasin&&(!doNotStop)) {
181 1 : std::string num; Tools::convert( basin, num );
182 1 : std::string str = "COMMITED TO BASIN " + num;
183 1 : ofile.addConstantField(str);
184 1 : ofile.printField();
185 1 : ofile.flush();
186 1 : plumed.stop();
187 : }
188 24 : }
189 :
190 : }
191 4839 : }
|
