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