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/ActionShortcut.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/ActionWithValue.h"
27 :
28 : //+PLUMEDOC MCOLVAR EUCLIDEAN_DISTANCE
29 : /*
30 : Calculate the euclidean distance between two vectors of arguments
31 :
32 : If we have two $n$-dimensional vectors $u$ and $v$ we can calculate the
33 : [Euclidean distance](https://en.wikipedia.org/wiki/Euclidean_distance) between the two points as
34 :
35 : $$
36 : d = \sqrt{ \sum_{i=1}^n (u_i - v_i)^2 }
37 : $$
38 :
39 : which can be expressed in matrix form as:
40 :
41 : $$
42 : d^2 = (u-v)^T (u-v)
43 : $$
44 :
45 : The input below shows an example where this is used to calculate the Euclidean distance
46 : between the instaneous values of some torsional angles and some reference values
47 : for these torsion. In this first example the input values are vectors:
48 :
49 : ```plumed
50 : c: CONSTANT VALUES=1,2,3
51 : d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
52 : dd: EUCLIDEAN_DISTANCE ARG1=c ARG2=d
53 : PRINT ARG=dd FILE=colvar
54 : ```
55 :
56 : while this second example does the same thing but uses scalars in input.
57 :
58 : ```plumed
59 : c1: CONSTANT VALUE=1
60 : d1: DISTANCE ATOMS=1,2
61 : c2: CONSTANT VALUE=2
62 : d2: DISTANCE ATOMS=3,4
63 : c3: CONSTANT VALUE=3
64 : d3: DISTANCE ATOMS=5,6
65 : dd: EUCLIDEAN_DISTANCE ARG1=c1,c2,c3 ARG2=d1,d2,d3
66 : PRINT ARG=dd FILE=colvar
67 : ```
68 :
69 : Lastly, note that if you want to calculate the square of the distance rather than the distance you can use
70 : the `SQUARED` flag as shown below:
71 :
72 : ```plumed
73 : c: CONSTANT VALUES=1,2,3
74 : d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
75 : dd: EUCLIDEAN_DISTANCE ARG1=c ARG2=d SQUARED
76 : PRINT ARG=dd FILE=colvar
77 : ```
78 :
79 : Calculating the square of the distance is slightly cheapter than computing the distance as you avoid taking the square root.
80 :
81 : ## Calculating multiple distances
82 :
83 : Suppose that we now have $m$ reference configurations we can define the following $m$ distances
84 : from these reference configurations:
85 :
86 : $$
87 : d_j^2 = (u-v_j)^T (u-v_j)
88 : $$
89 :
90 : Lets suppose that we put the $m$, $n$-dimensional $(u-v_j)$ vectors in this expression into a
91 : $n\times m$ matrix, $A$, by using the [DISPLACEMENT](DISPLACEMENT.md) command. It is then
92 : straightforward to show that the $d_j^2$ values in the above expression are then the diagonal
93 : elements of the matrix product $A^T A$.
94 :
95 : We can use this idea to calculate multiple EUCLIDEAN_DISTANCE values in the following inputs.
96 : This first example calculates the three distances between the instaneoues values of two torsions
97 : and three reference configurations.
98 :
99 : ```plumed
100 : ref_psi: CONSTANT VALUES=2.25,1.3,-1.5
101 : ref_phi: CONSTANT VALUES=-1.91,-0.6,2.4
102 :
103 : psi: TORSION ATOMS=1,2,3,4
104 : phi: TORSION ATOMS=13,14,15,16
105 :
106 : dd: EUCLIDEAN_DISTANCE ARG2=psi,phi ARG1=ref_psi,ref_phi
107 : PRINT ARG=dd FILE=colvar
108 : ```
109 :
110 : This section example calculates the three distances between a single reference value for the two
111 : torsions and three instances of this pair of torsions.
112 :
113 : ```plumed
114 : ref_psi: CONSTANT VALUES=2.25
115 : ref_phi: CONSTANT VALUES=-1.91
116 :
117 : psi: TORSION ATOMS1=1,2,3,4 ATOMS2=5,6,7,8 ATOMS3=9,10,11,12
118 : phi: TORSION ATOMS1=13,14,15,16 ATOMS2=17,18,19,20 ATOMS3=21,22,23,24
119 :
120 : dd: EUCLIDEAN_DISTANCE ARG1=psi,phi ARG2=ref_psi,ref_phi
121 : PRINT ARG=dd FILE=colvar
122 : ```
123 :
124 : This final example then computes three distances between three pairs of torsional angles and threee
125 : reference values for these three values.
126 :
127 : ```plumed
128 : ref_psi: CONSTANT VALUES=2.25,1.3,-1.5
129 : ref_phi: CONSTANT VALUES=-1.91,-0.6,2.4
130 :
131 : psi: TORSION ATOMS1=1,2,3,4 ATOMS2=5,6,7,8 ATOMS3=9,10,11,12
132 : phi: TORSION ATOMS1=13,14,15,16 ATOMS2=17,18,19,20 ATOMS3=21,22,23,24
133 :
134 : dd: EUCLIDEAN_DISTANCE ARG1=psi,phi ARG2=ref_psi,ref_phi
135 : PRINT ARG=dd FILE=colvar
136 : ```
137 :
138 : !!! note "scalars must be specified in ARG2"
139 :
140 : If you use a mixture of vectors are scalars when specifying the input to to this action the
141 : vectors should be passed using the ARG1 keyword and the scalars must be passed in the ARG2 keyword
142 : as is done in the example inputs above.
143 :
144 : */
145 : //+ENDPLUMEDOC
146 :
147 : namespace PLMD {
148 : namespace refdist {
149 :
150 : class EuclideanDistance : public ActionShortcut {
151 : public:
152 : static void registerKeywords( Keywords& keys );
153 : explicit EuclideanDistance(const ActionOptions&ao);
154 : };
155 :
156 : PLUMED_REGISTER_ACTION(EuclideanDistance,"EUCLIDEAN_DISTANCE")
157 :
158 35 : void EuclideanDistance::registerKeywords( Keywords& keys ) {
159 35 : ActionShortcut::registerKeywords(keys);
160 35 : keys.add("compulsory","ARG1","The poin that we are calculating the distance from");
161 35 : keys.add("compulsory","ARG2","The point that we are calculating the distance to");
162 35 : keys.addFlag("SQUARED",false,"The squared distance should be calculated");
163 70 : keys.setValueDescription("scalar/vector","the euclidean distances between the input vectors");
164 35 : keys.needsAction("DISPLACEMENT");
165 35 : keys.needsAction("CUSTOM");
166 35 : keys.needsAction("TRANSPOSE");
167 35 : keys.needsAction("MATRIX_PRODUCT_DIAGONAL");
168 35 : }
169 :
170 28 : EuclideanDistance::EuclideanDistance( const ActionOptions& ao):
171 : Action(ao),
172 28 : ActionShortcut(ao) {
173 : std::string arg1, arg2;
174 28 : parse("ARG1",arg1);
175 28 : parse("ARG2",arg2);
176 : // Vectors are in rows here
177 56 : readInputLine( getShortcutLabel() + "_diff: DISPLACEMENT ARG1=" + arg1 + " ARG2=" + arg2 );
178 : // Get the action that computes the differences
179 28 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_diff");
180 28 : plumed_assert( av );
181 : // Check if squared
182 : bool squared;
183 28 : parseFlag("SQUARED",squared);
184 28 : std::string olab = getShortcutLabel();
185 28 : if( !squared ) {
186 : olab += "_2";
187 : }
188 : // Deal with an annoying corner case when displacement has a single argument
189 28 : if( av->copyOutput(0)->getRank()==0 ) {
190 0 : readInputLine( olab + ": CUSTOM ARG=" + getShortcutLabel() + "_diff FUNC=x*x PERIODIC=NO");
191 : } else {
192 : // Notice that the vectors are in the columns here
193 56 : readInputLine( getShortcutLabel() + "_diffT: TRANSPOSE ARG=" + getShortcutLabel() + "_diff");
194 56 : readInputLine( olab + ": MATRIX_PRODUCT_DIAGONAL ARG=" + getShortcutLabel() + "_diff," + getShortcutLabel() + "_diffT");
195 : }
196 28 : if( !squared ) {
197 46 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
198 : }
199 28 : }
200 :
201 : }
202 : }
|
