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 : ## Calculating multiple distances
76 :
77 : Suppose that we now have $m$ reference configurations we can define the following $m$ distances
78 : from these reference configurations:
79 :
80 : $$
81 : d_j^2 = (u-v_j)^T a \odot (u-v_j)
82 : $$
83 :
84 : Lets suppose that we put the $m$, $n$-dimensional $(u-v_j)$ vectors in this expression into a
85 : $n\times m$ matrix, $A$, by using the [DISPLACEMENT](DISPLACEMENT.md) command. It is then
86 : straightforward to show that the $d_j^2$ values in the above expression are the diagonal
87 : elements of the matrix product $A^T K \odot A$, where $K$ is an $n \times $m$ matrix that contains
88 : $m$ copies of the inverse covariance matrix $a$ in its columns.
89 :
90 : We can use this idea to calculate multiple NORMALIZED_EUCLIDEAN_DISTANCE values in the following inputs.
91 : This first example calculates the three distances between the instaneoues values of two torsions
92 : and three reference configurations.
93 :
94 : ```plumed
95 : m: CONSTANT VALUES=0.1,0.2
96 : ref_psi: CONSTANT VALUES=2.25,1.3,-1.5
97 : ref_phi: CONSTANT VALUES=-1.91,-0.6,2.4
98 :
99 : psi: TORSION ATOMS=1,2,3,4
100 : phi: TORSION ATOMS=13,14,15,16
101 :
102 : dd: NORMALIZED_EUCLIDEAN_DISTANCE ARG1=psi,phi ARG2=ref_psi,ref_phi METRIC=m
103 : PRINT ARG=dd FILE=colvar
104 : ```
105 :
106 : This section example calculates the three distances between a single reference value for the two
107 : torsions and three instances of this pair of torsions.
108 :
109 : ```plumed
110 : m: CONSTANT VALUES=0.1,0.2
111 : ref_psi: CONSTANT VALUES=2.25
112 : ref_phi: CONSTANT VALUES=-1.91
113 :
114 : psi: TORSION ATOMS1=1,2,3,4 ATOMS2=5,6,7,8 ATOMS3=9,10,11,12
115 : phi: TORSION ATOMS1=13,14,15,16 ATOMS2=17,18,19,20 ATOMS3=21,22,23,24
116 :
117 : dd: NORMALIZED_EUCLIDEAN_DISTANCE ARG1=psi,phi ARG2=ref_psi,ref_phi METRIC=m
118 : PRINT ARG=dd FILE=colvar
119 : ```
120 :
121 : This final example then computes three distances between three pairs of torsional angles and threee
122 : reference values for these three values.
123 :
124 : ```plumed
125 : m: CONSTANT VALUES=0.1,0.2
126 : ref_psi: CONSTANT VALUES=2.25,1.3,-1.5
127 : ref_phi: CONSTANT VALUES=-1.91,-0.6,2.4
128 :
129 : psi: TORSION ATOMS1=1,2,3,4 ATOMS2=5,6,7,8 ATOMS3=9,10,11,12
130 : phi: TORSION ATOMS1=13,14,15,16 ATOMS2=17,18,19,20 ATOMS3=21,22,23,24
131 :
132 : dd: NORMALIZED_EUCLIDEAN_DISTANCE ARG1=psi,phi ARG2=ref_psi,ref_phi METRIC=m
133 : PRINT ARG=dd FILE=colvar
134 : ```
135 :
136 : */
137 : //+ENDPLUMEDOC
138 :
139 : namespace PLMD {
140 : namespace refdist {
141 :
142 : class NormalizedEuclideanDistance : public ActionShortcut {
143 : public:
144 : static void registerKeywords( Keywords& keys );
145 : explicit NormalizedEuclideanDistance(const ActionOptions&ao);
146 : };
147 :
148 : PLUMED_REGISTER_ACTION(NormalizedEuclideanDistance,"NORMALIZED_EUCLIDEAN_DISTANCE")
149 :
150 13 : void NormalizedEuclideanDistance::registerKeywords( Keywords& keys ) {
151 13 : ActionShortcut::registerKeywords(keys);
152 13 : keys.add("compulsory","ARG1","The poin that we are calculating the distance from");
153 13 : keys.add("compulsory","ARG2","The point that we are calculating the distance to");
154 13 : keys.add("compulsory","METRIC","The inverse covariance matrix that should be used when calculating the distance");
155 13 : keys.addFlag("SQUARED",false,"The squared distance should be calculated");
156 26 : keys.setValueDescription("scalar/vector","the normalized euclidean distances between the input vectors");
157 13 : keys.needsAction("DISPLACEMENT");
158 13 : keys.needsAction("CUSTOM");
159 13 : keys.needsAction("OUTER_PRODUCT");
160 13 : keys.needsAction("TRANSPOSE");
161 13 : keys.needsAction("MATRIX_PRODUCT_DIAGONAL");
162 13 : keys.needsAction("ONES");
163 13 : }
164 :
165 8 : NormalizedEuclideanDistance::NormalizedEuclideanDistance( const ActionOptions& ao):
166 : Action(ao),
167 8 : ActionShortcut(ao) {
168 : std::string arg1, arg2, metstr;
169 8 : parse("ARG1",arg1);
170 8 : parse("ARG2",arg2);
171 8 : parse("METRIC",metstr);
172 : // Vectors are in rows here
173 16 : readInputLine( getShortcutLabel() + "_diff: DISPLACEMENT ARG1=" + arg1 + " ARG2=" + arg2 );
174 : // Vectors are in columns here
175 16 : readInputLine( getShortcutLabel() + "_diffT: TRANSPOSE ARG=" + getShortcutLabel() + "_diff");
176 : // Get the action that computes the differences
177 8 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( getShortcutLabel() + "_diffT");
178 8 : plumed_assert( av );
179 : // If this is a matrix we need create a matrix to multiply by
180 8 : if( av->copyOutput(0)->getRank()==2 ) {
181 : // Create some ones
182 : std::string nones;
183 4 : Tools::convert( av->copyOutput(0)->getShape()[1], nones );
184 8 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + nones);
185 : // Now do some multiplication to create a matrix that can be multiplied by our "inverse variance" vector
186 4 : if( av->copyOutput(0)->getShape()[0]==1 ) {
187 4 : readInputLine( getShortcutLabel() + "_" + metstr + "T: CUSTOM ARG=" + metstr + "," + getShortcutLabel() + "_ones FUNC=x*y PERIODIC=NO");
188 4 : readInputLine( getShortcutLabel() + "_" + metstr + ": TRANSPOSE ARG=" + getShortcutLabel() + "_" + metstr + "T");
189 : } else {
190 4 : readInputLine( getShortcutLabel() + "_" + metstr + ": OUTER_PRODUCT ARG=" + metstr + "," + getShortcutLabel() + "_ones");
191 : }
192 8 : metstr = getShortcutLabel() + "_" + metstr;
193 : }
194 : // Now do the multiplication
195 16 : readInputLine( getShortcutLabel() + "_sdiff: CUSTOM ARG=" + metstr + "," + getShortcutLabel() +"_diffT FUNC=x*y PERIODIC=NO");
196 : bool squared;
197 8 : parseFlag("SQUARED",squared);
198 8 : std::string olab = getShortcutLabel();
199 8 : if( !squared ) {
200 : olab += "_2";
201 : }
202 16 : readInputLine( olab + ": MATRIX_PRODUCT_DIAGONAL ARG=" + getShortcutLabel() +"_diff," + getShortcutLabel() + "_sdiff");
203 8 : if( !squared ) {
204 10 : readInputLine( getShortcutLabel() + ": CUSTOM ARG=" + getShortcutLabel() + "_2 FUNC=sqrt(x) PERIODIC=NO");
205 : }
206 8 : }
207 :
208 : }
209 : }
|
