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 : */
81 : //+ENDPLUMEDOC
82 :
83 : namespace PLMD {
84 : namespace landmarks {
85 :
86 : class CollectFrames : public ActionShortcut {
87 : private:
88 : std::string fixArgumentName( const std::string& argin );
89 : public:
90 : static void registerKeywords( Keywords& keys );
91 : explicit CollectFrames( const ActionOptions& ao );
92 : };
93 :
94 : PLUMED_REGISTER_ACTION(CollectFrames,"COLLECT_FRAMES")
95 :
96 31 : void CollectFrames::registerKeywords( Keywords& keys ) {
97 31 : ActionShortcut::registerKeywords( keys );
98 62 : keys.addInputKeyword("optional","ARG","scalar/vector","the labels of the values whose time series you would like to collect for later analysis");
99 31 : keys.add("compulsory","STRIDE","1","the frequency with which data should be stored for analysis. By default data is collected on every step");
100 31 : keys.add("compulsory","CLEAR","0","the frequency with which data should all be deleted and restarted");
101 31 : keys.add("compulsory","ALIGN","OPTIMAL","if storing atoms how would you like the alignment to be done can be SIMPLE/OPTIMAL");
102 31 : keys.add("optional","ATOMS","list of atomic positions that you would like to collect and store for later analysis");
103 31 : keys.add("optional","LOGWEIGHTS","list of actions that calculates log weights that should be used to weight configurations when calculating averages");
104 62 : keys.addOutputComponent("data","default","matrix","the data that is being collected by this action");
105 62 : keys.addOutputComponent("logweights","default","vector","the logarithms of the weights of the data points");
106 31 : keys.needsAction("POSITION");
107 31 : keys.needsAction("CONCATENATE");
108 31 : keys.needsAction("MEAN");
109 31 : keys.needsAction("CUSTOM");
110 31 : keys.needsAction("CONCATENATE");
111 31 : keys.needsAction("COLLECT");
112 31 : keys.needsAction("TRANSPOSE");
113 31 : keys.needsAction("RMSD_VECTOR");
114 31 : keys.needsAction("COMBINE");
115 31 : keys.needsAction("VSTACK");
116 31 : keys.needsAction("CONSTANT");
117 31 : }
118 :
119 48 : std::string CollectFrames::fixArgumentName( const std::string& argin ) {
120 48 : std::string argout = argin;
121 48 : std::size_t dot=argin.find(".");
122 48 : if( dot!=std::string::npos ) {
123 8 : argout = argin.substr(0,dot) + "_" + argin.substr(dot+1);
124 : }
125 48 : return argout;
126 : }
127 :
128 19 : CollectFrames::CollectFrames( const ActionOptions& ao ):
129 : Action(ao),
130 19 : ActionShortcut(ao) {
131 : std::string stride, clearstride;
132 19 : parse("STRIDE",stride);
133 38 : parse("CLEAR",clearstride);
134 : std::vector<std::string> argn;
135 38 : parseVector("ARG",argn);
136 : std::vector<Value*> theargs;
137 19 : ActionWithArguments::interpretArgumentList( argn, plumed.getActionSet(), this, theargs );
138 : std::string indices;
139 38 : parse("ATOMS",indices);
140 19 : if( theargs.size()==0 && indices.length()==0 ) {
141 0 : error("no arguments or atoms were specified for collection");
142 : }
143 :
144 : // Create the values to collect the atomic positions
145 19 : if( indices.length()>0 ) {
146 : // Collect reference position
147 12 : readInputLine( getShortcutLabel() + "_getposx: POSITION ATOMS=" + indices );
148 : std::string align;
149 6 : parse("ALIGN",align);
150 12 : readInputLine( getShortcutLabel() + "_getpos: CONCATENATE ARG=" + getShortcutLabel() + "_getposx.x," + getShortcutLabel() + "_getposx.y," + getShortcutLabel() + "_getposx.z");
151 : // Find atomic center
152 12 : readInputLine( getShortcutLabel() + "_cposx: MEAN ARG=" + getShortcutLabel() + "_getposx.x PERIODIC=NO");
153 12 : readInputLine( getShortcutLabel() + "_cposy: MEAN ARG=" + getShortcutLabel() + "_getposx.y PERIODIC=NO");
154 12 : readInputLine( getShortcutLabel() + "_cposz: MEAN ARG=" + getShortcutLabel() + "_getposx.z PERIODIC=NO");
155 : // Subtract atomimc center
156 12 : readInputLine( getShortcutLabel() + "_refx: CUSTOM ARG=" + getShortcutLabel() + "_getposx.x," + getShortcutLabel() + "_cposx FUNC=x-y PERIODIC=NO");
157 12 : readInputLine( getShortcutLabel() + "_refy: CUSTOM ARG=" + getShortcutLabel() + "_getposx.y," + getShortcutLabel() + "_cposy FUNC=x-y PERIODIC=NO");
158 12 : readInputLine( getShortcutLabel() + "_refz: CUSTOM ARG=" + getShortcutLabel() + "_getposx.z," + getShortcutLabel() + "_cposz FUNC=x-y PERIODIC=NO");
159 12 : readInputLine( getShortcutLabel() + "_ref: CONCATENATE ARG=" + getShortcutLabel() + "_refx," + getShortcutLabel() + "_refy," + getShortcutLabel() + "_refz");
160 : // Store the reference position in a collect action
161 12 : readInputLine( getShortcutLabel() + "_refpos: COLLECT TYPE=matrix ARG=" + getShortcutLabel() + "_ref STRIDE=" + clearstride + " CLEAR=" + clearstride );
162 12 : readInputLine( getShortcutLabel() + "_refposT: TRANSPOSE ARG=" + getShortcutLabel() + "_refpos");
163 : // Calculate the RMSD between the instaneous position and the reference position
164 12 : readInputLine( getShortcutLabel() + "_rmsd: RMSD_VECTOR ARG=" + getShortcutLabel() + "_getpos," + getShortcutLabel() + "_refpos DISPLACEMENT SQUARED TYPE=" + align );
165 : // Add the reference position to the RMSD displacement
166 12 : readInputLine( getShortcutLabel() + "_fpos: COMBINE ARG=" + getShortcutLabel() + "_refposT," + getShortcutLabel() + "_rmsd.disp PERIODIC=NO");
167 : // Store the reference data
168 6 : std::string suffix = "_atomdata";
169 6 : if( theargs.size()==0 ) {
170 : suffix = "_data";
171 : }
172 12 : readInputLine( getShortcutLabel() + suffix + ": COLLECT TYPE=matrix ARG=" + getShortcutLabel() + "_fpos STRIDE=" + stride + " CLEAR=" + clearstride );
173 : }
174 :
175 : // Create all the collect actions for arguments
176 43 : for(unsigned i=0; i<theargs.size(); ++i) {
177 24 : if( theargs[i]->getNumberOfValues()!=theargs[0]->getNumberOfValues() ) {
178 0 : error("mismatch between number of arguments calculated by each collected argument");
179 : }
180 48 : readInputLine( getShortcutLabel() + "_" + fixArgumentName( theargs[i]->getName() ) + ": COLLECT ARG=" + theargs[i]->getName() + " STRIDE=" + stride + " CLEAR=" + clearstride );
181 : }
182 : // Make a list of collect actions
183 19 : if( theargs.size()>0 ) {
184 26 : std::string allcol = getShortcutLabel() + "_" + fixArgumentName( theargs[0]->getName() );
185 24 : for(unsigned i=1; i<theargs.size(); ++i) {
186 22 : allcol += "," + getShortcutLabel() + "_" + fixArgumentName( theargs[i]->getName() );
187 : }
188 : // And transfer everything to a matrix
189 13 : std::string suffix = "_argdata";
190 13 : if( indices.length()==0 ) {
191 : suffix = "_data";
192 : }
193 26 : readInputLine( getShortcutLabel() + suffix + ": VSTACK ARG=" + allcol );
194 : }
195 : // Merge all the collected data together into a single matrix
196 19 : if( theargs.size()>0 && indices.length()>0 ) {
197 0 : readInputLine( getShortcutLabel() + "_data: CONCATENATE MATRIX11=" + getShortcutLabel() + "_atomdata MATRIX12=" + getShortcutLabel() + "_argdata");
198 : }
199 :
200 : // Now get the logweights
201 : std::vector<std::string> logw;
202 38 : parseVector("LOGWEIGHTS",logw);
203 : std::vector<Value*> thew;
204 19 : if( logw.size()>0 ) {
205 0 : ActionWithArguments::interpretArgumentList( logw, plumed.getActionSet(), this, thew );
206 : }
207 19 : if( logw.size()>1 ) {
208 0 : error("maximum of one argument should be specified for logweights");
209 : }
210 :
211 19 : if( logw.size()==0 ) {
212 19 : std::string zeros="0";
213 19 : if( theargs.size()>0 ) {
214 13 : for(unsigned i=1; i<theargs[0]->getNumberOfValues(); ++i) {
215 : zeros += ",0";
216 : }
217 : }
218 38 : readInputLine( getShortcutLabel() + "_cweight: CONSTANT VALUE=" + zeros );
219 38 : readInputLine( getShortcutLabel() + "_logweights: COLLECT ARG=" + getShortcutLabel() + "_cweight STRIDE=" + stride + " CLEAR=" + clearstride );
220 : } else {
221 0 : if( theargs[0]->getNumberOfValues()!=thew[0]->getNumberOfValues() ) {
222 0 : error("mismatch between number of weights and number of collected arguments");
223 : }
224 0 : readInputLine( getShortcutLabel() + "_logweights: COLLECT ARG=" + thew[0]->getName() + " STRIDE=" + stride + " CLEAR=" + clearstride );
225 : }
226 : // 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
227 19 : readInputLine( getShortcutLabel() + "_one: CONSTANT VALUE=1");
228 38 : readInputLine( getShortcutLabel() + "_ones: COLLECT ARG=" + getShortcutLabel() + "_one STRIDE=" + stride + " CLEAR=" + clearstride );
229 57 : }
230 :
231 : }
232 : }
|
