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 : {
81 : private:
82 : bool ignore_time;
83 : bool ignore_forces;
84 : bool cloned_file;
85 : unsigned nlinesPerStep;
86 : std::string filename;
87 : /// Unique pointer with the same scope as ifile.
88 : std::unique_ptr<IFile> ifile_ptr;
89 : /// Pointer to input file.
90 : /// It is either pointing to the content of ifile_ptr
91 : /// or to the file it is cloned from.
92 : IFile* ifile;
93 : std::vector<std::unique_ptr<Value>> readvals;
94 : public:
95 : static void registerKeywords( Keywords& keys );
96 : explicit Read(const ActionOptions&);
97 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
98 : void prepare() override;
99 921429 : 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 : PLUMED_REGISTER_ACTION(Read,"READ")
109 :
110 95 : void Read::registerKeywords(Keywords& keys) {
111 95 : Action::registerKeywords(keys);
112 95 : ActionPilot::registerKeywords(keys);
113 95 : ActionWithValue::registerKeywords(keys);
114 190 : keys.add("compulsory","STRIDE","1","the frequency with which the file should be read.");
115 190 : 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.");
116 190 : keys.add("compulsory","VALUES","the values to read from the file");
117 190 : keys.add("compulsory","FILE","the name of the file from which to read these quantities");
118 190 : 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 190 : 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 95 : keys.remove("NUMERICAL_DERIVATIVES");
123 95 : keys.use("UPDATE_FROM");
124 95 : keys.use("UPDATE_UNTIL");
125 95 : ActionWithValue::useCustomisableComponents(keys);
126 95 : }
127 :
128 91 : Read::Read(const ActionOptions&ao):
129 : Action(ao),
130 : ActionPilot(ao),
131 : ActionWithValue(ao),
132 91 : ignore_time(false),
133 91 : ignore_forces(false),
134 91 : nlinesPerStep(1)
135 : {
136 : // Read the file name from the input line
137 91 : parse("FILE",filename);
138 : // Check if time is to be ignored
139 91 : parseFlag("IGNORE_TIME",ignore_time);
140 : // Check if forces are to be ignored
141 91 : parseFlag("IGNORE_FORCES",ignore_forces);
142 : // Open the file if it is not already opened
143 91 : cloned_file=false;
144 91 : std::vector<Read*> other_reads=plumed.getActionSet().select<Read*>();
145 141 : for(unsigned i=0; i<other_reads.size(); ++i) {
146 50 : if( other_reads[i]->getFilename()==filename ) {
147 47 : ifile=other_reads[i]->getFile();
148 47 : cloned_file=true;
149 : }
150 : }
151 91 : if( !cloned_file ) {
152 63 : ifile_ptr=Tools::make_unique<IFile>();
153 63 : ifile=ifile_ptr.get();
154 63 : if( !ifile->FileExist(filename) ) error("could not find file named " + filename);
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) log.printf(" only reading every %uth line of file %s\n",nlinesPerStep,filename.c_str() );
161 89 : else log.printf(" reading data from file %s\n",filename.c_str() );
162 : // Find out what we are reading
163 91 : std::vector<std::string> valread; parseVector("VALUES",valread);
164 :
165 91 : 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 91 : if( valread[0].find(".")!=std::string::npos ) {
169 11 : std::string label=valread[0].substr(0,dot);
170 11 : std::string name=valread[0].substr(dot+1);
171 11 : if( name=="*" ) {
172 2 : if( valread.size()>1 ) error("all values must be from the same Action when using READ");
173 : std::vector<std::string> fieldnames;
174 2 : ifile->scanFieldList( fieldnames );
175 16 : for(unsigned i=0; i<fieldnames.size(); ++i) {
176 14 : if( fieldnames[i].substr(0,dot)==label ) {
177 12 : readvals.emplace_back(Tools::make_unique<Value>(this, fieldnames[i], false) ); addComponentWithDerivatives( fieldnames[i].substr(dot+1) );
178 12 : if( ifile->FieldExist("min_" + fieldnames[i]) ) componentIsPeriodic( fieldnames[i].substr(dot+1), "-pi","pi" );
179 12 : else componentIsNotPeriodic( fieldnames[i].substr(dot+1) );
180 : }
181 : }
182 2 : } else {
183 9 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) ); addComponentWithDerivatives( name );
184 18 : if( ifile->FieldExist("min_" + valread[0]) ) componentIsPeriodic( valread[0].substr(dot+1), "-pi", "pi" );
185 18 : else componentIsNotPeriodic( valread[0].substr(dot+1) );
186 12 : for(unsigned i=1; i<valread.size(); ++i) {
187 6 : if( valread[i].substr(0,dot)!=label ) error("all values must be from the same Action when using READ");;
188 6 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[i], false) ); addComponentWithDerivatives( valread[i].substr(dot+1) );
189 6 : if( ifile->FieldExist("min_" + valread[i]) ) componentIsPeriodic( valread[i].substr(dot+1), "-pi", "pi" );
190 6 : else componentIsNotPeriodic( valread[i].substr(dot+1) );
191 : }
192 : }
193 : } else {
194 80 : if( valread.size()!=1 ) error("all values must be from the same Action when using READ");
195 80 : readvals.emplace_back(Tools::make_unique<Value>(this, valread[0], false) ); addValueWithDerivatives();
196 169 : if( ifile->FieldExist("min_" + valread[0]) ) setPeriodic( "-pi", "pi" );
197 71 : else setNotPeriodic();
198 80 : log.printf(" reading value %s and storing as %s\n",valread[0].c_str(),getLabel().c_str() );
199 : }
200 91 : checkRead();
201 182 : }
202 :
203 4 : std::string Read::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
204 4 : plumed_assert( !exists( getLabel() ) );
205 7 : for(unsigned i=0; i<readvals.size(); ++i) {
206 11 : if( readvals[i]->getName().find( cname )!=std::string::npos ) return "values from the column labelled " + readvals[i]->getName() + " in the file named " + filename;
207 : }
208 0 : plumed_error(); return "";
209 : }
210 :
211 50 : std::string Read::getFilename() const {
212 50 : return filename;
213 : }
214 :
215 47 : IFile* Read::getFile() {
216 47 : return ifile;
217 : }
218 :
219 0 : unsigned Read::getNumberOfDerivatives() {
220 0 : return 0;
221 : }
222 :
223 26 : void Read::turnOnDerivatives() {
224 26 : if( !ignore_forces ) error("cannot calculate derivatives for colvars that are read in from a file. If you are postprocessing and "
225 : "these forces do not matter add the flag IGNORE_FORCES to all READ actions");
226 26 : }
227 :
228 921429 : void Read::prepare() {
229 921429 : if( !cloned_file ) {
230 : double du_time;
231 403414 : if( !ifile->scanField("time",du_time) ) {
232 0 : error("Reached end of file " + filename + " before end of trajectory");
233 201707 : } else if( std::abs( du_time-getTime() )>getTimeStep() && !ignore_time ) {
234 0 : std::string str_dutime,str_ptime; Tools::convert(du_time,str_dutime); Tools::convert(getTime(),str_ptime);
235 0 : error("mismatched times in colvar files : colvar time=" + str_dutime + " plumed time=" + str_ptime + ". Add IGNORE_TIME to ignore error.");
236 : }
237 : }
238 921429 : }
239 :
240 921429 : void Read::calculate() {
241 : std::string smin, smax;
242 2019293 : for(unsigned i=0; i<readvals.size(); ++i) {
243 : // .get returns the raw pointer
244 : // ->get calls the Value::get() method
245 1097864 : ifile->scanField( readvals[i].get() );
246 1097864 : getPntrToComponent(i)->set( readvals[i]->get() );
247 1097864 : if( readvals[i]->isPeriodic() ) {
248 3309 : readvals[i]->getDomain( smin, smax );
249 3309 : getPntrToComponent(i)->setDomain( smin, smax );
250 : }
251 : }
252 921429 : }
253 :
254 921429 : void Read::update() {
255 921429 : if( !cloned_file ) {
256 403965 : for(unsigned i=0; i<nlinesPerStep; ++i) {
257 202258 : ifile->scanField(); double du_time;
258 404516 : if( !ifile->scanField("time",du_time) && !plumed.inputsAreActive() ) plumed.stop();
259 : }
260 : }
261 921429 : }
262 :
263 : }
264 : }
|
