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