Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 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 :
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 labelleled 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 quanties
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 : */
62 : //+ENDPLUMEDOC
63 :
64 : class Read :
65 : public ActionPilot,
66 : public ActionWithValue
67 : {
68 : private:
69 : bool ignore_time;
70 : bool ignore_forces;
71 : bool cloned_file;
72 : unsigned nlinesPerStep;
73 : std::string filename;
74 : IFile* ifile;
75 : std::vector<Value*> readvals;
76 : public:
77 : static void registerKeywords( Keywords& keys );
78 : explicit Read(const ActionOptions&);
79 : ~Read();
80 : void prepare();
81 2331 : void apply() {}
82 : void calculate();
83 : void update();
84 : std::string getFilename() const;
85 : IFile* getFile();
86 : unsigned getNumberOfDerivatives();
87 : void turnOnDerivatives();
88 : };
89 :
90 6489 : PLUMED_REGISTER_ACTION(Read,"READ")
91 :
92 38 : void Read::registerKeywords(Keywords& keys) {
93 38 : Action::registerKeywords(keys);
94 38 : ActionPilot::registerKeywords(keys);
95 38 : ActionWithValue::registerKeywords(keys);
96 190 : keys.add("compulsory","STRIDE","1","the frequency with which the file should be read.");
97 190 : keys.add("compulsory","EVERY","1","only read every ith line of the colvar file. This should be used if the colvar was written more frequently than the trajectory.");
98 152 : keys.add("compulsory","VALUES","the values to read from the file");
99 152 : keys.add("compulsory","FILE","the name of the file from which to read these quantities");
100 114 : keys.addFlag("IGNORE_TIME",false,"ignore the time in the colvar file. When this flag is not present read will be quite strict "
101 : "about the start time of the simulation and the stride between frames");
102 114 : 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 "
103 : "safely ignored if you are doing postprocessing that does not involve outputting forces");
104 76 : keys.remove("NUMERICAL_DERIVATIVES");
105 76 : keys.use("UPDATE_FROM");
106 76 : keys.use("UPDATE_UNTIL");
107 38 : ActionWithValue::useCustomisableComponents(keys);
108 38 : }
109 :
110 37 : Read::Read(const ActionOptions&ao):
111 : Action(ao),
112 : ActionPilot(ao),
113 : ActionWithValue(ao),
114 : ignore_time(false),
115 : ignore_forces(false),
116 74 : nlinesPerStep(1)
117 : {
118 : // Read the file name from the input line
119 74 : parse("FILE",filename);
120 : // Check if time is to be ignored
121 74 : parseFlag("IGNORE_TIME",ignore_time);
122 : // Check if forces are to be ignored
123 74 : parseFlag("IGNORE_FORCES",ignore_forces);
124 : // Open the file if it is not already opened
125 37 : cloned_file=false;
126 74 : std::vector<Read*> other_reads=plumed.getActionSet().select<Read*>();
127 80 : for(unsigned i=0; i<other_reads.size(); ++i) {
128 4 : if( other_reads[i]->getFilename()==filename ) {
129 2 : ifile=other_reads[i]->getFile();
130 2 : cloned_file=true;
131 : }
132 : }
133 37 : if( !cloned_file ) {
134 35 : ifile=new IFile();
135 35 : if( !ifile->FileExist(filename) ) error("could not find file named " + filename);
136 35 : ifile->link(*this);
137 35 : ifile->open(filename);
138 35 : ifile->allowIgnoredFields();
139 : }
140 74 : parse("EVERY",nlinesPerStep);
141 37 : if(nlinesPerStep>1) log.printf(" only reading every %uth line of file %s\n",nlinesPerStep,filename.c_str() );
142 74 : else log.printf(" reading data from file %s\n",filename.c_str() );
143 : // Find out what we are reading
144 111 : std::vector<std::string> valread; parseVector("VALUES",valread);
145 :
146 : std::size_t dot=valread[0].find_first_of('.');
147 37 : if( valread[0].find(".")!=std::string::npos ) {
148 3 : std::string label=valread[0].substr(0,dot);
149 3 : std::string name=valread[0].substr(dot+1);
150 3 : if( name=="*" ) {
151 1 : if( valread.size()>1 ) error("all values must be from the same Action when using READ");
152 1 : std::vector<std::string> fieldnames;
153 1 : ifile->scanFieldList( fieldnames );
154 23 : for(unsigned i=0; i<fieldnames.size(); ++i) {
155 14 : if( fieldnames[i].substr(0,dot)==label ) {
156 12 : readvals.push_back(new Value(this, fieldnames[i], false) ); addComponentWithDerivatives( fieldnames[i].substr(dot+1) );
157 9 : if( ifile->FieldExist("min_" + fieldnames[i]) ) componentIsPeriodic( fieldnames[i].substr(dot+1), "-pi","pi" );
158 6 : else componentIsNotPeriodic( fieldnames[i].substr(dot+1) );
159 : }
160 : }
161 : } else {
162 4 : readvals.push_back(new Value(this, valread[0], false) ); addComponentWithDerivatives( name );
163 6 : if( ifile->FieldExist("min_" + valread[0]) ) componentIsPeriodic( valread[0].substr(dot+1), "-pi", "pi" );
164 4 : else componentIsNotPeriodic( valread[0].substr(dot+1) );
165 7 : for(unsigned i=1; i<valread.size(); ++i) {
166 2 : if( valread[i].substr(0,dot)!=label ) error("all values must be from the same Action when using READ");;
167 4 : readvals.push_back(new Value(this, valread[i], false) ); addComponentWithDerivatives( valread[i].substr(dot+1) );
168 3 : if( ifile->FieldExist("min_" + valread[i]) ) componentIsPeriodic( valread[i].substr(dot+1), "-pi", "pi" );
169 2 : else componentIsNotPeriodic( valread[i].substr(dot+1) );
170 : }
171 : }
172 : } else {
173 34 : if( valread.size()!=1 ) error("all values must be from the same Action when using READ");
174 68 : readvals.push_back(new Value(this, valread[0], false) ); addValueWithDerivatives();
175 102 : if( ifile->FieldExist("min_" + valread[0]) ) setPeriodic( "-pi", "pi" );
176 34 : else setNotPeriodic();
177 68 : log.printf(" reading value %s and storing as %s\n",valread[0].c_str(),getLabel().c_str() );
178 : }
179 37 : checkRead();
180 37 : }
181 :
182 148 : Read::~Read() {
183 37 : if( !cloned_file ) { ifile->close(); delete ifile; }
184 194 : for(unsigned i=0; i<readvals.size(); ++i) delete readvals[i];
185 74 : }
186 :
187 2 : std::string Read::getFilename() const {
188 2 : return filename;
189 : }
190 :
191 2 : IFile* Read::getFile() {
192 2 : return ifile;
193 : }
194 :
195 0 : unsigned Read::getNumberOfDerivatives() {
196 0 : return 0;
197 : }
198 :
199 0 : void Read::turnOnDerivatives() {
200 0 : if( !ignore_forces ) error("cannot calculate derivatives for colvars that are read in from a file. If you are postprocessing and "
201 : "these forces do not matter add the flag IGNORE_FORCES to all READ actions");
202 0 : }
203 :
204 2331 : void Read::prepare() {
205 2331 : if( !cloned_file ) {
206 : double du_time;
207 2476 : if( !ifile->scanField("time",du_time) ) {
208 0 : error("Reached end of file " + filename + " before end of trajectory");
209 2476 : } else if( fabs( du_time-getTime() )>plumed.getAtoms().getTimeStep() && !ignore_time ) {
210 0 : std::string str_dutime,str_ptime; Tools::convert(du_time,str_dutime); Tools::convert(getTime(),str_ptime);
211 0 : error("mismatched times in colvar files : colvar time=" + str_dutime + " plumed time=" + str_ptime + ". Add IGNORE_TIME to ignore error.");
212 : }
213 : }
214 2331 : }
215 :
216 2331 : void Read::calculate() {
217 : std::string smin, smax;
218 15234 : for(unsigned i=0; i<readvals.size(); ++i) {
219 7048 : ifile->scanField( readvals[i] );
220 7048 : getPntrToComponent(i)->set( readvals[i]->get() );
221 3524 : if( readvals[i]->isPeriodic() ) {
222 0 : readvals[i]->getDomain( smin, smax );
223 0 : getPntrToComponent(i)->setDomain( smin, smax );
224 : }
225 : }
226 2331 : }
227 :
228 2331 : void Read::update() {
229 2331 : if( !cloned_file ) {
230 3714 : for(unsigned i=0; i<nlinesPerStep; ++i) {
231 1238 : ifile->scanField(); double du_time;
232 3714 : if( plumed.getAtoms().getNatoms()==0 && !ifile->scanField("time",du_time) ) plumed.stop();
233 : }
234 : }
235 2331 : }
236 :
237 : }
238 4839 : }
|
