Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020-2023 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 :
23 : #include "Function.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/IFile.h"
26 : #include <algorithm>
27 : #include <limits>
28 : #include <string>
29 :
30 : namespace PLMD {
31 : namespace function {
32 :
33 : //+PLUMEDOC FUNCTION FUNCPATHGENERAL
34 : /*
35 : This function calculates path collective variables (PCVs) using an arbitrary combination of collective variables.
36 :
37 : The method used to calculate the PCVs that is used in this method is described in the paper from the bibliography below.
38 :
39 : This variable computes the progress along a given set of frames that is provided in an input file ("s" component) and the distance from them ("z" component).
40 : The input file could be a colvar file generated with plumed driver on a trajectory containing the frames.
41 :
42 : The metric for the path collective variables takes the following form:
43 :
44 : $$
45 : R[X - X_i] = \sum_{j=1}^M c_j^2 (x_j - x_{i,j})^2\,.
46 : $$
47 :
48 : Here, the coefficients $c_j$ determine the relative weights of the collective variables $c_j$ in the metric.
49 : A value for the lambda coefficient also needs to be provided, typically chosen in such a way that it ensures a smooth variation of the "s" component.
50 :
51 : ## Examples
52 :
53 : This command calculates the PCVs using the values from the file COLVAR_TRAJ and the provided values for the lambda and the coefficients.
54 : Since the columns in the file were not specified, the first one will be ignored (assumed to correspond to the time) and the rest used.
55 :
56 : ```plumed
57 : FUNCPATHGENERAL ...
58 : LABEL=path
59 : LAMBDA=12.2
60 : REFERENCE=COLVAR_TRAJ
61 : COEFFICIENTS=0.3536,0.3536,0.3536,0.3536,0.7071
62 : ARG=d1,d2,d,t,drmsd
63 : ... FUNCPATHGENERAL
64 : ```
65 :
66 : The command below is a variation of the previous one, specifying a subset of the collective variables and using a neighbor list.
67 : The columns are zero-indexed.
68 : The neighbor list will include the 10 closest frames and will be recalculated every 20 steps.
69 :
70 : ```plumed
71 : FUNCPATHGENERAL ...
72 : LABEL=path
73 : LAMBDA=5.0
74 : REFERENCE=COLVAR_TRAJ
75 : COLUMNS=2,3,4
76 : COEFFICIENTS=0.3536,0.3536,0.3536
77 : ARG=d2,d,t
78 : NEIGH_SIZE=10
79 : NEIGH_STRIDE=20
80 : ... FUNCPATHGENERAL
81 : ```
82 :
83 : */
84 : //+ENDPLUMEDOC
85 :
86 : class FuncPathGeneral : public Function {
87 : double lambda;
88 : int neigh_size;
89 : double neigh_stride;
90 :
91 : std::vector<double> coefficients;
92 : std::vector< std::vector<double> > path_cv_values;
93 :
94 : // For faster calculation
95 : std::vector<double> expdists;
96 :
97 : // For calculating derivatives
98 : std::vector< std::vector<double> > numerators;
99 : std::vector<double> s_path_ders;
100 : std::vector<double> z_path_ders;
101 :
102 : // For handling periodicity
103 : std::vector<double> domains;
104 :
105 : std::string reference;
106 : std::vector<int> columns;
107 :
108 : std::vector< std::pair<int,double> > neighpair;
109 : std::vector <Value*> allArguments;
110 :
111 : // Methods
112 : void loadReference();
113 :
114 : struct pairordering {
115 : bool operator ()(std::pair<int, double> const& a, std::pair<int, double> const& b) {
116 0 : return (a).second < (b).second;
117 : }
118 : };
119 :
120 : public:
121 : explicit FuncPathGeneral(const ActionOptions&);
122 : // Active methods:
123 : virtual void calculate();
124 : virtual void prepare();
125 : static void registerKeywords(Keywords& keys);
126 : };
127 :
128 : PLUMED_REGISTER_ACTION(FuncPathGeneral, "FUNCPATHGENERAL")
129 :
130 1 : void FuncPathGeneral::loadReference() {
131 1 : IFile input;
132 1 : input.open(reference);
133 1 : if (!input) {
134 0 : plumed_merror("Could not open the reference file!");
135 : }
136 18 : while (input) {
137 17 : std::vector<std::string> strings;
138 17 : Tools::getParsedLine(input, strings);
139 17 : if (strings.empty()) {
140 : continue;
141 : }
142 16 : std::vector<double> colvarLine;
143 : double value;
144 16 : int max = columns.empty() ? strings.size() : columns.size();
145 128 : for (int i = 0; i < max; ++i) {
146 112 : int col = columns.empty() ? i : columns[i];
147 : // If no columns have been entered, ignore the first (time) and take the rest
148 112 : if (columns.empty() && i == 0) {
149 16 : continue;
150 : }
151 :
152 96 : if( col<0 || col>=static_cast<int>(strings.size()) ) {
153 0 : error("column " + std::to_string(col) + " not found in line of reference file " + reference);
154 : }
155 96 : Tools::convert(strings[col], value);
156 96 : colvarLine.push_back(value);
157 : }
158 16 : path_cv_values.push_back(colvarLine);
159 17 : }
160 1 : }
161 :
162 5 : void FuncPathGeneral::registerKeywords(Keywords& keys) {
163 5 : Function::registerKeywords(keys);
164 5 : keys.add("compulsory", "LAMBDA", "Lambda parameter required for smoothing");
165 5 : keys.add("compulsory", "COEFFICIENTS", "Coefficients to be assigned to the CVs");
166 5 : keys.add("compulsory", "REFERENCE", "Colvar file needed to provide the CV milestones");
167 5 : keys.add("optional", "COLUMNS", "List of columns in the reference colvar file specifying the CVs");
168 5 : keys.add("optional", "NEIGH_SIZE", "Size of the neighbor list");
169 5 : keys.add("optional", "NEIGH_STRIDE", "How often the neighbor list needs to be calculated in time units");
170 10 : keys.addOutputComponent("s", "default", "scalar","Position on the path");
171 10 : keys.addOutputComponent("z", "default", "scalar","Distance from the path");
172 5 : keys.addDOI("10.1021/acs.jctc.8b00563");
173 5 : }
174 :
175 1 : FuncPathGeneral::FuncPathGeneral(const ActionOptions&ao):
176 : Action(ao),
177 : Function(ao),
178 1 : neigh_size(-1),
179 1 : neigh_stride(-1.) {
180 1 : parse("LAMBDA", lambda);
181 1 : parse("NEIGH_SIZE", neigh_size);
182 1 : parse("NEIGH_STRIDE", neigh_stride);
183 1 : parse("REFERENCE", reference);
184 1 : parseVector("COEFFICIENTS", coefficients);
185 1 : parseVector("COLUMNS", columns);
186 1 : checkRead();
187 1 : log.printf(" lambda is %f\n", lambda);
188 1 : if (getNumberOfArguments() != coefficients.size()) {
189 0 : plumed_merror("The numbers of coefficients and CVs are different!");
190 : }
191 1 : if (!columns.empty()) {
192 0 : if (columns.size() != coefficients.size()) {
193 0 : plumed_merror("The numbers of coefficients and columns are different!");
194 : }
195 : }
196 1 : log.printf(" Consistency check completed! Your path cvs look good!\n");
197 1 : if( !(lambda>0) ) {
198 0 : error("LAMBDA should be positive");
199 : }
200 1 : if( neigh_size>0 && !(neigh_stride>0) ) {
201 0 : error("with NEIGH_SIZE you should also give a positive NEIGH_STRIDE");
202 : }
203 :
204 : // Load the reference colvar file
205 1 : loadReference();
206 1 : if( path_cv_values.empty() ) {
207 0 : error("no frames found in reference file " + reference);
208 : }
209 17 : for(const auto & frame : path_cv_values) {
210 16 : if( frame.size()<getNumberOfArguments() ) {
211 0 : error("each line of the reference file should contain at least one value for each argument");
212 : }
213 : }
214 :
215 : // Do some neighbour printout
216 1 : if (neigh_size > 0) {
217 0 : if (static_cast<unsigned>(neigh_size) > path_cv_values.size()) {
218 0 : log.printf(" List size required ( %d ) is too large: resizing to the maximum number of arg required: %d \n", neigh_size, getNumberOfArguments());
219 0 : neigh_size = path_cv_values.size();
220 : }
221 0 : log.printf(" Neighbour list enabled: \n");
222 0 : log.printf(" size : %d elements\n", neigh_size);
223 0 : log.printf(" stride : %f time \n", neigh_stride);
224 : } else {
225 1 : log.printf(" Neighbour list NOT enabled \n");
226 : }
227 :
228 2 : addComponentWithDerivatives("s");
229 2 : componentIsNotPeriodic("s");
230 2 : addComponentWithDerivatives("z");
231 2 : componentIsNotPeriodic("z");
232 :
233 : // Initialise vectors
234 1 : std::vector<double> temp (coefficients.size());
235 17 : for (unsigned i = 0; i < path_cv_values.size(); ++i) {
236 16 : numerators.push_back(temp);
237 16 : expdists.push_back(0.);
238 16 : s_path_ders.push_back(0.);
239 16 : z_path_ders.push_back(0.);
240 : }
241 :
242 : // Store the arguments
243 7 : for (unsigned i=0; i<getNumberOfArguments(); i++) {
244 6 : allArguments.push_back(getPntrToArgument(i));
245 : }
246 :
247 : // Get periodic domains, negative for not periodic, stores half the domain length (maximum difference)
248 7 : for (unsigned i = 0; i < allArguments.size(); ++i) {
249 6 : if (allArguments[i]->isPeriodic()) {
250 : double min_lim, max_lim;
251 0 : allArguments[i]->getDomain(min_lim, max_lim);
252 0 : domains.push_back((max_lim - min_lim) / 2);
253 : } else {
254 6 : domains.push_back(-1.);
255 : }
256 : }
257 1 : }
258 :
259 : // Calculator
260 175001 : void FuncPathGeneral::calculate() {
261 : double s_path = 0.;
262 : double partition = 0.;
263 : double tmp, value, diff, expdist, s_der, z_der;
264 : int ii;
265 :
266 : typedef std::vector< std::pair< int,double> >::iterator pairiter;
267 :
268 2975001 : for (pairiter it = neighpair.begin(); it != neighpair.end(); ++it) {
269 2800000 : (*it).second = 0.;
270 : }
271 :
272 175001 : if (neighpair.empty()) {
273 : // Resize at the first step
274 1 : neighpair.resize(path_cv_values.size());
275 17 : for (unsigned i = 0; i < path_cv_values.size(); ++i) {
276 16 : neighpair[i].first = i;
277 : }
278 : }
279 :
280 175001 : Value* val_s_path=getPntrToComponent("s");
281 175001 : Value* val_z_path=getPntrToComponent("z");
282 :
283 1225007 : for(unsigned j = 0; j < allArguments.size(); ++j) {
284 1050006 : value = allArguments[j]->get();
285 17850102 : for (pairiter it = neighpair.begin(); it != neighpair.end(); ++it) {
286 16800096 : diff = (value - path_cv_values[(*it).first][j]);
287 16800096 : if (domains[j] > 0) {
288 0 : if (diff > domains[j]) {
289 0 : diff -= 2 * domains[j];
290 : }
291 0 : if (diff < -domains[j]) {
292 0 : diff += 2 * domains[j];
293 : }
294 : }
295 33600192 : (*it).second += Tools::fastpow(coefficients[j] * diff, 2);
296 33600192 : numerators[(*it).first][j] = 2 * Tools::fastpow(coefficients[j], 2) * diff;
297 : }
298 : }
299 :
300 : // The weights are computed relative to the closest frame, exp(-lambda*(d_i-d_min)),
301 : // so that they do not all underflow to zero far from the path.
302 175001 : double dmin = neighpair[0].second;
303 2975017 : for (pairiter it = neighpair.begin(); it != neighpair.end(); ++it) {
304 3756614 : dmin = std::min(dmin, (*it).second);
305 : }
306 2975017 : for (pairiter it = neighpair.begin(); it != neighpair.end(); ++it) {
307 2800016 : expdist = std::exp(-lambda * ((*it).second - dmin));
308 2800016 : expdists[(*it).first] = expdist;
309 2800016 : s_path += ((*it).first + 1) * expdist;
310 2800016 : partition += expdist;
311 : }
312 : // the closest frame has weight one, so partition>=1; the check makes this explicit
313 : // (and keeps static analysis from flagging the logarithm below)
314 175001 : if (!(partition >= 1.0)) {
315 0 : error("the weights of the frames of the path are not finite: check LAMBDA and the arguments");
316 : return;
317 : }
318 :
319 175001 : s_path /= partition;
320 : val_s_path->set(s_path);
321 175001 : val_z_path->set(dmin - (1. / lambda) * std::log(partition));
322 :
323 : // Derivatives
324 2975017 : for (pairiter it = neighpair.begin(); it != neighpair.end(); ++it) {
325 2800016 : ii = (*it).first;
326 2800016 : tmp = lambda * expdists[ii] * (s_path - (ii + 1)) / partition;
327 2800016 : s_path_ders[ii] = tmp;
328 2800016 : z_path_ders[ii] = expdists[ii] / partition;
329 : }
330 1225007 : for (unsigned i = 0; i < coefficients.size(); ++i) {
331 : s_der = 0.;
332 : z_der = 0.;
333 17850102 : for (pairiter it = neighpair.begin(); it != neighpair.end(); ++it) {
334 16800096 : ii = (*it).first;
335 16800096 : s_der += s_path_ders[ii] * numerators[ii][i];
336 16800096 : z_der += z_path_ders[ii] * numerators[ii][i];
337 : }
338 : setDerivative(val_s_path, i, s_der);
339 : setDerivative(val_z_path, i, z_der);
340 : }
341 : }
342 :
343 : // Prepare the required arguments
344 175001 : void FuncPathGeneral::prepare() {
345 : // Neighbour list: rank and activate the chain for the next step
346 :
347 : // Neighbour list: if neigh_size < 0 never sort and keep the full vector
348 : // Neighbour list: if neigh_size > 0
349 : // if the size is full -> sort the vector and decide the dependencies for next step
350 : // if the size is not full -> check if next step will need the full dependency otherwise keep these dependencies
351 :
352 175001 : if (neigh_size > 0) {
353 0 : if (neighpair.size() == path_cv_values.size()) {
354 : // The complete round has been done: need to sort, shorten and give it a go
355 : // Sort the values
356 0 : std::sort(neighpair.begin(), neighpair.end(), pairordering());
357 : // Resize the effective list
358 0 : neighpair.resize(neigh_size);
359 0 : log.printf(" NEIGHBOUR LIST NOW INCLUDES INDICES: ");
360 0 : for (int i = 0; i < neigh_size; ++i) {
361 0 : log.printf(" %i ",neighpair[i].first);
362 : }
363 0 : log.printf(" \n");
364 : } else {
365 0 : const long long nstride = std::max(1LL, static_cast<long long>(neigh_stride / getTimeStep()));
366 0 : if (getStep() % nstride == 0) {
367 0 : log.printf(" Time %f : recalculating full neighbour list \n", getStep() * getTimeStep());
368 0 : neighpair.resize(path_cv_values.size());
369 0 : for (unsigned i = 0; i < path_cv_values.size(); ++i) {
370 0 : neighpair[i].first = i;
371 : }
372 : }
373 : }
374 : }
375 : // all the arguments are always needed, so the argument list never changes
376 175001 : }
377 :
378 : }
379 : }
|
