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 :
23 : #include "Analysis.h"
24 : #include "core/ActionSet.h"
25 : #include "core/ActionWithValue.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/Atoms.h"
28 : #include "tools/IFile.h"
29 : #include "reference/ReferenceConfiguration.h"
30 : #include "reference/ReferenceArguments.h"
31 : #include "reference/ReferenceAtoms.h"
32 : #include "reference/MetricRegister.h"
33 :
34 : namespace PLMD {
35 : namespace analysis {
36 :
37 5 : void Analysis::registerKeywords( Keywords& keys ) {
38 5 : vesselbase::ActionWithAveraging::registerKeywords( keys );
39 20 : keys.use("ARG"); keys.reset_style("ARG","optional");
40 20 : keys.add("atoms","ATOMS","the atoms whose positions we are tracking for the purpose of analysing the data");
41 25 : keys.add("compulsory","METRIC","EUCLIDEAN","how are we measuring the distances between configurations");
42 25 : keys.add("compulsory","RUN","0","the frequency with which to run the analysis algorithm. The default value of zero assumes you want to analyse the whole trajectory");
43 20 : keys.add("optional","FMT","the format that should be used in analysis output files");
44 15 : keys.addFlag("WRITE_CHECKPOINT",false,"write out a checkpoint so that the analysis can be restarted in a later run");
45 20 : keys.add("hidden","REUSE_DATA_FROM","eventually this will allow you to analyse the same set of data multiple times");
46 20 : keys.add("hidden","IGNORE_REWEIGHTING","this allows you to ignore any reweighting factors");
47 25 : keys.use("RESTART"); keys.use("UPDATE_FROM"); keys.use("UPDATE_UNTIL"); keys.remove("TOL");
48 5 : }
49 :
50 3 : Analysis::Analysis(const ActionOptions&ao):
51 : Action(ao),
52 : ActionWithAveraging(ao),
53 : nomemory(true),
54 : write_chq(false),
55 : reusing_data(false),
56 : ignore_reweight(false),
57 : idata(0),
58 : //firstAnalysisDone(false),
59 : //old_norm(0.0),
60 : ofmt("%f"),
61 : current_args(getNumberOfArguments()),
62 15 : argument_names(getNumberOfArguments())
63 : {
64 6 : parse("FMT",ofmt); // Read the format for output files
65 : // Make a vector containing all the argument names
66 11 : for(unsigned i=0; i<getNumberOfArguments(); ++i) argument_names[i]=getPntrToArgument(i)->getName();
67 : // Read in the metric style
68 6 : parse("METRIC",metricname); std::vector<AtomNumber> atom_numbers;
69 3 : ReferenceConfiguration* checkref=metricRegister().create<ReferenceConfiguration>( metricname );
70 : // Check if we should read atoms
71 3 : ReferenceAtoms* hasatoms=dynamic_cast<ReferenceAtoms*>( checkref );
72 3 : if( hasatoms ) {
73 2 : parseAtomList("ATOMS",atom_numbers); requestAtoms(atom_numbers);
74 1 : log.printf(" monitoring positions of atoms ");
75 68 : for(unsigned i=0; i<atom_numbers.size(); ++i) log.printf("%d ",atom_numbers[i].serial() );
76 1 : log.printf("\n");
77 : }
78 : // Check if we should read arguments
79 3 : ReferenceArguments* hasargs=dynamic_cast<ReferenceArguments*>( checkref );
80 4 : if( !hasargs && getNumberOfArguments()!=0 ) error("use of arguments with metric type " + metricname + " is invalid");
81 3 : if( hasatoms && hasargs ) error("currently dependencies break if you have both arguments and atoms");
82 : // And delte the fake reference we created
83 3 : delete checkref;
84 :
85 6 : std::string prev_analysis; parse("REUSE_DATA_FROM",prev_analysis);
86 3 : if( prev_analysis.length()>0 ) {
87 0 : reusing_data=true;
88 0 : mydatastash=plumed.getActionSet().selectWithLabel<Analysis*>( prev_analysis );
89 0 : if( !mydatastash ) error("could not find analysis action named " + prev_analysis );
90 0 : parseFlag("IGNORE_REWEIGHTING",ignore_reweight);
91 0 : if( ignore_reweight ) log.printf(" reusing data stored by %s but ignoring all reweighting\n",prev_analysis.c_str() );
92 0 : else log.printf(" reusing data stored by %s\n",prev_analysis.c_str() );
93 : } else {
94 6 : parse("RUN",freq);
95 3 : if( freq==0 ) {
96 2 : log.printf(" analyzing all data in trajectory\n");
97 : } else {
98 1 : if( freq%getStride()!=0 ) error("frequncy of running is not a multiple of the stride");
99 1 : log.printf(" running analysis every %u steps\n",freq);
100 1 : ndata=freq/getStride(); data.resize( ndata ); logweights.resize( ndata );
101 201 : for(unsigned i=0; i<ndata; ++i) {
102 200 : data[i]=metricRegister().create<ReferenceConfiguration>( metricname );
103 200 : data[i]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
104 : }
105 : }
106 6 : parseFlag("WRITE_CHECKPOINT",write_chq);
107 3 : if( write_chq ) {
108 0 : write_chq=false;
109 0 : warning("ignoring WRITE_CHECKPOINT flag because we are analyzing all data");
110 : }
111 :
112 : // We need no restart file if we are just collecting data and analyzing all of it
113 12 : std::string filename = getName() + "_" + getLabel() + ".chkpnt";
114 3 : if( write_chq ) rfile.link(*this);
115 6 : if( getRestart() ) {
116 0 : if( !write_chq ) warning("restarting without writing a checkpoint file is somewhat strange");
117 : // Read in data from input file
118 0 : readDataFromFile( filename );
119 : // Setup the restart file (append mode)
120 0 : if( write_chq ) rfile.open( filename.c_str() ); // In append mode automatically because of restart
121 : // Run the analysis if we stoped in the middle of it last time
122 0 : log.printf(" restarting analysis with %u points read from restart file\n",idata);
123 3 : } else if( write_chq ) {
124 : // Setup the restart file (delete any old one)
125 0 : rfile.open( filename.c_str() ); // In overwrite mode automatically because there is no restart
126 : }
127 3 : if( write_chq ) {
128 : //rfile.addConstantField("old_normalization");
129 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) rfile.setupPrintValue( getPntrToArgument(i) );
130 : }
131 : }
132 3 : }
133 :
134 0 : void Analysis::readDataFromFile( const std::string& filename ) {
135 0 : FILE* fp=fopen(filename.c_str(),"r");
136 0 : if(fp!=NULL) {
137 : double tstep, oldtstep;
138 : bool do_read=true, first=true;
139 0 : while (do_read) {
140 0 : PDB mypdb;
141 0 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
142 0 : if(do_read) {
143 0 : data[idata]->set( mypdb );
144 0 : data[idata]->parse("TIME",tstep);
145 0 : if( !first && ((tstep-oldtstep) - getStride()*plumed.getAtoms().getTimeStep())>plumed.getAtoms().getTimeStep() ) {
146 0 : error("frequency of data storage in " + filename + " is not equal to frequency of data storage plumed.dat file");
147 : }
148 0 : data[idata]->parse("LOG_WEIGHT",logweights[idata]);
149 : //data[idata]->parse("OLD_NORM",old_norm);
150 0 : data[idata]->checkRead();
151 0 : idata++; first=false; oldtstep=tstep;
152 : } else {
153 : break;
154 : }
155 : }
156 0 : fclose(fp);
157 : }
158 : // if(old_norm>0) firstAnalysisDone=true;
159 0 : }
160 :
161 4 : void Analysis::parseOutputFile( const std::string& key, std::string& filename ) {
162 4 : parse(key,filename);
163 4 : if(filename=="dont output") return;
164 :
165 8 : if( !getRestart() ) {
166 8 : OFile ofile; ofile.link(*this);
167 8 : ofile.setBackupString("analysis");
168 4 : ofile.backupAllFiles(filename);
169 : }
170 : }
171 :
172 1292 : void Analysis::accumulate() {
173 : // Don't store the first step (also don't store if we are getting data from elsewhere)
174 1292 : if( getStep()==0 || reusing_data ) return;
175 : // This is used when we have a full quota of data from the first run
176 1492 : if( freq>0 && idata==logweights.size() ) return;
177 : // Get the arguments ready to transfer to reference configuration
178 4276 : for(unsigned i=0; i<getNumberOfArguments(); ++i) current_args[i]=getArgument(i);
179 :
180 1292 : if( freq>0) {
181 : // Get the arguments and store them in a vector of vectors
182 800 : data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() );
183 400 : logweights[idata] = lweight;
184 : } else {
185 2184 : data.push_back( metricRegister().create<ReferenceConfiguration>( metricname ) );
186 : plumed_dbg_assert( data.size()==idata+1 );
187 2184 : data[idata]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
188 4368 : data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() );
189 1092 : logweights.push_back(lweight);
190 : }
191 :
192 : // Write data to checkpoint file
193 1292 : if( write_chq ) {
194 0 : rfile.rewind();
195 0 : data[idata]->print( rfile, getTime(), logweights[idata], atoms.getUnits().getLength()/0.1, 1.0 ); //old_norm );
196 0 : rfile.flush();
197 : }
198 : // Increment data counter
199 1292 : idata++;
200 : }
201 :
202 9 : Analysis::~Analysis() {
203 3582 : for(unsigned i=0; i<data.size(); ++i ) delete data[i];
204 3 : if( write_chq ) rfile.close();
205 3 : }
206 :
207 1292 : std::vector<double> Analysis::getMetric() const {
208 : // Add more exotic metrics in here -- FlexibleHill for instance
209 : std::vector<double> empty;
210 2584 : if( metricname=="EUCLIDEAN" ) {
211 746 : empty.resize( getNumberOfArguments(), 1.0 );
212 : }
213 1292 : return empty;
214 : }
215 :
216 2394732 : double Analysis::getWeight( const unsigned& idata ) const {
217 2394732 : if( !reusing_data ) {
218 : plumed_dbg_assert( idata<data.size() );
219 7184196 : return data[idata]->getWeight();
220 : } else {
221 0 : return mydatastash->getWeight(idata);
222 : }
223 : }
224 :
225 4 : void Analysis::finalizeWeights( const bool& ignore_weights ) {
226 : // Check that we have the correct ammount of data
227 8 : if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong. Am trying to run analysis but I don't have sufficient data");
228 :
229 : double norm=0; // Reset normalization constant
230 4 : if( ignore_weights ) {
231 0 : for(unsigned i=0; i<logweights.size(); ++i) {
232 0 : data[i]->setWeight(1.0); norm+=1.0;
233 : }
234 4 : } else if( nomemory ) {
235 : // Find the maximum weight
236 4 : double maxweight=logweights[0];
237 2580 : for(unsigned i=1; i<getNumberOfDataPoints(); ++i) {
238 2576 : if(logweights[i]>maxweight) maxweight=logweights[i];
239 : }
240 : // Calculate normalization constant
241 3884 : for(unsigned i=0; i<logweights.size(); ++i) {
242 1292 : norm+=exp( logweights[i]-maxweight );
243 : }
244 : // Calculate weights (no memory)
245 3884 : for(unsigned i=0; i<logweights.size(); ++i) {
246 2584 : data[i]->setWeight( exp( logweights[i]-maxweight ) );
247 : }
248 : // Calculate normalized weights (with memory)
249 : } else {
250 0 : plumed_merror("analysis can now only support block averages");
251 : // if( !firstAnalysisDone )
252 : // finalizeWeightsNoLogSums( 1.0 );
253 : // else finalizeWeightsNoLogSums( old_norm );
254 : }
255 4 : }
256 :
257 : // void Analysis::finalizeWeightsNoLogSums( const double& onorm ){
258 : // if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong. Am trying to run analysis but I don't have sufficient data");
259 : // // Calculate normalization constant
260 : // double norm=0; for(unsigned i=0;i<logweights.size();++i) norm+=exp( logweights[i] );
261 : // // Calculate weights (with memory)
262 : // for(unsigned i=0;i<logweights.size();++i) data[i]->setWeight( exp( logweights[i] ) / norm );
263 : // }
264 :
265 0 : void Analysis::getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const {
266 : plumed_dbg_assert( getNumberOfAtoms()==0 );
267 0 : if( !reusing_data ) {
268 : plumed_dbg_assert( idata<logweights.size() && point.size()==getNumberOfArguments() );
269 0 : for(unsigned i=0; i<point.size(); ++i) point[i]=data[idata]->getReferenceArgument(i);
270 0 : weight=data[idata]->getWeight();
271 : } else {
272 0 : return mydatastash->getDataPoint( idata, point, weight );
273 : }
274 : }
275 :
276 4 : void Analysis::runAnalysis() {
277 :
278 : // Note : could add multiple walkers here - simply read in the data from all
279 : // other walkers here if we are writing the check points.
280 :
281 : // Calculate the final weights from the log weights
282 4 : if( !reusing_data ) {
283 4 : finalizeWeights( ignore_reweight );
284 : } else {
285 0 : mydatastash->finalizeWeights( ignore_reweight );
286 : // norm=mydatastash->retrieveNorm();
287 : }
288 : // This ensures everything is set up to run the calculation
289 : // if( single_run ) setAnalysisStride( single_run, freq );
290 : // And run the analysis
291 4 : performAnalysis(); idata=0;
292 : // Update total normalization constant
293 : // old_norm+=norm; firstAnalysisDone=true;
294 :
295 4 : }
296 :
297 1292 : void Analysis::performOperations( const bool& from_update ) {
298 1292 : accumulate();
299 1292 : if( freq>0 ) {
300 200 : if( getStep()>0 && getStep()%freq==0 ) runAnalysis();
301 396 : else if( idata==logweights.size() ) error("something has gone wrong. Probably a wrong initial time on restart");
302 : }
303 1292 : }
304 :
305 0 : bool Analysis::getPeriodicityInformation(const unsigned& i, std::string& dmin, std::string& dmax) {
306 0 : bool isperiodic=getPntrToArgument(i)->isPeriodic();
307 0 : if(isperiodic) getPntrToArgument(i)->getDomain(dmin,dmax);
308 0 : return isperiodic;
309 : }
310 :
311 3 : void Analysis::runFinalJobs() {
312 3 : if( freq>0 ) return;
313 2 : if( getNumberOfDataPoints()==0 ) error("no data is available for analysis");
314 2 : runAnalysis();
315 : }
316 :
317 : }
318 4839 : }
|
