Skip to content
Snippets Groups Projects
Commit 61e9e629 authored by Gareth Tribello's avatar Gareth Tribello
Browse files

Added new functionality for output pdb in analysis

You can now use COLLECT_FRAMES to collect pdb data from the trajectory
and output a pdb file that contains the underlying atomic coordinates
in pdb format even if you are not using the rmsd distance to measure
the dissimilarity matrix for dimensionality reduction.  Also fixed
a number of bugs.
parent 32e150c5
No related branches found
No related tags found
No related merge requests found
Showing
with 2761 additions and 66 deletions
include ../../scripts/test.make
type=driver
# this is to test a different name
arg="--plumed plumed.dat --trajectory-stride 10 --timestep 0.005 --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt=%8.4f"
extra_files="../../trajectories/trajectory.xyz"
d1: DISTANCE ATOMS=1,2
d2: DISTANCE ATOMS=3,4
c1: COLLECT_FRAMES ATOMS=1-20 STRIDE=1
r1: EUCLIDEAN_DISSIMILARITIES FRAMES=c1 ARG=d1,d2
l1: LANDMARK_SELECT_STRIDE USE_OUTPUT_DATA_FROM=r1 NLANDMARKS=5
OUTPUT_ANALYSIS_DATA_TO_PDB USE_OUTPUT_DATA_FROM=l1 FILE=configs-stride.pdb FMT=%8.4f
c1: READ_ANALYSIS_FRAMES ATOMS=1-20 STRIDE=1 c1: COLLECT_FRAMES ATOMS=1-20 STRIDE=1
r1: READ_DISSIMILARITY_MATRIX TRAJ=c1 FILE=mymatrix.dat r1: READ_DISSIMILARITY_MATRIX FRAMES=c1 FILE=mymatrix.dat
PRINT_DISSIMILARITY_MATRIX USE_OUTPUT_DATA_FROM=r1 FILE=mymatrix_out.dat FMT=%8.4f PRINT_DISSIMILARITY_MATRIX USE_OUTPUT_DATA_FROM=r1 FILE=mymatrix_out.dat FMT=%8.4f
l1: LANDMARK_SELECT_STRIDE USE_OUTPUT_DATA_FROM=r1 NLANDMARKS=5 l1: LANDMARK_SELECT_STRIDE USE_OUTPUT_DATA_FROM=r1 NLANDMARKS=5
......
include ../scripts/module.make
include ../../scripts/test.make
This diff is collapsed.
#! FIELDS mds.1 mds.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.0302 -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
type=simplemd
This diff is collapsed.
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
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
#! FIELDS mds.1 mds.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.0410 0.0086
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.1212 -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
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
ff: COLLECT_FRAMES ATOMS=1-7 STRIDE=10 RUN=1000
oo: EUCLIDEAN_DISSIMILARITIES ARG=c1.moment-2,c1.moment-3 FRAMES=ff
CLASSICAL_MDS ...
USE_OUTPUT_DATA_FROM=oo
NLOW_DIM=2
LABEL=mds
... CLASSICAL_MDS
OUTPUT_ANALYSIS_DATA_TO_COLVAR USE_OUTPUT_DATA_FROM=mds FILE=list_embed FMT=%8.4f
OUTPUT_ANALYSIS_DATA_TO_PDB USE_OUTPUT_DATA_FROM=mds FILE=embed FMT=%8.4f
...@@ -52,7 +52,7 @@ mydata(NULL) ...@@ -52,7 +52,7 @@ mydata(NULL)
if( datastr.length()>0 ){ if( datastr.length()>0 ){
mydata=plumed.getActionSet().selectWithLabel<AnalysisBase*>( datastr ); mydata=plumed.getActionSet().selectWithLabel<AnalysisBase*>( datastr );
ReadAnalysisFrames* checkt = dynamic_cast<ReadAnalysisFrames*>( mydata ); ReadAnalysisFrames* checkt = dynamic_cast<ReadAnalysisFrames*>( mydata );
if( checkt ) error("READ_ANALYSIS_FRAMES should only be used in association with READ_DISSSIMILARITY_MATRIX"); if( checkt ) error("READ_ANALYSIS_FRAMES should only be used in input to the FRAMES keyword");
log.printf(" performing analysis on output from %s \n",datastr.c_str() ); log.printf(" performing analysis on output from %s \n",datastr.c_str() );
if( !mydata ) error("could not find analysis action named " + datastr ); if( !mydata ) error("could not find analysis action named " + datastr );
freq=mydata->freq; use_all_data=mydata->use_all_data; freq=mydata->freq; use_all_data=mydata->use_all_data;
......
...@@ -85,8 +85,6 @@ public: ...@@ -85,8 +85,6 @@ public:
virtual void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ; virtual void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ;
/// Get a reference configuration (in dimensionality reduction this returns the projection) /// Get a reference configuration (in dimensionality reduction this returns the projection)
virtual ReferenceConfiguration* getReferenceConfiguration( const unsigned& idata ); virtual ReferenceConfiguration* getReferenceConfiguration( const unsigned& idata );
/// Get the input configuration (in dimensionality reduction this returns the high dimensional configuration)
virtual ReferenceConfiguration* getInputReferenceConfiguration( const unsigned& idata );
/// This actually performs the analysis /// This actually performs the analysis
virtual void performAnalysis()=0; virtual void performAnalysis()=0;
/// These overwrite things from inherited classes (this is a bit of a fudge) /// These overwrite things from inherited classes (this is a bit of a fudge)
...@@ -170,11 +168,6 @@ ReferenceConfiguration* AnalysisBase::getReferenceConfiguration( const unsigned& ...@@ -170,11 +168,6 @@ ReferenceConfiguration* AnalysisBase::getReferenceConfiguration( const unsigned&
return mydata->getReferenceConfiguration( idata ); return mydata->getReferenceConfiguration( idata );
} }
inline
ReferenceConfiguration* AnalysisBase::getInputReferenceConfiguration( const unsigned& idata ){
return mydata->getInputReferenceConfiguration( idata );
}
} }
} }
......
...@@ -46,6 +46,9 @@ void AnalysisWithDataCollection::registerKeywords( Keywords& keys ){ ...@@ -46,6 +46,9 @@ void AnalysisWithDataCollection::registerKeywords( Keywords& keys ){
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("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 + RUN/USE_ALL_DATA. " 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 + RUN/USE_ALL_DATA. "
"For more information see \\ref analysisbas"); "For more information see \\ref analysisbas");
keys.reserve("optional","FRAMES","with this keyword you can specify the atomic configurations that you would like to store using "
"an COLLECT_FRAMES action. When the projection is output in dimensionality reduction you will then "
"print the underlying atoms and their projection");
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 + RUN/USE_ALL_DATA. " 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 + RUN/USE_ALL_DATA. "
"For more information see \\ref analysisbas"); "For more information see \\ref analysisbas");
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("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.");
...@@ -63,13 +66,25 @@ write_chq(false), ...@@ -63,13 +66,25 @@ write_chq(false),
rtemp(0), rtemp(0),
idata(0), idata(0),
firstAnalysisDone(false), firstAnalysisDone(false),
old_norm(0.0) old_norm(0.0),
myframes(NULL)
{ {
if( !mydata ){ if( !mydata ){
// Check for FRAMES keyword that allows us to read reference configurations from elsewhere
std::string instring; if( keywords.exists("FRAMES") ) parse("FRAMES",instring);
if( instring.length()>0 ){
if( mydata ) error("frames keyword is not compatible with USE_OUTPUT_DATA_FROM");
myframes=plumed.getActionSet().selectWithLabel<ReadAnalysisFrames*>( instring );
if( !myframes ) error( instring + " is not the name of an object that collects trajectory data");
log.printf(" frames are stored by object with label %s \n",instring.c_str() );
freq=myframes->freq; use_all_data=myframes->use_all_data; setStride( myframes->getStride() );
}
// Check if we are using the input data from another action // Check if we are using the input data from another action
std::string datastr; std::string datastr;
if( keywords.exists("REUSE_INPUT_DATA_FROM") ) parse("REUSE_INPUT_DATA_FROM",datastr); if( keywords.exists("REUSE_INPUT_DATA_FROM") ) parse("REUSE_INPUT_DATA_FROM",datastr);
if( datastr.length()>0 ) { if( datastr.length()>0 ) {
if( myframes ) error("FRAMES command is incompatible with REUSE_INPUT_DATA_FROM");
AnalysisWithDataCollection* checkd = plumed.getActionSet().selectWithLabel<AnalysisWithDataCollection*>( datastr ); 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"); if( !checkd) error("cannot reuse input data from action with label " + datastr + " as this does not store data");
ReadAnalysisFrames* checkt = dynamic_cast<ReadAnalysisFrames*>( checkd ); ReadAnalysisFrames* checkt = dynamic_cast<ReadAnalysisFrames*>( checkd );
...@@ -129,31 +144,49 @@ old_norm(0.0) ...@@ -129,31 +144,49 @@ old_norm(0.0)
metricname=""; metricname="";
} }
// Read in the information about how often to run the analysis (storage is read in in ActionPilot.cpp) bool dobias;
if( keywords.exists("USE_ALL_DATA") ){ if( !myframes ){
if( !keywords.exists("RUN") ) use_all_data=true; // Read in the information about how often to run the analysis (storage is read in in ActionPilot.cpp)
else parseFlag("USE_ALL_DATA",use_all_data); if( keywords.exists("USE_ALL_DATA") ){
} if( !keywords.exists("RUN") ) use_all_data=true;
if(!use_all_data){ else parseFlag("USE_ALL_DATA",use_all_data);
if( keywords.exists("RUN") ) parse("RUN",freq); }
// Setup everything given the ammount of data that we will have in each analysis if(!use_all_data){
if( freq%getStride()!= 0 ) error("Frequncy of running is not a multiple of the stride"); if( keywords.exists("RUN") ) parse("RUN",freq);
unsigned ndata=freq/getStride(); data.resize(ndata); logweights.resize( ndata ); // Setup everything given the ammount of data that we will have in each analysis
for(unsigned i=0;i<ndata;++i) data[i]=metricRegister().create<ReferenceConfiguration>( metricname ); if( freq%getStride()!= 0 ) error("Frequncy of running is not a multiple of the stride");
log.printf(" running analysis every %u steps\n",freq); // Check if we are doing block averaging
// Check if we are doing block averaging nomemory=false;
nomemory=false; if( keywords.exists("NOMEMORY") ) parseFlag("NOMEMORY",nomemory);
if( keywords.exists("NOMEMORY") ) parseFlag("NOMEMORY",nomemory); if(nomemory) log.printf(" doing block averaging and analysing each portion of trajectory separately\n");
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
if( keywords.exists("REWEIGHT_BIAS") ) parseFlag("REWEIGHT_BIAS",dobias);
// Reweighting for temperatures
rtemp=0;
if( keywords.exists("REWEIGHT_TEMP") ) 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;
if( keywords.exists("TEMP") ) 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");
} else { } else {
log.printf(" analysing all data in trajectory\n"); nomemory = myframes->nomemory; dobias = (myframes->biases.size()>0);
rtemp = myframes->rtemp; simtemp = myframes->simtemp;
} }
// Read in stuff for reweighting of trajectories // Setup bias reweighting
// Reweighting for biases
bool dobias;
if( keywords.exists("REWEIGHT_BIAS") ) parseFlag("REWEIGHT_BIAS",dobias);
if( dobias ){ if( dobias ){
std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
if( all.empty() ) error("your input file is not telling plumed to calculate anything"); if( all.empty() ) error("your input file is not telling plumed to calculate anything");
...@@ -171,20 +204,13 @@ old_norm(0.0) ...@@ -171,20 +204,13 @@ old_norm(0.0)
if( biases.empty() ) error("you are asking to reweight bias but there does not appear to be a bias acting on your system"); 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 ); requestArguments( arg );
} }
// Setup reference configuration objects to store data
// Reweighting for temperatures if( !use_all_data ){
rtemp=0; plumed_assert( !mydata );
if( keywords.exists("REWEIGHT_TEMP") ) parse("REWEIGHT_TEMP",rtemp); log.printf(" running analysis every %u steps\n",freq);
if( rtemp!=0 ){ unsigned ndata=freq/getStride(); data.resize(ndata); logweights.resize( ndata );
rtemp*=plumed.getAtoms().getKBoltzmann(); for(unsigned i=0;i<ndata;++i) data[i]=metricRegister().create<ReferenceConfiguration>( metricname );
log.printf(" reweighting simulation to probabilities at temperature %f\n",rtemp);
} }
// Now retrieve the temperature in the simulation
simtemp=0;
if( keywords.exists("TEMP") ) 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 // Check if a check point is required (this should be got rid of at some point when we have proper checkpointing) GAT
if( keywords.exists("WRITE_CHECKPOINT") ) parseFlag("WRITE_CHECKPOINT",write_chq); if( keywords.exists("WRITE_CHECKPOINT") ) parseFlag("WRITE_CHECKPOINT",write_chq);
...@@ -207,8 +233,8 @@ old_norm(0.0) ...@@ -207,8 +233,8 @@ old_norm(0.0)
for(unsigned i=0;i<getNumberOfArguments();++i) rfile.setupPrintValue( getPntrToArgument(i) ); 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 // 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(){ AnalysisWithDataCollection::~AnalysisWithDataCollection(){
...@@ -264,15 +290,14 @@ void AnalysisWithDataCollection::getDataPoint( const unsigned& idat, std::vector ...@@ -264,15 +290,14 @@ void AnalysisWithDataCollection::getDataPoint( const unsigned& idat, std::vector
} }
ReferenceConfiguration* AnalysisWithDataCollection::getReferenceConfiguration( const unsigned& idat ){ ReferenceConfiguration* AnalysisWithDataCollection::getReferenceConfiguration( const unsigned& idat ){
if( !mydata ){ plumed_dbg_assert( idat<data.size() ); return data[idat]; } if( !mydata ){
plumed_dbg_assert( idat<data.size() );
if( myframes ) return myframes->getReferenceConfiguration( idat );
return data[idat];
}
return AnalysisBase::getReferenceConfiguration( idat ); return AnalysisBase::getReferenceConfiguration( idat );
} }
ReferenceConfiguration* AnalysisWithDataCollection::getInputReferenceConfiguration( const unsigned& idat ){
if( !mydata ){ plumed_dbg_assert( idat<data.size() ); return data[idat]; }
return AnalysisBase::getInputReferenceConfiguration( idat );
}
void AnalysisWithDataCollection::update(){ void AnalysisWithDataCollection::update(){
if( mydata ){ AnalysisBase::update(); return; } 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)? // 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)?
......
...@@ -31,7 +31,11 @@ class ReferenceConfiguration; ...@@ -31,7 +31,11 @@ class ReferenceConfiguration;
namespace analysis { namespace analysis {
class ReadAnalysisFrames;
class AnalysisWithDataCollection : public AnalysisBase { class AnalysisWithDataCollection : public AnalysisBase {
friend class EuclideanDissimilarityMatrix;
friend class ReadDissimilarityMatrix;
private: private:
/// Are we treating each block of data separately /// Are we treating each block of data separately
bool nomemory; bool nomemory;
...@@ -53,6 +57,9 @@ private: ...@@ -53,6 +57,9 @@ private:
bool firstAnalysisDone; bool firstAnalysisDone;
/// The value of the old normalization constant /// The value of the old normalization constant
double norm, old_norm; double norm, old_norm;
/// This allows you to output trajectory data with the projections
/// in EUCLIDEAN_DISSIMLAIRITIES and READ_DISSIMILARITY_MATRIX
ReadAnalysisFrames* myframes;
/// Data is collected from the trajectory by passing it to this pdb. These pdb /// Data is collected from the trajectory by passing it to this pdb. These pdb
/// files are then read by the ReferenceConfigurations in data /// files are then read by the ReferenceConfigurations in data
PDB mypdb; PDB mypdb;
...@@ -93,8 +100,6 @@ public: ...@@ -93,8 +100,6 @@ public:
virtual void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ; virtual void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ;
/// Get a reference configuration (in dimensionality reduction this returns the projection) /// Get a reference configuration (in dimensionality reduction this returns the projection)
virtual ReferenceConfiguration* getReferenceConfiguration( const unsigned& idat ); virtual ReferenceConfiguration* getReferenceConfiguration( const unsigned& idat );
/// Get the underlying reference configuration (in dimensionality reduction this return the high dimensional point)
ReferenceConfiguration* getInputReferenceConfiguration( const unsigned& idat );
/// This ensures that the energy is stored if we are reweighting /// This ensures that the energy is stored if we are reweighting
void prepare(); void prepare();
/// This stores the data and calls the analysis to be performed /// This stores the data and calls the analysis to be performed
......
...@@ -55,7 +55,7 @@ PLUMED_REGISTER_ACTION(EuclideanDissimilarityMatrix,"EUCLIDEAN_DISSIMILARITIES") ...@@ -55,7 +55,7 @@ PLUMED_REGISTER_ACTION(EuclideanDissimilarityMatrix,"EUCLIDEAN_DISSIMILARITIES")
void EuclideanDissimilarityMatrix::registerKeywords( Keywords& keys ){ void EuclideanDissimilarityMatrix::registerKeywords( Keywords& keys ){
AnalysisWithDataCollection::registerKeywords( keys ); AnalysisWithDataCollection::registerKeywords( keys );
keys.reset_style("METRIC","atoms-1"); keys.reset_style("METRIC","atoms-1"); keys.use("FRAMES");
} }
EuclideanDissimilarityMatrix::EuclideanDissimilarityMatrix( const ActionOptions& ao ): EuclideanDissimilarityMatrix::EuclideanDissimilarityMatrix( const ActionOptions& ao ):
...@@ -72,8 +72,15 @@ void EuclideanDissimilarityMatrix::performAnalysis(){ ...@@ -72,8 +72,15 @@ void EuclideanDissimilarityMatrix::performAnalysis(){
double EuclideanDissimilarityMatrix::getDissimilarity( const unsigned& iframe, const unsigned& jframe ){ double EuclideanDissimilarityMatrix::getDissimilarity( const unsigned& iframe, const unsigned& jframe ){
plumed_dbg_assert( iframe<dissimilarities.nrows() && jframe<dissimilarities.ncols() ); plumed_dbg_assert( iframe<dissimilarities.nrows() && jframe<dissimilarities.ncols() );
if( dissimilarities(iframe,jframe)>0. ){ return dissimilarities(iframe,jframe); } if( dissimilarities(iframe,jframe)>0. ){ return dissimilarities(iframe,jframe); }
if( iframe!=jframe ){ if( iframe!=jframe ){
dissimilarities(iframe,jframe) = dissimilarities(jframe,iframe) = distance( getPbc(), getArguments(), getReferenceConfiguration(iframe), getReferenceConfiguration(jframe), true ); ReferenceConfiguration* myref1; ReferenceConfiguration* myref2;
if( mydata ){ myref1=AnalysisBase::getReferenceConfiguration(iframe); myref2=AnalysisBase::getReferenceConfiguration(jframe); }
else { myref1 = data[iframe]; myref2 = data[jframe]; }
if( myref1->getNumberOfProperties()>0 ){
dissimilarities(iframe,jframe) = dissimilarities(jframe,iframe) = property_distance( myref1, myref2, true );
} else {
dissimilarities(iframe,jframe) = dissimilarities(jframe,iframe) = distance( getPbc(), getArguments(), myref1, myref2, true );
}
return dissimilarities(iframe,jframe); return dissimilarities(iframe,jframe);
} }
return 0.0; return 0.0;
......
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
#include "core/ActionSet.h" #include "core/ActionSet.h"
#include "core/ActionRegister.h" #include "core/ActionRegister.h"
//+PLUMEDOC ANALYSIS READ_ANALYSIS_FRAMES //+PLUMEDOC ANALYSIS COLLECT_FRAMES
/* /*
This allows you to convert a trajectory and a dissimilarity matrix into a dissimilarity object This allows you to convert a trajectory and a dissimilarity matrix into a dissimilarity object
...@@ -36,21 +36,17 @@ This allows you to convert a trajectory and a dissimilarity matrix into a dissim ...@@ -36,21 +36,17 @@ This allows you to convert a trajectory and a dissimilarity matrix into a dissim
namespace PLMD { namespace PLMD {
namespace analysis { namespace analysis {
PLUMED_REGISTER_ACTION(ReadAnalysisFrames,"READ_ANALYSIS_FRAMES") PLUMED_REGISTER_ACTION(ReadAnalysisFrames,"COLLECT_FRAMES")
void ReadAnalysisFrames::registerKeywords( Keywords& keys ){ void ReadAnalysisFrames::registerKeywords( Keywords& keys ){
AnalysisWithDataCollection::registerKeywords( keys ); AnalysisWithDataCollection::registerKeywords( keys );
keys.reset_style("USE_ALL_DATA","hidden"); keys.remove("RUN"); keys.remove("REWEIGHT_BIAS"); keys.remove("REWEIGHT_TEMP"); keys.remove("ARG"); keys.remove("SERIAL"); keys.remove("USE_OUTPUT_DATA_FROM");
keys.remove("TEMP"); keys.remove("WRITE_CHECKPOINT"); keys.remove("NOMEMORY"); keys.remove("RESTART");
keys.remove("UPDATE_FROM"); keys.remove("UPDATE_UNTIL"); keys.remove("USE_OUTPUT_DATA_FROM");
keys.remove("ARG"); keys.remove("SERIAL"); keys.remove("REUSE_INPUT_DATA_FROM");
} }
ReadAnalysisFrames::ReadAnalysisFrames( const ActionOptions& ao ): ReadAnalysisFrames::ReadAnalysisFrames( const ActionOptions& ao ):
Action(ao), Action(ao),
AnalysisWithDataCollection(ao) AnalysisWithDataCollection(ao)
{ {
if( plumed.getActionSet().size()!=0 ) error("read analysis frames command must be at top of input file");
} }
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment