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