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