Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2023 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 generic {
30 :
31 : //+PLUMEDOC PRINTANALYSIS COMMITTOR
32 : /*
33 : Does a committor analysis.
34 :
35 : The following input monitors two torsional angles during a simulation,
36 : defines two basins (A and B) as a function of the two torsion angles and
37 : stops the simulation when it falls in one of the two. The log
38 : file shows the latest values for the CVs and the basin reached.
39 :
40 : ```plumed
41 : TORSION ATOMS=1,2,3,4 LABEL=r1
42 : TORSION ATOMS=2,3,4,5 LABEL=r2
43 : COMMITTOR ...
44 : ARG=r1,r2
45 : STRIDE=10
46 : BASIN_LL1=0.15,0.20
47 : BASIN_UL1=0.25,0.40
48 : BASIN_LL2=-0.25,-0.40
49 : BASIN_UL2=-0.15,-0.20
50 : ... COMMITTOR
51 : ```
52 :
53 : */
54 : //+ENDPLUMEDOC
55 :
56 : class Committor :
57 : public ActionPilot,
58 : public ActionWithArguments {
59 : private:
60 : std::string file;
61 : OFile ofile;
62 : std::string fmt;
63 : std::vector< std::vector<double> > lowerlimits;
64 : std::vector< std::vector<double> > upperlimits;
65 : unsigned nbasins;
66 : unsigned basin;
67 : bool doNotStop;
68 : public:
69 : static void registerKeywords( Keywords& keys );
70 : explicit Committor(const ActionOptions&ao);
71 : void calculate() override;
72 30 : void apply() override {}
73 : };
74 :
75 : PLUMED_REGISTER_ACTION(Committor,"COMMITTOR")
76 :
77 13 : void Committor::registerKeywords( Keywords& keys ) {
78 13 : Action::registerKeywords(keys);
79 13 : ActionPilot::registerKeywords(keys);
80 13 : ActionWithArguments::registerKeywords(keys);
81 26 : keys.addInputKeyword("compulsory","ARG","scalar","the labels of the values which is being used to define the committor surface");
82 13 : keys.add("numbered", "BASIN_LL","List of lower limits for basin #");
83 13 : keys.add("numbered", "BASIN_UL","List of upper limits for basin #");
84 26 : keys.reset_style("BASIN_LL","compulsory");
85 26 : keys.reset_style("BASIN_UL","compulsory");
86 13 : keys.add("compulsory","STRIDE","1","the frequency with which the CVs are analyzed");
87 13 : keys.add("optional","FILE","the name of the file on which to output the reached basin");
88 13 : keys.add("optional","FMT","the format that should be used to output real numbers");
89 13 : keys.addFlag("NOSTOP",false,"if true do not stop the simulation when reaching a basin but just keep track of it");
90 13 : }
91 :
92 9 : Committor::Committor(const ActionOptions&ao):
93 : Action(ao),
94 : ActionPilot(ao),
95 : ActionWithArguments(ao),
96 9 : fmt("%f"),
97 9 : nbasins(0),
98 9 : basin(0),
99 9 : doNotStop(false) {
100 9 : ofile.link(*this);
101 18 : parse("FILE",file);
102 9 : if(file.length()>0) {
103 2 : ofile.open(file);
104 2 : log.printf(" on file %s\n",file.c_str());
105 : } else {
106 7 : log.printf(" on plumed log file\n");
107 7 : ofile.link(log);
108 : }
109 9 : parse("FMT",fmt);
110 9 : fmt=" "+fmt;
111 9 : log.printf(" with format %s\n",fmt.c_str());
112 :
113 13 : for(unsigned b=1;; ++b ) {
114 22 : std::vector<double> tmpl, tmpu;
115 22 : parseNumberedVector("BASIN_LL", b, tmpl );
116 44 : parseNumberedVector("BASIN_UL", b, tmpu );
117 22 : if( tmpl.empty() && tmpu.empty() ) {
118 : break;
119 : }
120 13 : if( tmpl.size()!=getNumberOfArguments()) {
121 0 : error("Wrong number of values for BASIN_LL: they should be equal to the number of arguments");
122 : }
123 13 : if( tmpu.size()!=getNumberOfArguments()) {
124 0 : error("Wrong number of values for BASIN_UL: they should be equal to the number of arguments");
125 : }
126 13 : lowerlimits.push_back(tmpl);
127 13 : upperlimits.push_back(tmpu);
128 13 : nbasins=b;
129 22 : }
130 :
131 9 : parseFlag("NOSTOP", doNotStop);
132 :
133 9 : checkRead();
134 9 : if(nbasins==0) {
135 0 : error("no basin defined: use BASIN_LL1 and BASIN_UL1");
136 : }
137 :
138 :
139 22 : for(unsigned b=0; b<nbasins; b++) {
140 13 : log.printf(" BASIN %u definition:\n", b+1);
141 32 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
142 19 : if(lowerlimits[b][i]>upperlimits[b][i]) {
143 0 : error("COMMITTOR: UPPER bounds must always be greater than LOWER bounds");
144 : }
145 19 : log.printf(" %f - %f\n", lowerlimits[b][i], upperlimits[b][i]);
146 : }
147 13 : if(doNotStop) {
148 3 : log.printf(" COMMITOR will keep track of the visited basins without stopping the simulations\n");
149 : }
150 : }
151 :
152 20 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
153 11 : ofile.setupPrintValue( getPntrToArgument(i) );
154 : }
155 9 : }
156 :
157 30 : void Committor::calculate() {
158 30 : std::vector<unsigned> inbasin;
159 30 : inbasin.assign (nbasins,1);
160 :
161 : // check if current configuration belongs to a basin
162 106 : for(unsigned b=0; b<nbasins; ++b) {
163 221 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
164 145 : if(getArgument(i)>lowerlimits[b][i]&&getArgument(i)<upperlimits[b][i]) {
165 : inbasin[b]*=1;
166 : } else {
167 80 : inbasin[b]*=0;
168 : }
169 : }
170 : }
171 :
172 : // check in which basin we are if any and if this is the same or a new one
173 : bool inonebasin = false;
174 71 : for(unsigned b=0; b<nbasins; ++b) {
175 61 : if(inbasin[b]==1) {
176 20 : if(basin!=(b+1)) {
177 15 : basin = b+1;
178 15 : ofile.fmtField(" %f");
179 15 : ofile.printField("time",getTime());
180 38 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
181 23 : ofile.fmtField(fmt);
182 23 : ofile.printField( getPntrToArgument(i), getArgument(i) );
183 : }
184 15 : ofile.printField("basin", static_cast<int> (b+1));
185 15 : ofile.printField();
186 : }
187 : inonebasin = true;
188 : break;
189 : }
190 : }
191 : if(!inonebasin) {
192 10 : basin = 0;
193 : }
194 :
195 : // then check if the simulation should be stopped
196 20 : if(inonebasin&&(!doNotStop)) {
197 : std::string num;
198 8 : Tools::convert( basin, num );
199 8 : std::string str = "COMMITTED TO BASIN " + num;
200 8 : ofile.addConstantField(str);
201 8 : ofile.printField();
202 8 : ofile.flush();
203 8 : plumed.stop();
204 : }
205 30 : }
206 :
207 : }
208 : }
|
