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 41 : void CollectFrames::registerKeywords( Keywords& keys ) {
105 41 : ActionShortcut::registerKeywords( keys );
106 82 : keys.addInputKeyword("optional","ARG","scalar/vector","the labels of the values whose time series you would like to collect for later analysis");
107 41 : keys.add("compulsory","STRIDE","1","the frequency with which data should be stored for analysis. By default data is collected on every step");
108 41 : keys.add("compulsory","CLEAR","0","the frequency with which data should all be deleted and restarted");
109 41 : keys.add("compulsory","ALIGN","OPTIMAL","if storing atoms how would you like the alignment to be done can be SIMPLE/OPTIMAL");
110 41 : keys.add("optional","ATOMS","list of atomic positions that you would like to collect and store for later analysis");
111 41 : keys.add("optional","LOGWEIGHTS","list of actions that calculates log weights that should be used to weight configurations when calculating averages");
112 82 : keys.addOutputComponent("data","default","matrix","the data that is being collected by this action");
113 82 : keys.addOutputComponent("logweights","default","vector","the logarithms of the weights of the data points");
114 41 : keys.needsAction("POSITION");
115 41 : keys.needsAction("CONCATENATE");
116 41 : keys.needsAction("MEAN");
117 41 : keys.needsAction("CUSTOM");
118 41 : keys.needsAction("CONCATENATE");
119 41 : keys.needsAction("COLLECT");
120 41 : keys.needsAction("TRANSPOSE");
121 41 : keys.needsAction("RMSD_VECTOR");
122 41 : keys.needsAction("COMBINE");
123 41 : keys.needsAction("VSTACK");
124 41 : keys.needsAction("CONSTANT");
125 41 : }
126 :
127 56 : std::string CollectFrames::fixArgumentName( const std::string& argin ) {
128 56 : std::string argout = argin;
129 56 : std::size_t dot=argin.find(".");
130 56 : if( dot!=std::string::npos ) {
131 12 : argout = argin.substr(0,dot) + "_" + argin.substr(dot+1);
132 : }
133 56 : return argout;
134 : }
135 :
136 23 : CollectFrames::CollectFrames( const ActionOptions& ao ):
137 : Action(ao),
138 23 : ActionShortcut(ao) {
139 : std::string stride, clearstride;
140 23 : parse("STRIDE",stride);
141 23 : parse("CLEAR",clearstride);
142 23 : std::vector<std::string> argn;
143 23 : parseVector("ARG",argn);
144 23 : std::vector<Value*> theargs;
145 23 : ActionWithArguments::interpretArgumentList( argn, plumed.getActionSet(), this, theargs );
146 : std::string indices;
147 46 : parse("ATOMS",indices);
148 23 : 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 23 : 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 24 : 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 18 : readInputLine( getShortcutLabel() + "_refx: CUSTOM ARG=" + getShortcutLabel() + "_getposx.x," + getShortcutLabel() + "_cposx FUNC=x-y PERIODIC=NO");
165 18 : readInputLine( getShortcutLabel() + "_refy: CUSTOM ARG=" + getShortcutLabel() + "_getposx.y," + getShortcutLabel() + "_cposy FUNC=x-y PERIODIC=NO");
166 18 : readInputLine( getShortcutLabel() + "_refz: CUSTOM ARG=" + getShortcutLabel() + "_getposx.z," + getShortcutLabel() + "_cposz FUNC=x-y PERIODIC=NO");
167 24 : readInputLine( getShortcutLabel() + "_ref: CONCATENATE ARG=" + getShortcutLabel() + "_refx," + getShortcutLabel() + "_refy," + getShortcutLabel() + "_refz");
168 : // Store the reference position in a collect action
169 18 : 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 instantaneous position and the reference position
172 18 : readInputLine( getShortcutLabel() + "_rmsd: RMSD_VECTOR ARG=" + getShortcutLabel() + "_getpos," + getShortcutLabel() + "_refpos DISPLACEMENT SQUARED TYPE=" + align );
173 : // Add the reference position to the RMSD displacement
174 18 : 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 30 : readInputLine( getShortcutLabel() + suffix + ": COLLECT TYPE=matrix ARG=" + getShortcutLabel() + "_fpos STRIDE=" + stride + " CLEAR=" + clearstride );
181 : }
182 :
183 : // Create all the collect actions for arguments
184 51 : for(unsigned i=0; i<theargs.size(); ++i) {
185 28 : if( theargs[i]->getNumberOfValues()!=theargs[0]->getNumberOfValues() ) {
186 0 : error("mismatch between number of arguments calculated by each collected argument");
187 : }
188 140 : readInputLine( getShortcutLabel() + "_" + fixArgumentName( theargs[i]->getName() ) + ": COLLECT ARG=" + theargs[i]->getName() + " STRIDE=" + stride + " CLEAR=" + clearstride );
189 : }
190 : // Make a list of collect actions
191 23 : if( theargs.size()>0 ) {
192 34 : std::string allcol = getShortcutLabel() + "_" + fixArgumentName( theargs[0]->getName() );
193 28 : for(unsigned i=1; i<theargs.size(); ++i) {
194 33 : allcol += "," + getShortcutLabel() + "_" + fixArgumentName( theargs[i]->getName() );
195 : }
196 : // And transfer everything to a matrix
197 17 : std::string suffix = "_argdata";
198 17 : if( indices.length()==0 ) {
199 : suffix = "_data";
200 : }
201 51 : readInputLine( getShortcutLabel() + suffix + ": VSTACK ARG=" + allcol );
202 : }
203 : // Merge all the collected data together into a single matrix
204 23 : 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 23 : std::vector<std::string> logw;
210 23 : parseVector("LOGWEIGHTS",logw);
211 23 : std::vector<Value*> thew;
212 23 : if( logw.size()>0 ) {
213 0 : ActionWithArguments::interpretArgumentList( logw, plumed.getActionSet(), this, thew );
214 : }
215 23 : if( logw.size()>1 ) {
216 0 : error("maximum of one argument should be specified for logweights");
217 : }
218 :
219 23 : if( logw.size()==0 ) {
220 23 : std::string zeros="0";
221 23 : if( theargs.size()>0 ) {
222 17 : for(unsigned i=1; i<theargs[0]->getNumberOfValues(); ++i) {
223 : zeros += ",0";
224 : }
225 : }
226 23 : readInputLine( getShortcutLabel() + "_cweight: CONSTANT VALUE=" + zeros );
227 92 : 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 23 : readInputLine( getShortcutLabel() + "_one: CONSTANT VALUE=1");
236 92 : readInputLine( getShortcutLabel() + "_ones: COLLECT ARG=" + getShortcutLabel() + "_one STRIDE=" + stride + " CLEAR=" + clearstride );
237 69 : }
238 :
239 : }
240 : }
|
