Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 :
28 : //+PLUMEDOC ANALYSIS COLLECT_FRAMES
29 : /*
30 : Collect atomic positions or argument values from the trajectory for later analysis
31 :
32 : This shortcut uses [COLLECT](COLLECT.md) actions to collect data from the trajectory that is amenable for later analysis
33 : using PLUMED's landmark selection actions or dimensionality reduction methods. You can use this method to collect atomic position
34 : data as shown in the following example:
35 :
36 : ```plumed
37 : # This stores the positions of all the first 10 atoms in the system for later analysis
38 : cc: COLLECT_FRAMES ATOMS=1,2,3,4,5,6,7,8,9,10 ALIGN=OPTIMAL
39 :
40 : # This should output the atomic positions for the frames that were collected to a pdb file called traj.pdb
41 : DUMPPDB ATOMS=cc_data ATOM_INDICES=1,2,3,4,5,6,7,8,9,10 FILE=traj.pdb
42 : ```
43 :
44 : If you look at the expanded version of the shortcut above you can see how the position data is stored in a matrix. It is also worth
45 : noting that all the stored structured are aligned to the first frame before being stored. When we take the difference between the
46 : stored configurations we are thus excluding any rotation or the reference frame or translation of the center of mass that has taken place.
47 :
48 : Instead of storing the positions of the atoms you you can store a time series of argument values as shown below:
49 :
50 : ```plumed
51 : t1: TORSION ATOMS=1,2,3,4
52 : t2: TORSION ATOMS=5,6,7,8
53 : t3: TORSION ATOMS=9,10,11,12
54 :
55 : # This collects the three torsion angles
56 : cc: COLLECT_FRAMES ARG=t1,t2,t3 STRIDE=10 CLEAR=1000
57 : DUMPVECTOR ARG=cc.* FILE=timeseries STRIDE=1000
58 : ```
59 :
60 : Notice that the stored data is in a matrix once again here. Further note that we are storing the three torsions on every tenth step. In addition we also delete all the stored data on every 1000th step of trajectory.
61 : The time series for each consecutive 1000 step segment of trajectory will be output to separate files in the above input. These output files also contain weights for each of the stored frames. For the above input
62 : these weights are all equal to one. We can, however, pass the weight as an input argument. This functionality is useful if a bias is acting upon the system and we want to keep track of how much bias is acting upon
63 : each frame of the traejctory as has been done in the input below.
64 :
65 : ```plumed
66 : t1: TORSION ATOMS=1,2,3,4
67 : t2: TORSION ATOMS=5,6,7,8
68 : t3: TORSION ATOMS=9,10,11,12
69 :
70 : r: RESTRAINT ARG=t1 AT=pi/2 KAPPA=10
71 : bw: REWEIGHT_BIAS TEMP=300
72 :
73 : # This collects the three torsion angles
74 : cc: COLLECT_FRAMES ARG=t1,t2,t3 LOGWEIGHTS=bw
75 : DUMPVECTOR ARG=cc.* FILE=timeseries
76 : ```
77 :
78 : You can learn how to use COLLECT_FRAMES to perform dimensionality reduction calculations by working through [this tutorial](https://www.plumed-tutorials.org/lessons/21/006/data/DIMENSIONALITY.html)
79 :
80 : ## Calculating dissimilary matrices
81 :
82 : One of the simplest things that we can do if we have stored a set of trajectory frames using this action is we can calculate the dissimilarity between
83 : every pair of frames we have stored. This is the first step for many of the dimensionality reduction algorithms described in [the dimred module](module_dimred.md).
84 : To calculate these dissimilarities you use the [DISSIMILARITIES](DISSIMILARITIES.md) command.
85 :
86 : Notice that instead of reading in a trajectory and computing dissimilarities you can read in the dissimilarity matrix directly using [CONSTANT](CONSTANT.md) and then perform all analyses with PLUMED.
87 :
88 : */
89 : //+ENDPLUMEDOC
90 :
91 : namespace PLMD {
92 : namespace landmarks {
93 :
94 : class CollectFrames : public ActionShortcut {
95 : private:
96 : std::string fixArgumentName( const std::string& argin );
97 : public:
98 : static void registerKeywords( Keywords& keys );
99 : explicit CollectFrames( const ActionOptions& ao );
100 : };
101 :
102 : PLUMED_REGISTER_ACTION(CollectFrames,"COLLECT_FRAMES")
103 :
104 31 : void CollectFrames::registerKeywords( Keywords& keys ) {
105 31 : ActionShortcut::registerKeywords( keys );
106 62 : keys.addInputKeyword("optional","ARG","scalar/vector","the labels of the values whose time series you would like to collect for later analysis");
107 31 : keys.add("compulsory","STRIDE","1","the frequency with which data should be stored for analysis. By default data is collected on every step");
108 31 : keys.add("compulsory","CLEAR","0","the frequency with which data should all be deleted and restarted");
109 31 : keys.add("compulsory","ALIGN","OPTIMAL","if storing atoms how would you like the alignment to be done can be SIMPLE/OPTIMAL");
110 31 : keys.add("optional","ATOMS","list of atomic positions that you would like to collect and store for later analysis");
111 31 : keys.add("optional","LOGWEIGHTS","list of actions that calculates log weights that should be used to weight configurations when calculating averages");
112 62 : keys.addOutputComponent("data","default","matrix","the data that is being collected by this action");
113 62 : keys.addOutputComponent("logweights","default","vector","the logarithms of the weights of the data points");
114 31 : keys.needsAction("POSITION");
115 31 : keys.needsAction("CONCATENATE");
116 31 : keys.needsAction("MEAN");
117 31 : keys.needsAction("CUSTOM");
118 31 : keys.needsAction("CONCATENATE");
119 31 : keys.needsAction("COLLECT");
120 31 : keys.needsAction("TRANSPOSE");
121 31 : keys.needsAction("RMSD_VECTOR");
122 31 : keys.needsAction("COMBINE");
123 31 : keys.needsAction("VSTACK");
124 31 : keys.needsAction("CONSTANT");
125 31 : }
126 :
127 48 : std::string CollectFrames::fixArgumentName( const std::string& argin ) {
128 48 : std::string argout = argin;
129 48 : std::size_t dot=argin.find(".");
130 48 : if( dot!=std::string::npos ) {
131 8 : argout = argin.substr(0,dot) + "_" + argin.substr(dot+1);
132 : }
133 48 : return argout;
134 : }
135 :
136 19 : CollectFrames::CollectFrames( const ActionOptions& ao ):
137 : Action(ao),
138 19 : ActionShortcut(ao) {
139 : std::string stride, clearstride;
140 19 : parse("STRIDE",stride);
141 38 : parse("CLEAR",clearstride);
142 : std::vector<std::string> argn;
143 38 : parseVector("ARG",argn);
144 : std::vector<Value*> theargs;
145 19 : ActionWithArguments::interpretArgumentList( argn, plumed.getActionSet(), this, theargs );
146 : std::string indices;
147 38 : parse("ATOMS",indices);
148 19 : if( theargs.size()==0 && indices.length()==0 ) {
149 0 : error("no arguments or atoms were specified for collection");
150 : }
151 :
152 : // Create the values to collect the atomic positions
153 19 : if( indices.length()>0 ) {
154 : // Collect reference position
155 12 : readInputLine( getShortcutLabel() + "_getposx: POSITION ATOMS=" + indices );
156 : std::string align;
157 6 : parse("ALIGN",align);
158 12 : readInputLine( getShortcutLabel() + "_getpos: CONCATENATE ARG=" + getShortcutLabel() + "_getposx.x," + getShortcutLabel() + "_getposx.y," + getShortcutLabel() + "_getposx.z");
159 : // Find atomic center
160 12 : readInputLine( getShortcutLabel() + "_cposx: MEAN ARG=" + getShortcutLabel() + "_getposx.x PERIODIC=NO");
161 12 : readInputLine( getShortcutLabel() + "_cposy: MEAN ARG=" + getShortcutLabel() + "_getposx.y PERIODIC=NO");
162 12 : readInputLine( getShortcutLabel() + "_cposz: MEAN ARG=" + getShortcutLabel() + "_getposx.z PERIODIC=NO");
163 : // Subtract atomimc center
164 12 : readInputLine( getShortcutLabel() + "_refx: CUSTOM ARG=" + getShortcutLabel() + "_getposx.x," + getShortcutLabel() + "_cposx FUNC=x-y PERIODIC=NO");
165 12 : readInputLine( getShortcutLabel() + "_refy: CUSTOM ARG=" + getShortcutLabel() + "_getposx.y," + getShortcutLabel() + "_cposy FUNC=x-y PERIODIC=NO");
166 12 : readInputLine( getShortcutLabel() + "_refz: CUSTOM ARG=" + getShortcutLabel() + "_getposx.z," + getShortcutLabel() + "_cposz FUNC=x-y PERIODIC=NO");
167 12 : readInputLine( getShortcutLabel() + "_ref: CONCATENATE ARG=" + getShortcutLabel() + "_refx," + getShortcutLabel() + "_refy," + getShortcutLabel() + "_refz");
168 : // Store the reference position in a collect action
169 12 : readInputLine( getShortcutLabel() + "_refpos: COLLECT TYPE=matrix ARG=" + getShortcutLabel() + "_ref STRIDE=" + clearstride + " CLEAR=" + clearstride );
170 12 : readInputLine( getShortcutLabel() + "_refposT: TRANSPOSE ARG=" + getShortcutLabel() + "_refpos");
171 : // Calculate the RMSD between the instaneous position and the reference position
172 12 : readInputLine( getShortcutLabel() + "_rmsd: RMSD_VECTOR ARG=" + getShortcutLabel() + "_getpos," + getShortcutLabel() + "_refpos DISPLACEMENT SQUARED TYPE=" + align );
173 : // Add the reference position to the RMSD displacement
174 12 : readInputLine( getShortcutLabel() + "_fpos: COMBINE ARG=" + getShortcutLabel() + "_refposT," + getShortcutLabel() + "_rmsd.disp PERIODIC=NO");
175 : // Store the reference data
176 6 : std::string suffix = "_atomdata";
177 6 : if( theargs.size()==0 ) {
178 : suffix = "_data";
179 : }
180 12 : readInputLine( getShortcutLabel() + suffix + ": COLLECT TYPE=matrix ARG=" + getShortcutLabel() + "_fpos STRIDE=" + stride + " CLEAR=" + clearstride );
181 : }
182 :
183 : // Create all the collect actions for arguments
184 43 : for(unsigned i=0; i<theargs.size(); ++i) {
185 24 : if( theargs[i]->getNumberOfValues()!=theargs[0]->getNumberOfValues() ) {
186 0 : error("mismatch between number of arguments calculated by each collected argument");
187 : }
188 48 : readInputLine( getShortcutLabel() + "_" + fixArgumentName( theargs[i]->getName() ) + ": COLLECT ARG=" + theargs[i]->getName() + " STRIDE=" + stride + " CLEAR=" + clearstride );
189 : }
190 : // Make a list of collect actions
191 19 : if( theargs.size()>0 ) {
192 26 : std::string allcol = getShortcutLabel() + "_" + fixArgumentName( theargs[0]->getName() );
193 24 : for(unsigned i=1; i<theargs.size(); ++i) {
194 22 : allcol += "," + getShortcutLabel() + "_" + fixArgumentName( theargs[i]->getName() );
195 : }
196 : // And transfer everything to a matrix
197 13 : std::string suffix = "_argdata";
198 13 : if( indices.length()==0 ) {
199 : suffix = "_data";
200 : }
201 26 : readInputLine( getShortcutLabel() + suffix + ": VSTACK ARG=" + allcol );
202 : }
203 : // Merge all the collected data together into a single matrix
204 19 : if( theargs.size()>0 && indices.length()>0 ) {
205 0 : readInputLine( getShortcutLabel() + "_data: CONCATENATE MATRIX11=" + getShortcutLabel() + "_atomdata MATRIX12=" + getShortcutLabel() + "_argdata");
206 : }
207 :
208 : // Now get the logweights
209 : std::vector<std::string> logw;
210 38 : parseVector("LOGWEIGHTS",logw);
211 : std::vector<Value*> thew;
212 19 : if( logw.size()>0 ) {
213 0 : ActionWithArguments::interpretArgumentList( logw, plumed.getActionSet(), this, thew );
214 : }
215 19 : if( logw.size()>1 ) {
216 0 : error("maximum of one argument should be specified for logweights");
217 : }
218 :
219 19 : if( logw.size()==0 ) {
220 19 : std::string zeros="0";
221 19 : if( theargs.size()>0 ) {
222 13 : for(unsigned i=1; i<theargs[0]->getNumberOfValues(); ++i) {
223 : zeros += ",0";
224 : }
225 : }
226 38 : readInputLine( getShortcutLabel() + "_cweight: CONSTANT VALUE=" + zeros );
227 38 : readInputLine( getShortcutLabel() + "_logweights: COLLECT ARG=" + getShortcutLabel() + "_cweight STRIDE=" + stride + " CLEAR=" + clearstride );
228 : } else {
229 0 : if( theargs[0]->getNumberOfValues()!=thew[0]->getNumberOfValues() ) {
230 0 : error("mismatch between number of weights and number of collected arguments");
231 : }
232 0 : readInputLine( getShortcutLabel() + "_logweights: COLLECT ARG=" + thew[0]->getName() + " STRIDE=" + stride + " CLEAR=" + clearstride );
233 : }
234 : // And finally create a value that contains as many ones as there are data points (this is used if we want to do Classical MDS
235 19 : readInputLine( getShortcutLabel() + "_one: CONSTANT VALUE=1");
236 38 : readInputLine( getShortcutLabel() + "_ones: COLLECT ARG=" + getShortcutLabel() + "_one STRIDE=" + stride + " CLEAR=" + clearstride );
237 57 : }
238 :
239 : }
240 : }
|
