Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 : #include "core/ActionPilot.h"
23 : #include "core/ActionWithValue.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "tools/IFile.h"
28 : #include <memory>
29 :
30 : namespace PLMD {
31 : namespace generic {
32 :
33 : //+PLUMEDOC GENERIC READ
34 : /*
35 : Read quantities from a colvar file.
36 :
37 : This Action can be used with driver to read in a colvar file that was generated during
38 : an MD simulation. The following example shows how this works.
39 :
40 : ```plumed
41 : #SETTINGS INPUTFILES=regtest/basic/rt-fametad-1/Input-COLVAR
42 : sum_abs: READ VALUES=sum_abs FILE=regtest/basic/rt-fametad-1/Input-COLVAR IGNORE_FORCES
43 : PRINT ARG=sum_abs STRIDE=1 FILE=colvar
44 : ```
45 :
46 : The input file `Input-COLVAR` is a colvar file that was generated using a [PRINT](PRINT.md)
47 : command. The instruction `VALUES=sum_abs` tells PLUMED that we want to read the column headed
48 : `sum_abs` from that file. The value outputted from a read command is thus always a scalar.
49 :
50 : The `IGNORE_FORCES` flag tells PLUMED that any forces that are applied on
51 : this action can be safely ignored. If you try to run a simulation in which a bias acts upon a
52 : quantity that is read in from a file using a READ command PLUMED will crash with an error.
53 :
54 : ## Dealing with components
55 :
56 : The following example provides another example where the READ command is used:
57 :
58 : ```plumed
59 : #SETTINGS INPUTFILES=regtest/basic/rt41/input_colvar
60 : r1: READ VALUES=p2.X FILE=regtest/basic/rt41/input_colvar
61 : r2: READ VALUES=p3.* FILE=regtest/basic/rt41/input_colvar
62 : PRINT ARG=r1.X,r2.* FILE=colvar
63 : ```
64 :
65 : Notice that the READ command preseves the part of the name that appears after the ".". Consequently,
66 : when we read the value `p2.X` from the input file the output value the action with label `r1` calls the output
67 : value that contains this information `r1.X`. The READ command is implemented this way so that you can use the
68 : wildcard syntax that is illustrated in the command labelled `r2` in the above input.
69 :
70 : ## Reading COLVAR files and trajectories
71 :
72 : The example input below indicates a case where a trajectory is being post processed and where some COLVAR files
73 : that were generated when the MD simulation was run are read in.
74 :
75 : ```plumed
76 : #SETTINGS INPUTFILES=regtest/basic/rt19/input_colvar2
77 : # The distance here is being calculated from the trajectory
78 : d: DISTANCE ATOMS=1,2
79 : # This CV is being read in from a file that was output with the same frequency as frames
80 : # were output from the trajectory. Notice that you can read from zip files
81 : r1: READ VALUES=rmsd0 FILE=regtest/basic/rt19/input_colvar.gz
82 : # This CV is being read in from a file that was output with twice as frequency as frames
83 : r2: READ VALUES=rmsd0 FILE=regtest/basic/rt19/input_colvar2 EVERY=2 IGNORE_TIME
84 : # This outputs our three quantities
85 : PRINT ARG=d,r1,r2 FILE=colvar
86 : ```
87 :
88 : The driver command to run this script as follows:
89 :
90 : ````
91 : plumed driver --plumed plumed.dat --trajectory-stride 10 --timestep 0.005 --ixyz trajectory.xyz
92 : ````
93 :
94 : When you are doing analyses like these, which involve mixing using READ command and analysing a trajectory, PLUMED forces you
95 : to take care to ensure that everything in the input file was generated from the same step in the input trajectory. You must
96 : use the `--trajectory-stride` and `--timestep` commands when using driver above so PLUMED can correctly calculate the simulation
97 : time and compare it with the time stamps in any colvar files that are being read in using READ commands.
98 :
99 : If you want to turn off these checks either because you are confident that you have set things up correctly or if you are not
100 : mixing variables that have been calculated from a trajectory with variables that are being read from a file you can use the
101 : `IGNORE_TIME` flag. Notice also that you can use the `EVERY` flag to tell PLUMED to ignore parts of the COLVAR file if data
102 : has been output to that file more frequently that data has been output to the trajectory.
103 :
104 : */
105 : //+ENDPLUMEDOC
106 :
107 : class Read :
108 : public ActionPilot,
109 : public ActionWithValue {
110 : private:
111 : bool ignore_time;
112 : bool ignore_forces;
113 : bool cloned_file;
114 : unsigned nlinesPerStep;
115 : std::string filename;
116 : /// Unique pointer with the same scope as ifile.
117 : std::unique_ptr<IFile> ifile_ptr;
118 : /// Pointer to input file.
119 : /// It is either pointing to the content of ifile_ptr
120 : /// or to the file it is cloned from.
121 : IFile* ifile;
122 : std::vector<std::unique_ptr<Value>> readvals;
123 : public:
124 : static void registerKeywords( Keywords& keys );
125 : explicit Read(const ActionOptions&);
126 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
127 : void prepare() override;
128 921429 : void apply() override {}
129 : void calculate() override;
130 : void update() override;
131 : std::string getFilename() const;
132 : IFile* getFile();
133 : unsigned getNumberOfDerivatives() override;
134 : void turnOnDerivatives() override;
135 : };
136 :
137 : PLUMED_REGISTER_ACTION(Read,"READ")
138 :
139 95 : void Read::registerKeywords(Keywords& keys) {
140 95 : Action::registerKeywords(keys);
141 95 : ActionPilot::registerKeywords(keys);
142 95 : ActionWithValue::registerKeywords(keys);
143 95 : keys.add("compulsory","STRIDE","1","the frequency with which the file should be read.");
144 95 : keys.add("compulsory","EVERY","1","only read every nth line of the colvar file. This should be used if the colvar was written more frequently than the trajectory.");
145 95 : keys.add("compulsory","VALUES","the values to read from the file");
146 95 : keys.add("compulsory","FILE","the name of the file from which to read these quantities");
147 95 : keys.addFlag("IGNORE_TIME",false,"ignore the time in the colvar file. When this flag is not present read will be quite strict "
148 : "about the start time of the simulation and the stride between frames");
149 95 : keys.addFlag("IGNORE_FORCES",false,"use this flag if the forces added by any bias can be safely ignored. As an example forces can be "
150 : "safely ignored if you are doing post processing that does not involve outputting forces");
151 190 : keys.remove("NUMERICAL_DERIVATIVES");
152 95 : keys.use("UPDATE_FROM");
153 95 : keys.use("UPDATE_UNTIL");
154 95 : ActionWithValue::useCustomisableComponents(keys);
155 95 : }
156 :
157 91 : Read::Read(const ActionOptions&ao):
158 : Action(ao),
159 : ActionPilot(ao),
160 : ActionWithValue(ao),
161 91 : ignore_time(false),
162 91 : ignore_forces(false),
163 91 : nlinesPerStep(1) {
164 : // Read the file name from the input line
165 91 : parse("FILE",filename);
166 : // Check if time is to be ignored
167 91 : parseFlag("IGNORE_TIME",ignore_time);
168 : // Check if forces are to be ignored
169 91 : parseFlag("IGNORE_FORCES",ignore_forces);
170 : // Open the file if it is not already opened
171 91 : cloned_file=false;
172 91 : std::vector<Read*> other_reads=plumed.getActionSet().select<Read*>();
173 141 : for(unsigned i=0; i<other_reads.size(); ++i) {
174 50 : if( other_reads[i]->getFilename()==filename ) {
175 47 : ifile=other_reads[i]->getFile();
176 47 : cloned_file=true;
177 : }
178 : }
179 91 : if( !cloned_file ) {
180 63 : ifile_ptr=Tools::make_unique<IFile>();
181 63 : ifile=ifile_ptr.get();
182 63 : if( !ifile->FileExist(filename) ) {
183 0 : error("could not find file named " + filename);
184 : }
185 63 : ifile->link(*this);
186 63 : ifile->open(filename);
187 63 : ifile->allowIgnoredFields();
188 : }
189 91 : parse("EVERY",nlinesPerStep);
190 91 : if(nlinesPerStep>1) {
191 2 : log.printf(" only reading every %uth line of file %s\n",nlinesPerStep,filename.c_str() );
192 : } else {
193 89 : log.printf(" reading data from file %s\n",filename.c_str() );
194 : }
195 : // Find out what we are reading
196 : std::vector<std::string> valread;
197 91 : parseVector("VALUES",valread);
198 :
199 91 : if(nlinesPerStep>1 && cloned_file) {
200 0 : error("Opening a file multiple times and using EVERY is not allowed");
201 : }
202 :
203 : std::size_t dot=valread[0].find_first_of('.');
204 91 : if( valread[0].find(".")!=std::string::npos ) {
205 11 : std::string label=valread[0].substr(0,dot);
206 11 : std::string name=valread[0].substr(dot+1);
207 11 : if( name=="*" ) {
208 2 : if( valread.size()>1 ) {
209 0 : error("all values must be from the same Action when using READ");
210 : }
211 : std::vector<std::string> fieldnames;
212 2 : ifile->scanFieldList( fieldnames );
213 16 : for(unsigned i=0; i<fieldnames.size(); ++i) {
214 14 : if( fieldnames[i].substr(0,dot)==label ) {
215 6 : readvals.emplace_back(Tools::make_unique<Value>(this, fieldnames[i], false) );
216 12 : addComponentWithDerivatives( fieldnames[i].substr(dot+1) );
217 12 : if( ifile->FieldExist("min_" + fieldnames[i]) ) {
218 0 : componentIsPeriodic( fieldnames[i].substr(dot+1), "-pi","pi" );
219 : } else {
220 12 : componentIsNotPeriodic( fieldnames[i].substr(dot+1) );
221 : }
222 : }
223 : }
224 2 : } else {
225 9 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) );
226 9 : addComponentWithDerivatives( name );
227 18 : if( ifile->FieldExist("min_" + valread[0]) ) {
228 0 : componentIsPeriodic( valread[0].substr(dot+1), "-pi", "pi" );
229 : } else {
230 18 : componentIsNotPeriodic( valread[0].substr(dot+1) );
231 : }
232 12 : for(unsigned i=1; i<valread.size(); ++i) {
233 6 : if( valread[i].substr(0,dot)!=label ) {
234 0 : error("all values must be from the same Action when using READ");
235 : };
236 3 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[i], false) );
237 6 : addComponentWithDerivatives( valread[i].substr(dot+1) );
238 6 : if( ifile->FieldExist("min_" + valread[i]) ) {
239 0 : componentIsPeriodic( valread[i].substr(dot+1), "-pi", "pi" );
240 : } else {
241 6 : componentIsNotPeriodic( valread[i].substr(dot+1) );
242 : }
243 : }
244 : }
245 : } else {
246 80 : if( valread.size()!=1 ) {
247 0 : error("all values must be from the same Action when using READ");
248 : }
249 80 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) );
250 80 : addValueWithDerivatives();
251 160 : if( ifile->FieldExist("min_" + valread[0]) ) {
252 18 : setPeriodic( "-pi", "pi" );
253 : } else {
254 71 : setNotPeriodic();
255 : }
256 80 : log.printf(" reading value %s and storing as %s\n",valread[0].c_str(),getLabel().c_str() );
257 : }
258 91 : checkRead();
259 182 : }
260 :
261 4 : std::string Read::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
262 4 : plumed_assert( !exists( getLabel() ) );
263 7 : for(unsigned i=0; i<readvals.size(); ++i) {
264 7 : if( readvals[i]->getName().find( cname )!=std::string::npos ) {
265 8 : return "values from the column labelled " + readvals[i]->getName() + " in the file named " + filename;
266 : }
267 : }
268 0 : plumed_error();
269 : return "";
270 : }
271 :
272 50 : std::string Read::getFilename() const {
273 50 : return filename;
274 : }
275 :
276 47 : IFile* Read::getFile() {
277 47 : return ifile;
278 : }
279 :
280 0 : unsigned Read::getNumberOfDerivatives() {
281 0 : return 0;
282 : }
283 :
284 26 : void Read::turnOnDerivatives() {
285 26 : if( !ignore_forces )
286 0 : error("cannot calculate derivatives for colvars that are read in from a file. If you are postprocessing and "
287 : "these forces do not matter add the flag IGNORE_FORCES to all READ actions");
288 26 : }
289 :
290 921429 : void Read::prepare() {
291 921429 : if( !cloned_file ) {
292 : double du_time;
293 403414 : if( !ifile->scanField("time",du_time) ) {
294 0 : error("Reached end of file " + filename + " before end of trajectory");
295 201707 : } else if( std::abs( du_time-getTime() )>getTimeStep() && !ignore_time ) {
296 : std::string str_dutime,str_ptime;
297 0 : Tools::convert(du_time,str_dutime);
298 0 : Tools::convert(getTime(),str_ptime);
299 0 : error("mismatched times in colvar files : colvar time=" + str_dutime + " plumed time=" + str_ptime + ". Add IGNORE_TIME to ignore error.");
300 : }
301 : }
302 921429 : }
303 :
304 921429 : void Read::calculate() {
305 : std::string smin, smax;
306 2019293 : for(unsigned i=0; i<readvals.size(); ++i) {
307 : // .get returns the raw pointer
308 : // ->get calls the Value::get() method
309 1097864 : ifile->scanField( readvals[i].get() );
310 1097864 : getPntrToComponent(i)->set( readvals[i]->get() );
311 1097864 : if( readvals[i]->isPeriodic() ) {
312 3309 : readvals[i]->getDomain( smin, smax );
313 3309 : getPntrToComponent(i)->setDomain( smin, smax );
314 : }
315 : }
316 921429 : }
317 :
318 921429 : void Read::update() {
319 921429 : if( !cloned_file ) {
320 403965 : for(unsigned i=0; i<nlinesPerStep; ++i) {
321 202258 : ifile->scanField();
322 : double du_time;
323 404516 : if( !ifile->scanField("time",du_time) && !plumed.inputsAreActive() ) {
324 54 : plumed.stop();
325 : }
326 : }
327 : }
328 921429 : }
329 :
330 : }
331 : }
|
