diff --git a/developer-doc/AddingAnAnalysis.txt b/developer-doc/AddingAnAnalysis.txt new file mode 100644 index 0000000000000000000000000000000000000000..8ac30d91eb4ad417680738347aeb6c23d92fcae2 --- /dev/null +++ b/developer-doc/AddingAnAnalysis.txt @@ -0,0 +1,208 @@ +// This document is formatted for Doxygen +/** + +\page AddingAnAnalysis Implementing analysis methods in PLUMED + +Information on implementing methods that perform analysis on stored trajectory information eg dimensionality reduction + +Implementing methods for analysing trajectory data is more complex than implementing collective variables. +Consequently it is difficult to write a step by step guide like those we have written on implementing \ref AddingAColvar "colvars" or +\ref AddingAFunction "functions". Hence, this document tries to explain the things that we have considered and the way these +have been incorporated into the PLMD::analysis::AnalysisBase and PLMD::analysis::AnalysisWithDataCollection abstract base classes. +Hopefully this will provide some insight into our rationale in writing these parts of the code and will help you to understand how +any new analysis method can be implemented in the PLUMED code in a way that exploits those features that are already there. + +\section overview An overview of analysis methods in PLUMED + +There are two distinct ways in which one may wish to perform some form of analysis on a molecular dynamics trajectory. In the first method some quantity +is calculated for each of the atoms in the trajectory in each of the frames and this is then averaged over the whole trajectory. Velocity +autocorrelation functions or mean squared displacements are examples of forms of analysis of this type. The methods implemented in PLMD::analysis are +not of this type. These methods are designed to collect set of snapshots of the trajectory and to perform some analysis of these snapshots. +These trajectory snapshots might be the values of a particular set of collective variables for each of the frames, they might be the +instantaneous positions of the atoms in each frame or they might be some combination of the above. The assumption then when running one of these +analysis methods is that a representation (or snapshot) will be collected intermittently from the trajectory and then once a sufficiently large + collection of these snapshots are collected they will be analysed. + +\section antraj Analysis on the fly + +It is important to remember that PLUMED is primarily a code for performing biased molecular dynamics. The code's original purpose was +to be a plugin that could be easily added to a number of different molecular dynamics engines that allowed additional forces to be +incorporated when integrating the equations of motion. The command line tools that allows one to analyse trajectories during post +processing were added later as an afterthought. This consideration is particularly important when considering analysis algorithms +in this code because most analysis codes that are used in the community read in the trajectory and to do all their calculations +during post processing. The analysis functions that have been implemented in PLUMED can be used to post-process trajectories - you +simply make use of the command line tool driver - however, they can also be used to analyse trajectories on the fly. We believe this +is useful for a number of reasons: + +- Computers are becoming more powerful so it is possible to run simulations for much longer. At the same time, however, hard drives +space is at a premium and it is thus difficult to store these large trajectories for post-processing. If trajectories can be analysed +on the fly this presents less of a problem. +- A number of free energy methods (e.g. Gaussian mixture umbrella sampling, adaptive umbrella sampling and reconnaissance metadynamics) +work by performing a sophistacated analysis of the trajectory and then using the result from this analysis to design a bias for the +dynamics. +- Analysis methods implemented in PLUMED can take advantage of the many different collective variables that we have implemented in +this code and are thus extremely flexible implementations of these techniques. +- Analysis methods implemented in PLUMED make use of the PLUMED input syntax, which hopefully allows users already familiar with +PLUMED to get to grips with using these tools more rapidly. + +\section genphil General Phillosopy + +The aim in the PLMD::analysis and PLMD::dimred modules is to write code that is very flexible and that allows the user a great deal of +flexibility in the input. For example we would like to be able to write the code so that a user can collect data from the trajectory and +then at analysis time they can: + +- Select a subset of landmark points from the stored data +- Generate projections of these landmark points using sketch-map +- Project the remaining non-landmark points using the sketch-map projection generated and construct a histogram as a function of the sketch-map coordinates. + +Furthermore, we would like to be able to do all the above with a minimum of syntax for the user and with a minimum amount of code to maintain. +This is why the analysis class is structured in the way it is. The general idea is that one PLMD::analysis::AnalysisWithDataCollection object +collects the trajectory data as the simulation progresses (for the example above this would be an object of type PLMD::analysis::EuclideanDissimilarityMatrix). +Then when it is time to analyse a chain of analysis objects are run on the data collected by the PLMD::analysis::AnalysisWithDataCollection. +There are thus two types of analysis actions: + +- Those that inherit from PLMD::analysis::AnalysisWithDataCollection - these can collect data from a trajectory +- Those that inherit from PLMD::analysis::AnalysisBase - these cannot collect data from a trajectory. They get their input from another PLMD::analysis::AnalysisBase Action + +If you look at the code most of the routines in the PLMD::analysis and PLMD::dimred modules do not inherit from PLMD::analysis::AnalysisWithDataCollection. +In fact the only ones where this is an option that users really see in the manual are PLMD::analysis::Histogram and PLMD::analysis::EuclideanDissimilarityMatrix +(many of the other analysis actions that inherit from PLMD::analysis::AnalysisWithDataCollection only do so because this allows us to write straightforward +regression tests for this part of the code). The remaining analysis actions inherit from PLMD::analysis::AnalysisBase because in general they require some +dissimilarity information in order to function and because this information is calculated within PLMD::analysis::EuclideanDissimilarityMatrix. +Consequently, if you are writing a new analysis class it should probably not inherit from PLMD::analysis::AnalysisWithDataCollection. + +\section storing Storing trajectory data for analysis + +As discussed in the previous section the methods in PLMD::analysis all work by storing a snapshot of the trajectory every \f$N\f$ steps +and by then running an analysis method every \f$M\f$ steps, where \f$M\f$ is some substantial multiple of \f$N\f$. This intermittent +storing of trajectory frames and occasional analysis of the trajectory is all looked after within the PLMD::analysis::AnalysisWithDataCollection +abstract base class. Users can then set of a chain of analysis actions on the stored data by using actions that inherit from +PLMD::analysis::AnalysisBase. Any class inheriting from PLMD::analysis::AnalysisBase must have a method within it named performAnalysis(), +which will actually perform the analysis on the stored trajectory frames. When implementing a new analysis method the majority of your +development time will probably be spent implementing some part of this performAnalysis method. + +\section reweight Reweighting trajectories + +As discussed in previous sections PLUMED is primarily a code for doing biased molecular dynamics simulations. This bias is used to force rare events +to occur in the short time scales that are accessible within a molecular dynamics simulation. Oftentimes when analysing a trajectory we would like +to remove the contribution of the bias and to reweight the simulation so as to get the underlying free energy surface. When performing any analysis of +the trajectory one may similarly wish to remove the effect of the bias and to consider what weight each sampled point would have had if it had been sampled +in accordance with the underlying canonical probability distribution function. This process of reweighting points - of ascribing a weight to each snapshot +that discounts the effect of the simulation bias - is again looked after within PLMD::analysis::AnalysisWithDataCollection. If you wish to take these +weights into account in your analysis method you should use the method PLMD::analysis::AnalysisBase::getWeight to access them. Obviously, if you have +no simulation bias on the system then each point will have a weight of one and this will be the weight returned by PLMD::analysis::AnalysisBase::getWeight. + +\section output Outputting data files + +The fact that PLMD::analysis::AnalysisWithDataCollection can be used to run trajectory analysis in either post processing or on the fly during a trajectory means +that this class must also look after a number of things. For example one might wish to perform multiple analyses of the trajectory +during a simulation. Obviously, you do not want to overwrite the output file from your first analysis when you perform the second +analysis of the trajectory. In addition, you do not want to overwrite files from earlier runs if you choose to rerun your analysis +in a directory where you had already run an earlier calculation. For these reasons whenever you wish to read in the name of an output file +you should use the following code to make sure that any old files are backed up on restart: + +\verbatim +if( !getRestart() ){ OFile ofile; ofile.link(*this); ofile.setBackupString("analysis"); ofile.backupAllFiles(fname); } +\endverbatim + +where fname is the name of your output file. On top of this when you open an output file in your analysis method you should use the following +set of commands: + +\verbatim +OFile gfile; gfile.link(*this); +gfile.setBackupString("analysis"); +gfile.open( ofilename.c_str() ); +\endverbatim + +The second line ensures that files are named analysis.0.ofilename, analysis.1.ofilename and so on. Having said that it is probably best +to not write routines to output data in analysis classes and to simply ensure that you can pass any output data from your method to the +PLMD::analysis::OutputColvarFile and PLMD::analysis::OutputPDBFile methods. If you have done everything properly these classes should be +able to interact with your new class through methods that are in PLMD::analysis::AnalysisBase. + +\section metrics Calculating dissimilarity matrices + +One of the most basic ways in which you can analyse a collection of trajectory snapshots is to calculate the matrix of dissimilarities between each of the pairs +of trajectories frames. In plumed this is done in by the class PLMD::analysis::EuclideanDissimilarityMatrix. Notice as well that this class makes full use +of the PLMD::reference module that is discussed in \ref AddingAMetric and so this class alone allows you to calculate the dissimilairity between any pair of +trajectory frames in a wide variety of different ways. In addition, you can use PLMD::analysis::ReadDissimilarityMatrix to get the dissimilarities from an +input file. There should thus be little need to implement new objects for calculating dissimilarity matrices -- if you really feel you need something other than +what is already there or that is not implementable by \ref AddingAMetric then you are doing a calculation that is very unusual. + +\section landmarks Landmark selection algorithms + +Many analyses methods scale poorly with the number of trajectory frames that you wish to analyse. This happens in part because you need to compute the matrix +of pairwise disimiarities (a square matrix in which the number of columns is equal to the number of trajectory frames) but also because you then have to +do some algebra involving this matrix. To alleviate these problems a common strategy is to perform the analysis on a set of so-called landmark frames and to +then project the non-landmark snapshots from the trajectory using some out-of-sample extension of your analysis algorithm. Classes that inherit from +PLMD::analysis::LandmarkSelectionBase are implementations of the various landmark selection algorithms that are commonly employed. If your favourite landmark +selection algorithm is not there you may choose to implement a new landmark selection algorithm by writing a new PLMD::Action that inherits from this class. +Within this new class you need to only define a registerKeywords function, a constructor and a method that actually selects the landmarks that will be a function +that must be called selectLandmarks. Once you are satisfied that you want frame \f$k\f$ in your landmark set then you can select this using +PLMD::analysis::LandmarkSelectionBase::selectFrame. Please note that you can get the number of landmarks to select by calling: + +\verbatim +getNumberOfDataPoints() +\endverbatim + +If you would like to get the total number of frames from which you can get your subset of landmarks you should call: + +\verbatim +mydata->getNumberOfDataPoints() +\endverbatim + +Lastly, if you are using a method, which like PLMD::analysis::FarthestPointSampling uses the distances between your input points in some way, you should +add something akin to: + +\verbatim +if( !dissimilaritiesWereSet() ) error("dissimilarities have not been calcualted in input action"); +\endverbatim + +in your constructor as this will ensure that users are told what they are doing wrong if they do not specify how to calculate distances between points in the +input. To then get the dissimilirity between input point \f$i\f$ and input point \f$j\f$ use: + +\verbatim +mydata->getDissimilarity( landmarks[i], k ); +\endverbatim + +Calling PLMD::analysis::AnalysisBase::getDissimilarity will give you the distance between a pair of landmarks, which is not what you need. + +\section dimred Dimensionality reduction + +The aim when writing any large code such as PLUMED is to minise the number of lines of code as fewer lines of code means fewer bugs on average. +Hence, as explained in other sections of this developer manual, all the object oriented programming, inheritance and polymorphism. Given this +consider how we would go about implementing a library of dimensionality reduction algorithms. In LLE, ISOMAP, sketch-map or MDS the aim is to +generate a low-dimensional projection of some set of high-dimensional data points. For all these methods we can use the same code to to store +the high and low dimensional points and to output this data to a file. In fact the only things that differ in these various different methods are +the ways in which the dissimilarities between the high-dimensional points are calculated and the manner in which the low-dimensional projections +are generated. We have already discussed how PLUMED calculates matrices of dissimilarities between points using PLMD::analysis::EuclideanDissimilarityMatrix +and how one can think about introduce new methods of calculating dissimilarities. Priting to files meanwhile can be looked after by PLMD::analysis::OutputColvarFile +and PLMD::analysis::OutputPDBFile. Furthermore, if dimensionality reduction classes are written properly it is even possible to pass the projections +generated by them to PLMD::analysis::Histogram and to output histograms and free energy surfaces as a function of the low-dimensional coordinates. +As such in PLMD::dimred::DimensionalityReductionBase and its daughter classes we are thus solely concerned with calculating projections of data points. +The dissimilarities between the input high dimensional frames will always be calculated by some method akin to PLMD::analysis::EuclideanDissimilarityMatrix. +PLMD::dimred::DimensionalityReductionBase inherits from PLMD::analysis::AnalysisBase and expects another PLMD::analysis::AnalysisBase object as input. +Furthermore, with the exception of PLMD::dimred::ClassicalMultiDimensionalScaling the input PLMD::AnalysisBase must be an +PLMD::dimred::DimensionalityReductionBase as initial guesses must be suggested when using an interative optimization algorithm such as PLMD::dimred::SMACOF. + +Much of the work of dimensionality reduction done in the base PLMD::dimred::DimensionalityReductionBase class. When implementing any new +dimensionality reduction algorithm your focus will be on writing a PLMD::dimred::calculateProjections routines. This function will take as +input the matrix of pairwise dissimilarities between points (be they landmark points or otherwise) and is expected to return a matrix +containing the projections of the high-dimensional data points. Depending on the stress function you minimise to find projections you +may also have to implement PLMD::dimred::DimensionalityReductionBase::setTargetDistance and PLMD::dimred::DimensionalityReductionBase::calculateStress +functions. This is necessary with PLMD::dimred::SketchMapBase for example because sketch-map uses transformed versions of the dissimilarities +and distances in its stress function. These two routines are used within PLMD::dimred::ProjectNonLandmarkPoints which will product optimal projections of +points that were not selected as landmarks. + +\section cluster Clustering trajectory data + +There are currently no clustering methods implemented in the PLMD::analysis module. This section is thus here to explain how I (Gareth Tribello) imagined one +might go about implementing these methods. Again there are many commonalities between methods such as kmeans, Gaussian mixture models and so on, which should +be thought about when constructing an abstract base class. Furthermore, this abstract base class should (like PLMD::analysis::DimensionalityReductionBase) +inherit from PLMD::analysis::AnalysisBase and be implemented in a way that allows one to exploit the use of landmark selection algorithms that inherit from +PLMD::analysis::LandmarkSelectionBase and the measure actions such as PLMD::analysis::EuclideanDissimilarityMatrix. It should also be able to work with the +weights you get by reweiting the trajectories with the bias and so on. I don't think that you should have this inheriting from +PLMD::analysis::AnalysisWithDataCollection as I believe you want the clustering to work with projections of data generated by dimensionality reduction algorithms. +Obviously, if you are thinking of adding methods to cluster trajectory frames within PLUMED please feel free to get in touch with me (gareth.tribello\@gmail.com). +I will be more than happy to discuss these ideas with you. + +*/ diff --git a/developer-doc/Doxyfile b/developer-doc/Doxyfile index dcf741eb8b07197109ddb9ba8cd8268bf2abe79a..136576e65c3071242c35e15458ac1857b64bb124 100644 --- a/developer-doc/Doxyfile +++ b/developer-doc/Doxyfile @@ -742,7 +742,8 @@ INPUT = ./Intro.txt \ ./AddingAColvar.txt \ ./AddingAFunction.txt \ ./AddingACLTool.txt \ - ./AddingAMultColvar.txt + ./AddingAMultColvar.txt \ + ./AddingAnAnalysis.txt # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses diff --git a/regtest/analysis/rt-mds/Makefile b/regtest/analysis/rt-mds/Makefile deleted file mode 100644 index 3703b27cea227aa053fb6d1d73f861e4384dbcee..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/Makefile +++ /dev/null @@ -1 +0,0 @@ -include ../../scripts/test.make diff --git a/regtest/analysis/rt-mds/analysis.0.embed.reference b/regtest/analysis/rt-mds/analysis.0.embed.reference deleted file mode 100644 index 2e7d4808e54541161da11dea5fa411b231f08a1e..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/analysis.0.embed.reference +++ /dev/null @@ -1,401 +0,0 @@ -DESCRIPTION: results from classical mds analysis performed at time 5.0000 -REMARK WEIGHT=0.0054 CLASSICAL_MDS.1=0.0873 CLASSICAL_MDS.2=0.0013 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7551 c1.moment-3=1.3322 -END -REMARK WEIGHT=0.0285 CLASSICAL_MDS.1=0.0535 CLASSICAL_MDS.2=0.0044 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7434 c1.moment-3=1.3003 -END -REMARK WEIGHT=0.0359 CLASSICAL_MDS.1=0.0150 CLASSICAL_MDS.2=0.0044 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7333 c1.moment-3=1.2632 -END -REMARK WEIGHT=0.0136 CLASSICAL_MDS.1=-0.0136 CLASSICAL_MDS.2=0.0019 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7283 c1.moment-3=1.2349 -END -REMARK WEIGHT=0.0047 CLASSICAL_MDS.1=-0.0308 CLASSICAL_MDS.2=-0.0012 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7268 c1.moment-3=1.2175 -END -REMARK WEIGHT=0.0042 CLASSICAL_MDS.1=-0.0355 CLASSICAL_MDS.2=-0.0018 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7262 c1.moment-3=1.2128 -END -REMARK WEIGHT=0.0289 CLASSICAL_MDS.1=-0.0090 CLASSICAL_MDS.2=-0.0015 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7327 c1.moment-3=1.2385 -END -REMARK WEIGHT=0.0481 CLASSICAL_MDS.1=0.0145 CLASSICAL_MDS.2=0.0004 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7371 c1.moment-3=1.2616 -END -REMARK WEIGHT=0.0430 CLASSICAL_MDS.1=0.0191 CLASSICAL_MDS.2=0.0022 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7366 c1.moment-3=1.2666 -END -REMARK WEIGHT=0.1867 CLASSICAL_MDS.1=0.0278 CLASSICAL_MDS.2=0.0015 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7394 c1.moment-3=1.2748 -END -REMARK WEIGHT=0.0735 CLASSICAL_MDS.1=0.0091 CLASSICAL_MDS.2=-0.0004 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7365 c1.moment-3=1.2562 -END -REMARK WEIGHT=0.0012 CLASSICAL_MDS.1=-0.0359 CLASSICAL_MDS.2=-0.0012 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7255 c1.moment-3=1.2126 -END -REMARK WEIGHT=0.0006 CLASSICAL_MDS.1=-0.0540 CLASSICAL_MDS.2=0.0004 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7193 c1.moment-3=1.1955 -END -REMARK WEIGHT=0.0082 CLASSICAL_MDS.1=-0.0322 CLASSICAL_MDS.2=-0.0034 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7286 c1.moment-3=1.2156 -END -REMARK WEIGHT=0.0036 CLASSICAL_MDS.1=-0.0205 CLASSICAL_MDS.2=-0.0150 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7429 c1.moment-3=1.2239 -END -REMARK WEIGHT=0.0005 CLASSICAL_MDS.1=-0.0282 CLASSICAL_MDS.2=-0.0289 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7542 c1.moment-3=1.2128 -END -REMARK WEIGHT=0.0002 CLASSICAL_MDS.1=-0.0530 CLASSICAL_MDS.2=-0.0358 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7544 c1.moment-3=1.1871 -END -REMARK WEIGHT=0.0002 CLASSICAL_MDS.1=-0.0769 CLASSICAL_MDS.2=-0.0312 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7438 c1.moment-3=1.1652 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.0893 CLASSICAL_MDS.2=-0.0259 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7355 c1.moment-3=1.1546 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.0697 CLASSICAL_MDS.2=-0.0246 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7393 c1.moment-3=1.1739 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.0301 CLASSICAL_MDS.2=-0.0196 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7448 c1.moment-3=1.2133 -END -REMARK WEIGHT=0.0019 CLASSICAL_MDS.1=0.0093 CLASSICAL_MDS.2=-0.0072 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7430 c1.moment-3=1.2547 -END -REMARK WEIGHT=0.0160 CLASSICAL_MDS.1=0.0274 CLASSICAL_MDS.2=0.0009 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7400 c1.moment-3=1.2742 -END -REMARK WEIGHT=0.0740 CLASSICAL_MDS.1=0.0301 CLASSICAL_MDS.2=0.0016 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7400 c1.moment-3=1.2770 -END -REMARK WEIGHT=0.0150 CLASSICAL_MDS.1=0.0265 CLASSICAL_MDS.2=-0.0033 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7438 c1.moment-3=1.2723 -END -REMARK WEIGHT=0.0029 CLASSICAL_MDS.1=0.0433 CLASSICAL_MDS.2=-0.0053 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7501 c1.moment-3=1.2879 -END -REMARK WEIGHT=0.0061 CLASSICAL_MDS.1=0.0514 CLASSICAL_MDS.2=-0.0019 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7489 c1.moment-3=1.2967 -END -REMARK WEIGHT=0.0328 CLASSICAL_MDS.1=0.0222 CLASSICAL_MDS.2=0.0031 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7365 c1.moment-3=1.2698 -END -REMARK WEIGHT=0.0299 CLASSICAL_MDS.1=-0.0129 CLASSICAL_MDS.2=0.0061 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7245 c1.moment-3=1.2366 -END -REMARK WEIGHT=0.0072 CLASSICAL_MDS.1=-0.0322 CLASSICAL_MDS.2=0.0084 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7172 c1.moment-3=1.2186 -END -REMARK WEIGHT=0.0039 CLASSICAL_MDS.1=-0.0355 CLASSICAL_MDS.2=0.0080 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7168 c1.moment-3=1.2153 -END -REMARK WEIGHT=0.0045 CLASSICAL_MDS.1=-0.0337 CLASSICAL_MDS.2=0.0067 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7185 c1.moment-3=1.2168 -END -REMARK WEIGHT=0.0059 CLASSICAL_MDS.1=-0.0308 CLASSICAL_MDS.2=0.0047 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7211 c1.moment-3=1.2191 -END -REMARK WEIGHT=0.0125 CLASSICAL_MDS.1=-0.0253 CLASSICAL_MDS.2=0.0016 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7255 c1.moment-3=1.2235 -END -REMARK WEIGHT=0.0144 CLASSICAL_MDS.1=-0.0132 CLASSICAL_MDS.2=-0.0021 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7323 c1.moment-3=1.2343 -END -REMARK WEIGHT=0.0307 CLASSICAL_MDS.1=0.0040 CLASSICAL_MDS.2=-0.0031 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7377 c1.moment-3=1.2506 -END -REMARK WEIGHT=0.0437 CLASSICAL_MDS.1=0.0107 CLASSICAL_MDS.2=0.0011 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7354 c1.moment-3=1.2582 -END -REMARK WEIGHT=0.1212 CLASSICAL_MDS.1=0.0299 CLASSICAL_MDS.2=0.0060 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7357 c1.moment-3=1.2780 -END -REMARK WEIGHT=0.8415 CLASSICAL_MDS.1=0.0540 CLASSICAL_MDS.2=0.0070 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7410 c1.moment-3=1.3015 -END -REMARK WEIGHT=0.6449 CLASSICAL_MDS.1=0.0650 CLASSICAL_MDS.2=0.0039 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7468 c1.moment-3=1.3114 -END -REMARK WEIGHT=0.0450 CLASSICAL_MDS.1=0.0415 CLASSICAL_MDS.2=0.0025 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7421 c1.moment-3=1.2882 -END -REMARK WEIGHT=0.0093 CLASSICAL_MDS.1=0.0074 CLASSICAL_MDS.2=0.0033 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7325 c1.moment-3=1.2555 -END -REMARK WEIGHT=0.0134 CLASSICAL_MDS.1=-0.0185 CLASSICAL_MDS.2=0.0059 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7232 c1.moment-3=1.2312 -END -REMARK WEIGHT=0.0068 CLASSICAL_MDS.1=-0.0413 CLASSICAL_MDS.2=0.0073 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7158 c1.moment-3=1.2096 -END -REMARK WEIGHT=0.0035 CLASSICAL_MDS.1=-0.0548 CLASSICAL_MDS.2=0.0074 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7123 c1.moment-3=1.1966 -END -REMARK WEIGHT=0.0043 CLASSICAL_MDS.1=-0.0385 CLASSICAL_MDS.2=0.0041 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7197 c1.moment-3=1.2114 -END -REMARK WEIGHT=0.0021 CLASSICAL_MDS.1=-0.0383 CLASSICAL_MDS.2=-0.0006 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7243 c1.moment-3=1.2104 -END -REMARK WEIGHT=0.0007 CLASSICAL_MDS.1=-0.0584 CLASSICAL_MDS.2=-0.0008 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7192 c1.moment-3=1.1910 -END -REMARK WEIGHT=0.0034 CLASSICAL_MDS.1=-0.0584 CLASSICAL_MDS.2=0.0023 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7162 c1.moment-3=1.1918 -END -REMARK WEIGHT=0.0039 CLASSICAL_MDS.1=-0.0553 CLASSICAL_MDS.2=0.0090 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7105 c1.moment-3=1.1965 -END -REMARK WEIGHT=0.0004 CLASSICAL_MDS.1=-0.0772 CLASSICAL_MDS.2=0.0146 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6995 c1.moment-3=1.1768 -END -REMARK WEIGHT=0.0006 CLASSICAL_MDS.1=-0.0757 CLASSICAL_MDS.2=0.0157 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6988 c1.moment-3=1.1785 -END -REMARK WEIGHT=0.0031 CLASSICAL_MDS.1=-0.0644 CLASSICAL_MDS.2=0.0144 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7030 c1.moment-3=1.1891 -END -REMARK WEIGHT=0.0033 CLASSICAL_MDS.1=-0.0816 CLASSICAL_MDS.2=0.0141 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6988 c1.moment-3=1.1724 -END -REMARK WEIGHT=0.0016 CLASSICAL_MDS.1=-0.1045 CLASSICAL_MDS.2=0.0148 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6922 c1.moment-3=1.1504 -END -REMARK WEIGHT=0.0022 CLASSICAL_MDS.1=-0.0867 CLASSICAL_MDS.2=0.0111 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7004 c1.moment-3=1.1667 -END -REMARK WEIGHT=0.0035 CLASSICAL_MDS.1=-0.0365 CLASSICAL_MDS.2=0.0014 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7228 c1.moment-3=1.2126 -END -REMARK WEIGHT=0.0011 CLASSICAL_MDS.1=0.0113 CLASSICAL_MDS.2=-0.0067 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7431 c1.moment-3=1.2567 -END -REMARK WEIGHT=0.0038 CLASSICAL_MDS.1=0.0272 CLASSICAL_MDS.2=-0.0074 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7479 c1.moment-3=1.2718 -END -REMARK WEIGHT=0.0211 CLASSICAL_MDS.1=0.0381 CLASSICAL_MDS.2=-0.0012 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7448 c1.moment-3=1.2840 -END -REMARK WEIGHT=0.0278 CLASSICAL_MDS.1=0.0548 CLASSICAL_MDS.2=0.0021 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7459 c1.moment-3=1.3010 -END -REMARK WEIGHT=0.0246 CLASSICAL_MDS.1=0.0470 CLASSICAL_MDS.2=0.0012 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7447 c1.moment-3=1.2933 -END -REMARK WEIGHT=0.0024 CLASSICAL_MDS.1=0.0227 CLASSICAL_MDS.2=-0.0019 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7414 c1.moment-3=1.2690 -END -REMARK WEIGHT=0.0146 CLASSICAL_MDS.1=0.0132 CLASSICAL_MDS.2=-0.0022 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7392 c1.moment-3=1.2597 -END -REMARK WEIGHT=0.0126 CLASSICAL_MDS.1=0.0040 CLASSICAL_MDS.2=0.0004 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7343 c1.moment-3=1.2515 -END -REMARK WEIGHT=0.0031 CLASSICAL_MDS.1=-0.0045 CLASSICAL_MDS.2=0.0008 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7317 c1.moment-3=1.2434 -END -REMARK WEIGHT=0.0019 CLASSICAL_MDS.1=-0.0145 CLASSICAL_MDS.2=0.0006 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7293 c1.moment-3=1.2337 -END -REMARK WEIGHT=0.0042 CLASSICAL_MDS.1=-0.0116 CLASSICAL_MDS.2=0.0025 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7282 c1.moment-3=1.2370 -END -REMARK WEIGHT=0.0093 CLASSICAL_MDS.1=-0.0031 CLASSICAL_MDS.2=0.0044 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7286 c1.moment-3=1.2457 -END -REMARK WEIGHT=0.0099 CLASSICAL_MDS.1=-0.0148 CLASSICAL_MDS.2=0.0086 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7216 c1.moment-3=1.2355 -END -REMARK WEIGHT=0.0052 CLASSICAL_MDS.1=-0.0318 CLASSICAL_MDS.2=0.0109 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7149 c1.moment-3=1.2196 -END -REMARK WEIGHT=0.0030 CLASSICAL_MDS.1=-0.0166 CLASSICAL_MDS.2=0.0078 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7218 c1.moment-3=1.2336 -END -REMARK WEIGHT=0.0136 CLASSICAL_MDS.1=0.0154 CLASSICAL_MDS.2=0.0036 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7342 c1.moment-3=1.2633 -END -REMARK WEIGHT=0.0938 CLASSICAL_MDS.1=0.0451 CLASSICAL_MDS.2=0.0018 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7437 c1.moment-3=1.2915 -END -REMARK WEIGHT=0.1720 CLASSICAL_MDS.1=0.0364 CLASSICAL_MDS.2=0.0021 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7411 c1.moment-3=1.2832 -END -REMARK WEIGHT=0.0154 CLASSICAL_MDS.1=0.0125 CLASSICAL_MDS.2=0.0036 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7335 c1.moment-3=1.2606 -END -REMARK WEIGHT=0.0938 CLASSICAL_MDS.1=0.0330 CLASSICAL_MDS.2=0.0045 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7379 c1.moment-3=1.2806 -END -REMARK WEIGHT=1.0000 CLASSICAL_MDS.1=0.0605 CLASSICAL_MDS.2=0.0048 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7448 c1.moment-3=1.3072 -END -REMARK WEIGHT=0.3337 CLASSICAL_MDS.1=0.0417 CLASSICAL_MDS.2=0.0041 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7406 c1.moment-3=1.2889 -END -REMARK WEIGHT=0.0764 CLASSICAL_MDS.1=0.0115 CLASSICAL_MDS.2=0.0011 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7356 c1.moment-3=1.2590 -END -REMARK WEIGHT=0.1257 CLASSICAL_MDS.1=0.0257 CLASSICAL_MDS.2=-0.0018 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7421 c1.moment-3=1.2719 -END -REMARK WEIGHT=0.1215 CLASSICAL_MDS.1=0.0451 CLASSICAL_MDS.2=-0.0029 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7483 c1.moment-3=1.2904 -END -REMARK WEIGHT=0.0290 CLASSICAL_MDS.1=0.0387 CLASSICAL_MDS.2=-0.0032 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7469 c1.moment-3=1.2841 -END -REMARK WEIGHT=0.0162 CLASSICAL_MDS.1=0.0254 CLASSICAL_MDS.2=-0.0017 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7420 c1.moment-3=1.2716 -END -REMARK WEIGHT=0.0639 CLASSICAL_MDS.1=0.0311 CLASSICAL_MDS.2=-0.0001 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7418 c1.moment-3=1.2775 -END -REMARK WEIGHT=0.2186 CLASSICAL_MDS.1=0.0448 CLASSICAL_MDS.2=0.0006 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7448 c1.moment-3=1.2909 -END -REMARK WEIGHT=0.3093 CLASSICAL_MDS.1=0.0499 CLASSICAL_MDS.2=0.0020 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7448 c1.moment-3=1.2962 -END -REMARK WEIGHT=0.5713 CLASSICAL_MDS.1=0.0543 CLASSICAL_MDS.2=0.0037 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7443 c1.moment-3=1.3009 -END -REMARK WEIGHT=0.2377 CLASSICAL_MDS.1=0.0486 CLASSICAL_MDS.2=0.0054 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7412 c1.moment-3=1.2959 -END -REMARK WEIGHT=0.0425 CLASSICAL_MDS.1=0.0459 CLASSICAL_MDS.2=0.0051 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7407 c1.moment-3=1.2932 -END -REMARK WEIGHT=0.1714 CLASSICAL_MDS.1=0.0715 CLASSICAL_MDS.2=-0.0001 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7524 c1.moment-3=1.3166 -END -REMARK WEIGHT=0.0831 CLASSICAL_MDS.1=0.0764 CLASSICAL_MDS.2=-0.0063 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7596 c1.moment-3=1.3196 -END -REMARK WEIGHT=0.0194 CLASSICAL_MDS.1=0.0603 CLASSICAL_MDS.2=-0.0119 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7609 c1.moment-3=1.3027 -END -REMARK WEIGHT=0.0149 CLASSICAL_MDS.1=0.0378 CLASSICAL_MDS.2=-0.0153 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7583 c1.moment-3=1.2800 -END -REMARK WEIGHT=0.0217 CLASSICAL_MDS.1=0.0214 CLASSICAL_MDS.2=-0.0124 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7512 c1.moment-3=1.2649 -END -REMARK WEIGHT=0.0127 CLASSICAL_MDS.1=0.0063 CLASSICAL_MDS.2=-0.0091 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7441 c1.moment-3=1.2513 -END -REMARK WEIGHT=0.0063 CLASSICAL_MDS.1=-0.0095 CLASSICAL_MDS.2=-0.0016 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7327 c1.moment-3=1.2380 -END -REMARK WEIGHT=0.0067 CLASSICAL_MDS.1=-0.0259 CLASSICAL_MDS.2=0.0036 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7234 c1.moment-3=1.2235 -END -REMARK WEIGHT=0.0085 CLASSICAL_MDS.1=-0.0234 CLASSICAL_MDS.2=0.0061 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7217 c1.moment-3=1.2265 -END -REMARK WEIGHT=0.0069 CLASSICAL_MDS.1=0.0012 CLASSICAL_MDS.2=0.0060 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7282 c1.moment-3=1.2502 -END diff --git a/regtest/analysis/rt-mds/analysis.0.list_embed.reference b/regtest/analysis/rt-mds/analysis.0.list_embed.reference deleted file mode 100644 index a5fdcd999b09c16b758fa6803ca3388b314be2a8..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/analysis.0.list_embed.reference +++ /dev/null @@ -1,101 +0,0 @@ -#! FIELDS @17.1 @17.2 - 0.0873 0.0013 - 0.0535 0.0044 - 0.0150 0.0044 - -0.0136 0.0019 - -0.0308 -0.0012 - -0.0355 -0.0018 - -0.0090 -0.0015 - 0.0145 0.0004 - 0.0191 0.0022 - 0.0278 0.0015 - 0.0091 -0.0004 - -0.0359 -0.0012 - -0.0540 0.0004 - -0.0322 -0.0034 - -0.0205 -0.0150 - -0.0282 -0.0289 - -0.0530 -0.0358 - -0.0769 -0.0312 - -0.0893 -0.0259 - -0.0697 -0.0246 - -0.0301 -0.0196 - 0.0093 -0.0072 - 0.0274 0.0009 - 0.0301 0.0016 - 0.0265 -0.0033 - 0.0433 -0.0053 - 0.0514 -0.0019 - 0.0222 0.0031 - -0.0129 0.0061 - -0.0322 0.0084 - -0.0355 0.0080 - -0.0337 0.0067 - -0.0308 0.0047 - -0.0253 0.0016 - -0.0132 -0.0021 - 0.0040 -0.0031 - 0.0107 0.0011 - 0.0299 0.0060 - 0.0540 0.0070 - 0.0650 0.0039 - 0.0415 0.0025 - 0.0074 0.0033 - -0.0185 0.0059 - -0.0413 0.0073 - -0.0548 0.0074 - -0.0385 0.0041 - -0.0383 -0.0006 - -0.0584 -0.0008 - -0.0584 0.0023 - -0.0553 0.0090 - -0.0772 0.0146 - -0.0757 0.0157 - -0.0644 0.0144 - -0.0816 0.0141 - -0.1045 0.0148 - -0.0867 0.0111 - -0.0365 0.0014 - 0.0113 -0.0067 - 0.0272 -0.0074 - 0.0381 -0.0012 - 0.0548 0.0021 - 0.0470 0.0012 - 0.0227 -0.0019 - 0.0132 -0.0022 - 0.0040 0.0004 - -0.0045 0.0008 - -0.0145 0.0006 - -0.0116 0.0025 - -0.0031 0.0044 - -0.0148 0.0086 - -0.0318 0.0109 - -0.0166 0.0078 - 0.0154 0.0036 - 0.0451 0.0018 - 0.0364 0.0021 - 0.0125 0.0036 - 0.0330 0.0045 - 0.0605 0.0048 - 0.0417 0.0041 - 0.0115 0.0011 - 0.0257 -0.0018 - 0.0451 -0.0029 - 0.0387 -0.0032 - 0.0254 -0.0017 - 0.0311 -0.0001 - 0.0448 0.0006 - 0.0499 0.0020 - 0.0543 0.0037 - 0.0486 0.0054 - 0.0459 0.0051 - 0.0715 -0.0001 - 0.0764 -0.0063 - 0.0603 -0.0119 - 0.0378 -0.0153 - 0.0214 -0.0124 - 0.0063 -0.0091 - -0.0095 -0.0016 - -0.0259 0.0036 - -0.0234 0.0061 - 0.0012 0.0060 diff --git a/regtest/analysis/rt-mds/config b/regtest/analysis/rt-mds/config deleted file mode 100644 index 57d886c44c9b5b615e6e31bf117a02928408814b..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/config +++ /dev/null @@ -1,2 +0,0 @@ -type=simplemd - diff --git a/regtest/analysis/rt-mds/embed.reference b/regtest/analysis/rt-mds/embed.reference deleted file mode 100644 index b9a6a8a057f789ab3b966abd3e229839cdeca138..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/embed.reference +++ /dev/null @@ -1,401 +0,0 @@ -DESCRIPTION: results from classical mds analysis performed at time 10.0000 -REMARK WEIGHT=0.0073 CLASSICAL_MDS.1=0.0594 CLASSICAL_MDS.2=0.0053 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7340 c1.moment-3=1.2672 -END -REMARK WEIGHT=0.0102 CLASSICAL_MDS.1=0.0683 CLASSICAL_MDS.2=0.0034 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7371 c1.moment-3=1.2758 -END -REMARK WEIGHT=0.0032 CLASSICAL_MDS.1=0.0705 CLASSICAL_MDS.2=0.0006 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7401 c1.moment-3=1.2777 -END -REMARK WEIGHT=0.0085 CLASSICAL_MDS.1=0.0729 CLASSICAL_MDS.2=-0.0019 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7428 c1.moment-3=1.2796 -END -REMARK WEIGHT=0.0190 CLASSICAL_MDS.1=0.0874 CLASSICAL_MDS.2=-0.0069 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7497 c1.moment-3=1.2934 -END -REMARK WEIGHT=0.0162 CLASSICAL_MDS.1=0.0961 CLASSICAL_MDS.2=-0.0095 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7533 c1.moment-3=1.3018 -END -REMARK WEIGHT=0.0202 CLASSICAL_MDS.1=0.0867 CLASSICAL_MDS.2=-0.0062 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7489 c1.moment-3=1.2928 -END -REMARK WEIGHT=0.0262 CLASSICAL_MDS.1=0.0719 CLASSICAL_MDS.2=-0.0008 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7417 c1.moment-3=1.2788 -END -REMARK WEIGHT=0.1019 CLASSICAL_MDS.1=0.0575 CLASSICAL_MDS.2=0.0031 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7360 c1.moment-3=1.2651 -END -REMARK WEIGHT=0.0543 CLASSICAL_MDS.1=0.0552 CLASSICAL_MDS.2=0.0034 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7354 c1.moment-3=1.2628 -END -REMARK WEIGHT=0.0187 CLASSICAL_MDS.1=0.0682 CLASSICAL_MDS.2=-0.0002 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7405 c1.moment-3=1.2753 -END -REMARK WEIGHT=0.0250 CLASSICAL_MDS.1=0.0781 CLASSICAL_MDS.2=-0.0024 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7440 c1.moment-3=1.2847 -END -REMARK WEIGHT=0.0086 CLASSICAL_MDS.1=0.0785 CLASSICAL_MDS.2=-0.0025 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7441 c1.moment-3=1.2852 -END -REMARK WEIGHT=0.0031 CLASSICAL_MDS.1=0.0763 CLASSICAL_MDS.2=-0.0027 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7441 c1.moment-3=1.2829 -END -REMARK WEIGHT=0.0113 CLASSICAL_MDS.1=0.0755 CLASSICAL_MDS.2=-0.0041 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7454 c1.moment-3=1.2820 -END -REMARK WEIGHT=0.0598 CLASSICAL_MDS.1=0.0792 CLASSICAL_MDS.2=-0.0069 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7486 c1.moment-3=1.2853 -END -REMARK WEIGHT=0.0825 CLASSICAL_MDS.1=0.0778 CLASSICAL_MDS.2=-0.0062 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7477 c1.moment-3=1.2840 -END -REMARK WEIGHT=0.0659 CLASSICAL_MDS.1=0.0727 CLASSICAL_MDS.2=-0.0025 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7434 c1.moment-3=1.2795 -END -REMARK WEIGHT=0.1233 CLASSICAL_MDS.1=0.0726 CLASSICAL_MDS.2=0.0014 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7396 c1.moment-3=1.2798 -END -REMARK WEIGHT=0.2603 CLASSICAL_MDS.1=0.0756 CLASSICAL_MDS.2=0.0030 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7384 c1.moment-3=1.2830 -END -REMARK WEIGHT=0.1474 CLASSICAL_MDS.1=0.0821 CLASSICAL_MDS.2=0.0017 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7405 c1.moment-3=1.2893 -END -REMARK WEIGHT=0.0116 CLASSICAL_MDS.1=0.0816 CLASSICAL_MDS.2=0.0001 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7420 c1.moment-3=1.2886 -END -REMARK WEIGHT=0.1441 CLASSICAL_MDS.1=0.0774 CLASSICAL_MDS.2=0.0008 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7407 c1.moment-3=1.2845 -END -REMARK WEIGHT=0.2022 CLASSICAL_MDS.1=0.0766 CLASSICAL_MDS.2=0.0008 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7407 c1.moment-3=1.2837 -END -REMARK WEIGHT=0.1476 CLASSICAL_MDS.1=0.0838 CLASSICAL_MDS.2=-0.0005 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7429 c1.moment-3=1.2906 -END -REMARK WEIGHT=0.2285 CLASSICAL_MDS.1=0.0931 CLASSICAL_MDS.2=-0.0016 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7452 c1.moment-3=1.2997 -END -REMARK WEIGHT=0.4499 CLASSICAL_MDS.1=0.0936 CLASSICAL_MDS.2=-0.0010 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7446 c1.moment-3=1.3004 -END -REMARK WEIGHT=0.4344 CLASSICAL_MDS.1=0.0923 CLASSICAL_MDS.2=-0.0009 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7443 c1.moment-3=1.2990 -END -REMARK WEIGHT=0.3822 CLASSICAL_MDS.1=0.0906 CLASSICAL_MDS.2=-0.0008 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7440 c1.moment-3=1.2974 -END -REMARK WEIGHT=0.2334 CLASSICAL_MDS.1=0.0706 CLASSICAL_MDS.2=0.0025 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7382 c1.moment-3=1.2780 -END -REMARK WEIGHT=0.1038 CLASSICAL_MDS.1=0.0664 CLASSICAL_MDS.2=0.0050 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7352 c1.moment-3=1.2741 -END -REMARK WEIGHT=0.4459 CLASSICAL_MDS.1=0.0901 CLASSICAL_MDS.2=0.0019 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7412 c1.moment-3=1.2972 -END -REMARK WEIGHT=0.7086 CLASSICAL_MDS.1=0.1059 CLASSICAL_MDS.2=-0.0015 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7467 c1.moment-3=1.3125 -END -REMARK WEIGHT=0.1776 CLASSICAL_MDS.1=0.0917 CLASSICAL_MDS.2=0.0006 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7428 c1.moment-3=1.2986 -END -REMARK WEIGHT=0.1777 CLASSICAL_MDS.1=0.0876 CLASSICAL_MDS.2=0.0021 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7407 c1.moment-3=1.2948 -END -REMARK WEIGHT=0.1789 CLASSICAL_MDS.1=0.1008 CLASSICAL_MDS.2=-0.0001 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7446 c1.moment-3=1.3076 -END -REMARK WEIGHT=0.2332 CLASSICAL_MDS.1=0.0955 CLASSICAL_MDS.2=-0.0014 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7452 c1.moment-3=1.3022 -END -REMARK WEIGHT=0.1479 CLASSICAL_MDS.1=0.0858 CLASSICAL_MDS.2=-0.0018 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7444 c1.moment-3=1.2925 -END -REMARK WEIGHT=0.2978 CLASSICAL_MDS.1=0.0894 CLASSICAL_MDS.2=-0.0005 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7436 c1.moment-3=1.2963 -END -REMARK WEIGHT=0.7063 CLASSICAL_MDS.1=0.0922 CLASSICAL_MDS.2=0.0014 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7420 c1.moment-3=1.2992 -END -REMARK WEIGHT=1.0000 CLASSICAL_MDS.1=0.0917 CLASSICAL_MDS.2=0.0010 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7423 c1.moment-3=1.2987 -END -REMARK WEIGHT=0.2296 CLASSICAL_MDS.1=0.0925 CLASSICAL_MDS.2=-0.0043 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7477 c1.moment-3=1.2989 -END -REMARK WEIGHT=0.0617 CLASSICAL_MDS.1=0.0967 CLASSICAL_MDS.2=-0.0102 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7541 c1.moment-3=1.3023 -END -REMARK WEIGHT=0.3089 CLASSICAL_MDS.1=0.1073 CLASSICAL_MDS.2=-0.0096 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7548 c1.moment-3=1.3129 -END -REMARK WEIGHT=0.2276 CLASSICAL_MDS.1=0.1065 CLASSICAL_MDS.2=-0.0047 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7499 c1.moment-3=1.3126 -END -REMARK WEIGHT=0.0747 CLASSICAL_MDS.1=0.0935 CLASSICAL_MDS.2=-0.0003 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7439 c1.moment-3=1.3003 -END -REMARK WEIGHT=0.0472 CLASSICAL_MDS.1=0.0668 CLASSICAL_MDS.2=0.0049 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7353 c1.moment-3=1.2745 -END -REMARK WEIGHT=0.0408 CLASSICAL_MDS.1=0.0411 CLASSICAL_MDS.2=0.0087 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7284 c1.moment-3=1.2494 -END -REMARK WEIGHT=0.0419 CLASSICAL_MDS.1=0.0264 CLASSICAL_MDS.2=0.0110 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7242 c1.moment-3=1.2352 -END -REMARK WEIGHT=0.0119 CLASSICAL_MDS.1=0.0254 CLASSICAL_MDS.2=0.0125 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7226 c1.moment-3=1.2344 -END -REMARK WEIGHT=0.0125 CLASSICAL_MDS.1=0.0503 CLASSICAL_MDS.2=0.0085 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7297 c1.moment-3=1.2586 -END -REMARK WEIGHT=0.0217 CLASSICAL_MDS.1=0.0855 CLASSICAL_MDS.2=0.0007 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7419 c1.moment-3=1.2925 -END -REMARK WEIGHT=0.0392 CLASSICAL_MDS.1=0.0975 CLASSICAL_MDS.2=-0.0028 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7469 c1.moment-3=1.3039 -END -REMARK WEIGHT=0.2898 CLASSICAL_MDS.1=0.0983 CLASSICAL_MDS.2=-0.0035 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7477 c1.moment-3=1.3047 -END -REMARK WEIGHT=0.4917 CLASSICAL_MDS.1=0.0997 CLASSICAL_MDS.2=-0.0044 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7487 c1.moment-3=1.3059 -END -REMARK WEIGHT=0.0181 CLASSICAL_MDS.1=0.0877 CLASSICAL_MDS.2=-0.0021 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7449 c1.moment-3=1.2944 -END -REMARK WEIGHT=0.0013 CLASSICAL_MDS.1=0.0501 CLASSICAL_MDS.2=0.0053 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7328 c1.moment-3=1.2580 -END -REMARK WEIGHT=0.0036 CLASSICAL_MDS.1=0.0104 CLASSICAL_MDS.2=0.0140 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7192 c1.moment-3=1.2197 -END -REMARK WEIGHT=0.0076 CLASSICAL_MDS.1=-0.0199 CLASSICAL_MDS.2=0.0221 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7074 c1.moment-3=1.1907 -END -REMARK WEIGHT=0.0033 CLASSICAL_MDS.1=-0.0416 CLASSICAL_MDS.2=0.0277 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6990 c1.moment-3=1.1698 -END -REMARK WEIGHT=0.0002 CLASSICAL_MDS.1=-0.0688 CLASSICAL_MDS.2=0.0338 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6895 c1.moment-3=1.1436 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.0858 CLASSICAL_MDS.2=0.0354 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6858 c1.moment-3=1.1269 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.0897 CLASSICAL_MDS.2=0.0313 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6894 c1.moment-3=1.1226 -END -REMARK WEIGHT=0.0002 CLASSICAL_MDS.1=-0.0920 CLASSICAL_MDS.2=0.0254 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6950 c1.moment-3=1.1195 -END -REMARK WEIGHT=0.0004 CLASSICAL_MDS.1=-0.1168 CLASSICAL_MDS.2=0.0227 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6945 c1.moment-3=1.0946 -END -REMARK WEIGHT=0.0007 CLASSICAL_MDS.1=-0.1447 CLASSICAL_MDS.2=0.0227 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6910 c1.moment-3=1.0669 -END -REMARK WEIGHT=0.0031 CLASSICAL_MDS.1=-0.1636 CLASSICAL_MDS.2=0.0221 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6892 c1.moment-3=1.0481 -END -REMARK WEIGHT=0.0095 CLASSICAL_MDS.1=-0.1461 CLASSICAL_MDS.2=0.0182 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6953 c1.moment-3=1.0650 -END -REMARK WEIGHT=0.0030 CLASSICAL_MDS.1=-0.1212 CLASSICAL_MDS.2=0.0139 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7027 c1.moment-3=1.0891 -END -REMARK WEIGHT=0.0040 CLASSICAL_MDS.1=-0.0775 CLASSICAL_MDS.2=0.0056 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7164 c1.moment-3=1.1314 -END -REMARK WEIGHT=0.0032 CLASSICAL_MDS.1=-0.0008 CLASSICAL_MDS.2=-0.0041 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7358 c1.moment-3=1.2063 -END -REMARK WEIGHT=0.0057 CLASSICAL_MDS.1=0.0469 CLASSICAL_MDS.2=-0.0038 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7415 c1.moment-3=1.2537 -END -REMARK WEIGHT=0.0638 CLASSICAL_MDS.1=0.0771 CLASSICAL_MDS.2=-0.0004 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7419 c1.moment-3=1.2840 -END -REMARK WEIGHT=0.3314 CLASSICAL_MDS.1=0.0917 CLASSICAL_MDS.2=0.0005 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7428 c1.moment-3=1.2986 -END -REMARK WEIGHT=0.1319 CLASSICAL_MDS.1=0.0656 CLASSICAL_MDS.2=0.0042 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7359 c1.moment-3=1.2732 -END -REMARK WEIGHT=0.0139 CLASSICAL_MDS.1=-0.0034 CLASSICAL_MDS.2=0.0115 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7200 c1.moment-3=1.2057 -END -REMARK WEIGHT=0.0053 CLASSICAL_MDS.1=-0.0901 CLASSICAL_MDS.2=0.0141 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7064 c1.moment-3=1.1200 -END -REMARK WEIGHT=0.0040 CLASSICAL_MDS.1=-0.1784 CLASSICAL_MDS.2=0.0069 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7024 c1.moment-3=1.0315 -END -REMARK WEIGHT=0.0020 CLASSICAL_MDS.1=-0.2253 CLASSICAL_MDS.2=-0.0077 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7109 c1.moment-3=0.9832 -END -REMARK WEIGHT=0.0006 CLASSICAL_MDS.1=-0.2372 CLASSICAL_MDS.2=-0.0248 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7264 c1.moment-3=0.9692 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.2415 CLASSICAL_MDS.2=-0.0448 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7456 c1.moment-3=0.9624 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.2562 CLASSICAL_MDS.2=-0.0539 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7529 c1.moment-3=0.9466 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.2901 CLASSICAL_MDS.2=-0.0433 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7381 c1.moment-3=0.9144 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.3173 CLASSICAL_MDS.2=-0.0215 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7130 c1.moment-3=0.8902 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.3117 CLASSICAL_MDS.2=0.0096 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6829 c1.moment-3=0.8996 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.2787 CLASSICAL_MDS.2=0.0317 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6651 c1.moment-3=0.9352 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.2319 CLASSICAL_MDS.2=0.0372 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6656 c1.moment-3=0.9823 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.1897 CLASSICAL_MDS.2=0.0305 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6776 c1.moment-3=1.0233 -END -REMARK WEIGHT=0.0005 CLASSICAL_MDS.1=-0.1592 CLASSICAL_MDS.2=0.0142 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.6975 c1.moment-3=1.0515 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.1211 CLASSICAL_MDS.2=-0.0049 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7213 c1.moment-3=1.0868 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.0960 CLASSICAL_MDS.2=-0.0150 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7345 c1.moment-3=1.1105 -END -REMARK WEIGHT=0.0010 CLASSICAL_MDS.1=-0.0842 CLASSICAL_MDS.2=-0.0162 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7372 c1.moment-3=1.1221 -END -REMARK WEIGHT=0.0005 CLASSICAL_MDS.1=-0.0813 CLASSICAL_MDS.2=-0.0187 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7401 c1.moment-3=1.1246 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.0823 CLASSICAL_MDS.2=-0.0193 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7406 c1.moment-3=1.1235 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.0773 CLASSICAL_MDS.2=-0.0179 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7398 c1.moment-3=1.1287 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.0621 CLASSICAL_MDS.2=-0.0202 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7440 c1.moment-3=1.1435 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.0449 CLASSICAL_MDS.2=-0.0215 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7475 c1.moment-3=1.1604 -END -REMARK WEIGHT=0.0000 CLASSICAL_MDS.1=-0.0269 CLASSICAL_MDS.2=-0.0307 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7589 c1.moment-3=1.1770 -END -REMARK WEIGHT=0.0001 CLASSICAL_MDS.1=-0.0075 CLASSICAL_MDS.2=-0.0368 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7673 c1.moment-3=1.1955 -END -REMARK WEIGHT=0.0015 CLASSICAL_MDS.1=0.0266 CLASSICAL_MDS.2=-0.0342 -REMARK ARG=c1.moment-2,c1.moment-3 -REMARK c1.moment-2=0.7690 c1.moment-3=1.2297 -END diff --git a/regtest/analysis/rt-mds/in b/regtest/analysis/rt-mds/in deleted file mode 100644 index 030417936fdc6c7f1f4d0edff090893b1485f13d..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/in +++ /dev/null @@ -1,11 +0,0 @@ -inputfile input.xyz -outputfile output.xyz -temperature 0.2 -tstep 0.005 -friction 1 -forcecutoff 2.5 -listcutoff 3.0 -ndim 2 -nstep 2000 -nconfig 1000 trajectory.xyz -nstat 1000 energies.dat diff --git a/regtest/analysis/rt-mds/input.xyz b/regtest/analysis/rt-mds/input.xyz deleted file mode 100644 index 9e1e59174c1ce7a8d31fa8220d64ea05a9edbc14..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/input.xyz +++ /dev/null @@ -1,9 +0,0 @@ -7 -100. 100. 100. -Ar 7.3933470660 -2.6986483924 0.0000000000 -Ar 7.8226765198 -0.7390907295 0.0000000000 -Ar 7.1014969839 -1.6164766614 0.0000000000 -Ar 8.2357184242 -1.7097824975 0.0000000000 -Ar 6.7372520842 -0.5111536183 0.0000000000 -Ar 6.3777119489 -2.4640437401 0.0000000000 -Ar 5.9900631495 -1.3385375043 0.0000000000 diff --git a/regtest/analysis/rt-mds/list_embed.reference b/regtest/analysis/rt-mds/list_embed.reference deleted file mode 100644 index 6f40ac6e8c9694f04ba5de8a9359facec1e04047..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/list_embed.reference +++ /dev/null @@ -1,101 +0,0 @@ -#! FIELDS @17.1 @17.2 - 0.0594 0.0053 - 0.0683 0.0034 - 0.0705 0.0006 - 0.0729 -0.0019 - 0.0874 -0.0069 - 0.0961 -0.0095 - 0.0867 -0.0062 - 0.0719 -0.0008 - 0.0575 0.0031 - 0.0552 0.0034 - 0.0682 -0.0002 - 0.0781 -0.0024 - 0.0785 -0.0025 - 0.0763 -0.0027 - 0.0755 -0.0041 - 0.0792 -0.0069 - 0.0778 -0.0062 - 0.0727 -0.0025 - 0.0726 0.0014 - 0.0756 0.0030 - 0.0821 0.0017 - 0.0816 0.0001 - 0.0774 0.0008 - 0.0766 0.0008 - 0.0838 -0.0005 - 0.0931 -0.0016 - 0.0936 -0.0010 - 0.0923 -0.0009 - 0.0906 -0.0008 - 0.0706 0.0025 - 0.0664 0.0050 - 0.0901 0.0019 - 0.1059 -0.0015 - 0.0917 0.0006 - 0.0876 0.0021 - 0.1008 -0.0001 - 0.0955 -0.0014 - 0.0858 -0.0018 - 0.0894 -0.0005 - 0.0922 0.0014 - 0.0917 0.0010 - 0.0925 -0.0043 - 0.0967 -0.0102 - 0.1073 -0.0096 - 0.1065 -0.0047 - 0.0935 -0.0003 - 0.0668 0.0049 - 0.0411 0.0087 - 0.0264 0.0110 - 0.0254 0.0125 - 0.0503 0.0085 - 0.0855 0.0007 - 0.0975 -0.0028 - 0.0983 -0.0035 - 0.0997 -0.0044 - 0.0877 -0.0021 - 0.0501 0.0053 - 0.0104 0.0140 - -0.0199 0.0221 - -0.0416 0.0277 - -0.0688 0.0338 - -0.0858 0.0354 - -0.0897 0.0313 - -0.0920 0.0254 - -0.1168 0.0227 - -0.1447 0.0227 - -0.1636 0.0221 - -0.1461 0.0182 - -0.1212 0.0139 - -0.0775 0.0056 - -0.0008 -0.0041 - 0.0469 -0.0038 - 0.0771 -0.0004 - 0.0917 0.0005 - 0.0656 0.0042 - -0.0034 0.0115 - -0.0901 0.0141 - -0.1784 0.0069 - -0.2253 -0.0077 - -0.2372 -0.0248 - -0.2415 -0.0448 - -0.2562 -0.0539 - -0.2901 -0.0433 - -0.3173 -0.0215 - -0.3117 0.0096 - -0.2787 0.0317 - -0.2319 0.0372 - -0.1897 0.0305 - -0.1592 0.0142 - -0.1211 -0.0049 - -0.0960 -0.0150 - -0.0842 -0.0162 - -0.0813 -0.0187 - -0.0823 -0.0193 - -0.0773 -0.0179 - -0.0621 -0.0202 - -0.0449 -0.0215 - -0.0269 -0.0307 - -0.0075 -0.0368 - 0.0266 -0.0342 diff --git a/regtest/analysis/rt-mds/plumed.dat b/regtest/analysis/rt-mds/plumed.dat deleted file mode 100755 index 2d8d22c09ce508dc3313b4ebcface5d989294bea..0000000000000000000000000000000000000000 --- a/regtest/analysis/rt-mds/plumed.dat +++ /dev/null @@ -1,28 +0,0 @@ -UNITS NATURAL -COM ATOMS=1-7 LABEL=com -DISTANCE ATOMS=1,com LABEL=d1 -UPPER_WALLS ARG=d1 AT=2.0 KAPPA=100. -DISTANCE ATOMS=2,com LABEL=d2 -UPPER_WALLS ARG=d2 AT=2.0 KAPPA=100. -DISTANCE ATOMS=3,com LABEL=d3 -UPPER_WALLS ARG=d3 AT=2.0 KAPPA=100. -DISTANCE ATOMS=4,com LABEL=d4 -UPPER_WALLS ARG=d4 AT=2.0 KAPPA=100. -DISTANCE ATOMS=5,com LABEL=d5 -UPPER_WALLS ARG=d5 AT=2.0 KAPPA=100. -DISTANCE ATOMS=6,com LABEL=d6 -UPPER_WALLS ARG=d6 AT=2.0 KAPPA=100. -DISTANCE ATOMS=7,com LABEL=d7 -UPPER_WALLS ARG=d7 AT=2.0 KAPPA=100. - -COORDINATIONNUMBER SPECIES=1-7 MOMENTS=2-3 SWITCH={RATIONAL R_0=1.5 NN=8 MM=16} LABEL=c1 -CLASSICAL_MDS ... - ARG=c1.moment-2,c1.moment-3 - REWEIGHT_TEMP=0.1 TEMP=0.2 - STRIDE=10 - RUN=1000 - NLOW_DIM=2 - FMT=%8.4f - OUTPUT_FILE=list_embed - EMBEDDING_OFILE=embed -... CLASSICAL_MDS diff --git a/src/analysis/Analysis.cpp b/src/analysis/Analysis.cpp deleted file mode 100644 index 86c8c1cb59ed620f05de0dd87df744bf794b174e..0000000000000000000000000000000000000000 --- a/src/analysis/Analysis.cpp +++ /dev/null @@ -1,427 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2012-2015 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ - -#include "Analysis.h" -#include "core/ActionSet.h" -#include "core/ActionWithValue.h" -#include "core/PlumedMain.h" -#include "core/Atoms.h" -#include "tools/IFile.h" -#include "reference/ReferenceConfiguration.h" -#include "reference/ReferenceArguments.h" -#include "reference/ReferenceAtoms.h" -#include "reference/MetricRegister.h" - -namespace PLMD { -namespace analysis { - -//+PLUMEDOC INTERNAL reweighting -/* -Calculate free energies from a biassed/higher temperature trajectory. - -We can use our knowledge of the Boltzmann distribution in the cannonical ensemble to reweight the data -contained in trajectories. Using this procedure we can take trajectory at temperature \f$T_1\f$ and use it to -extract probabilities at a different temperature, \f$T_2\f$, using: - -\f[ -P(s',t) = \frac{ \sum_{t'}^t \delta( s(x) - s' ) \exp\left( +( \left[\frac{1}{T_1} - \frac{1}{T_2}\right] \frac{U(x,t')}{k_B} \right) }{ \sum_{t'}^t \exp\left( +\left[\frac{1}{T_1} - \frac{1}{T_2}\right] \frac{U(x,t')}{k_B} \right) } -\f] - -where \f$U(x,t')\f$ is the potential energy of the system. Alternatively, if a static or pseudo-static bias \f$V(x,t')\f$ is acting on -the system we can remove this bias and get the unbiased probability distribution using: - -\f[ -P(s',t) = \frac{ \sum_{t'}^t \delta( s(x) - s' ) \exp\left( +\frac{V(x,t')}{k_B T} \right) }{ \sum_t'^t \exp\left( +\frac{V(x,t')}{k_B T} \right) } -\f] - -*/ -//+ENDPLUMEDOC - -void Analysis::registerKeywords( Keywords& keys ){ - Action::registerKeywords( keys ); - ActionPilot::registerKeywords( keys ); - ActionAtomistic::registerKeywords( keys ); - ActionWithArguments::registerKeywords( keys ); - keys.use("ARG"); keys.reset_style("ARG","optional"); - keys.add("atoms","ATOMS","the atoms whose positions we are tracking for the purpose of analysing the data"); - keys.add("compulsory","METRIC","EUCLIDEAN","how are we measuring the distances between configurations"); - keys.add("compulsory","STRIDE","1","the frequency with which data should be stored for analysis"); - keys.addFlag("USE_ALL_DATA",false,"use the data from the entire trajectory to perform the analysis"); - keys.add("compulsory","RUN","the frequency with which to run the analysis algorithm. This is not required if you specify USE_ALL_DATA"); - keys.add("optional","FMT","the format that should be used in analysis output files"); - keys.addFlag("REWEIGHT_BIAS",false,"reweight the data using all the biases acting on the dynamics. For more information see \\ref reweighting."); - keys.add("optional","TEMP","the system temperature. This is required if you are reweighting or doing free energies."); - keys.add("optional","REWEIGHT_TEMP","reweight data from a trajectory at one temperature and output the probability " - "distribution at a second temperature. For more information see \\ref reweighting. " - "This is not possible during postprocessing."); - keys.addFlag("WRITE_CHECKPOINT",false,"write out a checkpoint so that the analysis can be restarted in a later run"); - keys.add("hidden","REUSE_DATA_FROM","eventually this will allow you to analyse the same set of data multiple times"); - keys.add("hidden","IGNORE_REWEIGHTING","this allows you to ignore any reweighting factors"); - keys.reserveFlag("NOMEMORY",false,"analyse each block of data separately"); - keys.use("RESTART"); - keys.use("UPDATE_FROM"); - keys.use("UPDATE_UNTIL"); - ActionWithVessel::registerKeywords( keys ); keys.remove("TOL"); -} - -Analysis::Analysis(const ActionOptions&ao): -Action(ao), -ActionPilot(ao), -ActionAtomistic(ao), -ActionWithArguments(ao), -ActionWithVessel(ao), -single_run(false), -nomemory(true), -write_chq(false), -reusing_data(false), -ignore_reweight(false), -needeng(false), -idata(0), -firstAnalysisDone(false), -old_norm(0.0), -ofmt("%f"), -current_args(getNumberOfArguments()), -argument_names(getNumberOfArguments()) -{ - parse("FMT",ofmt); // Read the format for output files - // Make a vector containing all the argument names - for(unsigned i=0;i<getNumberOfArguments();++i) argument_names[i]=getPntrToArgument(i)->getName(); - // Read in the metric style - parse("METRIC",metricname); std::vector<AtomNumber> atom_numbers; - ReferenceConfiguration* checkref=metricRegister().create<ReferenceConfiguration>( metricname ); - // Check if we should read atoms - ReferenceAtoms* hasatoms=dynamic_cast<ReferenceAtoms*>( checkref ); - if( hasatoms ){ - parseAtomList("ATOMS",atom_numbers); requestAtoms(atom_numbers); - log.printf(" monitoring positions of atoms "); - for(unsigned i=0;i<atom_numbers.size();++i) log.printf("%d ",atom_numbers[i].serial() ); - log.printf("\n"); - } - // Check if we should read arguments - ReferenceArguments* hasargs=dynamic_cast<ReferenceArguments*>( checkref ); - if( !hasargs && getNumberOfArguments()!=0 ) error("use of arguments with metric type " + metricname + " is invalid"); - if( hasatoms && hasargs ) error("currently dependencies break if you have both arguments and atoms"); - // And delte the fake reference we created - delete checkref; - - std::string prev_analysis; parse("REUSE_DATA_FROM",prev_analysis); - if( prev_analysis.length()>0 ){ - reusing_data=true; - mydatastash=plumed.getActionSet().selectWithLabel<Analysis*>( prev_analysis ); - if( !mydatastash ) error("could not find analysis action named " + prev_analysis ); - parseFlag("IGNORE_REWEIGHTING",ignore_reweight); - if( ignore_reweight ) log.printf(" reusing data stored by %s but ignoring all reweighting\n",prev_analysis.c_str() ); - else log.printf(" reusing data stored by %s\n",prev_analysis.c_str() ); - } else { - if( keywords.exists("REWEIGHT_BIAS") ){ - bool dobias; parseFlag("REWEIGHT_BIAS",dobias); - if( dobias ){ - std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); - if( all.empty() ) error("your input file is not telling plumed to calculate anything"); - std::vector<Value*> arg( getArguments() ); - log.printf(" reweigting using the following biases "); - for(unsigned j=0;j<all.size();j++){ - std::string flab; flab=all[j]->getLabel() + ".bias"; - if( all[j]->exists(flab) ){ - biases.push_back( all[j]->copyOutput(flab) ); - arg.push_back( all[j]->copyOutput(flab) ); - log.printf(" %s",flab.c_str()); - } - } - log.printf("\n"); - if( biases.empty() ) error("you are asking to reweight bias but there does not appear to be a bias acting on your system"); - requestArguments( arg ); - } - } - - rtemp=0; - if( keywords.exists("REWEIGHT_TEMP") ) parse("REWEIGHT_TEMP",rtemp); - if( rtemp!=0 ){ - needeng=true; - log.printf(" reweighting simulation to probabilities at temperature %f\n",rtemp); - rtemp*=plumed.getAtoms().getKBoltzmann(); - } - simtemp=0.; parse("TEMP",simtemp); - if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann(); - else simtemp=plumed.getAtoms().getKbT(); - - if( rtemp>0 || !biases.empty() ){ - if(simtemp==0) error("The MD engine does not pass the temperature to plumed so you have to specify it using TEMP"); - } - - parseFlag("USE_ALL_DATA",single_run); - if( !single_run ){ - parse("RUN",freq ); - log.printf(" running analysis every %u steps\n",freq); - if( freq%getStride()!= 0 ) error("Frequncy of running is not a multiple of the stride"); - ndata=freq/getStride(); - data.resize( ndata ); - for(unsigned i=0;i<ndata;++i){ - data[i]=metricRegister().create<ReferenceConfiguration>( metricname ); - data[i]->setNamesAndAtomNumbers( atom_numbers, argument_names ); - } - logweights.resize( ndata ); - } else { - log.printf(" analyzing all data in trajectory\n"); - } - if( keywords.exists("NOMEMORY") ){ nomemory=false; parseFlag("NOMEMORY",nomemory); } - if(nomemory) log.printf(" doing a separate analysis for each block of data\n"); - parseFlag("WRITE_CHECKPOINT",write_chq); - if( write_chq && single_run ){ - write_chq=false; - warning("ignoring WRITE_CHECKPOINT flag because we are analyzing all data"); - } - - // We need no restart file if we are just collecting data and analyzing all of it - std::string filename = getName() + "_" + getLabel() + ".chkpnt"; - if( write_chq ) rfile.link(*this); - if( getRestart() ){ - if( single_run ) error("cannot restart histogram when using the USE_ALL_DATA option"); - if( !write_chq ) warning("restarting without writing a checkpoint file is somewhat strange"); - // Read in data from input file - readDataFromFile( filename ); - // Setup the restart file (append mode) - if( write_chq ) rfile.open( filename.c_str() ); // In append mode automatically because of restart - // Run the analysis if we stoped in the middle of it last time - log.printf(" restarting analysis with %u points read from restart file\n",idata); - } else if( write_chq ){ - // Setup the restart file (delete any old one) - rfile.open( filename.c_str() ); // In overwrite mode automatically because there is no restart - } - if( write_chq ){ - rfile.addConstantField("old_normalization"); - for(unsigned i=0;i<getNumberOfArguments();++i) rfile.setupPrintValue( getPntrToArgument(i) ); - } - } -} - -void Analysis::readDataFromFile( const std::string& filename ){ - FILE* fp=fopen(filename.c_str(),"r"); double tstep, oldtstep; - if(fp!=NULL){ - bool do_read=true, first=true; - while (do_read) { - PDB mypdb; - do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()); - if(do_read){ - data[idata]->set( mypdb ); - data[idata]->parse("TIME",tstep); - if( !first && ((tstep-oldtstep) - getStride()*plumed.getAtoms().getTimeStep())>plumed.getAtoms().getTimeStep() ){ - error("frequency of data storage in " + filename + " is not equal to frequency of data storage plumed.dat file"); - } - data[idata]->parse("LOG_WEIGHT",logweights[idata]); - data[idata]->parse("OLD_NORM",old_norm); - data[idata]->checkRead(); - idata++; first=false; oldtstep=tstep; - } else{ - break; - } - } - fclose(fp); - } - if(old_norm>0) firstAnalysisDone=true; -} - -void Analysis::parseOutputFile( const std::string& key, std::string& filename ){ - parse(key,filename); - if(filename=="dont output") return; - - if( !getRestart() ){ - OFile ofile; ofile.link(*this); - ofile.setBackupString("analysis"); - ofile.backupAllFiles(filename); - } -} - -void Analysis::prepare(){ - if(needeng) plumed.getAtoms().setCollectEnergy(true); -} - -void Analysis::calculate(){ -// do nothing -} - -void Analysis::accumulate(){ - // Don't store the first step (also don't store if we are getting data from elsewhere) - if( (!single_run && getStep()==0) || reusing_data ) return; - // This is used when we have a full quota of data from the first run - if( !single_run && idata==logweights.size() ) return; - - // Retrieve the bias - double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get(); - - double ww=0; - if(needeng){ - double energy=plumed.getAtoms().getEnergy()+bias; - // Reweighting because of temperature difference - ww=-( (1.0/rtemp) - (1.0/simtemp) )*(energy+bias); - } - // Reweighting because of biases - if( !biases.empty() ) ww += bias/ simtemp; - - // Get the arguments ready to transfer to reference configuration - for(unsigned i=0;i<getNumberOfArguments();++i) current_args[i]=getArgument(i); - - if(single_run){ - data.push_back( metricRegister().create<ReferenceConfiguration>( metricname ) ); - plumed_dbg_assert( data.size()==idata+1 ); - data[idata]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names ); - data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() ); - logweights.push_back(ww); - } else { - // Get the arguments and store them in a vector of vectors - data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() ); - logweights[idata] = ww; - } - - // Write data to checkpoint file - if( write_chq ){ - rfile.rewind(); - data[idata]->print( rfile, getTime(), logweights[idata], old_norm ); - rfile.flush(); - } - // Increment data counter - idata++; -} - -Analysis::~Analysis(){ - for(unsigned i=0;i<data.size();++i ) delete data[i]; - if( write_chq ) rfile.close(); -} - -std::vector<double> Analysis::getMetric() const { - // Add more exotic metrics in here -- FlexibleHill for instance - std::vector<double> empty; - if( metricname=="EUCLIDEAN" ){ - empty.resize( getNumberOfArguments(), 1.0 ); - } - return empty; -} - -double Analysis::getWeight( const unsigned& idata ) const { - if( !reusing_data ){ - plumed_dbg_assert( idata<data.size() ); - return data[idata]->getWeight(); - } else { - return mydatastash->getWeight(idata); - } -} - -void Analysis::finalizeWeights( const bool& ignore_weights ){ - // Check that we have the correct ammount of data - if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong. Am trying to run analysis but I don't have sufficient data"); - - norm=0; // Reset normalization constant - if( ignore_weights ){ - for(unsigned i=0;i<logweights.size();++i){ - data[i]->setWeight(1.0); norm+=1.0; - } - } else if( nomemory ){ - // Find the maximum weight - double maxweight=logweights[0]; - for(unsigned i=1;i<getNumberOfDataPoints();++i){ - if(logweights[i]>maxweight) maxweight=logweights[i]; - } - // Calculate normalization constant - for(unsigned i=0;i<logweights.size();++i){ - norm+=exp( logweights[i]-maxweight ); - } - // Calculate weights (no memory) - for(unsigned i=0;i<logweights.size();++i){ - data[i]->setWeight( exp( logweights[i]-maxweight ) ); - } - // Calculate normalized weights (with memory) - } else { - // Calculate normalization constant - for(unsigned i=0;i<logweights.size();++i){ - norm+=exp( logweights[i] ); - } - if( !firstAnalysisDone ) old_norm=1.0; - // Calculate weights (with memory) - for(unsigned i=0;i<logweights.size();++i){ - data[i]->setWeight( exp( logweights[i] ) / old_norm ); - } - if( !firstAnalysisDone ) old_norm=0.0; - } - -} - -void Analysis::getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const { - plumed_dbg_assert( getNumberOfAtoms()==0 ); - if( !reusing_data ){ - plumed_dbg_assert( idata<logweights.size() && point.size()==getNumberOfArguments() ); - for(unsigned i=0;i<point.size();++i) point[i]=data[idata]->getReferenceArgument(i); - weight=data[idata]->getWeight(); - } else { - return mydatastash->getDataPoint( idata, point, weight ); - } -} - -void Analysis::runAnalysis(){ - - // Note : could add multiple walkers here - simply read in the data from all - // other walkers here if we are writing the check points. - - // Calculate the final weights from the log weights - if( !reusing_data ){ - finalizeWeights( ignore_reweight ); - } else { - mydatastash->finalizeWeights( ignore_reweight ); - norm=mydatastash->retrieveNorm(); - } - // And run the analysis - performAnalysis(); idata=0; - // Update total normalization constant - old_norm+=norm; firstAnalysisDone=true; - -} - -double Analysis::getNormalization() const { - if( nomemory || !firstAnalysisDone ) return norm; - return ( 1. + norm/old_norm ); -} - -double Analysis::getTemp() const { - return simtemp; -} - -void Analysis::update(){ - accumulate(); - if( !single_run ){ - if( getStep()>0 && getStep()%freq==0 ) runAnalysis(); - else if( idata==logweights.size() ) error("something has gone wrong. Probably a wrong initial time on restart"); - } -} - -bool Analysis::getPeriodicityInformation(const unsigned& i, std::string& dmin, std::string& dmax){ - bool isperiodic=getPntrToArgument(i)->isPeriodic(); - if(isperiodic) getPntrToArgument(i)->getDomain(dmin,dmax); - return isperiodic; -} - -void Analysis::runFinalJobs() { - if( !single_run ) return; - if( getNumberOfDataPoints()==0 ) error("no data is available for analysis"); - runAnalysis(); -} - -} -} diff --git a/src/analysis/Analysis.h b/src/analysis/Analysis.h deleted file mode 100644 index 33a01f9726c1be0dd4253b4e9f7468c9c98563f7..0000000000000000000000000000000000000000 --- a/src/analysis/Analysis.h +++ /dev/null @@ -1,216 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2012-2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifndef __PLUMED_analysis_Analysis_h -#define __PLUMED_analysis_Analysis_h - -#include "core/ActionPilot.h" -#include "core/ActionAtomistic.h" -#include "core/ActionWithArguments.h" -#include "vesselbase/ActionWithVessel.h" - -#define PLUMED_ANALYSIS_INIT(ao) Action(ao),Analysis(ao) - -namespace PLMD { - -class ReferenceConfiguration; - -namespace analysis { - -/** -\ingroup INHERIT -This is the abstract base class to use for implementing new methods for analyzing the trajectory, within it there -is information as to how to go about implementing a new analysis method. - -*/ - -class Analysis : - public ActionPilot, - public ActionAtomistic, - public ActionWithArguments, - public vesselbase::ActionWithVessel - { -private: -/// Are we running only once for the whole trajectory - bool single_run; -/// Are we treating each block of data separately - bool nomemory; -/// Are we writing a checkpoint file - bool write_chq; -/// Are we reusing data stored by another analysis action - bool reusing_data; -/// If we are reusing data are we ignoring the reweighting in that data - bool ignore_reweight; -/// The Analysis action that we are reusing data from - Analysis* mydatastash; -/// The frequency with which we are performing analysis - unsigned freq; -/// Number of data point we need to run analysis - unsigned ndata; -/// The temperature at which we are running the calculation - double simtemp; -/// The temperature at which we want the histogram - double rtemp; -/// Do we need the energy (are we reweighting at a different temperature) - bool needeng; -/// The biases we are using in reweighting and the args we store them separately - std::vector<Value*> biases; -/// The piece of data we are inserting - unsigned idata; -/// The weights of all the data points - std::vector<double> logweights; -/// Have we analyzed the data for the first time - bool firstAnalysisDone; -/// The value of the old normalization constant - double norm, old_norm; -/// The format to use in output files - std::string ofmt; -/// Tempory vector to store values of arguments - std::vector<double> current_args; -/// List of argument names - std::vector<std::string> argument_names; -/// The type of metric we are using to measure distances - std::string metricname; -/// The checkpoint file - OFile rfile; -/// Read in data from a file - void readDataFromFile( const std::string& filename ); -/// This retrieves the value of norm from the analysis action. -/// It is only used to transfer data from one analysis action to -/// another. You should never need to use it. If you think you need it -/// you probably need getNormalization() - double retrieveNorm() const ; -/// Get the metric if we are using malonobius distance and flexible hill - std::vector<double> getMetric() const ; -protected: -/// This is used to read in output file names for analysis methods. When -/// this method is used and the calculation is not restarted old analysis -/// files are backed up. - void parseOutputFile( const std::string& key, std::string& filename ); -/// The data we are going to analyze - std::vector<ReferenceConfiguration*> data; -/// Get the name of the metric we are using to measure distances - std::string getMetricName() const ; -/// Return the number of arguments (this overwrites the one in ActionWithArguments) - unsigned getNumberOfArguments() const; -/// Return the number of data points - unsigned getNumberOfDataPoints() const; -/// Return the weight of the ith point - double getWeight( const unsigned& idata ) const ; -/// Retrieve the ith point - void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ; -/// Returns true if argument i is periodic together with the domain - bool getPeriodicityInformation(const unsigned& i, std::string& dmin, std::string& dmax); -/// Return the normalization constant - double getNormalization() const; -/// Return the set temperature (N.B. k_B T is what is returned by this function) - double getTemp () const; -/// Are we analyzing each data block separately (if we are not this also returns the old normalization ) - bool usingMemory() const; -/// Convert the stored log weights to proper weights - void finalizeWeights( const bool& ignore_weights ); -/// Overwrite ActionWithArguments getArguments() so that we don't return -/// the bias - std::vector<Value*> getArguments(); -/// Return the format to use for numbers in output files - std::string getOutputFormat() const ; -public: - static void registerKeywords( Keywords& keys ); - Analysis(const ActionOptions&); - ~Analysis(); - void prepare(); - void calculate(); - void update(); - void accumulate(); - virtual void performAnalysis()=0; - void apply(){} - void runFinalJobs(); - void runAnalysis(); - void lockRequests(); - void unlockRequests(); - void calculateNumericalDerivatives( ActionWithValue* a=NULL ){ plumed_error(); } - bool isPeriodic(){ plumed_error(); return false; } - unsigned getNumberOfDerivatives(){ plumed_error(); return 0; } -}; - -inline -std::string Analysis::getMetricName() const { - return metricname; -} - -inline -void Analysis::lockRequests(){ - ActionAtomistic::lockRequests(); - ActionWithArguments::lockRequests(); -} - -inline -void Analysis::unlockRequests(){ - ActionAtomistic::unlockRequests(); - ActionWithArguments::unlockRequests(); -} - -inline -unsigned Analysis::getNumberOfDataPoints() const { - if( !reusing_data ){ - plumed_dbg_assert( data.size()==logweights.size() ); - return data.size(); - } else { - return mydatastash->getNumberOfDataPoints(); - } -} - -inline -bool Analysis::usingMemory() const { - if( !reusing_data ){ - return !nomemory; - } else { - return mydatastash->usingMemory(); - } -} - -inline -unsigned Analysis::getNumberOfArguments() const { - unsigned nargs=ActionWithArguments::getNumberOfArguments(); - return nargs - biases.size(); -} - -inline -double Analysis::retrieveNorm() const { - return norm; -} - -inline -std::vector<Value*> Analysis::getArguments(){ - std::vector<Value*> arg_vals( ActionWithArguments::getArguments() ); - for(unsigned i=0;i<biases.size();++i) arg_vals.erase(arg_vals.end()-1); - return arg_vals; -} - -inline -std::string Analysis::getOutputFormat() const { - return ofmt; -} - -} -} - -#endif diff --git a/src/analysis/AnalysisBase.cpp b/src/analysis/AnalysisBase.cpp new file mode 100644 index 0000000000000000000000000000000000000000..26fb94e17540d6d70c2d41ec4858ffb8198f015a --- /dev/null +++ b/src/analysis/AnalysisBase.cpp @@ -0,0 +1,79 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2012-2015 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed-code.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see <http://www.gnu.org/licenses/>. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "AnalysisBase.h" +#include "core/PlumedMain.h" +#include "core/ActionSet.h" + +namespace PLMD { +namespace analysis { + +void AnalysisBase::registerKeywords( Keywords& keys ){ + Action::registerKeywords( keys ); + ActionPilot::registerKeywords( keys ); + ActionAtomistic::registerKeywords( keys ); + ActionWithArguments::registerKeywords( keys ); + ActionWithVessel::registerKeywords( keys ); keys.remove("TOL"); keys.remove("LOWMEM"); + keys.add("atoms-2","USE_OUTPUT_DATA_FROM","use the ouput of the analysis performed by this object as input to your new analysis object"); +} + +AnalysisBase::AnalysisBase(const ActionOptions&ao): +Action(ao), +ActionPilot(ao), +ActionAtomistic(ao), +ActionWithArguments(ao), +ActionWithVessel(ao), +mydata(NULL) +{ + // We have an if statement here so that this doesn't break with READ_DISSIMILARITIES + std::string datastr; if( keywords.exists("USE_OUTPUT_DATA_FROM") ) parse("USE_OUTPUT_DATA_FROM",datastr); + if( keywords.exists("USE_OUTPUT_DATA_FROM") && + datastr.length()==0 && + !keywords.exists("REUSE_INPUT_DATA_FROM") ) error("input analysis action was not specified use USE_OUTPUT_DATA_FROM"); + if( datastr.length()>0 ){ + mydata=plumed.getActionSet().selectWithLabel<AnalysisBase*>( datastr ); + log.printf(" performing analysis on output from %s \n",datastr.c_str() ); + if( !mydata ) error("could not find analysis action named " + datastr ); + freq=mydata->freq; use_all_data=mydata->use_all_data; + if( !use_all_data ) setStride( freq ); + } +} + +void AnalysisBase::update(){ + // Do nothing if we are just analysing at the end of the calculation + if( use_all_data ) return ; + // Check that we are on step + plumed_dbg_assert( getStep()%freq==0 ); + // And do the analysis + if( getStep()>0 ) performAnalysis(); +} + +void AnalysisBase::runFinalJobs(){ + // Nothing to do if we are not analysing all the data in the trajectory + if( !use_all_data ) return; + // Erm ... user has done something weird + if( getNumberOfDataPoints()==0 ) error("no data is available for analysis"); + // And do the analysis + if( use_all_data ) performAnalysis(); +} + +} +} diff --git a/src/analysis/AnalysisBase.h b/src/analysis/AnalysisBase.h new file mode 100644 index 0000000000000000000000000000000000000000..cfa1a613a1d21ece74eaf3c2a38b91cfeae1c298 --- /dev/null +++ b/src/analysis/AnalysisBase.h @@ -0,0 +1,166 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2012-2014 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed-code.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see <http://www.gnu.org/licenses/>. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_analysis_AnalysisBase_h +#define __PLUMED_analysis_AnalysisBase_h + +#include "core/ActionPilot.h" +#include "core/ActionAtomistic.h" +#include "core/ActionWithArguments.h" +#include "vesselbase/ActionWithVessel.h" + +namespace PLMD { + +class ReferenceConfiguration; + +namespace analysis { + +/** +\ingroup INHERIT +This is the abstract base class to use for implementing new methods for analyzing the trajectory. You can find +\ref AddingAnAnalysis "information" on how to use it to implement new analysis methods here. + +*/ + +class AnalysisBase : + public ActionPilot, + public ActionAtomistic, + public ActionWithArguments, + public vesselbase::ActionWithVessel + { +friend class ReadDissimilarityMatrix; +friend class AnalysisWithDataCollection; +private: +/// Just run the analysis a single time + bool use_all_data; +/// The frequency with which we are performing analysis + unsigned freq; +protected: +/// The Analysis action that we are reusing data from + AnalysisBase* mydata; +public: + static void registerKeywords( Keywords& keys ); + AnalysisBase(const ActionOptions&); +/// These are required because we inherit from both ActionAtomistic and ActionWithArguments + void lockRequests(); + void unlockRequests(); +/// Return the number of data points + virtual unsigned getNumberOfDataPoints() const ; +/// Return the index of the data point in the base class + virtual unsigned getDataPointIndexInBase( const unsigned& idata ) const ; +/// Return the weight of the ith point + virtual double getWeight( const unsigned& idata ) const ; +/// Are we using memory in this calculation this affects the weights of points + virtual bool usingMemory() const ; +/// Return the normalisation constant for the calculation + virtual double getNormalization() const ; +/// Ensures that dissimilarities were set somewhere + virtual bool dissimilaritiesWereSet() const ; +/// Get the squared dissimilarity between two reference configurations + virtual double getDissimilarity( const unsigned& i, const unsigned& j ); +/// This returns the label of the object that contains the base data + virtual std::string getBaseDataLabel() const ; +/// Get the ith data point + virtual void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ; +/// Get a reference configuration + virtual ReferenceConfiguration* getReferenceConfiguration( const unsigned& idata ); +/// This actually performs the analysis + virtual void performAnalysis()=0; +/// These overwrite things from inherited classes (this is a bit of a fudge) + bool isPeriodic(){ plumed_error(); return false; } + unsigned getNumberOfDerivatives(){ plumed_error(); return 0; } + void calculateNumericalDerivatives( ActionWithValue* a=NULL ){ plumed_error(); } +/// Calculate and apply do nothing all analysis is done during update step + void calculate(){} + void apply(){} +/// This will call the analysis to be performed + virtual void update(); +/// This calls the analysis to be performed in the final step of the calculation +/// i.e. when use_all_data is true + virtual void runFinalJobs(); +}; + +inline +void AnalysisBase::lockRequests(){ + ActionAtomistic::lockRequests(); + ActionWithArguments::lockRequests(); +} + +inline +void AnalysisBase::unlockRequests(){ + ActionAtomistic::unlockRequests(); + ActionWithArguments::unlockRequests(); +} + +inline +unsigned AnalysisBase::getNumberOfDataPoints() const { + return mydata->getNumberOfDataPoints(); +} + +inline +unsigned AnalysisBase::getDataPointIndexInBase( const unsigned& idata ) const { + return mydata->getDataPointIndexInBase( idata ); +} + +inline +double AnalysisBase::getWeight( const unsigned& idata ) const { + return mydata->getWeight( idata ); +} + +inline +bool AnalysisBase::usingMemory() const { + return mydata->usingMemory(); +} + +inline +double AnalysisBase::getNormalization() const { + return mydata->getNormalization(); +} + +inline +bool AnalysisBase::dissimilaritiesWereSet() const { + return mydata->dissimilaritiesWereSet(); +} + +inline +double AnalysisBase::getDissimilarity( const unsigned& i, const unsigned& j ){ + return mydata->getDissimilarity( i, j ); +} + +inline +std::string AnalysisBase::getBaseDataLabel() const { + return mydata->getBaseDataLabel(); +} + +inline +void AnalysisBase::getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const { + mydata->getDataPoint( idata, point, weight ); +} + +inline +ReferenceConfiguration* AnalysisBase::getReferenceConfiguration( const unsigned& idata ){ + return mydata->getReferenceConfiguration( idata ); +} + +} +} + +#endif diff --git a/src/analysis/AnalysisWithDataCollection.cpp b/src/analysis/AnalysisWithDataCollection.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e19409dbdcc4219b0fcf8b8177bf4e001fb4ea83 --- /dev/null +++ b/src/analysis/AnalysisWithDataCollection.cpp @@ -0,0 +1,343 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2012-2015 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed-code.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see <http://www.gnu.org/licenses/>. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "AnalysisWithDataCollection.h" +#include "reference/ReferenceConfiguration.h" +#include "reference/ReferenceArguments.h" +#include "reference/ReferenceAtoms.h" +#include "reference/MetricRegister.h" +#include "core/PlumedMain.h" +#include "core/ActionSet.h" +#include "core/Atoms.h" +#include "tools/IFile.h" + +namespace PLMD { +namespace analysis { + +void AnalysisWithDataCollection::registerKeywords( Keywords& keys ){ + AnalysisBase::registerKeywords( keys ); + keys.use("ARG"); keys.reset_style("ARG","atoms-1"); + keys.add("atoms-1","ATOMS","the atoms whose positions we are tracking for the purpose of analysing the data"); + keys.add("hidden","METRIC","how are we measuring the distances between configurations. If you have only arguments this will by default be the euclidean metric. You must specify a metric if you are analysing atoms"); + keys.add("atoms-1","STRIDE","the frequency with which data should be stored for analysis. By default data is collected on every step"); + keys.add("atoms-1","RUN","the frequency with which to run the analysis algorithms."); + keys.addFlag("USE_ALL_DATA",false,"just analyse all the data in the trajectory. This option should be used in tandem with ATOMS/ARG + STRIDE"); + keys.addFlag("REWEIGHT_BIAS",false,"reweight the data using all the biases acting on the dynamics. This option must be used in tandem with ATOMS/ARG + STRIDE. For more information see \\ref reweighting"); + keys.add("optional","REWEIGHT_TEMP","reweight data from a trajectory at one temperature and output the probability distribution at a second temperature. This option must be used in tandem with ATOMS/ARG + STRIDE. For more information see \\ref reweighting"); + keys.add("optional","TEMP","the system temperature. This is required if you are reweighting (REWEIGHT_BIAS/REWEIGHT_TEMP) or if you are calculating free energies. You are not required to specify the temperature if this is passed by the underlying MD code."); + keys.add("atoms-3","REUSE_INPUT_DATA_FROM","do a second form of analysis on the data stored by a previous analysis object"); + keys.addFlag("WRITE_CHECKPOINT",false,"write out a checkpoint so that the analysis can be restarted in a later run. This option must be used in tandem with ATOMS/ARG + STRIDE."); + keys.addFlag("NOMEMORY",false,"do a block averaging i.e. analyse each block of data separately. This option must be used in tandem with ATOMS/ARG + STRIDE."); + keys.use("RESTART"); keys.use("UPDATE_FROM"); keys.use("UPDATE_UNTIL"); +} + +AnalysisWithDataCollection::AnalysisWithDataCollection( const ActionOptions& ao ): +Action(ao), +AnalysisBase(ao), +nomemory(false), +write_chq(false), +rtemp(0), +idata(0), +firstAnalysisDone(false), +old_norm(0.0) +{ + if( !mydata ){ + // Check if we are using the input data from another action + std::string datastr; parse("REUSE_INPUT_DATA_FROM",datastr); + if( datastr.length()>0 ) { + AnalysisWithDataCollection* checkd = plumed.getActionSet().selectWithLabel<AnalysisWithDataCollection*>( datastr ); + if( !checkd) error("cannot reuse input data from action with label " + datastr + " as this does not store data"); + mydata=dynamic_cast<AnalysisBase*>( checkd ); + log.printf(" performing analysis on input data stored in from %s \n",datastr.c_str() ); + freq=mydata->freq; use_all_data=mydata->use_all_data; + if( !use_all_data ) setStride( freq ); + // If we are not using input data from elsewhere or output data from another object then + // we must collect data from the trajectory + } else { + // Get information on numbers of atoms and argument names + std::vector<AtomNumber> atom_numbers; std::vector<std::string> argument_names( getNumberOfArguments() ); + for(unsigned i=0;i<getNumberOfArguments();++i) argument_names[i]=getPntrToArgument(i)->getName(); + + // Read in information on the metric that is being used in this analysis object + parse("METRIC",metricname); + if( metricname.length()==0 ) metricname="EUCLIDEAN"; + ReferenceConfiguration* checkref=metricRegister().create<ReferenceConfiguration>( metricname ); + // Check if we should read atoms + ReferenceAtoms* hasatoms=dynamic_cast<ReferenceAtoms*>( checkref ); + if( hasatoms ){ + parseAtomList("ATOMS",atom_numbers); requestAtoms(atom_numbers); + if( atom_numbers.size()==0 ) error("no atom positions have been specified in input"); + log.printf(" monitoring positions of atoms "); + for(unsigned i=0;i<atom_numbers.size();++i) log.printf("%d ",atom_numbers[i].serial() ); + log.printf("\n"); + } + // Check if we should read arguments + ReferenceArguments* hasargs=dynamic_cast<ReferenceArguments*>( checkref ); + if( !hasargs && getNumberOfArguments()!=0 ) error("use of arguments with metric type " + metricname + " is invalid"); + if( hasargs && getNumberOfArguments()==0 ) error("no arguments have been specified in input"); + if( hasatoms && hasargs ) error("currently dependencies break if you have both arguments and atoms"); // Not sure if this is really a problem anymore + // And delte the fake reference we created + delete checkref; + log.printf(" storing data as %s type reference objects \n",metricname.c_str() ); + + // Read in the information about how often to run the analysis (storage is read in in ActionPilot.cpp) + parseFlag("USE_ALL_DATA",use_all_data); + if(!use_all_data){ + parse("RUN",freq); + // Setup everything given the ammount of data that we will have in each analysis + if( freq%getStride()!= 0 ) error("Frequncy of running is not a multiple of the stride"); + unsigned ndata=freq/getStride(); data.resize(ndata); logweights.resize( ndata ); + for(unsigned i=0;i<ndata;++i){ + data[i]=metricRegister().create<ReferenceConfiguration>( metricname ); + data[i]->setNamesAndAtomNumbers( atom_numbers, argument_names ); + } + log.printf(" running analysis every %u steps\n",freq); + // Check if we are doing block averaging + parseFlag("NOMEMORY",nomemory); + if(nomemory) log.printf(" doing block averaging and analysing each portion of trajectory separately\n"); + } else { + log.printf(" analysing all data in trajectory\n"); + } + + // Read in stuff for reweighting of trajectories + + // Reweighting for biases + bool dobias; parseFlag("REWEIGHT_BIAS",dobias); + if( dobias ){ + std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); + if( all.empty() ) error("your input file is not telling plumed to calculate anything"); + std::vector<Value*> arg( getArguments() ); + log.printf(" reweigting using the following biases "); + for(unsigned j=0;j<all.size();j++){ + std::string flab; flab=all[j]->getLabel() + ".bias"; + if( all[j]->exists(flab) ){ + biases.push_back( all[j]->copyOutput(flab) ); + arg.push_back( all[j]->copyOutput(flab) ); + log.printf(" %s",flab.c_str()); + } + } + log.printf("\n"); + if( biases.empty() ) error("you are asking to reweight bias but there does not appear to be a bias acting on your system"); + requestArguments( arg ); + } + + // Reweighting for temperatures + rtemp=0; parse("REWEIGHT_TEMP",rtemp); + if( rtemp!=0 ){ + rtemp*=plumed.getAtoms().getKBoltzmann(); + log.printf(" reweighting simulation to probabilities at temperature %f\n",rtemp); + } + // Now retrieve the temperature in the simulation + simtemp=0; parse("TEMP",simtemp); + if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann(); + else simtemp=plumed.getAtoms().getKbT(); + if(simtemp==0 && (rtemp!=0 || !biases.empty()) ) error("The MD engine does not pass the temperature to plumed so you have to specify it using TEMP"); + + // Check if a check point is required (this should be got rid of at some point when we have proper checkpointing) GAT + parseFlag("WRITE_CHECKPOINT",write_chq); + std::string filename = getName() + "_" + getLabel() + ".chkpnt"; + if( write_chq ) rfile.link(*this); + if( getRestart() ){ + if( !write_chq ) warning("restarting without writing a checkpoint file is somewhat strange"); + // Read in data from input file + readCheckPointFile( filename ); + // Setup the restart file (append mode) + if( write_chq ) rfile.open( filename.c_str() ); // In append mode automatically because of restart + // Run the analysis if we stoped in the middle of it last time + log.printf(" restarting analysis with %u points read from restart file\n",idata); + } else if( write_chq ){ + // Setup the restart file (delete any old one) + rfile.open( filename.c_str() ); // In overwrite mode automatically because there is no restart + } + if( write_chq ){ + rfile.addConstantField("old_normalization"); + for(unsigned i=0;i<getNumberOfArguments();++i) rfile.setupPrintValue( getPntrToArgument(i) ); + } + // This is the end of the stuff for the checkpoint file - hopefully get rid of all this in not too distant future GAT + } + } +} + +AnalysisWithDataCollection::~AnalysisWithDataCollection(){ + for(unsigned i=0;i<data.size();++i ) delete data[i]; + if( write_chq ) rfile.close(); +} + + +void AnalysisWithDataCollection::readCheckPointFile( const std::string& filename ){ + FILE* fp=fopen(filename.c_str(),"r"); double tstep, oldtstep; + if(fp!=NULL){ + bool do_read=true, first=true; + while (do_read) { + PDB mypdb; + do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()); + if(do_read){ + data[idata]->set( mypdb ); + data[idata]->parse("TIME",tstep); + if( !first && ((tstep-oldtstep) - getStride()*plumed.getAtoms().getTimeStep())>plumed.getAtoms().getTimeStep() ){ + error("frequency of data storage in " + filename + " is not equal to frequency of data storage plumed.dat file"); + } + data[idata]->parse("LOG_WEIGHT",logweights[idata]); + data[idata]->parse("OLD_NORM",old_norm); + data[idata]->checkRead(); + idata++; first=false; oldtstep=tstep; + } else{ + break; + } + } + fclose(fp); + } + if(old_norm>0) firstAnalysisDone=true; +} + +void AnalysisWithDataCollection::prepare(){ + if(rtemp!=0) plumed.getAtoms().setCollectEnergy(true); +} + +double AnalysisWithDataCollection::getWeight( const unsigned& idat ) const { + if( !mydata ){ plumed_dbg_assert( idat<data.size() ); return data[idat]->getWeight(); } + return AnalysisBase::getWeight( idat ); +} + +void AnalysisWithDataCollection::getDataPoint( const unsigned& idat, std::vector<double>& point, double& weight ) const { + if( !mydata ){ + plumed_dbg_assert( idat<data.size() && point.size()==getNumberOfArguments() ); + weight = data[idat]->getWeight(); + for(unsigned i=0;i<point.size();++i) point[i]=data[idat]->getReferenceArgument(i); + } else { + AnalysisBase::getDataPoint( idat, point, weight ); + } +} + +ReferenceConfiguration* AnalysisWithDataCollection::getReferenceConfiguration( const unsigned& idat ){ + if( !mydata ){ plumed_dbg_assert( idat<data.size() ); return data[idat]; } + return AnalysisBase::getReferenceConfiguration( idat ); +} + +void AnalysisWithDataCollection::update(){ + if( mydata ){ AnalysisBase::update(); return; } + + // Ignore first bit of data if we are not using all data - this is a weird choice - I am not sure I understand GAT (perhaps this should be changed)? + if( !use_all_data && getStep()==0 ) return ; + + // This bit does the collection of data + if( idata<logweights.size() || use_all_data ){ + // Retrieve the bias + double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get(); + + double ww=0; + if(rtemp!=0){ + double energy=plumed.getAtoms().getEnergy()+bias; + // Reweighting because of temperature difference + ww=-( (1.0/rtemp) - (1.0/simtemp) )*(energy+bias); + } + // Reweighting because of biases + if( !biases.empty() ) ww += bias/simtemp; + + // Get the arguments ready to transfer to reference configuration + std::vector<double> current_args( getNumberOfArguments() ); + for(unsigned i=0;i<getNumberOfArguments();++i) current_args[i]=getArgument(i); + // Could add stuff for fancy metrics here eventually but for now unecessary + std::vector<double> mymetric( getNumberOfArguments(), 1.0 ); + + if(use_all_data){ + data.push_back( metricRegister().create<ReferenceConfiguration>( metricname ) ); + plumed_dbg_assert( data.size()==idata+1 ); + std::vector<std::string> argument_names( getNumberOfArguments() ); + for(unsigned i=0;i<getNumberOfArguments();++i) argument_names[i] = getPntrToArgument(i)->getName(); + data[idata]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names ); + data[idata]->setReferenceConfig( getPositions(), current_args, mymetric ); + logweights.push_back(ww); + } else { + // Get the arguments and store them in a vector of vectors + data[idata]->setReferenceConfig( getPositions(), current_args, mymetric ); + logweights[idata] = ww; + } + + // Write data to checkpoint file + if( write_chq ){ + rfile.rewind(); + if( plumed.getAtoms().usingNaturalUnits() ) data[idata]->print( 1.0, rfile, getTime(), logweights[idata], old_norm ); + else data[idata]->print( atoms.getUnits().getLength()/0.1, rfile, getTime(), logweights[idata], old_norm ); + rfile.flush(); + } + // Increment data counter + idata++; + } + + // This makes sure analysis is done when it should be done + if( !use_all_data ){ + if( getStep()>0 && getStep()%freq==0 ) runAnalysis(); + else if( idata==logweights.size() ) error("something has gone wrong. Probably a wrong initial time on restart"); + } +} + +void AnalysisWithDataCollection::runFinalJobs(){ + if( mydata ){ AnalysisBase::runFinalJobs(); return; } + if( use_all_data ) runAnalysis(); +} + +void AnalysisWithDataCollection::runAnalysis(){ + // Note : could add multiple walkers here - simply read in the data from all + // other walkers here if we are writing the check points. + + // Ensure weights are not finalized if we are using readDissimilarityMatrix + if( logweights.size()>0 ){ + // Check for errors + if( getNumberOfDataPoints()==0 ) error("no data is available for analysis"); + if( idata!=logweights.size() ) error("something has gone wrong. Am trying to run analysis but I don't have sufficient data"); + norm=0; // Reset normalization constant + if( biases.empty() && rtemp==0 ){ + for(unsigned i=0;i<logweights.size();++i){ + data[i]->setWeight(1.0); norm+=1.0; + } + } else if( nomemory ){ + // Find the maximum weight + double maxweight=logweights[0]; + for(unsigned i=1;i<getNumberOfDataPoints();++i){ + if(logweights[i]>maxweight) maxweight=logweights[i]; + } + // Calculate normalization constant + for(unsigned i=0;i<logweights.size();++i){ + norm+=exp( logweights[i]-maxweight ); + } + // Calculate weights (no memory) + for(unsigned i=0;i<logweights.size();++i) data[i]->setWeight( exp( logweights[i]-maxweight ) ); + // Calculate normalized weights (with memory) + } else { + // Calculate normalization constant + for(unsigned i=0;i<logweights.size();++i){ + norm+=exp( logweights[i] ); + } + if( !firstAnalysisDone ) old_norm=1.0; + // Calculate weights (with memory) + for(unsigned i=0;i<logweights.size();++i) data[i]->setWeight( exp( logweights[i] ) / old_norm ); + if( !firstAnalysisDone ) old_norm=0.0; + } + } + // And run the analysis + performAnalysis(); idata=0; + // Update total normalization constant + old_norm+=norm; firstAnalysisDone=true; +} + +} +} diff --git a/src/analysis/AnalysisWithDataCollection.h b/src/analysis/AnalysisWithDataCollection.h new file mode 100644 index 0000000000000000000000000000000000000000..2d4dd9bc6b56ba3ea0339a3dea43c5e57f2dbe66 --- /dev/null +++ b/src/analysis/AnalysisWithDataCollection.h @@ -0,0 +1,148 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2012-2014 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed-code.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see <http://www.gnu.org/licenses/>. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_analysis_AnalysisWithDataCollection_h +#define __PLUMED_analysis_AnalysisWithDataCollection_h + +#include "AnalysisBase.h" + +namespace PLMD { + +class ReferenceConfiguration; + +namespace analysis { + +class AnalysisWithDataCollection : public AnalysisBase { +private: +/// Are we treating each block of data separately + bool nomemory; +/// Are we writing a checkpoint file + bool write_chq; +/// The temperature at which we are running the calculation + double simtemp; +/// The temperature at which we want the histogram + double rtemp; +/// The biases we are using in reweighting and the args we store them separately + std::vector<Value*> biases; +/// The piece of data we are inserting + unsigned idata; +/// The data we are going to analyze + std::vector<ReferenceConfiguration*> data; +/// The weights of all the data points + std::vector<double> logweights; +/// Have we analyzed the data for the first time + bool firstAnalysisDone; +/// The value of the old normalization constant + double norm, old_norm; +/// List of argument names + std::vector<std::string> argument_names; +/// The type of metric we are using to measure distances + std::string metricname; +/// The checkpoint file --- really I would like to get rid of this and have some universal mechanism and a single file GT + OFile rfile; +/// Read the checkpoint file + void readCheckPointFile( const std::string& filename ); +/// Perform the analysis -- we have a funciton as it is called from both runFinalJobs() and upate() + void runAnalysis(); +protected: +/// Return the temperature (used by Histogram) + double getTemp() const { return simtemp; } +public: + static void registerKeywords( Keywords& keys ); + AnalysisWithDataCollection(const ActionOptions&); + ~AnalysisWithDataCollection(); +/// Return the number of arguments (this overwrites the one in ActionWithArguments) + unsigned getNumberOfArguments() const; +/// Return the number of data points + virtual unsigned getNumberOfDataPoints() const ; +/// Return the index of the data point in the base class + virtual unsigned getDataPointIndexInBase( const unsigned& idata ) const ; +/// Return the weight of the ith point + virtual double getWeight( const unsigned& idata ) const ; +/// Are we using memory in this calculation this affects the weights + virtual bool usingMemory() const ; +/// Return the normalisation constant for the calculation + virtual double getNormalization() const ; +/// By default dissimilarities are not set - they are only set in dissimilarity objects + virtual bool dissimilaritiesWereSet() const ; +/// This returns the label of the object that contains the base data + std::string getBaseDataLabel() const ; +/// Get the ith data point + virtual void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ; +/// Get a reference configuration + virtual ReferenceConfiguration* getReferenceConfiguration( const unsigned& idat ); +/// This ensures that the energy is stored if we are reweighting + void prepare(); +/// This stores the data and calls the analysis to be performed + void update(); +/// This does the analysis if it is to be done in the final step of the calculation + void runFinalJobs(); +}; + +inline +unsigned AnalysisWithDataCollection::getNumberOfDataPoints() const { + if( !mydata ) return data.size(); + return AnalysisBase::getNumberOfDataPoints(); +} + +inline +unsigned AnalysisWithDataCollection::getDataPointIndexInBase( const unsigned& idata ) const { + if( !mydata ) return idata; + return AnalysisBase::getDataPointIndexInBase( idata ); +} + +inline +bool AnalysisWithDataCollection::usingMemory() const { + if( !mydata ) return !nomemory; + return AnalysisBase::usingMemory(); +} + +inline +bool AnalysisWithDataCollection::dissimilaritiesWereSet() const { + if( !mydata ) return false; + return AnalysisBase::dissimilaritiesWereSet(); +} + +inline +double AnalysisWithDataCollection::getNormalization() const { + if( !mydata ){ + if( nomemory || !firstAnalysisDone ) return norm; + return ( 1. + norm/old_norm ); + } + return AnalysisBase::getNormalization(); +} + +inline +std::string AnalysisWithDataCollection::getBaseDataLabel() const { + if( !mydata ) return getLabel(); + return AnalysisBase::getBaseDataLabel(); +} + +inline +unsigned AnalysisWithDataCollection::getNumberOfArguments() const { + unsigned nargs=ActionWithArguments::getNumberOfArguments(); + return nargs - biases.size(); +} + +} +} + +#endif diff --git a/src/analysis/AnalysisWithLandmarks.cpp b/src/analysis/AnalysisWithLandmarks.cpp deleted file mode 100644 index 4a4548b884a11ec9317ba71c904d336be0426fbc..0000000000000000000000000000000000000000 --- a/src/analysis/AnalysisWithLandmarks.cpp +++ /dev/null @@ -1,79 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2013,2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "AnalysisWithLandmarks.h" -#include "LandmarkRegister.h" -#include "LandmarkSelectionBase.h" - -//+PLUMEDOC INTERNAL landmarkselection -/* -This is currently a filler page. - -Just use LANDMARKS=ALL. More complex versions will appear in later versions. - -*/ -//+ENDPLUMEDOC - -namespace PLMD { -namespace analysis { - -void AnalysisWithLandmarks::registerKeywords( Keywords& keys ){ - Analysis::registerKeywords( keys ); - keys.add("compulsory","LANDMARKS","ALL","only use a subset of the data that was collected. " - "For more information on the landmark selection algorithms that are available in " - "plumed see \\ref landmarkselection."); -} - -AnalysisWithLandmarks::AnalysisWithLandmarks( const ActionOptions& ao): -Action(ao), -Analysis(ao), -data_to_analyze(NULL) -{ - std::string linput; parse("LANDMARKS",linput); - std::vector<std::string> words=Tools::getWords(linput); - landmarkSelector=landmarkRegister().create( LandmarkSelectionOptions(words,this) ); - log.printf(" %s\n", landmarkSelector->description().c_str() ); -} - -AnalysisWithLandmarks::~AnalysisWithLandmarks(){ - delete landmarkSelector; -} - -void AnalysisWithLandmarks::setDataToAnalyze( MultiReferenceBase* mydata ){ - data_to_analyze=mydata; -} - -unsigned AnalysisWithLandmarks::getNumberOfLandmarks() const { - return landmarkSelector->getNumberOfLandmarks(); -} - -void AnalysisWithLandmarks::performAnalysis(){ - plumed_assert( data_to_analyze ); - landmarkSelector->selectLandmarks( data_to_analyze ); - analyzeLandmarks(); -} - -void AnalysisWithLandmarks::performTask( const unsigned& taskIndex, const unsigned& current, MultiValue& myvals ) const { - plumed_merror("Should not be here"); -} - -} -} diff --git a/src/analysis/AnalysisWithLandmarks.h b/src/analysis/AnalysisWithLandmarks.h deleted file mode 100644 index 86426fc7961c612b86af1e9c0bf6ad25cd96991b..0000000000000000000000000000000000000000 --- a/src/analysis/AnalysisWithLandmarks.h +++ /dev/null @@ -1,61 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2012-2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifndef __PLUMED_analysis_AnalysisWithLandmarks_h -#define __PLUMED_analysis_AnalysisWithLandmarks_h - -#include "Analysis.h" - -namespace PLMD { - -class MultiReferenceBase; - -namespace analysis { - -class LandmarkSelectionBase; - -class AnalysisWithLandmarks : public Analysis { -friend class LandmarkSelectionBase; -friend class CopyAllFrames; -private: -/// This object selects landmarks from the data - LandmarkSelectionBase* landmarkSelector; -/// A pointer to the data we are analyzing - MultiReferenceBase* data_to_analyze; -protected: -/// Set the data that needs to be analyzed - void setDataToAnalyze( MultiReferenceBase* mydata ); -/// Return the number of landmarks we are selecting - unsigned getNumberOfLandmarks() const ; -public: - static void registerKeywords( Keywords& keys ); - AnalysisWithLandmarks( const ActionOptions& ); - ~AnalysisWithLandmarks(); -/// Do the analysis - void performAnalysis(); - virtual void analyzeLandmarks()=0; -/// This does nothing - void performTask( const unsigned& , const unsigned& , MultiValue& ) const ; -}; - -} -} -#endif diff --git a/src/analysis/ClassicalMultiDimensionalScaling.cpp b/src/analysis/ClassicalMultiDimensionalScaling.cpp deleted file mode 100644 index f64e5211d5a1f1a4c2be805f11ac494c526516bf..0000000000000000000000000000000000000000 --- a/src/analysis/ClassicalMultiDimensionalScaling.cpp +++ /dev/null @@ -1,248 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2013,2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "AnalysisWithLandmarks.h" -#include "ClassicalScaling.h" -#include "reference/PointWiseMapping.h" -#include "core/ActionRegister.h" - -namespace PLMD { -namespace analysis { - -//+PLUMEDOC ANALYSIS CLASSICAL_MDS -/* -Create a low-dimensional projection of a trajectory using the classical multidimensional -scaling algorithm. - -Multidimensional scaling (MDS) is similar to what is done when you make a map. You start with distances -between London, Belfast, Paris and Dublin and then you try to arrange points on a piece of paper so that the (suitably scaled) -distances between the points in your map representing each of those cities are related to the true distances between the cities. -Stating this more mathematically MDS endeavors to find an <a href="http://en.wikipedia.org/wiki/Isometry">isometry</a> -between points distributed in a high-dimensional space and a set of points distributed in a low-dimensional plane. -In other words, if we have \f$M\f$ \f$D\f$-dimensional points, \f$\mathbf{X}\f$, -and we can calculate dissimilarities between pairs them, \f$D_{ij}\f$, we can, with an MDS calculation, try to create \f$M\f$ projections, -\f$\mathbf{x}\f$, of the high dimensionality points in a \f$d\f$-dimensional linear space by trying to arrange the projections so that the -Euclidean distances between pairs of them, \f$d_{ij}\f$, resemble the dissimilarities between the high dimensional points. In short we minimize: - -\f[ -\chi^2 = \sum_{i \ne j} \left( D_{ij} - d_{ij} \right)^2 -\f] - -where \f$D_{ij}\f$ is the distance between point \f$X^{i}\f$ and point \f$X^{j}\f$ and \f$d_{ij}\f$ is the distance between the projection -of \f$X^{i}\f$, \f$x^i\f$, and the projection of \f$X^{j}\f$, \f$x^j\f$. A tutorial on this approach can be used to analyse simulations -can be found in the tutorial \ref belfast-3 and in the following <a href="https://www.youtube.com/watch?v=ofC2qz0_9_A&feature=youtu.be" > short video.</a> - -\par Examples - -The following command instructs plumed to construct a classical multidimensional scaling projection of a trajectory. -The RMSD distance between atoms 1-256 have moved is used to measure the distances in the high-dimensional space. - -\verbatim -CLASSICAL_MDS ... - ATOMS=1-256 - METRIC=OPTIMAL-FAST - USE_ALL_DATA - NLOW_DIM=2 - OUTPUT_FILE=rmsd-embed -... CLASSICAL_MDS -\endverbatim - -The following section is for people who are interested in how this method works in detail. A solid understanding of this material is -not necessary to use MDS. - -\section dim-sec Method of optimisation - -The stress function can be minimized using a standard optimization algorithm such as conjugate gradients or steepest descent. -However, it is more common to do this minimization using a technique known as classical scaling. Classical scaling works by -recognizing that each of the distances $D_{ij}$ in the above sum can be written as: - -\f[ -D_{ij}^2 = \sum_{\alpha} (X^i_\alpha - X^j_\alpha)^2 = \sum_\alpha (X^i_\alpha)^2 + (X^j_\alpha)^2 - 2X^i_\alpha X^j_\alpha -\f] - -We can use this expression and matrix algebra to calculate multiple distances at once. For instance if we have three points, -\f$\mathbf{X}\f$, we can write distances between them as: - -\f{eqnarray*}{ -D^2(\mathbf{X}) &=& \left[ \begin{array}{ccc} -0 & d_{12}^2 & d_{13}^2 \\ -d_{12}^2 & 0 & d_{23}^2 \\ -d_{13}^2 & d_{23}^2 & 0 -\end{array}\right] \\ -&=& -\sum_\alpha \left[ \begin{array}{ccc} -(X^1_\alpha)^2 & (X^1_\alpha)^2 & (X^1_\alpha)^2 \\ -(X^2_\alpha)^2 & (X^2_\alpha)^2 & (X^2_\alpha)^2 \\ -(X^3_\alpha)^2 & (X^3_\alpha)^2 & (X^3_\alpha)^2 \\ -\end{array}\right] - + \sum_\alpha \left[ \begin{array}{ccc} -(X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\ -(X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\ -(X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\ -\end{array}\right] -- 2 \sum_\alpha \left[ \begin{array}{ccc} -X^1_\alpha X^1_\alpha & X^1_\alpha X^2_\alpha & X^1_\alpha X^3_\alpha \\ -X^2_\alpha X^1_\alpha & X^2_\alpha X^2_\alpha & X^2_\alpha X^3_\alpha \\ -X^1_\alpha X^3_\alpha & X^3_\alpha X^2_\alpha & X^3_\alpha X^3_\alpha -\end{array}\right] \nonumber \\ -&=& \mathbf{c 1^T} + \mathbf{1 c^T} - 2 \sum_\alpha \mathbf{x}_a \mathbf{x}^T_a = \mathbf{c 1^T} + \mathbf{1 c^T} - 2\mathbf{X X^T} -\f} - -This last equation can be extended to situations when we have more than three points. In it \f$\mathbf{X}\f$ is a matrix that has -one high-dimensional point on each of its rows and \f$\mathbf{X^T}\f$ is its transpose. \f$\mathbf{1}\f$ is an \f$M \times 1\f$ vector -of ones and \f$\mathbf{c}\f$ is a vector with components given by: - -\f[ -c_i = \sum_\alpha (x_\alpha^i)^2 -\f] - -These quantities are the diagonal elements of \f$\mathbf{X X^T}\f$, which is a dot product or Gram Matrix that contains the -dot product of the vector \f$X_i\f$ with the vector \f$X_j\f$ in element \f$i,j\f$. - -In classical scaling we introduce a centering matrix \f$\mathbf{J}\f$ that is given by: - -\f[ -\mathbf{J} = \mathbf{I} - \frac{1}{M} \mathbf{11^T} -\f] - -where \f$\mathbf{I}\f$ is the identity. Multiplying the equations above from the front and back by this matrix and a factor of a \f$-\frac{1}{2}\f$ gives: - -\f{eqnarray*}{ - -\frac{1}{2} \mathbf{J} \mathbf{D}^2(\mathbf{X}) \mathbf{J} &=& -\frac{1}{2}\mathbf{J}( \mathbf{c 1^T} + \mathbf{1 c^T} - 2\mathbf{X X^T})\mathbf{J} \\ - &=& -\frac{1}{2}\mathbf{J c 1^T J} - \frac{1}{2} \mathbf{J 1 c^T J} + \frac{1}{2} \mathbf{J}(2\mathbf{X X^T})\mathbf{J} \\ - &=& \mathbf{ J X X^T J } = \mathbf{X X^T } \label{eqn:scaling} -\f} - -The fist two terms in this expression disappear because \f$\mathbf{1^T J}=\mathbf{J 1} =\mathbf{0}\f$, where \f$\mathbf{0}\f$ -is a matrix containing all zeros. In the final step meanwhile we use the fact that the matrix of squared distances will not -change when we translate all the points. We can thus assume that the mean value, \f$\mu\f$, for each of the components, \f$\alpha\f$: -\f[ -\mu_\alpha = \frac{1}{M} \sum_{i=1}^N \mathbf{X}^i_\alpha -\f] -is equal to 0 so the columns of \f$\mathbf{X}\f$ add up to 0. This in turn means that each of the columns of -\f$\mathbf{X X^T}\f$ adds up to zero, which is what allows us to write \f$\mathbf{ J X X^T J } = \mathbf{X X^T }\f$. - -The matrix of squared distances is symmetric and positive-definite we can thus use the spectral decomposition to decompose it as: - -\f[ -\Phi= \mathbf{V} \Lambda \mathbf{V}^T -\f] - -Furthermore, because the matrix we are diagonalizing, \f$\mathbf{X X^T}\f$, is the product of a matrix and its transpose -we can use this decomposition to write: - -\f[ -\mathbf{X} =\mathbf{V} \Lambda^\frac{1}{2} -\f] - -Much as in PCA there are generally a small number of large eigenvalues in \f$\Lambda\f$ and many small eigenvalues. -We can safely use only the large eigenvalues and their corresponding eigenvectors to express the relationship between -the coordinates \f$\mathbf{X}\f$. This gives us our set of low-dimensional projections. - -This derivation makes a number of assumptions about the how the low dimensional points should best be arranged to minimise -the stress. If you use an interative optimization algorithm such as SMACOF you may thus be able to find a better -(lower-stress) projection of the points. For more details on the assumptions made -see <a href="http://quest4rigor.com/tag/multidimensional-scaling/"> this website.</a> -*/ -//+ENDPLUMEDOC - -class ClassicalMultiDimensionalScaling : public AnalysisWithLandmarks { -private: - unsigned nlow; - std::string ofilename; - std::string efilename; - PointWiseMapping* myembedding; -public: - static void registerKeywords( Keywords& keys ); - ClassicalMultiDimensionalScaling( const ActionOptions& ao ); - ~ClassicalMultiDimensionalScaling(); - void analyzeLandmarks(); -}; - -PLUMED_REGISTER_ACTION(ClassicalMultiDimensionalScaling,"CLASSICAL_MDS") - -void ClassicalMultiDimensionalScaling::registerKeywords( Keywords& keys ){ - AnalysisWithLandmarks::registerKeywords( keys ); - keys.add("compulsory","NLOW_DIM","number of low-dimensional coordinates required"); - keys.add("compulsory","OUTPUT_FILE","file on which to output the final embedding coordinates"); - keys.add("compulsory","EMBEDDING_OFILE","dont output","file on which to output the embedding in plumed input format"); -} - -ClassicalMultiDimensionalScaling::ClassicalMultiDimensionalScaling( const ActionOptions& ao ): -Action(ao), -AnalysisWithLandmarks(ao) -{ - myembedding = new PointWiseMapping( getMetricName(), false ); - setDataToAnalyze( dynamic_cast<MultiReferenceBase*>(myembedding) ); - - parse("NLOW_DIM",nlow); - if( nlow<1 ) error("dimensionality of low dimensional space must be at least one"); - std::vector<std::string> propnames( nlow ); std::string num; - for(unsigned i=0;i<propnames.size();++i){ - Tools::convert(i+1,num); std::string lab=getLabel(); - if(lab.find("@")!=std::string::npos) propnames[i]=getName() + "." + num; - else propnames[i]=getLabel() + "." + num; - } - myembedding->setPropertyNames( propnames, false ); - - parseOutputFile("EMBEDDING_OFILE",efilename); - parseOutputFile("OUTPUT_FILE",ofilename); -} - -ClassicalMultiDimensionalScaling::~ClassicalMultiDimensionalScaling(){ - delete myembedding; -} - -void ClassicalMultiDimensionalScaling::analyzeLandmarks(){ - // Calculate all pairwise diatances - myembedding->calculateAllDistances( getPbc(), getArguments(), comm, myembedding->modifyDmat(), true ); - - // Run multidimensional scaling - ClassicalScaling::run( myembedding ); - - // Output the embedding as long lists of data -// std::string gfname=saveResultsFromPreviousAnalyses( ofilename ); - OFile gfile; gfile.link(*this); - gfile.setBackupString("analysis"); - gfile.fmtField(getOutputFormat()+" "); - gfile.open( ofilename.c_str() ); - - // Print embedding coordinates - for(unsigned i=0;i<myembedding->getNumberOfReferenceFrames();++i){ - for(unsigned j=0;j<nlow;++j){ - std::string num; Tools::convert(j+1,num); - gfile.printField( getLabel() + "." + num , myembedding->getProjectionCoordinate(i,j) ); - } - gfile.printField(); - } - gfile.close(); - - // Output the embedding in plumed format - if( efilename!="dont output"){ - OFile afile; afile.link(*this); afile.setBackupString("analysis"); - afile.open( efilename.c_str() ); - myembedding->print( "classical mds", getTime(), afile, getOutputFormat() ); - afile.close(); - } -} - -} -} diff --git a/src/analysis/ClassicalScaling.cpp b/src/analysis/ClassicalScaling.cpp deleted file mode 100644 index 051ffbd92012435ed371a95297ac6a09ed8f77a6..0000000000000000000000000000000000000000 --- a/src/analysis/ClassicalScaling.cpp +++ /dev/null @@ -1,56 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2012-2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "ClassicalScaling.h" -#include "reference/PointWiseMapping.h" - -namespace PLMD { -namespace analysis { - -void ClassicalScaling::run( PointWiseMapping* mymap ){ - // Retrieve the distances from the dimensionality reduction object - double half=(-0.5); Matrix<double> distances( half*mymap->modifyDmat() ); - - // Apply centering transtion - unsigned n=distances.nrows(); double sum; - // First HM - for(unsigned i=0;i<n;++i){ - sum=0; for(unsigned j=0;j<n;++j) sum+=distances(i,j); - for(unsigned j=0;j<n;++j) distances(i,j) -= sum/n; - } - // Now (HM)H - for(unsigned i=0;i<n;++i){ - sum=0; for(unsigned j=0;j<n;++j) sum+=distances(j,i); - for(unsigned j=0;j<n;++j) distances(j,i) -= sum/n; - } - - // Diagonalize matrix - std::vector<double> eigval(n); Matrix<double> eigvec(n,n); - diagMat( distances, eigval, eigvec ); - - // Pass final projections to map object - for(unsigned i=0;i<n;++i){ - for(unsigned j=0;j<mymap->getNumberOfProperties();++j) mymap->setProjectionCoordinate( i, j, sqrt(eigval[n-1-j])*eigvec(n-1-j,i) ); - } -} - -} -} diff --git a/src/analysis/ClassicalScaling.h b/src/analysis/ClassicalScaling.h deleted file mode 100644 index 90cd03561fc4a867be22f21773d48cab3f1f7e62..0000000000000000000000000000000000000000 --- a/src/analysis/ClassicalScaling.h +++ /dev/null @@ -1,40 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2012-2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifndef __PLUMED_analysis_ClassicalScaling_h -#define __PLUMED_analysis_ClassicalScaling_h - -#include <vector> - -namespace PLMD { - -class PointWiseMapping; - -namespace analysis { - -class ClassicalScaling { -public: - static void run( PointWiseMapping* mymap ); -}; - -} -} -#endif diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp index 761d8a26a9dc12c2a22c5ce66cdc4ef6c3c49d0a..c0f8cc961994f0184429415c6ee69c97409242a0 100644 --- a/src/analysis/Histogram.cpp +++ b/src/analysis/Histogram.cpp @@ -19,7 +19,7 @@ You should have received a copy of the GNU Lesser General Public License along with plumed. If not, see <http://www.gnu.org/licenses/>. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "Analysis.h" +#include "AnalysisWithDataCollection.h" #include "core/PlumedMain.h" #include "core/ActionRegister.h" #include "tools/Grid.h" @@ -111,11 +111,12 @@ HISTOGRAM ... */ //+ENDPLUMEDOC -class Histogram : public Analysis { +class Histogram : public AnalysisWithDataCollection { private: std::vector<std::string> gmin, gmax; std::vector<double> point, bw; std::vector<unsigned> gbin; + std::string fmt; std::string gridfname; std::string kerneltype; bool fenergy; @@ -130,8 +131,8 @@ public: PLUMED_REGISTER_ACTION(Histogram,"HISTOGRAM") void Histogram::registerKeywords( Keywords& keys ){ - Analysis::registerKeywords( keys ); keys.reset_style("METRIC","hidden"); - keys.remove("ATOMS"); keys.reset_style("ARG","compulsory"); + AnalysisWithDataCollection::registerKeywords( keys ); + keys.remove("ATOMS"); keys.add("compulsory","GRID_MIN","the lower bounds for the grid"); keys.add("compulsory","GRID_MAX","the upper bounds for the grid"); keys.add("optional","GRID_BIN","the number of bins for the grid"); @@ -141,16 +142,20 @@ void Histogram::registerKeywords( Keywords& keys ){ keys.add("optional","BANDWIDTH","the bandwdith for kernel density estimation"); keys.addFlag("FREE-ENERGY",false,"Set to TRUE if you want a FREE ENERGY instead of a probabilty density (you need to set TEMP)."); keys.addFlag("UNNORMALIZED",false,"Set to TRUE if you don't want histogram to be normalized or free energy to be shifted."); + keys.add("optional","FMT","the format that should be used in the output file"); keys.add("compulsory","GRID_WFILE","histogram","the file on which to write the grid"); - keys.use("NOMEMORY"); } Histogram::Histogram(const ActionOptions&ao): -PLUMED_ANALYSIS_INIT(ao), +Action(ao), +AnalysisWithDataCollection(ao), point(getNumberOfArguments()), +fmt("%f"), fenergy(false), unnormalized(false) { + // Read output format + parse("FMT",fmt); // Read stuff for Grid parseVector("GRID_MIN",gmin); if(gmin.size()!=getNumberOfArguments()) error("Wrong number of values for GRID_MIN: they should be equal to the number of arguments"); @@ -177,7 +182,8 @@ unnormalized(false) unsigned n=((b-a)/gspacing[i])+1; if(gbin[i]<n) gbin[i]=n; } - parseOutputFile("GRID_WFILE",gridfname); + parse("GRID_WFILE",gridfname); + if( !getRestart() ){ OFile ofile; ofile.link(*this); ofile.setBackupString("analysis"); ofile.backupAllFiles(gridfname); } // Read the type of kernel we are using parse("KERNEL",kerneltype); @@ -223,8 +229,9 @@ void Histogram::performAnalysis(){ pmin.resize(getNumberOfArguments()); pmax.resize(getNumberOfArguments()); for(unsigned i=0;i<getNumberOfArguments();++i){ - pbc.push_back( getPeriodicityInformation(i,dmin,dmax) ); + pbc.push_back( getPntrToArgument(i)->isPeriodic() ); if(pbc[i]){ + getPntrToArgument(i)->getDomain(dmin,dmax); Tools::convert(dmin,gmin[i]); Tools::convert(dmax,gmax[i]); Tools::convert(dmin,pmin[i]); @@ -232,16 +239,22 @@ void Histogram::performAnalysis(){ } } - Grid* gg; IFile oldf; oldf.link(*this); + Grid* gg; IFile oldf; oldf.link(*this); + + // Retrieve the arguments remembering to delete all biases + std::vector<Value*> arg_vals( ActionWithArguments::getArguments() ); + unsigned nbias = arg_vals.size() - getNumberOfArguments(); + for(unsigned i=0;i<nbias;++i) arg_vals.erase(arg_vals.end()-1); + if( usingMemory() && oldf.FileExist(gridfname) ){ oldf.open(gridfname); - gg = Grid::create( "probs", getArguments(), oldf, gmin, gmax, gbin, false, false, false ); + gg = Grid::create( "probs", arg_vals, oldf, gmin, gmax, gbin, false, false, false ); oldf.close(); } else { - gg = new Grid( "probs", getArguments(), gmin, gmax, gbin,false,false); + gg = new Grid( "probs", arg_vals, gmin, gmax, gbin,false,false); } // Set output format for grid - gg->setOutputFmt( getOutputFormat() ); + gg->setOutputFmt( fmt ); // Now build the histogram double weight; std::vector<double> point( getNumberOfArguments() ); diff --git a/src/analysis/LandmarkRegister.cpp b/src/analysis/LandmarkRegister.cpp deleted file mode 100644 index 375315d7d038c7829df3741b98a07c6ab4ca6f04..0000000000000000000000000000000000000000 --- a/src/analysis/LandmarkRegister.cpp +++ /dev/null @@ -1,69 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2011-2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "LandmarkRegister.h" -#include <iostream> - -namespace PLMD{ -namespace analysis{ - -LandmarkRegister::~LandmarkRegister(){ - if(m.size()>0){ - std::string names=""; - for(std::map<std::string,creator_pointer>::iterator p=m.begin();p!=m.end();++p) names+=p->first+" "; - std::cerr<<"WARNING: ReferenceConfiguration "+ names +" has not been properly unregistered. This might lead to memory leak!!\n"; - } -} - -LandmarkRegister& landmarkRegister(){ - static LandmarkRegister ans; - return ans; -} - -void LandmarkRegister::remove(creator_pointer f){ - for(std::map<std::string,creator_pointer>::iterator p=m.begin();p!=m.end();++p){ - if((*p).second==f){ - m.erase(p); break; - } - } -} - -void LandmarkRegister::add( std::string type, creator_pointer f ){ - plumed_massert(m.count(type)==0,"type has already been registered"); - m.insert(std::pair<std::string,creator_pointer>(type,f)); -} - -bool LandmarkRegister::check(std::string type){ - if( m.count(type)>0 ) return true; - return false; -} - -LandmarkSelectionBase* LandmarkRegister::create( const LandmarkSelectionOptions& lo ){ - LandmarkSelectionBase* lselect; - if( check(lo.words[0]) ){ - lselect=m[lo.words[0]](lo); - lselect->checkRead(); - } else lselect=NULL; - return lselect; -} - -} -} diff --git a/src/analysis/LandmarkRegister.h b/src/analysis/LandmarkRegister.h deleted file mode 100644 index ef7b3cfb650350a3315279bfdcf41f398f54a1b3..0000000000000000000000000000000000000000 --- a/src/analysis/LandmarkRegister.h +++ /dev/null @@ -1,68 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2013,2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifndef __PLUMED_analysis_LandmarkRegister_h -#define __PLUMED_analysis_LandmarkRegister_h - -#include <string> -#include <cstring> -#include <vector> -#include <map> -#include "LandmarkSelectionBase.h" - -namespace PLMD{ - -class PDB; - -namespace analysis{ - -class LandmarkRegister{ -private: -/// Pointer to a function which, given the type for a ReferenceConfiguration, creates it - typedef LandmarkSelectionBase*(*creator_pointer)(const LandmarkSelectionOptions&); -/// The set of possible landmark selection algorithms we can work with - std::map<std::string,creator_pointer> m; -public: -/// The destructor - ~LandmarkRegister(); -/// Add a new landmark selection style to the register of landmark selectors - void add( std::string type, creator_pointer ); -/// Remove a landmark selection style from the register of metrics - void remove(creator_pointer f); -/// Verify if a landmark selection style is present in the register - bool check(std::string type); -/// Create a landmark selection object - LandmarkSelectionBase* create( const LandmarkSelectionOptions& lo ); -}; - -LandmarkRegister& landmarkRegister(); - -#define PLUMED_REGISTER_LANDMARKS(classname,type) \ - static class classname##RegisterMe{ \ - static LandmarkSelectionBase * create(const LandmarkSelectionOptions&lo){return new classname(lo);} \ - public: \ - classname##RegisterMe(){landmarkRegister().add(type,create);}; \ - ~classname##RegisterMe(){landmarkRegister().remove(create);}; \ - } classname##RegisterMeObject; - -} -} -#endif diff --git a/src/analysis/LandmarkSelectionBase.cpp b/src/analysis/LandmarkSelectionBase.cpp deleted file mode 100644 index 505d91dac04a163dfafe48f90b2307624b1f7a3c..0000000000000000000000000000000000000000 --- a/src/analysis/LandmarkSelectionBase.cpp +++ /dev/null @@ -1,116 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2013-2015 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "LandmarkSelectionBase.h" -#include "reference/MultiReferenceBase.h" - -namespace PLMD { -namespace analysis { - -LandmarkSelectionOptions::LandmarkSelectionOptions( const std::vector<std::string>& input, AnalysisWithLandmarks* myanalysis ): -words(input), -action(myanalysis) -{ -} - -LandmarkSelectionBase::LandmarkSelectionBase( const LandmarkSelectionOptions& lo ): -style(lo.words[0]), -input(lo.words), -action(lo.action) -{ - input.erase( input.begin() ); - if( style=="ALL" ){ - novoronoi=true; - } else { - parse("N",nlandmarks); - parseFlag("NOVORONOI",novoronoi); - } - parseFlag("IGNORE_WEIGHTS",noweights); -} - -LandmarkSelectionBase::~LandmarkSelectionBase(){ -} - -void LandmarkSelectionBase::parseFlag(const std::string& key, bool& t){ - Tools::parseFlag(input,key,t); -} - -void LandmarkSelectionBase::checkRead() const { - if(!input.empty()){ - std::string msg="cannot understand the following words from landmark selection input : "; - for(unsigned i=0;i<input.size();++i) msg = msg + input[i] + ", "; - plumed_merror(msg); - } -} - -std::string LandmarkSelectionBase::description(){ - std::ostringstream ostr; - if( style=="ALL"){ - ostr<<"using all data"; - } else { - ostr<<"selecting "<<nlandmarks<<" using "<<style<<" algorithm to analyze\n"; - ostr<<" "<<rest_of_description()<<"\n"; - if(noweights) ostr<<" ignoring all reweighting of data during landmark selection\n"; - if(novoronoi) ostr<<" voronoi weights will not be ascribed to points\n"; - } - return ostr.str(); -} - -double LandmarkSelectionBase::getWeightOfFrame( const unsigned& iframe ){ - if(noweights) return 1.0; - return action->getWeight(iframe); -} -double LandmarkSelectionBase::getDistanceBetweenFrames( const unsigned& iframe, const unsigned& jframe ){ - return distance( action->getPbc(), action->getArguments(), action->data[iframe], action->data[jframe], false ); -} - -void LandmarkSelectionBase::selectFrame( const unsigned& iframe, MultiReferenceBase* myframes){ - plumed_assert( myframes->getNumberOfReferenceFrames()<nlandmarks ); - myframes->copyFrame( action->data[iframe] ); -} - -void LandmarkSelectionBase::selectLandmarks( MultiReferenceBase* myframes ){ - // Select landmarks - myframes->clearFrames(); select( myframes ); - plumed_assert( myframes->getNumberOfReferenceFrames()==nlandmarks ); - - // Now calculate voronoi weights - if( !novoronoi ){ - unsigned rank=action->comm.Get_rank(); - unsigned size=action->comm.Get_size(); - std::vector<double> weights( nlandmarks, 0.0 ); - for(unsigned i=rank;i<action->data.size();i+=size){ - unsigned closest=0; - double mindist=distance( action->getPbc(), action->getArguments(), action->data[i], myframes->getFrame(0), false ); - for(unsigned j=1;j<nlandmarks;++j){ - double dist=distance( action->getPbc(), action->getArguments(), action->data[i], myframes->getFrame(j), false ); - if( dist<mindist ){ mindist=dist; closest=j; } - } - weights[closest] += getWeightOfFrame(i); - } - action->comm.Sum( &weights[0], weights.size() ); - myframes->setWeights( weights ); - } -} - -} -} - diff --git a/src/analysis/LandmarkSelectionBase.h b/src/analysis/LandmarkSelectionBase.h deleted file mode 100644 index 123dd14cfe8d72eb9424b4fed2b869277b7262f5..0000000000000000000000000000000000000000 --- a/src/analysis/LandmarkSelectionBase.h +++ /dev/null @@ -1,110 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2013-2015 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#ifndef __PLUMED_analysis_LandmarkSelectionBase_h -#define __PLUMED_analysis_LandmarkSelectionBase_h - -#include "reference/ReferenceConfiguration.h" -#include "AnalysisWithLandmarks.h" - -namespace PLMD { -namespace analysis { - -class LandmarkSelectionOptions{ -friend class LandmarkRegister; -friend class LandmarkSelectionBase; -private: - std::vector<std::string> words; - AnalysisWithLandmarks* action; -public: - LandmarkSelectionOptions( const std::vector<std::string>& input, AnalysisWithLandmarks* myanalysis ); -}; - -class LandmarkSelectionBase { -friend class AnalysisWithLandmarks; -friend class CopyAllFrames; -private: -/// Name of the method we are using for landmark selection - std::string style; -/// The number of landmarks we are selecting - unsigned nlandmarks; -/// The input to the landmark selection object - std::vector<std::string> input; -/// A pointer to the underlying action - AnalysisWithLandmarks* action; -/// How do we treat weights - bool novoronoi, noweights; -protected: -/// Return the numbe of landmarks - unsigned getNumberOfLandmarks() const ; -/// Return the communicator - Communicator& getCommunicator(); -/// Read a keywords from the input - template <class T> - void parse(const std::string& ,T& ); -/// Read a flag from the input - void parseFlag(const std::string& key, bool& t); -/// Get the number of frames in the underlying action - unsigned getNumberOfFrames() const; -/// Get the weight of the ith frame - double getWeightOfFrame( const unsigned& ); -/// Calculate the distance between the ith and jth frames - double getDistanceBetweenFrames( const unsigned& , const unsigned& ); -/// Transfer frame i in the underlying action to the object we are going to analyze - void selectFrame( const unsigned& , MultiReferenceBase* ); -public: - LandmarkSelectionBase( const LandmarkSelectionOptions& lo ); - virtual ~LandmarkSelectionBase(); -/// Check everything was read in - void checkRead() const ; -/// Return a description of the landmark selection protocol - std::string description(); -/// Overwrite this to have a more descriptive output - virtual std::string rest_of_description(){ return ""; }; -/// Actually do landmark selection - void selectLandmarks( MultiReferenceBase* ); - virtual void select( MultiReferenceBase* )=0; -}; - -inline -unsigned LandmarkSelectionBase::getNumberOfLandmarks() const { - return nlandmarks; -} - -inline -Communicator& LandmarkSelectionBase::getCommunicator(){ - return action->comm; -} - -inline -unsigned LandmarkSelectionBase::getNumberOfFrames() const { - return action->getNumberOfDataPoints(); -} - -template <class T> -void LandmarkSelectionBase::parse( const std::string& key, T& t ){ - bool found=Tools::parse(input,key,t); - if(!found) plumed_merror("landmark seleciton style " + style + " requires " + key + " keyword"); -} - -} -} -#endif diff --git a/src/analysis/SelectAllFrames.cpp b/src/analysis/SelectAllFrames.cpp deleted file mode 100644 index b0f7512643bcd3b4814b74b0912276340d9c49d7..0000000000000000000000000000000000000000 --- a/src/analysis/SelectAllFrames.cpp +++ /dev/null @@ -1,47 +0,0 @@ -/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - Copyright (c) 2012-2014 The plumed team - (see the PEOPLE file at the root of the distribution for a list of names) - - See http://www.plumed-code.org for more information. - - This file is part of plumed, version 2. - - plumed is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - plumed is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with plumed. If not, see <http://www.gnu.org/licenses/>. -+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ -#include "LandmarkSelectionBase.h" -#include "LandmarkRegister.h" - -namespace PLMD { -namespace analysis { - -class CopyAllFrames : public LandmarkSelectionBase { -public: - CopyAllFrames( const LandmarkSelectionOptions& lo ); - void select( MultiReferenceBase* ); -}; - -PLUMED_REGISTER_LANDMARKS(CopyAllFrames,"ALL") - -CopyAllFrames::CopyAllFrames( const LandmarkSelectionOptions& lo ): -LandmarkSelectionBase(lo) -{ -} - -void CopyAllFrames::select( MultiReferenceBase* myframes ){ - nlandmarks = action->getNumberOfDataPoints(); - for(unsigned i=0;i<getNumberOfFrames();++i) selectFrame( i, myframes ); -} - -} -} diff --git a/src/core/ActionPilot.cpp b/src/core/ActionPilot.cpp index 34b53adba45229e6a74c79ef628a9cbcd12a1089..ab3db6a3c292178d12ff13ca78a364c6a2a1df27 100644 --- a/src/core/ActionPilot.cpp +++ b/src/core/ActionPilot.cpp @@ -31,8 +31,14 @@ ActionPilot::ActionPilot(const ActionOptions&ao): Action(ao), stride(1) { - parse("STRIDE",stride); - log.printf(" with stride %d\n",stride); + if( keywords.exists("STRIDE") ){ + parse("STRIDE",stride); + log.printf(" with stride %d\n",stride); + } +} + +void ActionPilot::setStride( const unsigned& s ){ + stride=s; } bool ActionPilot::onStep()const{ diff --git a/src/core/ActionPilot.h b/src/core/ActionPilot.h index 551c73d0b901504bf5a2c5233096406a60996e31..bc44b5da17bf405271ce62ef6ebeb11a2ea625d9 100644 --- a/src/core/ActionPilot.h +++ b/src/core/ActionPilot.h @@ -42,6 +42,7 @@ class ActionPilot: int stride; // multiple time step protected: int getStride()const; + void setStride( const unsigned& s ); public: ActionPilot(const ActionOptions&); /// Create the keywords for actionPilot diff --git a/src/reference/PointWiseMapping.cpp b/src/reference/PointWiseMapping.cpp index 034f7fd772d5e26cd6d5118b2607293cb96e2bba..30b4ae1d3f4bede1a7448aacce33a80cbeb401c8 100644 --- a/src/reference/PointWiseMapping.cpp +++ b/src/reference/PointWiseMapping.cpp @@ -77,24 +77,24 @@ unsigned PointWiseMapping::getPropertyIndex( const std::string& name ) const { return 0; } -void PointWiseMapping::print( const std::string& method, const double & time, OFile& afile, const std::string& fmt ){ - std::string descr2, descr="DESCRIPTION: results from %s analysis performed at time " + fmt +"\n"; - afile.printf(descr.c_str(), method.c_str(), time ); - if(fmt.find("-")!=std::string::npos){ - descr="REMARK WEIGHT=" + fmt + " %s=" + fmt + " "; descr2="%s=" + fmt; - } else { - // This ensures numbers are left justified (i.e. next to the equals sign - std::size_t psign=fmt.find("%"); - plumed_assert( psign!=std::string::npos ); - descr="REMARK WEIGHT=%-" + fmt.substr(psign+1) + " %s=%-" + fmt.substr(psign+1) + " "; - descr2="%s=%-" + fmt.substr(psign+1); - } - for(unsigned i=0;i<frames.size();++i){ - afile.printf(descr.c_str(), frames[i]->getWeight(), property[0].c_str(), low_dim[i][0] ); - for(unsigned j=1;j<property.size();++j) afile.printf(descr2.c_str(), property[j].c_str(), low_dim[i][j]); - afile.printf("\n"); - frames[i]->print( afile, fmt ); - } -} +// void PointWiseMapping::print( const std::string& method, const double & time, OFile& afile, const std::string& fmt ){ +// std::string descr2, descr="DESCRIPTION: results from %s analysis performed at time " + fmt +"\n"; +// afile.printf(descr.c_str(), method.c_str(), time ); +// if(fmt.find("-")!=std::string::npos){ +// descr="REMARK WEIGHT=" + fmt + " %s=" + fmt + " "; descr2="%s=" + fmt; +// } else { +// // This ensures numbers are left justified (i.e. next to the equals sign +// std::size_t psign=fmt.find("%"); +// plumed_assert( psign!=std::string::npos ); +// descr="REMARK WEIGHT=%-" + fmt.substr(psign+1) + " %s=%-" + fmt.substr(psign+1) + " "; +// descr2="%s=%-" + fmt.substr(psign+1); +// } +// for(unsigned i=0;i<frames.size();++i){ +// afile.printf(descr.c_str(), frames[i]->getWeight(), property[0].c_str(), low_dim[i][0] ); +// for(unsigned j=1;j<property.size();++j) afile.printf(descr2.c_str(), property[j].c_str(), low_dim[i][j]); +// afile.printf("\n"); +// frames[i]->print( afile, fmt ); +// } +// } } diff --git a/src/reference/PointWiseMapping.h b/src/reference/PointWiseMapping.h index b5358f5133cc70af63f38879635c67536e7ec14f..950d949e6b4e1124524a0c6c7c2b90a42a0360f2 100644 --- a/src/reference/PointWiseMapping.h +++ b/src/reference/PointWiseMapping.h @@ -71,7 +71,7 @@ public: /// Get a pointer to the matrix of pairwise distances Matrix<double>& modifyDmat(); /// Print out the low dimensional mapping - void print( const std::string& method, const double & time, OFile& afile, const std::string& fmt ); +// void print( const std::string& method, const double & time, OFile& afile, const std::string& fmt ); /// Get the low dimensional embedding coordinate double getProjectionCoordinate( const unsigned& iframe, const unsigned& jcoord ) const ; /// Set the value of the projection coordinate diff --git a/src/reference/ReferenceAtoms.cpp b/src/reference/ReferenceAtoms.cpp index 1fecbaef9b4bf7f5d62db12952d69b3af550597a..c28dc682e8f07b6cdb25c1a9e18bdc50068761cf 100644 --- a/src/reference/ReferenceAtoms.cpp +++ b/src/reference/ReferenceAtoms.cpp @@ -50,11 +50,11 @@ void ReferenceAtoms::setAtomNumbers( const std::vector<AtomNumber>& numbers ){ } } -void ReferenceAtoms::printAtoms( OFile& ofile ) const { +void ReferenceAtoms::printAtoms( const double& lunits, OFile& ofile ) const { for(unsigned i=0;i<reference_atoms.size();++i){ ofile.printf("ATOM %4d X RES %4u %8.3f %8.3f %8.3f %6.2f %6.2f\n", indices[i].serial(), i, - reference_atoms[i][0], reference_atoms[i][1], reference_atoms[i][2], + lunits*reference_atoms[i][0], lunits*reference_atoms[i][1], lunits*reference_atoms[i][2], align[i], displace[i] ); } } diff --git a/src/reference/ReferenceAtoms.h b/src/reference/ReferenceAtoms.h index b8ce2d36ad2b6a4ad8a74f3790005013958a4a3d..5e7e0c12e79cef387520c281df99d9c228f678dc 100644 --- a/src/reference/ReferenceAtoms.h +++ b/src/reference/ReferenceAtoms.h @@ -97,7 +97,7 @@ public: /// Set the positions of the reference atoms virtual void setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in )=0; /// Print the atomic positions - void printAtoms( OFile& ofile ) const ; + void printAtoms( const double& lunits, OFile& ofile ) const ; /// Return all atom indexes const std::vector<AtomNumber>& getAbsoluteIndexes(); /// This returns how many atoms there should be diff --git a/src/reference/ReferenceConfiguration.cpp b/src/reference/ReferenceConfiguration.cpp index dafbc0f4670643bc3cc546db619117eca34cbd3c..530013ef2039928db35a0e7fd58ba6c010f7db0b 100644 --- a/src/reference/ReferenceConfiguration.cpp +++ b/src/reference/ReferenceConfiguration.cpp @@ -156,16 +156,16 @@ double ReferenceConfiguration::calculate( const std::vector<Vector>& pos, const // virialWasSet=ref->virialWasSet; virial=ref->virial; // } -void ReferenceConfiguration::print( OFile& ofile, const double& time, const double& weight, const double& old_norm ){ +void ReferenceConfiguration::print( const double& lunits, OFile& ofile, const double& time, const double& weight, const double& old_norm ){ ofile.printf("REMARK TIME=%f LOG_WEIGHT=%f OLD_NORM=%f\n",time, weight, old_norm ); - print( ofile, "%f" ); // HARD CODED FORMAT HERE AS THIS IS FOR CHECKPOINT FILE + print( lunits, ofile, "%f" ); // HARD CODED FORMAT HERE AS THIS IS FOR CHECKPOINT FILE } -void ReferenceConfiguration::print( OFile& ofile, const std::string& fmt ){ +void ReferenceConfiguration::print( const double& lunits, OFile& ofile, const std::string& fmt ){ ReferenceArguments* args=dynamic_cast<ReferenceArguments*>(this); if(args) args->printArguments( ofile, fmt ); ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>(this); - if(atoms) atoms->printAtoms( ofile ); + if(atoms) atoms->printAtoms( lunits, ofile ); ofile.printf("END\n"); } diff --git a/src/reference/ReferenceConfiguration.h b/src/reference/ReferenceConfiguration.h index 37509736bc5dcdd5ff302100f9cd19cd89f7dbae..6864a19efa5537004b52050c1f8c9ce84e1f53e9 100644 --- a/src/reference/ReferenceConfiguration.h +++ b/src/reference/ReferenceConfiguration.h @@ -142,8 +142,8 @@ public: /// Set the reference structure (perhaps should also pass the pbc and align and displace ) void setReferenceConfig( const std::vector<Vector>& pos, const std::vector<double>& arg, const std::vector<double>& metric ); /// Print a pdb file containing the reference configuration - void print( OFile& ofile, const double& time, const double& weight, const double& old_norm ); - void print( OFile& ofile, const std::string& fmt ); + void print( const double& lunits, OFile& ofile, const double& time, const double& weight, const double& old_norm ); + void print( const double& lunits, OFile& ofile, const std::string& fmt ); /// Get one of the referene arguments virtual double getReferenceArgument( const unsigned& i ) const { plumed_error(); return 0.0; } /// These are overwritten in ReferenceArguments and ReferenceAtoms but are required here