c34dc772157c64063627759370a8a936a2d46c13
[namd.git] / colvars / src / colvarcomp_coordnums.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 <cmath>
11
12 #include "colvarmodule.h"
13 #include "colvarparse.h"
14 #include "colvaratoms.h"
15 #include "colvarvalue.h"
16 #include "colvar.h"
17 #include "colvarcomp.h"
18
19
20
21
22 template<bool calculate_gradients>
23 cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
24                                                int const &en,
25                                                int const &ed,
26                                                cvm::atom &A1,
27                                                cvm::atom &A2)
28 {
29   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
30   cvm::real const l2 = diff.norm2()/(r0*r0);
31
32   // Assume en and ed are even integers, and avoid sqrt in the following
33   int const en2 = en/2;
34   int const ed2 = ed/2;
35
36   cvm::real const xn = cvm::integer_power(l2, en2);
37   cvm::real const xd = cvm::integer_power(l2, ed2);
38   cvm::real const func = (1.0-xn)/(1.0-xd);
39
40   if (calculate_gradients) {
41     cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
42     cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
43     A1.grad += (-1.0)*dFdl2*dl2dx;
44     A2.grad +=        dFdl2*dl2dx;
45   }
46
47   return func;
48 }
49
50
51 template<bool calculate_gradients>
52 cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
53                                                int const &en,
54                                                int const &ed,
55                                                cvm::atom &A1,
56                                                cvm::atom &A2)
57 {
58   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
59   cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
60   cvm::real const l2 = scal_diff.norm2();
61
62   // Assume en and ed are even integers, and avoid sqrt in the following
63   int const en2 = en/2;
64   int const ed2 = ed/2;
65
66   cvm::real const xn = cvm::integer_power(l2, en2);
67   cvm::real const xd = cvm::integer_power(l2, ed2);
68   cvm::real const func = (1.0-xn)/(1.0-xd);
69
70   if (calculate_gradients) {
71     cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
72     cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
73                              (2.0/(r0_vec.y*r0_vec.y))*diff.y,
74                              (2.0/(r0_vec.z*r0_vec.z))*diff.z);
75     A1.grad += (-1.0)*dFdl2*dl2dx;
76     A2.grad +=        dFdl2*dl2dx;
77   }
78   return func;
79 }
80
81
82 colvar::coordnum::coordnum(std::string const &conf)
83   : cvc(conf), b_anisotropic(false), b_group2_center_only(false)
84 {
85   function_type = "coordnum";
86   x.type(colvarvalue::type_scalar);
87
88   group1 = parse_group(conf, "group1");
89   group2 = parse_group(conf, "group2");
90
91   if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
92     cvm::error("Error: group1 and group2 share a common atom (number: " +
93       cvm::to_str(atom_number) + ")\n");
94     return;
95   }
96
97   if (group1->b_dummy) {
98     cvm::error("Error: only group2 is allowed to be a dummy atom\n");
99     return;
100   }
101
102   bool const b_isotropic = get_keyval(conf, "cutoff", r0,
103                                       cvm::real(4.0 * cvm::unit_angstrom()));
104
105   if (get_keyval(conf, "cutoff3", r0_vec, cvm::rvector(4.0 * cvm::unit_angstrom(),
106                                                        4.0 * cvm::unit_angstrom(),
107                                                        4.0 * cvm::unit_angstrom()))) {
108     if (b_isotropic) {
109       cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
110                  INPUT_ERROR);
111       return;
112     }
113
114     b_anisotropic = true;
115     // remove meaningless negative signs
116     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
117     if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
118     if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
119   }
120
121   get_keyval(conf, "expNumer", en, 6);
122   get_keyval(conf, "expDenom", ed, 12);
123
124   if ( (en%2) || (ed%2) ) {
125     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
126                INPUT_ERROR);
127   }
128
129   if ( (en <= 0) || (ed <= 0) ) {
130     cvm::error("Error: negative exponent(s) provided.\n",
131                INPUT_ERROR);
132   }
133
134   if (!is_enabled(f_cvc_pbc_minimum_image)) {
135     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
136   }
137
138   get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
139 }
140
141
142 colvar::coordnum::coordnum()
143   : b_anisotropic(false), b_group2_center_only(false)
144 {
145   function_type = "coordnum";
146   x.type(colvarvalue::type_scalar);
147 }
148
149
150 void colvar::coordnum::calc_value()
151 {
152   x.real_value = 0.0;
153
154   if (b_group2_center_only) {
155
156     // create a fake atom to hold the group2 com coordinates
157     cvm::atom group2_com_atom;
158     group2_com_atom.pos = group2->center_of_mass();
159
160     if (b_anisotropic) {
161       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
162         x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);
163     } else {
164       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
165         x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);
166     }
167
168   } else {
169
170     if (b_anisotropic) {
171       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
172         for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
173           x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);
174         }
175     } else {
176       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
177         for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
178           x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);
179         }
180     }
181   }
182 }
183
184
185 void colvar::coordnum::calc_gradients()
186 {
187   if (b_group2_center_only) {
188
189     // create a fake atom to hold the group2 com coordinates
190     cvm::atom group2_com_atom;
191     group2_com_atom.pos = group2->center_of_mass();
192
193
194     if (b_anisotropic) {
195       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
196         switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);
197     } else {
198       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
199         switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
200     }
201
202     group2->set_weighted_gradient(group2_com_atom.grad);
203
204   } else {
205
206     if (b_anisotropic) {
207       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
208         for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
209           switching_function<true>(r0_vec, en, ed, *ai1, *ai2);
210         }
211     } else {
212       for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
213         for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
214           switching_function<true>(r0, en, ed, *ai1, *ai2);
215         }
216     }
217   }
218 }
219
220
221 void colvar::coordnum::apply_force(colvarvalue const &force)
222 {
223   if (!group1->noforce)
224     group1->apply_colvar_force(force.real_value);
225
226   if (!group2->noforce)
227     group2->apply_colvar_force(force.real_value);
228 }
229
230
231 simple_scalar_dist_functions(coordnum)
232
233
234
235 // h_bond member functions
236
237 colvar::h_bond::h_bond(std::string const &conf)
238   : cvc(conf)
239 {
240   if (cvm::debug())
241     cvm::log("Initializing h_bond object.\n");
242
243   function_type = "h_bond";
244   x.type(colvarvalue::type_scalar);
245
246   int a_num, d_num;
247   get_keyval(conf, "acceptor", a_num, -1);
248   get_keyval(conf, "donor",    d_num, -1);
249
250   if ( (a_num == -1) || (d_num == -1) ) {
251     cvm::error("Error: either acceptor or donor undefined.\n");
252     return;
253   }
254
255   cvm::atom acceptor = cvm::atom(a_num);
256   cvm::atom donor    = cvm::atom(d_num);
257   register_atom_group(new cvm::atom_group);
258   atom_groups[0]->add_atom(acceptor);
259   atom_groups[0]->add_atom(donor);
260
261   get_keyval(conf, "cutoff",   r0, (3.3 * cvm::unit_angstrom()));
262   get_keyval(conf, "expNumer", en, 6);
263   get_keyval(conf, "expDenom", ed, 8);
264
265   if ( (en%2) || (ed%2) ) {
266     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
267                INPUT_ERROR);
268   }
269
270   if ( (en <= 0) || (ed <= 0) ) {
271     cvm::error("Error: negative exponent(s) provided.\n",
272                INPUT_ERROR);
273   }
274
275   if (cvm::debug())
276     cvm::log("Done initializing h_bond object.\n");
277 }
278
279
280 colvar::h_bond::h_bond(cvm::atom const &acceptor,
281                        cvm::atom const &donor,
282                        cvm::real r0_i, int en_i, int ed_i)
283   : r0(r0_i), en(en_i), ed(ed_i)
284 {
285   function_type = "h_bond";
286   x.type(colvarvalue::type_scalar);
287
288   register_atom_group(new cvm::atom_group);
289   atom_groups[0]->add_atom(acceptor);
290   atom_groups[0]->add_atom(donor);
291 }
292
293
294 colvar::h_bond::h_bond()
295   : cvc()
296 {
297   function_type = "h_bond";
298   x.type(colvarvalue::type_scalar);
299 }
300
301
302 colvar::h_bond::~h_bond()
303 {
304   delete atom_groups[0];
305 }
306
307
308 void colvar::h_bond::calc_value()
309 {
310   x.real_value = colvar::coordnum::switching_function<false>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
311 }
312
313
314 void colvar::h_bond::calc_gradients()
315 {
316   colvar::coordnum::switching_function<true>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
317 }
318
319
320 void colvar::h_bond::apply_force(colvarvalue const &force)
321 {
322   (atom_groups[0])->apply_colvar_force(force);
323 }
324
325
326 simple_scalar_dist_functions(h_bond)
327
328
329
330 colvar::selfcoordnum::selfcoordnum(std::string const &conf)
331   : cvc(conf)
332 {
333   function_type = "selfcoordnum";
334   x.type(colvarvalue::type_scalar);
335
336   group1 = parse_group(conf, "group1");
337
338   get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom()));
339   get_keyval(conf, "expNumer", en, 6);
340   get_keyval(conf, "expDenom", ed, 12);
341
342
343   if ( (en%2) || (ed%2) ) {
344     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
345                INPUT_ERROR);
346   }
347
348   if ( (en <= 0) || (ed <= 0) ) {
349     cvm::error("Error: negative exponent(s) provided.\n",
350                INPUT_ERROR);
351   }
352
353   if (!is_enabled(f_cvc_pbc_minimum_image)) {
354     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
355   }
356 }
357
358
359 colvar::selfcoordnum::selfcoordnum()
360 {
361   function_type = "selfcoordnum";
362   x.type(colvarvalue::type_scalar);
363 }
364
365
366 void colvar::selfcoordnum::calc_value()
367 {
368   x.real_value = 0.0;
369   for (size_t i = 0; i < group1->size() - 1; i++) {
370     for (size_t j = i + 1; j < group1->size(); j++) {
371       x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, (*group1)[i], (*group1)[j]);
372     }
373   }
374 }
375
376
377 void colvar::selfcoordnum::calc_gradients()
378 {
379   for (size_t i = 0; i < group1->size() - 1; i++) {
380     for (size_t j = i + 1; j < group1->size(); j++) {
381       colvar::coordnum::switching_function<true>(r0, en, ed, (*group1)[i], (*group1)[j]);
382     }
383   }
384 }
385
386 void colvar::selfcoordnum::apply_force(colvarvalue const &force)
387 {
388   if (!group1->noforce) {
389     group1->apply_colvar_force(force.real_value);
390   }
391 }
392
393
394 simple_scalar_dist_functions(selfcoordnum)
395
396
397 // groupcoordnum member functions
398 colvar::groupcoordnum::groupcoordnum(std::string const &conf)
399   : distance(conf), b_anisotropic(false)
400 {
401   function_type = "groupcoordnum";
402   x.type(colvarvalue::type_scalar);
403
404   // group1 and group2 are already initialized by distance()
405   if (group1->b_dummy || group2->b_dummy) {
406     cvm::error("Error: neither group can be a dummy atom\n");
407     return;
408   }
409
410   bool const b_scale = get_keyval(conf, "cutoff", r0,
411                                   cvm::real(4.0 * cvm::unit_angstrom()));
412
413   if (get_keyval(conf, "cutoff3", r0_vec,
414                  cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
415
416     if (b_scale) {
417       cvm::error("Error: cannot specify \"scale\" and "
418                        "\"scale3\" at the same time.\n");
419       return;
420     }
421     b_anisotropic = true;
422     // remove meaningless negative signs
423     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
424     if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
425     if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
426   }
427
428   get_keyval(conf, "expNumer", en, 6);
429   get_keyval(conf, "expDenom", ed, 12);
430
431   if ( (en%2) || (ed%2) ) {
432     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
433                INPUT_ERROR);
434   }
435
436   if ( (en <= 0) || (ed <= 0) ) {
437     cvm::error("Error: negative exponent(s) provided.\n",
438                INPUT_ERROR);
439   }
440
441   if (!is_enabled(f_cvc_pbc_minimum_image)) {
442     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
443   }
444
445 }
446
447
448 colvar::groupcoordnum::groupcoordnum()
449   : b_anisotropic(false)
450 {
451   function_type = "groupcoordnum";
452   x.type(colvarvalue::type_scalar);
453 }
454
455
456 template<bool calculate_gradients>
457 cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
458                                                     int const &en,
459                                                     int const &ed,
460                                                     cvm::atom &A1,
461                                                     cvm::atom &A2)
462 {
463   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
464   cvm::real const l2 = diff.norm2()/(r0*r0);
465
466   // Assume en and ed are even integers, and avoid sqrt in the following
467   int const en2 = en/2;
468   int const ed2 = ed/2;
469
470   cvm::real const xn = cvm::integer_power(l2, en2);
471   cvm::real const xd = cvm::integer_power(l2, ed2);
472   cvm::real const func = (1.0-xn)/(1.0-xd);
473
474   if (calculate_gradients) {
475     cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
476     cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
477     A1.grad += (-1.0)*dFdl2*dl2dx;
478     A2.grad +=        dFdl2*dl2dx;
479   }
480
481   return func;
482 }
483
484
485 #if 0  // AMG: I don't think there's any reason to support anisotropic,
486        //      and I don't have those flags below in calc_value, but
487        //      if I need them, I'll also need to uncomment this method
488 template<bool calculate_gradients>
489 cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
490                                                     int const &en,
491                                                     int const &ed,
492                                                     cvm::atom &A1,
493                                                     cvm::atom &A2)
494 {
495   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
496   cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
497   cvm::real const l2 = scal_diff.norm2();
498
499   // Assume en and ed are even integers, and avoid sqrt in the following
500   int const en2 = en/2;
501   int const ed2 = ed/2;
502
503   cvm::real const xn = cvm::integer_power(l2, en2);
504   cvm::real const xd = cvm::integer_power(l2, ed2);
505   cvm::real const func = (1.0-xn)/(1.0-xd);
506
507   if (calculate_gradients) {
508     cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
509     cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
510                              (2.0/(r0_vec.y*r0_vec.y))*diff.y,
511                              (2.0/(r0_vec.z*r0_vec.z))*diff.z);
512     A1.grad += (-1.0)*dFdl2*dl2dx;
513     A2.grad +=        dFdl2*dl2dx;
514   }
515   return func;
516 }
517 #endif
518
519
520 void colvar::groupcoordnum::calc_value()
521 {
522
523   // create fake atoms to hold the com coordinates
524   cvm::atom group1_com_atom;
525   cvm::atom group2_com_atom;
526   group1_com_atom.pos = group1->center_of_mass();
527   group2_com_atom.pos = group2->center_of_mass();
528
529   x.real_value = coordnum::switching_function<false>(r0, en, ed,
530                                                      group1_com_atom, group2_com_atom);
531 }
532
533
534 void colvar::groupcoordnum::calc_gradients()
535 {
536   cvm::atom group1_com_atom;
537   cvm::atom group2_com_atom;
538   group1_com_atom.pos = group1->center_of_mass();
539   group2_com_atom.pos = group2->center_of_mass();
540
541   coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
542   group1->set_weighted_gradient(group1_com_atom.grad);
543   group2->set_weighted_gradient(group2_com_atom.grad);
544
545 }
546
547
548 void colvar::groupcoordnum::apply_force(colvarvalue const &force)
549 {
550   if (!group1->noforce)
551     group1->apply_colvar_force(force.real_value);
552
553   if (!group2->noforce)
554     group2->apply_colvar_force(force.real_value);
555 }
556
557
558 simple_scalar_dist_functions(groupcoordnum)