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