Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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/ActionRegister.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/ActionSet.h"
25 : #include "core/ActionShortcut.h"
26 : #include "core/ActionWithValue.h"
27 :
28 : //+PLUMEDOC MCOLVAR DISPLACEMENT
29 : /*
30 : Calculate the displacement vector between the pair of input vectors
31 :
32 : This shortcut can be used to calculate the vector of displacements between two input vectors as shown below.
33 :
34 : ```plumed
35 : c: CONSTANT VALUES=1,2,3
36 : d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
37 : dd: DISPLACEMENT ARG1=c ARG2=d
38 : PRINT ARG=dd FILE=colvar
39 : ```
40 :
41 : The output here, `dd`, is a $1 \times 3$ matrix for reasons that will become clear later in this documentation.
42 : Notice that we can obtain the same result by specifying the input vectors here as two sets of three scalars as shown
43 : below:
44 :
45 : ```plumed
46 : c1: CONSTANT VALUE=1
47 : d1: DISTANCE ATOMS=1,2
48 : c2: CONSTANT VALUE=2
49 : d2: DISTANCE ATOMS=3,4
50 : c3: CONSTANT VALUE=3
51 : d3: DISTANCE ATOMS=5,6
52 : dd: DISPLACEMENT ARG1=c1,c2,c3 ARG2=d1,d2,d3
53 : PRINT ARG=dd FILE=colvar
54 : ```
55 :
56 : The DISPLACEMENT command that has been introduced in the above inputs is primarily used within the [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md),
57 : [NORMALIZED_EUCLIDEAN_DISTANCE](NORMALIZED_EUCLIDEAN_DISTANCE.md) and [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md) shortcuts. If the $1 \times N$ matrix
58 : of displacements that that we obtainfrom these commands is, $D$, these three actions calculate
59 :
60 : $$
61 : d = D M D^T
62 : $$
63 :
64 : The $N \times N$ matrix $M$ here is the identity if you are using [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md), a diagonal matrix if you are using
65 : [NORMALIZED_EUCLIDEAN_DISTANCE](NORMALIZED_EUCLIDEAN_DISTANCE.md) and a full matrix if you are computing the [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md).
66 :
67 : ## Calculating multiple displacement vectors
68 :
69 : The reason the output of DISPLACEMENT is a $1 \times 3$ matrix here becomes clearer once we consider the following input:
70 :
71 : ```plumed
72 : ref_psi: CONSTANT VALUES=2.25
73 : ref_phi: CONSTANT VALUES=-1.91
74 :
75 : psi: TORSION ATOMS1=1,2,3,4 ATOMS2=5,6,7,8 ATOMS3=9,10,11,12
76 : phi: TORSION ATOMS1=13,14,15,16 ATOMS2=17,18,19,20 ATOMS3=21,22,23,24
77 :
78 : dd: DISPLACEMENT ARG1=psi,phi ARG2=ref_psi,ref_phi
79 : PRINT ARG=dd FILE=colvar
80 : ```
81 :
82 : The output from the input above is a $3\times 2$ matrix. The rows of this matrix run over the 3 different
83 : torsion values that have been specified in the `psi` and `phi` commands. The first column of the matrix
84 : contains the differences between each of the instantaneous `psi` aingles and the reference value for this
85 : angle, while the second columns contains the differences between the `phi` angles and the reference.
86 :
87 : In other words, we can calculate multiple displacement vectors at once as each row of the final output matrix will
88 : contain a vector of displacements between two vectors. Notice that we can use a similar input to calculate the
89 : differences between the instantaneous value of a pair of torsions and 3 reference values as shown below:
90 :
91 : ```plumed
92 : ref_psi: CONSTANT VALUES=2.25,1.3,-1.5
93 : ref_phi: CONSTANT VALUES=-1.91,-0.6,2.4
94 :
95 : psi: TORSION ATOMS=1,2,3,4
96 : phi: TORSION ATOMS=13,14,15,16
97 :
98 : dd: DISPLACEMENT ARG1=psi,phi ARG2=ref_psi,ref_phi
99 : PRINT ARG=dd FILE=colvar
100 : ```
101 :
102 : The output here will again be a $3\times 2$ matrix with each of the three rows holding a vector of displacements
103 : between the 2 instananeous values and one of the three sets of reference values.
104 :
105 : Lastly, we can use two sets of vectors in the input to DISPLACEMENT as shown below:
106 :
107 : ```plumed
108 : ref_psi: CONSTANT VALUES=2.25,1.3,-1.5
109 : ref_phi: CONSTANT VALUES=-1.91,-0.6,2.4
110 :
111 : psi: TORSION ATOMS1=1,2,3,4 ATOMS2=5,6,7,8 ATOMS3=9,10,11,12
112 : phi: TORSION ATOMS1=13,14,15,16 ATOMS2=17,18,19,20 ATOMS3=21,22,23,24
113 :
114 : dd: DISPLACEMENT ARG1=psi,phi ARG2=ref_psi,ref_phi
115 : PRINT ARG=dd FILE=colvar
116 : ```
117 :
118 : The output here is still a $3 \times 2$ matrix. Now, however, each of the three instantaneous angles we have calculated
119 : has its own set of reference values. A different pair of instaneous and reference values is used to calculate each element
120 : of the resulting matrix.
121 :
122 : DISPLACEMENT actions that compute $M\times N$ matrices, $D$, are used within the [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md),
123 : [NORMALIZED_EUCLIDEAN_DISTANCE](NORMALIZED_EUCLIDEAN_DISTANCE.md) and [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md) shortcuts.
124 : Doing so is useful as if you take the diagonal elements of a product of matrices that is similar to the product of vectors and matrices that we introduced
125 : earlier:
126 :
127 : $$
128 : d = D M D^T
129 : $$
130 :
131 : you can calculate $M$ values for the [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md), [NORMALIZED_EUCLIDEAN_DISTANCE](NORMALIZED_EUCLIDEAN_DISTANCE.md)
132 : and [MAHALANOBIS_DISTANCE](MAHALANOBIS_DISTANCE.md).
133 :
134 : */
135 : //+ENDPLUMEDOC
136 :
137 : namespace PLMD {
138 : namespace refdist {
139 :
140 : class Displacement : public ActionShortcut {
141 : public:
142 : static std::string fixArgumentDot( const std::string& argin );
143 : static void registerKeywords( Keywords& keys );
144 : Value* getValueWithLabel( const std::string& name ) const ;
145 : explicit Displacement(const ActionOptions&ao);
146 : };
147 :
148 : PLUMED_REGISTER_ACTION(Displacement,"DISPLACEMENT")
149 :
150 104 : void Displacement::registerKeywords( Keywords& keys ) {
151 104 : ActionShortcut::registerKeywords(keys);
152 104 : keys.add("compulsory","ARG1","The point that we are calculating the distance from");
153 104 : keys.add("compulsory","ARG2","The point that we are calculating the distance to");
154 208 : keys.setValueDescription("vector/matrix","the differences between the input arguments");
155 104 : keys.needsAction("DIFFERENCE");
156 104 : keys.needsAction("TRANSPOSE");
157 104 : keys.needsAction("VSTACK");
158 104 : }
159 :
160 354 : std::string Displacement::fixArgumentDot( const std::string& argin ) {
161 354 : std::string argout = argin;
162 354 : std::size_t dot=argin.find(".");
163 354 : if( dot!=std::string::npos ) {
164 16 : argout = argin.substr(0,dot) + "_" + argin.substr(dot+1);
165 : }
166 354 : return argout;
167 : }
168 :
169 52 : Displacement::Displacement( const ActionOptions& ao):
170 : Action(ao),
171 52 : ActionShortcut(ao) {
172 : // Read in argument names
173 : std::vector<std::string> arg1f, arg2f;
174 52 : parseVector("ARG1",arg1f);
175 104 : parseVector("ARG2",arg2f);
176 : // Check if one of the input arguments is a reference cluster
177 52 : if( arg1f.size()!=arg2f.size() ) {
178 0 : error("number of arguments specified to ARG1 should be same as number for ARG2");
179 : }
180 :
181 52 : Value* val1=getValueWithLabel( arg1f[0] );
182 52 : if( arg1f.size()==1 && val1->getRank()!=0 ) {
183 7 : Value* val2=getValueWithLabel( arg2f[0] );
184 7 : if( val1->getNumberOfValues()==val2->getNumberOfValues() ) {
185 10 : readInputLine( getShortcutLabel() + "_" + fixArgumentDot(arg1f[0]) + "_diff: DIFFERENCE ARG=" + arg1f[0] + "," + arg2f[0] );
186 10 : readInputLine( getShortcutLabel() + ": TRANSPOSE ARG=" + getShortcutLabel() + "_" + fixArgumentDot(arg1f[0]) + "_diff");
187 : } else {
188 4 : readInputLine( getShortcutLabel() + ": DIFFERENCE ARG=" + arg1f[0] + "," + arg2f[0] );
189 : }
190 : } else {
191 217 : for(unsigned i=0; i<arg1f.size(); ++i) {
192 344 : readInputLine( getShortcutLabel() + "_" + fixArgumentDot(arg1f[i]) + "_diff: DIFFERENCE ARG=" + arg1f[i] + "," + arg2f[i] );
193 : }
194 90 : std::string argdat = "ARG=" + getShortcutLabel() + "_" + fixArgumentDot(arg1f[0]) + "_diff";
195 172 : for(unsigned i=1; i<arg1f.size(); ++i) {
196 254 : argdat += "," + getShortcutLabel() + "_" + fixArgumentDot(arg1f[i]) + "_diff";
197 : }
198 90 : readInputLine( getShortcutLabel() + ": VSTACK " + argdat );
199 : }
200 52 : }
201 :
202 59 : Value* Displacement::getValueWithLabel( const std::string& name ) const {
203 59 : std::size_t dot=name.find(".");
204 59 : std::string sname = name.substr(0,dot);
205 59 : ActionWithValue* vv=plumed.getActionSet().selectWithLabel<ActionWithValue*>( sname );
206 59 : if( !vv ) {
207 0 : error("cannot find value with name " + name );
208 : }
209 59 : if( dot==std::string::npos ) {
210 57 : return vv->copyOutput(0);
211 : }
212 2 : if( !vv->exists(name) ) {
213 0 : error("cannot find value with name " + name );
214 : }
215 2 : return vv->copyOutput( name );
216 : }
217 :
218 : }
219 : }
|
