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
39 :
40 : \par Description of components
41 :
42 : The READ command will read those fields that are labelled with the text string given to the
43 : VALUE keyword. It will also read in any fields that are labeled with the text string
44 : given to the VALUE keyword followed by a dot and a further string. If a single Value is read in
45 : this value can be referenced using the label of the Action. Alternatively, if multiple quantities
46 : are read in, they can be referenced elsewhere in the input by using the label for the Action
47 : followed by a dot and the character string that appeared after the dot in the title of the field.
48 :
49 : \par Examples
50 :
51 : This input reads in data from a file called input_colvar.data that was generated
52 : in a calculation that involved PLUMED. The first command reads in the data from the
53 : column headed phi1 while the second reads in the data from the column headed phi2.
54 :
55 : \plumedfile
56 : rphi1: READ FILE=input_colvar.data VALUES=phi1
57 : rphi2: READ FILE=input_colvar.data VALUES=phi2
58 : PRINT ARG=rphi1,rphi2 STRIDE=500 FILE=output_colvar.data
59 : \endplumedfile
60 :
61 : The file input_colvar.data is just a normal colvar file as shown below
62 :
63 : \auxfile{input_colvar.data}
64 : #! FIELDS time phi psi metad.bias metad.rbias metad.rct
65 : #! SET min_phi -pi
66 : #! SET max_phi pi
67 : #! SET min_psi -pi
68 : #! SET max_psi pi
69 : 0.000000 -1.2379 0.8942 0.0000 0.0000 0.0000
70 : 1.000000 -1.4839 1.0482 0.0000 0.0000 0.0089
71 : 2.000000 -1.3243 0.6055 0.0753 0.0664 0.0184
72 : \endauxfile
73 :
74 : */
75 : //+ENDPLUMEDOC
76 :
77 : class Read :
78 : public ActionPilot,
79 : public ActionWithValue {
80 : private:
81 : bool ignore_time;
82 : bool ignore_forces;
83 : bool cloned_file;
84 : unsigned nlinesPerStep;
85 : std::string filename;
86 : /// Unique pointer with the same scope as ifile.
87 : std::unique_ptr<IFile> ifile_ptr;
88 : /// Pointer to input file.
89 : /// It is either pointing to the content of ifile_ptr
90 : /// or to the file it is cloned from.
91 : IFile* ifile;
92 : std::vector<std::unique_ptr<Value>> readvals;
93 : public:
94 : static void registerKeywords( Keywords& keys );
95 : explicit Read(const ActionOptions&);
96 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
97 : void prepare() override;
98 921429 : void apply() override {}
99 : void calculate() override;
100 : void update() override;
101 : std::string getFilename() const;
102 : IFile* getFile();
103 : unsigned getNumberOfDerivatives() override;
104 : void turnOnDerivatives() override;
105 : };
106 :
107 : PLUMED_REGISTER_ACTION(Read,"READ")
108 :
109 97 : void Read::registerKeywords(Keywords& keys) {
110 97 : Action::registerKeywords(keys);
111 97 : ActionPilot::registerKeywords(keys);
112 97 : ActionWithValue::registerKeywords(keys);
113 194 : keys.add("compulsory","STRIDE","1","the frequency with which the file should be read.");
114 194 : 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.");
115 194 : keys.add("compulsory","VALUES","the values to read from the file");
116 194 : keys.add("compulsory","FILE","the name of the file from which to read these quantities");
117 194 : keys.addFlag("IGNORE_TIME",false,"ignore the time in the colvar file. When this flag is not present read will be quite strict "
118 : "about the start time of the simulation and the stride between frames");
119 194 : 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 "
120 : "safely ignored if you are doing post processing that does not involve outputting forces");
121 97 : keys.remove("NUMERICAL_DERIVATIVES");
122 97 : keys.use("UPDATE_FROM");
123 97 : keys.use("UPDATE_UNTIL");
124 97 : ActionWithValue::useCustomisableComponents(keys);
125 97 : }
126 :
127 91 : Read::Read(const ActionOptions&ao):
128 : Action(ao),
129 : ActionPilot(ao),
130 : ActionWithValue(ao),
131 91 : ignore_time(false),
132 91 : ignore_forces(false),
133 91 : nlinesPerStep(1) {
134 : // Read the file name from the input line
135 91 : parse("FILE",filename);
136 : // Check if time is to be ignored
137 91 : parseFlag("IGNORE_TIME",ignore_time);
138 : // Check if forces are to be ignored
139 91 : parseFlag("IGNORE_FORCES",ignore_forces);
140 : // Open the file if it is not already opened
141 91 : cloned_file=false;
142 91 : std::vector<Read*> other_reads=plumed.getActionSet().select<Read*>();
143 141 : for(unsigned i=0; i<other_reads.size(); ++i) {
144 50 : if( other_reads[i]->getFilename()==filename ) {
145 47 : ifile=other_reads[i]->getFile();
146 47 : cloned_file=true;
147 : }
148 : }
149 91 : if( !cloned_file ) {
150 63 : ifile_ptr=Tools::make_unique<IFile>();
151 63 : ifile=ifile_ptr.get();
152 63 : if( !ifile->FileExist(filename) ) {
153 0 : error("could not find file named " + filename);
154 : }
155 63 : ifile->link(*this);
156 63 : ifile->open(filename);
157 63 : ifile->allowIgnoredFields();
158 : }
159 91 : parse("EVERY",nlinesPerStep);
160 91 : if(nlinesPerStep>1) {
161 2 : log.printf(" only reading every %uth line of file %s\n",nlinesPerStep,filename.c_str() );
162 : } else {
163 89 : log.printf(" reading data from file %s\n",filename.c_str() );
164 : }
165 : // Find out what we are reading
166 : std::vector<std::string> valread;
167 91 : parseVector("VALUES",valread);
168 :
169 91 : if(nlinesPerStep>1 && cloned_file) {
170 0 : error("Opening a file multiple times and using EVERY is not allowed");
171 : }
172 :
173 : std::size_t dot=valread[0].find_first_of('.');
174 91 : if( valread[0].find(".")!=std::string::npos ) {
175 11 : std::string label=valread[0].substr(0,dot);
176 11 : std::string name=valread[0].substr(dot+1);
177 11 : if( name=="*" ) {
178 2 : if( valread.size()>1 ) {
179 0 : error("all values must be from the same Action when using READ");
180 : }
181 : std::vector<std::string> fieldnames;
182 2 : ifile->scanFieldList( fieldnames );
183 16 : for(unsigned i=0; i<fieldnames.size(); ++i) {
184 14 : if( fieldnames[i].substr(0,dot)==label ) {
185 6 : readvals.emplace_back(Tools::make_unique<Value>(this, fieldnames[i], false) );
186 12 : addComponentWithDerivatives( fieldnames[i].substr(dot+1) );
187 12 : if( ifile->FieldExist("min_" + fieldnames[i]) ) {
188 0 : componentIsPeriodic( fieldnames[i].substr(dot+1), "-pi","pi" );
189 : } else {
190 12 : componentIsNotPeriodic( fieldnames[i].substr(dot+1) );
191 : }
192 : }
193 : }
194 2 : } else {
195 9 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) );
196 9 : addComponentWithDerivatives( name );
197 18 : if( ifile->FieldExist("min_" + valread[0]) ) {
198 0 : componentIsPeriodic( valread[0].substr(dot+1), "-pi", "pi" );
199 : } else {
200 18 : componentIsNotPeriodic( valread[0].substr(dot+1) );
201 : }
202 12 : for(unsigned i=1; i<valread.size(); ++i) {
203 6 : if( valread[i].substr(0,dot)!=label ) {
204 0 : error("all values must be from the same Action when using READ");
205 : };
206 3 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[i], false) );
207 6 : addComponentWithDerivatives( valread[i].substr(dot+1) );
208 6 : if( ifile->FieldExist("min_" + valread[i]) ) {
209 0 : componentIsPeriodic( valread[i].substr(dot+1), "-pi", "pi" );
210 : } else {
211 6 : componentIsNotPeriodic( valread[i].substr(dot+1) );
212 : }
213 : }
214 : }
215 : } else {
216 80 : if( valread.size()!=1 ) {
217 0 : error("all values must be from the same Action when using READ");
218 : }
219 80 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) );
220 80 : addValueWithDerivatives();
221 160 : if( ifile->FieldExist("min_" + valread[0]) ) {
222 18 : setPeriodic( "-pi", "pi" );
223 : } else {
224 71 : setNotPeriodic();
225 : }
226 80 : log.printf(" reading value %s and storing as %s\n",valread[0].c_str(),getLabel().c_str() );
227 : }
228 91 : checkRead();
229 182 : }
230 :
231 4 : std::string Read::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
232 4 : plumed_assert( !exists( getLabel() ) );
233 7 : for(unsigned i=0; i<readvals.size(); ++i) {
234 7 : if( readvals[i]->getName().find( cname )!=std::string::npos ) {
235 8 : return "values from the column labelled " + readvals[i]->getName() + " in the file named " + filename;
236 : }
237 : }
238 0 : plumed_error();
239 : return "";
240 : }
241 :
242 50 : std::string Read::getFilename() const {
243 50 : return filename;
244 : }
245 :
246 47 : IFile* Read::getFile() {
247 47 : return ifile;
248 : }
249 :
250 0 : unsigned Read::getNumberOfDerivatives() {
251 0 : return 0;
252 : }
253 :
254 26 : void Read::turnOnDerivatives() {
255 26 : if( !ignore_forces )
256 0 : error("cannot calculate derivatives for colvars that are read in from a file. If you are postprocessing and "
257 : "these forces do not matter add the flag IGNORE_FORCES to all READ actions");
258 26 : }
259 :
260 921429 : void Read::prepare() {
261 921429 : if( !cloned_file ) {
262 : double du_time;
263 403414 : if( !ifile->scanField("time",du_time) ) {
264 0 : error("Reached end of file " + filename + " before end of trajectory");
265 201707 : } else if( std::abs( du_time-getTime() )>getTimeStep() && !ignore_time ) {
266 : std::string str_dutime,str_ptime;
267 0 : Tools::convert(du_time,str_dutime);
268 0 : Tools::convert(getTime(),str_ptime);
269 0 : error("mismatched times in colvar files : colvar time=" + str_dutime + " plumed time=" + str_ptime + ". Add IGNORE_TIME to ignore error.");
270 : }
271 : }
272 921429 : }
273 :
274 921429 : void Read::calculate() {
275 : std::string smin, smax;
276 2019293 : for(unsigned i=0; i<readvals.size(); ++i) {
277 : // .get returns the raw pointer
278 : // ->get calls the Value::get() method
279 1097864 : ifile->scanField( readvals[i].get() );
280 1097864 : getPntrToComponent(i)->set( readvals[i]->get() );
281 1097864 : if( readvals[i]->isPeriodic() ) {
282 3309 : readvals[i]->getDomain( smin, smax );
283 3309 : getPntrToComponent(i)->setDomain( smin, smax );
284 : }
285 : }
286 921429 : }
287 :
288 921429 : void Read::update() {
289 921429 : if( !cloned_file ) {
290 403965 : for(unsigned i=0; i<nlinesPerStep; ++i) {
291 202258 : ifile->scanField();
292 : double du_time;
293 404516 : if( !ifile->scanField("time",du_time) && !plumed.inputsAreActive() ) {
294 54 : plumed.stop();
295 : }
296 : }
297 : }
298 921429 : }
299 :
300 : }
301 : }
|
