Update Colvars to version 2018-05-15
[namd.git] / colvars / src / colvarcomp.cpp
1 // -*- c++ -*-
2
3 // This file is part of the Collective Variables module (Colvars).
4 // The original version of Colvars and its updates are located at:
5 // https://github.com/colvars/colvars
6 // Please update all Colvars source files before making any changes.
7 // If you wish to distribute your changes, please submit them to the
8 // Colvars repository at GitHub.
9
10 #include "colvarmodule.h"
11 #include "colvarvalue.h"
12 #include "colvar.h"
13 #include "colvarcomp.h"
14
15
16
17 colvar::cvc::cvc()
18   : sup_coeff(1.0),
19     sup_np(1),
20     b_periodic(false),
21     b_try_scalable(true)
22 {
23   init_cvc_requires();
24   sup_coeff = 1.0;
25   period = 0.0;
26   wrap_center = 0.0;
27 }
28
29
30 colvar::cvc::cvc(std::string const &conf)
31   : sup_coeff(1.0),
32     sup_np(1),
33     b_periodic(false),
34     b_try_scalable(true)
35 {
36   init_cvc_requires();
37   sup_coeff = 1.0;
38   period = 0.0;
39   wrap_center = 0.0;
40   init(conf);
41 }
42
43
44 int colvar::cvc::init(std::string const &conf)
45 {
46   if (cvm::debug())
47     cvm::log("Initializing cvc base object.\n");
48
49   get_keyval(conf, "name", this->name, this->name);
50   if (name.size() > 0) {
51     // Temporary description until child object is initialized
52     description = "cvc " + name;
53   } else {
54     description = "uninitialized cvc";
55   }
56
57   get_keyval(conf, "componentCoeff", sup_coeff, sup_coeff);
58   get_keyval(conf, "componentExp", sup_np, sup_np);
59
60   get_keyval(conf, "period", period, period);
61   get_keyval(conf, "wrapAround", wrap_center, wrap_center);
62
63   get_keyval_feature(dynamic_cast<colvarparse *>(this), conf, "debugGradients",
64                      f_cvc_debug_gradient, false, parse_silent);
65
66   bool b_no_PBC = !is_enabled(f_cvc_pbc_minimum_image); // Enabled by default
67   get_keyval(conf, "forceNoPBC", b_no_PBC, b_no_PBC);
68   if (b_no_PBC) {
69     disable(f_cvc_pbc_minimum_image);
70   } else {
71     enable(f_cvc_pbc_minimum_image);
72   }
73
74   // Attempt scalable calculations when in parallel? (By default yes, if available)
75   get_keyval(conf, "scalable", b_try_scalable, b_try_scalable);
76
77   if (cvm::debug())
78     cvm::log("Done initializing cvc base object.\n");
79
80   return cvm::get_error();
81 }
82
83
84 int colvar::cvc::init_total_force_params(std::string const &conf)
85 {
86   if (cvm::get_error()) return COLVARS_ERROR;
87
88   if (get_keyval_feature(this, conf, "oneSiteSystemForce",
89                          f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
90     cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: "
91              "please use \"oneSiteTotalForce\" instead.\n");
92   }
93   if (get_keyval_feature(this, conf, "oneSiteTotalForce",
94                          f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
95     cvm::log("Computing total force on group 1 only");
96   }
97
98   if (! is_enabled(f_cvc_one_site_total_force)) {
99     // check whether any of the other atom groups is dummy
100     std::vector<cvm::atom_group *>::iterator agi = atom_groups.begin();
101     agi++;
102     for ( ; agi != atom_groups.end(); agi++) {
103       if ((*agi)->b_dummy) {
104         provide(f_cvc_inv_gradient, false);
105         provide(f_cvc_Jacobian, false);
106       }
107     }
108   }
109
110   return COLVARS_OK;
111 }
112
113
114 cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
115                                           char const *group_key,
116                                           bool optional)
117 {
118   cvm::atom_group *group = NULL;
119   std::string group_conf;
120
121   if (key_lookup(conf, group_key, &group_conf)) {
122     group = new cvm::atom_group(group_key);
123
124     if (b_try_scalable) {
125       if (is_available(f_cvc_scalable_com)
126           && is_enabled(f_cvc_com_based)
127           && !is_enabled(f_cvc_debug_gradient)) {
128         enable(f_cvc_scalable_com);
129         enable(f_cvc_scalable);
130         // The CVC makes the feature available;
131         // the atom group will enable it unless it needs to compute a rotational fit
132         group->provide(f_ag_scalable_com);
133       }
134
135       // TODO check for other types of parallelism here
136     }
137
138     if (group_conf.size() == 0) {
139       cvm::error("Error: atom group \""+group->key+
140                  "\" is set, but has no definition.\n",
141                  INPUT_ERROR);
142       return group;
143     }
144
145     cvm::increase_depth();
146     if (group->parse(group_conf) == COLVARS_OK) {
147       register_atom_group(group);
148     }
149     group->check_keywords(group_conf, group_key);
150     if (cvm::get_error()) {
151       cvm::error("Error parsing definition for atom group \""+
152                  std::string(group_key)+"\"\n.", INPUT_ERROR);
153     }
154     cvm::decrease_depth();
155
156   } else {
157     if (! optional) {
158       cvm::error("Error: definition for atom group \""+
159                  std::string(group_key)+"\" not found.\n");
160     }
161   }
162
163   return group;
164 }
165
166
167 int colvar::cvc::setup()
168 {
169   description = "cvc " + name;
170   return COLVARS_OK;
171 }
172
173
174 colvar::cvc::~cvc()
175 {
176   free_children_deps();
177   remove_all_children();
178   for (size_t i = 0; i < atom_groups.size(); i++) {
179     if (atom_groups[i] != NULL) delete atom_groups[i];
180   }
181 }
182
183 void colvar::cvc::read_data()
184 {
185   size_t ig;
186   for (ig = 0; ig < atom_groups.size(); ig++) {
187     cvm::atom_group &atoms = *(atom_groups[ig]);
188     atoms.reset_atoms_data();
189     atoms.read_positions();
190     atoms.calc_required_properties();
191     // each atom group will take care of its own fitting_group, if defined
192   }
193
194 ////  Don't try to get atom velocities, as no back-end currently implements it
195 //   if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) {
196 //     for (i = 0; i < cvcs.size(); i++) {
197 //       for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
198 //         cvcs[i]->atom_groups[ig]->read_velocities();
199 //       }
200 //     }
201 //   }
202 }
203
204
205 void colvar::cvc::calc_force_invgrads()
206 {
207   cvm::error("Error: calculation of inverse gradients is not implemented "
208              "for colvar components of type \""+function_type+"\".\n",
209              COLVARS_NOT_IMPLEMENTED);
210 }
211
212
213 void colvar::cvc::calc_Jacobian_derivative()
214 {
215   cvm::error("Error: calculation of inverse gradients is not implemented "
216              "for colvar components of type \""+function_type+"\".\n",
217              COLVARS_NOT_IMPLEMENTED);
218 }
219
220
221 void colvar::cvc::calc_fit_gradients()
222 {
223   for (size_t ig = 0; ig < atom_groups.size(); ig++) {
224     atom_groups[ig]->calc_fit_gradients();
225   }
226 }
227
228
229 void colvar::cvc::debug_gradients()
230 {
231   // this function should work for any scalar cvc:
232   // the only difference will be the name of the atom group (here, "group")
233   // NOTE: this assumes that groups for this cvc are non-overlapping,
234   // since atom coordinates are modified only within the current group
235
236   cvm::log("Debugging gradients for " + description);
237
238   for (size_t ig = 0; ig < atom_groups.size(); ig++) {
239     cvm::atom_group *group = atom_groups[ig];
240     if (group->b_dummy) continue;
241
242     cvm::rotation const rot_0 = group->rot;
243     cvm::rotation const rot_inv = group->rot.inverse();
244
245     cvm::real x_0 = x.real_value;
246     if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_0 = x[0];
247
248     // cvm::log("gradients     = "+cvm::to_str (gradients)+"\n");
249
250     cvm::atom_group *group_for_fit = group->fitting_group ? group->fitting_group : group;
251     cvm::atom_pos fit_gradient_sum, gradient_sum;
252
253     // print the values of the fit gradients
254     if (group->b_rotate || group->b_center) {
255       if (group->is_enabled(f_ag_fit_gradients)) {
256         size_t j;
257
258         // fit_gradients are in the simulation frame: we should print them in the rotated frame
259         cvm::log("Fit gradients:\n");
260         for (j = 0; j < group_for_fit->fit_gradients.size(); j++) {
261           cvm::log((group->fitting_group ? std::string("refPosGroup") : group->key) +
262                   "[" + cvm::to_str(j) + "] = " +
263                   (group->b_rotate ?
264                     cvm::to_str(rot_0.rotate(group_for_fit->fit_gradients[j])) :
265                     cvm::to_str(group_for_fit->fit_gradients[j])));
266         }
267       }
268     }
269
270     // debug the gradients
271     for (size_t ia = 0; ia < group->size(); ia++) {
272
273       // tests are best conducted in the unrotated (simulation) frame
274       cvm::rvector const atom_grad = (group->b_rotate ?
275                                       rot_inv.rotate((*group)[ia].grad) :
276                                       (*group)[ia].grad);
277       gradient_sum += atom_grad;
278
279       for (size_t id = 0; id < 3; id++) {
280         // (re)read original positions
281         group->read_positions();
282         // change one coordinate
283         (*group)[ia].pos[id] += cvm::debug_gradients_step_size;
284         group->calc_required_properties();
285         calc_value();
286         cvm::real x_1 = x.real_value;
287         if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0];
288         cvm::log("Atom "+cvm::to_str(ia)+", component "+cvm::to_str(id)+":\n");
289         cvm::log("dx(actual) = "+cvm::to_str(x_1 - x_0,
290                               21, 14)+"\n");
291         cvm::real const dx_pred = (group->fit_gradients.size()) ?
292           (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) :
293           (cvm::debug_gradients_step_size * atom_grad[id]);
294         cvm::log("dx(interp) = "+cvm::to_str(dx_pred,
295                               21, 14)+"\n");
296         cvm::log("|dx(actual) - dx(interp)|/|dx(actual)| = "+
297                   cvm::to_str(std::fabs(x_1 - x_0 - dx_pred) /
298                               std::fabs(x_1 - x_0), 12, 5)+"\n");
299       }
300     }
301
302     if ((group->is_enabled(f_ag_fit_gradients)) && (group->fitting_group != NULL)) {
303       cvm::atom_group *ref_group = group->fitting_group;
304       group->read_positions();
305       group->calc_required_properties();
306
307       for (size_t ia = 0; ia < ref_group->size(); ia++) {
308
309         // fit gradients are in the unrotated (simulation) frame
310         cvm::rvector const atom_grad = ref_group->fit_gradients[ia];
311         fit_gradient_sum += atom_grad;
312
313         for (size_t id = 0; id < 3; id++) {
314           // (re)read original positions
315           group->read_positions();
316           ref_group->read_positions();
317           // change one coordinate
318           (*ref_group)[ia].pos[id] += cvm::debug_gradients_step_size;
319           group->calc_required_properties();
320           calc_value();
321
322           cvm::real const x_1 = x.real_value;
323           cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n");
324           cvm::log("dx(actual) = "+cvm::to_str (x_1 - x_0,
325                                 21, 14)+"\n");
326
327           cvm::real const dx_pred = cvm::debug_gradients_step_size * atom_grad[id];
328
329           cvm::log("dx(interp) = "+cvm::to_str (dx_pred,
330                                 21, 14)+"\n");
331           cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
332                     cvm::to_str(std::fabs (x_1 - x_0 - dx_pred) /
333                                 std::fabs (x_1 - x_0),
334                                 12, 5)+
335                     ".\n");
336         }
337       }
338     }
339
340     cvm::log("Gradient sum: " +  cvm::to_str(gradient_sum) +
341           "  Fit gradient sum: " + cvm::to_str(fit_gradient_sum) +
342           "  Total " + cvm::to_str(gradient_sum + fit_gradient_sum));
343   }
344   return;
345 }
346
347
348 cvm::real colvar::cvc::dist2(colvarvalue const &x1,
349                              colvarvalue const &x2) const
350 {
351   return x1.dist2(x2);
352 }
353
354
355 colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
356                                      colvarvalue const &x2) const
357 {
358   return x1.dist2_grad(x2);
359 }
360
361
362 colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
363                                      colvarvalue const &x2) const
364 {
365   return x2.dist2_grad(x1);
366 }
367
368
369 void colvar::cvc::wrap(colvarvalue &x) const
370 {
371   return;
372 }
373
374
375 // Static members
376
377 std::vector<colvardeps::feature *> colvar::cvc::cvc_features;