Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016,2017 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/ActionWithVector.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/ActionSet.h"
25 : #include "PathProjectionCalculator.h"
26 :
27 : //+PLUMEDOC COLVAR GEOMETRIC_PATH
28 : /*
29 : Distance along and from a path calculated using geometric formulas
30 :
31 : The Path Collective Variables developed by Branduardi and that are described in the first paper that is cited below allow one
32 : to compute the progress along a high-dimensional path and the distance from the high-dimensional path. The method introduced in that
33 : paper is implemented in the shortcut [PATH](PATH.md). This action provides an implementation of the alternative method for calculating
34 : the position on and distance from the path that was proposed by Diaz Leines and Ensing in the second paper that is cited below. In their
35 : method, the progress along the path $s$ is calculated using:
36 :
37 : $$
38 : s = i_2 + \textrm{sign}(i_2-i_1) \frac{ \sqrt{( \mathbf{v}_1\cdot\mathbf{v}_2 )^2 - |\mathbf{v}_3|^2(|\mathbf{v}_1|^2 - |\mathbf{v}_2|^2) } }{2|\mathbf{v}_3|^2} - \frac{\mathbf{v}_1\cdot\mathbf{v}_3 - |\mathbf{v}_3|^2}{2|\mathbf{v}_3|^2}
39 : $$
40 :
41 : where $\mathbf{v}_1$ and $\mathbf{v}_3$ are the vectors connecting the current position to the closest and second closest node of the path,
42 : respectfully and $i_1$ and $i_2$ are the projections of the closest and second closest frames of the path. $\mathbf{v}_2$, meanwhile, is the
43 : vector connecting the closest frame to the second closest frame. The distance from the path, $z$ is calculated using:
44 :
45 : $$
46 : z = \sqrt{ \left[ |\mathbf{v}_1|^2 - |\mathbf{v}_2| \left( \frac{ \sqrt{( \mathbf{v}_1\cdot\mathbf{v}_2 )^2 - |\mathbf{v}_3|^2(|\mathbf{v}_1|^2 - |\mathbf{v}_2|^2) } }{2|\mathbf{v}_3|^2} - \frac{\mathbf{v}_1\cdot\mathbf{v}_3 - |\mathbf{v}_3|^2}{2|\mathbf{v}_3|^2} \right) \right]^2 }
47 : $$
48 :
49 : The symbols here are as they were for $s$.
50 :
51 : Please note that the shortcut [GPATH](GPATH.md) provides an interface to this action that is less flexible than using this action directly but that is easier to use.
52 :
53 : ## Examples
54 :
55 : The example input below shows how to use this action
56 :
57 : ```plumed
58 : #SETTINGS INPUTFILES=regtest/mapping/rt-adapt/mypath.pdb
59 : # Calculate the instantaneous position in the high dimensional space
60 : d1: DISTANCE ATOMS=1,2 COMPONENTS
61 : # Read in the refence coordinates from the path
62 : ref_d1x: PDB2CONSTANT REFERENCE=regtest/mapping/rt-adapt/mypath.pdb ARG=d1.x
63 : ref_d1y: PDB2CONSTANT REFERENCE=regtest/mapping/rt-adapt/mypath.pdb ARG=d1.y
64 : path: CONSTANT VALUES=-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5
65 : # Now calculate the differences between our instantaneous position and the milestones on the high dimensional space
66 : d1_dataP: DISPLACEMENT ARG2=d1.x,d1.y ARG1=ref_d1x,ref_d1y
67 : # We need the vectors pointing from the reference positions to the instantaneous positions here which is
68 : # minus the vector calculated above. Remember we can only specify the vectors in ARG1 for the DISPLACEMENT command
69 : d1_data: CUSTOM ARG=d1_dataP FUNC=-x PERIODIC=NO
70 : # And finally the path CV
71 : pp: GEOMETRIC_PATH ARG=d1_data METRIC={DIFFERENCE} REFERENCE=ref_d1x,ref_d1y PROPERTY=path
72 : PRINT ARG=d1.x,d1.y,pp.* FILE=colvar
73 : ```
74 :
75 : The curved path here is defined using a series of points in a two dimensional space. The input above includes actions ([PDB2CONSTANT](PDB2CONSTANT.md)) that set up
76 : constant values that hold the values of the arguments for each of the milestones on the path. A [DISPLACEMENT](DISPLACEMENT.md) shortcut is then used to calculate the vectors
77 : connecting each of the path milestones to the instantenous coordinates of the system in the high dimensional space. Notice, how the METRIC keyword is used in the input
78 : to GEOMETRIC_PATH. This keyword is necessary as the vectors connecting adjacent milestones on the path appear in the equations for $s$ and $z$ above. The METRIC keyword tells
79 : this action how distances between frames should be calculated. There are essentially two options that we have tested here. You can use [DIFFERENCE](DIFFERENCE.md) as has been done above. In that
80 : case the assumption is the path milestones are vectors of coordinates and the difference $\mathbf{v}$ between two of them, $\mathbf{m}_1$ and $\mathbf{m}_2$, can be calculated as:
81 :
82 : $$
83 : \mathbf{v} = \mathbf{m}_1 - \mathbf{m}_2
84 : $$
85 :
86 : with periodic variables treated appropriately. Alternatively, the path milestones might be the positions of atoms and you may want to calculate differences between
87 : them using an [RMSD](RMSD.md) action as has been done in the example below.
88 :
89 : ```plumed
90 : ##SETTINGS INPUTFILES=regtest/trajectories/path_msd/all.pdb
91 : # Positions for all the frames in the path
92 : path: CONSTANT VALUES=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42
93 : # Calculate the vectors connecting the instantaneous position to the path milesones
94 : vecs: RMSD DISPLACEMENT SQUARED REFERENCE=regtest/trajectories/path_msd/all.pdb TYPE=OPTIMAL
95 : p1b: GEOMETRIC_PATH ARG=vecs.disp REFERENCE=vecs_ref METRIC={RMSD DISPLACEMENT TYPE=OPTIMAL ALIGN=1,1,1,1,1,1,1,1,1,1,1,1,1 DISPLACE=1,1,1,1,1,1,1,1,1,1,1,1,1} PROPERTY=path METRIC_COMPONENT=disp
96 : PRINT ARG=p1b.* FILE=colvar_b STRIDE=1
97 : ```
98 :
99 : Here an [RMSD](RMSD.md) action with label `vecs` is used to calculate the vector connecting the instantaneous position and each of the path milestones. These vectors thus describe
100 : how far each atom has been displaced in going from one structure to the other in a way that neglects translation of the center of mass and rotation of the reference frame. When we
101 : compute the formula aboves, which include the vectors that connect various milestone configurations on our path we must compute the vector connecting the milestones using the same
102 : approach. The input associated with the METRIC keyword is thus the same input as was used in the RMSD action that was used to calculate the vectors connecting the milestones to the
103 : instantaneous configuration. Furthermore, we must also use METRIC_COMPONENT here to specify that it is the `disp` component of the [RMSD](RMSD.md) action that contains the vector of
104 : interest.
105 :
106 : It is worth noting that the METRIC keyword in the inputs above can accept inputs for one or multiple PLUMED actions. There is thus flexibility to implement more complicated variants on the geometric path
107 : implemented here by using the input for different actions as input to the METRIC keyword. If you are interested in using these features you can look at the code within `src/mapping/PathProjectionCalculator.h
108 : and `src/mapping/PathProjectionCalculator.cpp` to see how this method called a separate PLUMED instance to calculate the vectors connecting the various milesones in the path.
109 :
110 :
111 : */
112 : //+ENDPLUMEDOC
113 :
114 : namespace PLMD {
115 : namespace mapping {
116 :
117 : class GeometricPath : public ActionWithVector {
118 : private:
119 : PathProjectionCalculator path_projector;
120 : public:
121 : static void registerKeywords(Keywords& keys);
122 : explicit GeometricPath(const ActionOptions&);
123 : unsigned getNumberOfDerivatives() override ;
124 : void calculate() override ;
125 : };
126 :
127 : PLUMED_REGISTER_ACTION(GeometricPath,"GEOMETRIC_PATH")
128 :
129 22 : void GeometricPath::registerKeywords(Keywords& keys) {
130 22 : ActionWithVector::registerKeywords(keys);
131 22 : PathProjectionCalculator::registerKeywords(keys);
132 44 : keys.addInputKeyword("compulsory","PROPERTY","vector","the label of a value that contains the coordinates we are projecting these points onto");
133 44 : keys.addOutputComponent("s","default","scalar","the position on the path");
134 44 : keys.addOutputComponent("z","default","scalar","the distance from the path");
135 22 : keys.addDOI("10.1063/1.2432340");
136 22 : keys.addDOI("10.1103/PhysRevLett.109.020601");
137 22 : keys.needsAction("GEOMETRIC_PATH");
138 22 : keys.needsAction("PDB2CONSTANT");
139 22 : }
140 :
141 9 : GeometricPath::GeometricPath(const ActionOptions&ao):
142 : Action(ao),
143 : ActionWithVector(ao),
144 9 : path_projector(this) {
145 : // Get the coordinates in the low dimensional space
146 9 : std::vector<std::string> pcoord;
147 9 : parseVector("PROPERTY", pcoord );
148 9 : std::vector<Value*> theprop;
149 9 : ActionWithArguments::interpretArgumentList( pcoord, plumed.getActionSet(), this, theprop );
150 9 : if( theprop.size()!=1 ) {
151 0 : error("did not find property to project on");
152 : }
153 9 : if( theprop[0]->getNumberOfValues()!=getPntrToArgument(0)->getShape()[0] ) {
154 0 : error("mismatch between number of frames and property of interest");
155 : }
156 9 : log.printf(" projecting onto : %s \n", theprop[0]->getName().c_str() );
157 9 : std::vector<Value*> args( getArguments() );
158 9 : args.push_back( theprop[0] );
159 9 : requestArguments( args );
160 : // Create the values to store the output
161 18 : addComponentWithDerivatives("s");
162 18 : componentIsNotPeriodic("s");
163 18 : addComponentWithDerivatives("z");
164 18 : componentIsNotPeriodic("z");
165 9 : }
166 :
167 4413 : unsigned GeometricPath::getNumberOfDerivatives() {
168 4413 : return getPntrToArgument(0)->getShape()[0]*getPntrToArgument(0)->getShape()[1] + getPntrToArgument(1)->getShape()[0];
169 : }
170 :
171 4395 : void GeometricPath::calculate() {
172 4395 : unsigned k=0, iclose1=0, iclose2=0;
173 : double v1v1=0, v3v3=0;
174 4395 : unsigned nrows = getPntrToArgument(0)->getShape()[0];
175 4395 : unsigned ncols = getPntrToArgument(0)->getShape()[1];
176 160555 : for(unsigned i=0; i<nrows; ++i) {
177 : double dist = 0;
178 4959540 : for(unsigned j=0; j<ncols; ++j) {
179 4803380 : double tmp = getPntrToArgument(0)->get(k);
180 4803380 : dist += tmp*tmp;
181 4803380 : k++;
182 : }
183 156160 : if( i==0 ) {
184 : v1v1 = dist;
185 4395 : iclose1 = 0;
186 151765 : } else if( dist<v1v1 ) {
187 : v3v3=v1v1;
188 : v1v1=dist;
189 60282 : iclose2=iclose1;
190 60282 : iclose1=i;
191 91483 : } else if( i==1 ) {
192 : v3v3=dist;
193 489 : iclose2=1;
194 90994 : } else if( dist<v3v3 ) {
195 : v3v3=dist;
196 2848 : iclose2=i;
197 : }
198 : }
199 : // And find third closest point
200 4395 : int isign = iclose1 - iclose2;
201 : if( isign>1 ) {
202 : isign=1;
203 4395 : } else if( isign<-1 ) {
204 : isign=-1;
205 : }
206 4395 : int iclose3 = iclose1 + isign;
207 4395 : unsigned ifrom=iclose1, ito=iclose3;
208 4395 : if( iclose3<0 || static_cast<unsigned>(iclose3)>=nrows ) {
209 345 : ifrom=iclose2;
210 345 : ito=iclose1;
211 : }
212 :
213 : // And calculate projection of vector connecting current point to closest frame on vector connecting nearest two frames
214 4395 : std::vector<double> displace;
215 4395 : path_projector.getDisplaceVector( ifrom, ito, displace );
216 : double v2v2=0, v1v2=0;
217 4395 : k=ncols*iclose1;
218 151861 : for(unsigned i=0; i<displace.size(); ++i) {
219 147466 : v2v2 += displace[i]*displace[i];
220 147466 : v1v2 += displace[i]*getPntrToArgument(0)->get(k+i);
221 : }
222 :
223 : // This computes s value
224 4395 : double spacing = getPntrToArgument(1)->get(iclose1) - getPntrToArgument(1)->get(iclose2);
225 4395 : double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
226 4395 : double dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
227 4395 : double path_s = getPntrToArgument(1)->get(iclose1) + spacing * dx;
228 4395 : Value* sp = getPntrToComponent(0);
229 : sp->set( path_s );
230 4395 : if( !doNotCalculateDerivatives() ) {
231 45318 : for(unsigned i=0; i<ncols; ++i) {
232 43680 : sp->addDerivative( ncols*iclose1 + i, 0.5*spacing*(v1v2*displace[i]/v2v2 - getPntrToArgument(0)->get(ncols*iclose1 + i))/root + 0.5*spacing*displace[i]/v2v2 );
233 43680 : sp->addDerivative( ncols*iclose2 + i, 0.5*spacing*getPntrToArgument(0)->get(ncols*iclose2 + i)/root );
234 : }
235 : }
236 :
237 : // This computes z value
238 4395 : path_projector.getDisplaceVector( iclose2, iclose1, displace );
239 : double v4v4=0, proj=0;
240 4395 : k=ncols*iclose1;
241 151861 : for(unsigned i=0; i<displace.size(); ++i) {
242 147466 : v4v4 += displace[i]*displace[i];
243 147466 : proj += displace[i]*getPntrToArgument(0)->get(k+i);
244 : }
245 4395 : double path_z = v1v1 + dx*dx*v4v4 - 2*dx*proj;
246 4395 : path_z = sqrt(path_z);
247 4395 : Value* zp = getPntrToComponent(1);
248 : zp->set( path_z );
249 4395 : if( !doNotCalculateDerivatives() ) {
250 45318 : for(unsigned i=0; i<ncols; ++i) {
251 43680 : zp->addDerivative( ncols*iclose1 + i, (1/path_z)*(getPntrToArgument(0)->get(ncols*iclose1 + i) +
252 43680 : (v4v4*dx-proj)*sp->getDerivative(ncols*iclose1 + i)/spacing -
253 43680 : dx*displace[i]) );
254 43680 : zp->addDerivative( ncols*iclose2 + i, (v4v4*dx-proj)*sp->getDerivative(ncols*iclose2 + i)/(path_z*spacing) );
255 : }
256 : }
257 4395 : }
258 :
259 : }
260 : }
|
