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 : # We start reading this CV from this file after we have read the first 100 ps of the trajectory
85 : # and stop reading from it after we have read the first 800 ps of the trajectory
86 : r3: READ VALUES=rmsd0 FILE=regtest/basic/rt19/input_colvar2 UPDATE_FROM=100 UPDATE_UNTIL=800 IGNORE_TIME
87 : # This outputs our the three quantities that are determined for every step
88 : PRINT ARG=d,r1,r2 FILE=colvar
89 : # This outputs the data that we only have for the 700 ps of the trajectory after the first 100 ps
90 : PRINT ARG=d,r1,r2,r3 FILE=colvar2 UPDATE_FROM=100 UPDATE_UNTIL=800
91 : ```
92 :
93 : The driver command to run this script as follows:
94 :
95 : ````
96 : plumed driver --plumed plumed.dat --trajectory-stride 10 --timestep 0.005 --ixyz trajectory.xyz
97 : ````
98 :
99 : When you are doing analyses like these, which involve mixing using READ command and analysing a trajectory, PLUMED forces you
100 : to take care to ensure that everything in the input file was generated from the same step in the input trajectory. You must
101 : use the `--trajectory-stride` and `--timestep` commands when using driver above so PLUMED can correctly calculate the simulation
102 : time and compare it with the time stamps in any colvar files that are being read in using READ commands.
103 :
104 : 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
105 : mixing variables that have been calculated from a trajectory with variables that are being read from a file you can use the
106 : `IGNORE_TIME` flag. Notice also that you can use the `EVERY` flag to tell PLUMED to ignore parts of the COLVAR file if data
107 : has been output to that file more frequently that data has been output to the trajectory.
108 :
109 : */
110 : //+ENDPLUMEDOC
111 :
112 : class Read :
113 : public ActionPilot,
114 : public ActionWithValue {
115 : private:
116 : bool ignore_time;
117 : bool ignore_forces;
118 : bool cloned_file;
119 : unsigned nlinesPerStep;
120 : std::string filename;
121 : /// Unique pointer with the same scope as ifile.
122 : std::unique_ptr<IFile> ifile_ptr;
123 : /// Pointer to input file.
124 : /// It is either pointing to the content of ifile_ptr
125 : /// or to the file it is cloned from.
126 : IFile* ifile;
127 : std::vector<std::unique_ptr<Value>> readvals;
128 : public:
129 : static void registerKeywords( Keywords& keys );
130 : explicit Read(const ActionOptions&);
131 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
132 : void prepare() override;
133 921430 : void apply() override {}
134 : void calculate() override;
135 : void update() override;
136 : std::string getFilename() const;
137 : IFile* getFile();
138 : unsigned getNumberOfDerivatives() override;
139 : void turnOnDerivatives() override;
140 : };
141 :
142 : PLUMED_REGISTER_ACTION(Read,"READ")
143 :
144 98 : void Read::registerKeywords(Keywords& keys) {
145 98 : Action::registerKeywords(keys);
146 98 : ActionPilot::registerKeywords(keys);
147 98 : ActionWithValue::registerKeywords(keys);
148 98 : keys.add("compulsory","STRIDE","1","the frequency with which the file should be read.");
149 98 : 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.");
150 98 : keys.add("compulsory","VALUES","the values to read from the file");
151 98 : keys.add("compulsory","FILE","the name of the file from which to read these quantities");
152 98 : keys.addFlag("IGNORE_TIME",false,"ignore the time in the colvar file. When this flag is not present read will be quite strict "
153 : "about the start time of the simulation and the stride between frames");
154 98 : 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 "
155 : "safely ignored if you are doing post processing that does not involve outputting forces");
156 196 : keys.remove("NUMERICAL_DERIVATIVES");
157 98 : keys.use("UPDATE_FROM");
158 98 : keys.use("UPDATE_UNTIL");
159 98 : ActionWithValue::useCustomisableComponents(keys);
160 98 : }
161 :
162 92 : Read::Read(const ActionOptions&ao):
163 : Action(ao),
164 : ActionPilot(ao),
165 : ActionWithValue(ao),
166 92 : ignore_time(false),
167 92 : ignore_forces(false),
168 92 : nlinesPerStep(1) {
169 : // Read the file name from the input line
170 92 : parse("FILE",filename);
171 : // Check if time is to be ignored
172 92 : parseFlag("IGNORE_TIME",ignore_time);
173 : // Check if forces are to be ignored
174 92 : parseFlag("IGNORE_FORCES",ignore_forces);
175 : // Open the file if it is not already opened
176 92 : cloned_file=false;
177 92 : std::vector<Read*> other_reads=plumed.getActionSet().select<Read*>();
178 142 : for(unsigned i=0; i<other_reads.size(); ++i) {
179 50 : if( other_reads[i]->getFilename()==filename ) {
180 47 : ifile=other_reads[i]->getFile();
181 47 : cloned_file=true;
182 : }
183 : }
184 92 : if( !cloned_file ) {
185 64 : ifile_ptr=Tools::make_unique<IFile>();
186 64 : ifile=ifile_ptr.get();
187 64 : if( !ifile->FileExist(filename) ) {
188 0 : error("could not find file named " + filename);
189 : }
190 64 : ifile->link(*this);
191 64 : ifile->open(filename);
192 64 : ifile->allowIgnoredFields();
193 : }
194 92 : parse("EVERY",nlinesPerStep);
195 92 : if(nlinesPerStep>1) {
196 2 : log.printf(" only reading every %uth line of file %s\n",nlinesPerStep,filename.c_str() );
197 : } else {
198 90 : log.printf(" reading data from file %s\n",filename.c_str() );
199 : }
200 : // Find out what we are reading
201 92 : std::vector<std::string> valread;
202 92 : parseVector("VALUES",valread);
203 :
204 92 : if(nlinesPerStep>1 && cloned_file) {
205 0 : error("Opening a file multiple times and using EVERY is not allowed");
206 : }
207 :
208 : std::size_t dot=valread[0].find_first_of('.');
209 92 : if( valread[0].find(".")!=std::string::npos ) {
210 11 : std::string labelVal=valread[0].substr(0,dot);
211 11 : std::string nameVal=valread[0].substr(dot+1);
212 11 : if( nameVal=="*" ) {
213 2 : if( valread.size()>1 ) {
214 0 : error("all values must be from the same Action when using READ");
215 : }
216 2 : std::vector<std::string> fieldnames;
217 2 : ifile->scanFieldList( fieldnames );
218 16 : for(unsigned i=0; i<fieldnames.size(); ++i) {
219 14 : if( fieldnames[i].substr(0,dot)==labelVal ) {
220 6 : readvals.emplace_back(Tools::make_unique<Value>(this, fieldnames[i], false) );
221 12 : addComponentWithDerivatives( fieldnames[i].substr(dot+1) );
222 12 : if( ifile->FieldExist("min_" + fieldnames[i]) ) {
223 0 : componentIsPeriodic( fieldnames[i].substr(dot+1), "-pi","pi" );
224 : } else {
225 12 : componentIsNotPeriodic( fieldnames[i].substr(dot+1) );
226 : }
227 : }
228 : }
229 2 : } else {
230 9 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) );
231 9 : addComponentWithDerivatives( nameVal );
232 18 : if( ifile->FieldExist("min_" + valread[0]) ) {
233 0 : componentIsPeriodic( valread[0].substr(dot+1), "-pi", "pi" );
234 : } else {
235 18 : componentIsNotPeriodic( valread[0].substr(dot+1) );
236 : }
237 12 : for(unsigned i=1; i<valread.size(); ++i) {
238 3 : if( valread[i].substr(0,dot)!=labelVal ) {
239 0 : error("all values must be from the same Action when using READ");
240 : };
241 3 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[i], false) );
242 6 : addComponentWithDerivatives( valread[i].substr(dot+1) );
243 6 : if( ifile->FieldExist("min_" + valread[i]) ) {
244 0 : componentIsPeriodic( valread[i].substr(dot+1), "-pi", "pi" );
245 : } else {
246 6 : componentIsNotPeriodic( valread[i].substr(dot+1) );
247 : }
248 : }
249 : }
250 : } else {
251 81 : if( valread.size()!=1 ) {
252 0 : error("all values must be from the same Action when using READ");
253 : }
254 81 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) );
255 81 : addValueWithDerivatives();
256 162 : if( ifile->FieldExist("min_" + valread[0]) ) {
257 18 : setPeriodic( "-pi", "pi" );
258 : } else {
259 72 : setNotPeriodic();
260 : }
261 81 : log.printf(" reading value %s and storing as %s\n",valread[0].c_str(),getLabel().c_str() );
262 : }
263 92 : checkRead();
264 184 : }
265 :
266 4 : std::string Read::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
267 4 : plumed_assert( !exists( getLabel() ) );
268 7 : for(unsigned i=0; i<readvals.size(); ++i) {
269 7 : if( readvals[i]->getName().find( cname )!=std::string::npos ) {
270 12 : return "values from the column labelled " + readvals[i]->getName() + " in the file named " + filename;
271 : }
272 : }
273 0 : plumed_error();
274 : return "";
275 : }
276 :
277 50 : std::string Read::getFilename() const {
278 50 : return filename;
279 : }
280 :
281 47 : IFile* Read::getFile() {
282 47 : return ifile;
283 : }
284 :
285 0 : unsigned Read::getNumberOfDerivatives() {
286 0 : return 0;
287 : }
288 :
289 26 : void Read::turnOnDerivatives() {
290 26 : if( !ignore_forces )
291 0 : error("cannot calculate derivatives for colvars that are read in from a file. If you are postprocessing and "
292 : "these forces do not matter add the flag IGNORE_FORCES to all READ actions");
293 26 : }
294 :
295 921430 : void Read::prepare() {
296 921430 : if( !cloned_file ) {
297 : double du_time;
298 403420 : if( !ifile->scanField("time",du_time) ) {
299 0 : error("Reached end of file " + filename + " before end of trajectory");
300 201710 : } else if( std::abs( du_time-getTime() )>getTimeStep() && !ignore_time ) {
301 : std::string str_dutime,str_ptime;
302 0 : Tools::convert(du_time,str_dutime);
303 0 : Tools::convert(getTime(),str_ptime);
304 0 : error("mismatched times in colvar files : colvar time=" + str_dutime + " plumed time=" + str_ptime + ". Add IGNORE_TIME to ignore error.");
305 : }
306 : }
307 921430 : }
308 :
309 921430 : void Read::calculate() {
310 : std::string smin, smax;
311 2019295 : for(unsigned i=0; i<readvals.size(); ++i) {
312 : // .get returns the raw pointer
313 : // ->get calls the Value::get() method
314 1097865 : ifile->scanField( readvals[i].get() );
315 1097865 : getPntrToComponent(i)->set( readvals[i]->get() );
316 1097865 : if( readvals[i]->isPeriodic() ) {
317 3309 : readvals[i]->getDomain( smin, smax );
318 3309 : getPntrToComponent(i)->setDomain( smin, smax );
319 : }
320 : }
321 921430 : }
322 :
323 921430 : void Read::update() {
324 921430 : if( !cloned_file ) {
325 403971 : for(unsigned i=0; i<nlinesPerStep; ++i) {
326 202261 : ifile->scanField();
327 : double du_time;
328 404522 : if( !ifile->scanField("time",du_time) && !plumed.inputsAreActive() ) {
329 55 : plumed.stop();
330 : }
331 : }
332 : }
333 921430 : }
334 :
335 : }
336 : }
|
