Update Colvars to version 2018-05-15
[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 template<int flags>
22 cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
23                                                cvm::rvector const &r0_vec,
24                                                int en,
25                                                int ed,
26                                                cvm::atom &A1,
27                                                cvm::atom &A2,
28                                                bool **pairlist_elem,
29                                                cvm::real pairlist_tol)
30 {
31   if ((flags & ef_use_pairlist) && !(flags & ef_rebuild_pairlist)) {
32     bool const within = **pairlist_elem;
33     (*pairlist_elem)++;
34     if (!within) {
35       return 0.0;
36     }
37   }
38
39   cvm::rvector const r0sq_vec(r0_vec.x*r0_vec.x,
40                               r0_vec.y*r0_vec.y,
41                               r0_vec.z*r0_vec.z);
42
43   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
44
45   cvm::rvector const scal_diff(diff.x/((flags & ef_anisotropic) ?
46                                        r0_vec.x : r0),
47                                diff.y/((flags & ef_anisotropic) ?
48                                        r0_vec.y : r0),
49                                diff.z/((flags & ef_anisotropic) ?
50                                        r0_vec.z : r0));
51   cvm::real const l2 = scal_diff.norm2();
52
53   // Assume en and ed are even integers, and avoid sqrt in the following
54   int const en2 = en/2;
55   int const ed2 = ed/2;
56
57   cvm::real const xn = cvm::integer_power(l2, en2);
58   cvm::real const xd = cvm::integer_power(l2, ed2);
59   cvm::real const func = (1.0-xn)/(1.0-xd);
60
61   if (flags & ef_gradients) {
62     cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) -
63                                             func*ed2*(xd/l2))*(-1.0);
64     cvm::rvector const dl2dx((2.0/((flags & ef_anisotropic) ? r0sq_vec.x :
65                                    r0*r0)) * diff.x,
66                              (2.0/((flags & ef_anisotropic) ? r0sq_vec.y :
67                                    r0*r0)) * diff.y,
68                              (2.0/((flags & ef_anisotropic) ? r0sq_vec.z :
69                                    r0*r0)) * diff.z);
70     A1.grad += (-1.0)*dFdl2*dl2dx;
71     A2.grad +=        dFdl2*dl2dx;
72   }
73
74   if (flags & ef_rebuild_pairlist) {
75     **pairlist_elem = (func > pairlist_tol) ? true : false;
76     (*pairlist_elem)++;
77   }
78
79   return func;
80 }
81
82
83 colvar::coordnum::coordnum(std::string const &conf)
84   : cvc(conf), b_anisotropic(false), group2_center(NULL)
85
86 {
87   function_type = "coordnum";
88   x.type(colvarvalue::type_scalar);
89
90   group1 = parse_group(conf, "group1");
91   group2 = parse_group(conf, "group2");
92
93   if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
94     cvm::error("Error: group1 and group2 share a common atom (number: " +
95                cvm::to_str(atom_number) + ")\n", INPUT_ERROR);
96     return;
97   }
98
99   if (group1->b_dummy) {
100     cvm::error("Error: only group2 is allowed to be a dummy atom\n",
101                INPUT_ERROR);
102     return;
103   }
104
105   bool const b_isotropic = get_keyval(conf, "cutoff", r0,
106                                       cvm::real(4.0 * cvm::unit_angstrom()));
107
108   if (get_keyval(conf, "cutoff3", r0_vec,
109                  cvm::rvector(4.0 * cvm::unit_angstrom(),
110                               4.0 * cvm::unit_angstrom(),
111                               4.0 * cvm::unit_angstrom()))) {
112     if (b_isotropic) {
113       cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" "
114                  "at the same time.\n",
115                  INPUT_ERROR);
116       return;
117     }
118
119     b_anisotropic = true;
120     // remove meaningless negative signs
121     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
122     if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
123     if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
124   }
125
126   get_keyval(conf, "expNumer", en, 6);
127   get_keyval(conf, "expDenom", ed, 12);
128
129   if ( (en%2) || (ed%2) ) {
130     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
131                INPUT_ERROR);
132   }
133
134   if ( (en <= 0) || (ed <= 0) ) {
135     cvm::error("Error: negative exponent(s) provided.\n",
136                INPUT_ERROR);
137   }
138
139   if (!is_enabled(f_cvc_pbc_minimum_image)) {
140     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
141   }
142
143   bool b_group2_center_only = false;
144   get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
145   if (b_group2_center_only) {
146     if (!group2_center) {
147       group2_center = new cvm::atom_group();
148       group2_center->add_atom(cvm::atom());
149     }
150   }
151
152   get_keyval(conf, "tolerance", tolerance, 0.0);
153   if (tolerance > 0) {
154     get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
155     if ( ! (pairlist_freq > 0) ) {
156       cvm::error("Error: non-positive pairlistfrequency provided.\n",
157                  INPUT_ERROR);
158       return; // and do not allocate the pairlists below
159     }
160     if (b_group2_center_only) {
161       pairlist = new bool[group1->size()];
162     }
163     else {
164       pairlist = new bool[group1->size() * group2->size()];
165     }
166   }
167
168 }
169
170
171 colvar::coordnum::coordnum()
172   : b_anisotropic(false), group2_center(NULL)
173 {
174   function_type = "coordnum";
175   x.type(colvarvalue::type_scalar);
176 }
177
178
179 colvar::coordnum::~coordnum()
180 {
181   if (pairlist != NULL) {
182     delete [] pairlist;
183   }
184   if (group2_center != NULL) {
185     delete group2_center;
186   }
187 }
188
189
190 template<int compute_flags> int colvar::coordnum::compute_coordnum()
191 {
192   if (group2_center) {
193     (*group2_center)[0].pos = group2->center_of_mass();
194     group2_center->calc_required_properties();
195   }
196   cvm::atom_group *group2p = group2_center ? group2_center : group2;
197
198   bool const use_pairlist = (pairlist != NULL);
199   bool const rebuild_pairlist = (pairlist != NULL) &&
200     (cvm::step_relative() % pairlist_freq == 0);
201
202   bool *pairlist_elem = use_pairlist ? pairlist : NULL;
203   cvm::atom_iter ai1 = group1->begin(), ai2 = group2p->begin();
204   cvm::atom_iter const ai1_end = group1->end();
205   cvm::atom_iter const ai2_end = group2p->end();
206
207   if (b_anisotropic) {
208
209     if (use_pairlist) {
210
211       if (rebuild_pairlist) {
212
213         int const flags = compute_flags | ef_anisotropic | ef_use_pairlist |
214           ef_rebuild_pairlist;
215         for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
216           for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
217             x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
218                                                       *ai1, *ai2,
219                                                       &pairlist_elem,
220                                                       tolerance);
221           }
222         }
223
224       } else {
225
226         int const flags = compute_flags | ef_anisotropic | ef_use_pairlist;
227         for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
228           for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
229             x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
230                                                       *ai1, *ai2,
231                                                       &pairlist_elem,
232                                                       tolerance);
233           }
234         }
235       }
236
237     } else { // if (use_pairlist) {
238
239       int const flags = compute_flags | ef_anisotropic;
240       for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
241         for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
242           x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
243                                                     *ai1, *ai2,
244                                                     NULL, 0.0);
245         }
246       }
247     }
248
249   } else {
250
251     if (use_pairlist) {
252
253       if (rebuild_pairlist) {
254
255         int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
256         for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
257           for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
258             x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
259                                                       *ai1, *ai2,
260                                                       &pairlist_elem,
261                                                       tolerance);
262           }
263         }
264
265       } else {
266
267         int const flags = compute_flags | ef_use_pairlist;
268         for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
269           for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
270             x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
271                                                       *ai1, *ai2,
272                                                       &pairlist_elem,
273                                                       tolerance);
274           }
275         }
276       }
277
278     } else { // if (use_pairlist) {
279
280       int const flags = compute_flags;
281       for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
282         for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
283           x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
284                                                     *ai1, *ai2,
285                                                     NULL, 0.0);
286         }
287       }
288     }
289   }
290
291   if (compute_flags & ef_gradients) {
292     if (group2_center) {
293       group2->set_weighted_gradient((*group2_center)[0].grad);
294     }
295   }
296
297   return COLVARS_OK;
298 }
299
300
301 void colvar::coordnum::calc_value()
302 {
303   x.real_value = 0.0;
304   if (is_enabled(f_cvc_gradient)) {
305     compute_coordnum<ef_gradients>();
306   } else {
307     compute_coordnum<ef_null>();
308   }
309 }
310
311
312 void colvar::coordnum::calc_gradients()
313 {
314   // Gradients are computed by calc_value() if f_cvc_gradients is enabled
315 }
316
317
318 void colvar::coordnum::apply_force(colvarvalue const &force)
319 {
320   if (!group1->noforce)
321     group1->apply_colvar_force(force.real_value);
322
323   if (!group2->noforce)
324     group2->apply_colvar_force(force.real_value);
325 }
326
327
328 simple_scalar_dist_functions(coordnum)
329
330
331
332 // h_bond member functions
333
334 colvar::h_bond::h_bond(std::string const &conf)
335   : cvc(conf)
336 {
337   if (cvm::debug())
338     cvm::log("Initializing h_bond object.\n");
339
340   function_type = "h_bond";
341   x.type(colvarvalue::type_scalar);
342
343   int a_num, d_num;
344   get_keyval(conf, "acceptor", a_num, -1);
345   get_keyval(conf, "donor",    d_num, -1);
346
347   if ( (a_num == -1) || (d_num == -1) ) {
348     cvm::error("Error: either acceptor or donor undefined.\n");
349     return;
350   }
351
352   cvm::atom acceptor = cvm::atom(a_num);
353   cvm::atom donor    = cvm::atom(d_num);
354   register_atom_group(new cvm::atom_group);
355   atom_groups[0]->add_atom(acceptor);
356   atom_groups[0]->add_atom(donor);
357
358   get_keyval(conf, "cutoff",   r0, (3.3 * cvm::unit_angstrom()));
359   get_keyval(conf, "expNumer", en, 6);
360   get_keyval(conf, "expDenom", ed, 8);
361
362   if ( (en%2) || (ed%2) ) {
363     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
364                INPUT_ERROR);
365   }
366
367   if ( (en <= 0) || (ed <= 0) ) {
368     cvm::error("Error: negative exponent(s) provided.\n",
369                INPUT_ERROR);
370   }
371
372   if (cvm::debug())
373     cvm::log("Done initializing h_bond object.\n");
374 }
375
376
377 colvar::h_bond::h_bond(cvm::atom const &acceptor,
378                        cvm::atom const &donor,
379                        cvm::real r0_i, int en_i, int ed_i)
380   : r0(r0_i), en(en_i), ed(ed_i)
381 {
382   function_type = "h_bond";
383   x.type(colvarvalue::type_scalar);
384
385   register_atom_group(new cvm::atom_group);
386   atom_groups[0]->add_atom(acceptor);
387   atom_groups[0]->add_atom(donor);
388 }
389
390
391 colvar::h_bond::h_bond()
392   : cvc()
393 {
394   function_type = "h_bond";
395   x.type(colvarvalue::type_scalar);
396 }
397
398
399 colvar::h_bond::~h_bond()
400 {
401   delete atom_groups[0];
402 }
403
404
405 void colvar::h_bond::calc_value()
406 {
407   int const flags = coordnum::ef_null;
408   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
409   x.real_value =
410     coordnum::switching_function<flags>(r0, r0_vec, en, ed,
411                                         (*atom_groups[0])[0],
412                                         (*atom_groups[0])[1],
413                                         NULL, 0.0);
414 }
415
416
417 void colvar::h_bond::calc_gradients()
418 {
419   int const flags = coordnum::ef_gradients;
420   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
421   coordnum::switching_function<flags>(r0, r0_vec, en, ed,
422                                       (*atom_groups[0])[0],
423                                       (*atom_groups[0])[1],
424                                       NULL, 0.0);
425 }
426
427
428 void colvar::h_bond::apply_force(colvarvalue const &force)
429 {
430   (atom_groups[0])->apply_colvar_force(force);
431 }
432
433
434 simple_scalar_dist_functions(h_bond)
435
436
437
438 colvar::selfcoordnum::selfcoordnum(std::string const &conf)
439   : cvc(conf)
440 {
441   function_type = "selfcoordnum";
442   x.type(colvarvalue::type_scalar);
443
444   group1 = parse_group(conf, "group1");
445
446   get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom()));
447   get_keyval(conf, "expNumer", en, 6);
448   get_keyval(conf, "expDenom", ed, 12);
449
450
451   if ( (en%2) || (ed%2) ) {
452     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
453                INPUT_ERROR);
454   }
455
456   if ( (en <= 0) || (ed <= 0) ) {
457     cvm::error("Error: negative exponent(s) provided.\n",
458                INPUT_ERROR);
459   }
460
461   if (!is_enabled(f_cvc_pbc_minimum_image)) {
462     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
463   }
464
465   get_keyval(conf, "tolerance", tolerance, 0.0);
466   if (tolerance > 0) {
467     get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
468     if ( ! (pairlist_freq > 0) ) {
469       cvm::error("Error: non-positive pairlistfrequency provided.\n",
470                  INPUT_ERROR);
471       return;
472     }
473     pairlist = new bool[(group1->size()-1) * (group1->size()-1)];
474   }
475 }
476
477
478 colvar::selfcoordnum::selfcoordnum()
479 {
480   function_type = "selfcoordnum";
481   x.type(colvarvalue::type_scalar);
482 }
483
484
485 colvar::selfcoordnum::~selfcoordnum()
486 {
487   if (pairlist != NULL) {
488     delete [] pairlist;
489   }
490 }
491
492
493 template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
494 {
495   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
496
497   bool const use_pairlist = (pairlist != NULL);
498   bool const rebuild_pairlist = (pairlist != NULL) &&
499     (cvm::step_relative() % pairlist_freq == 0);
500
501   bool *pairlist_elem = use_pairlist ? pairlist : NULL;
502   size_t i = 0, j = 0;
503   size_t const n = group1->size();
504
505   // Always isotropic (TODO: enable the ellipsoid?)
506
507   if (use_pairlist) {
508
509     if (rebuild_pairlist) {
510       int const flags = compute_flags | coordnum::ef_use_pairlist |
511         coordnum::ef_rebuild_pairlist;
512       for (i = 0; i < n - 1; i++) {
513         for (j = i + 1; j < n; j++) {
514           x.real_value +=
515             coordnum::switching_function<flags>(r0, r0_vec, en, ed,
516                                                 (*group1)[i],
517                                                 (*group1)[j],
518                                                 &pairlist_elem,
519                                                 tolerance);
520         }
521       }
522     } else {
523       int const flags = compute_flags | coordnum::ef_use_pairlist;
524       for (i = 0; i < n - 1; i++) {
525         for (j = i + 1; j < n; j++) {
526           x.real_value +=
527             coordnum::switching_function<flags>(r0, r0_vec, en, ed,
528                                                 (*group1)[i],
529                                                 (*group1)[j],
530                                                 &pairlist_elem,
531                                                 tolerance);
532         }
533       }
534     }
535
536   } else { // if (use_pairlist) {
537
538     int const flags = compute_flags | coordnum::ef_null;
539     for (i = 0; i < n - 1; i++) {
540       for (j = i + 1; j < n; j++) {
541         x.real_value +=
542           coordnum::switching_function<flags>(r0, r0_vec, en, ed,
543                                               (*group1)[i],
544                                               (*group1)[j],
545                                               &pairlist_elem,
546                                               tolerance);
547       }
548     }
549   }
550
551   return COLVARS_OK;
552 }
553
554
555 void colvar::selfcoordnum::calc_value()
556 {
557   x.real_value = 0.0;
558   if (is_enabled(f_cvc_gradient)) {
559     compute_selfcoordnum<coordnum::ef_gradients>();
560   } else {
561     compute_selfcoordnum<coordnum::ef_null>();
562   }
563 }
564
565
566 void colvar::selfcoordnum::calc_gradients()
567 {
568   // Gradients are computed by calc_value() if f_cvc_gradients is enabled
569 }
570
571
572 void colvar::selfcoordnum::apply_force(colvarvalue const &force)
573 {
574   if (!group1->noforce) {
575     group1->apply_colvar_force(force.real_value);
576   }
577 }
578
579
580 simple_scalar_dist_functions(selfcoordnum)
581
582
583
584 // groupcoordnum member functions
585 colvar::groupcoordnum::groupcoordnum(std::string const &conf)
586   : distance(conf), b_anisotropic(false)
587 {
588   function_type = "groupcoordnum";
589   x.type(colvarvalue::type_scalar);
590
591   // group1 and group2 are already initialized by distance()
592   if (group1->b_dummy || group2->b_dummy) {
593     cvm::error("Error: neither group can be a dummy atom\n");
594     return;
595   }
596
597   bool const b_scale = get_keyval(conf, "cutoff", r0,
598                                   cvm::real(4.0 * cvm::unit_angstrom()));
599
600   if (get_keyval(conf, "cutoff3", r0_vec,
601                  cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
602
603     if (b_scale) {
604       cvm::error("Error: cannot specify \"scale\" and "
605                        "\"scale3\" at the same time.\n");
606       return;
607     }
608     b_anisotropic = true;
609     // remove meaningless negative signs
610     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
611     if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
612     if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
613   }
614
615   get_keyval(conf, "expNumer", en, 6);
616   get_keyval(conf, "expDenom", ed, 12);
617
618   if ( (en%2) || (ed%2) ) {
619     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
620                INPUT_ERROR);
621   }
622
623   if ( (en <= 0) || (ed <= 0) ) {
624     cvm::error("Error: negative exponent(s) provided.\n",
625                INPUT_ERROR);
626   }
627
628   if (!is_enabled(f_cvc_pbc_minimum_image)) {
629     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
630   }
631
632 }
633
634
635 colvar::groupcoordnum::groupcoordnum()
636   : b_anisotropic(false)
637 {
638   function_type = "groupcoordnum";
639   x.type(colvarvalue::type_scalar);
640 }
641
642
643 void colvar::groupcoordnum::calc_value()
644 {
645   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
646
647   // create fake atoms to hold the com coordinates
648   cvm::atom group1_com_atom;
649   cvm::atom group2_com_atom;
650   group1_com_atom.pos = group1->center_of_mass();
651   group2_com_atom.pos = group2->center_of_mass();
652   if (b_anisotropic) {
653     int const flags = coordnum::ef_anisotropic;
654     x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
655                                                        group1_com_atom,
656                                                        group2_com_atom,
657                                                        NULL, 0.0);
658   } else {
659     int const flags = coordnum::ef_null;
660     x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
661                                                        group1_com_atom,
662                                                        group2_com_atom,
663                                                        NULL, 0.0);
664   }
665 }
666
667
668 void colvar::groupcoordnum::calc_gradients()
669 {
670   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
671
672   cvm::atom group1_com_atom;
673   cvm::atom group2_com_atom;
674   group1_com_atom.pos = group1->center_of_mass();
675   group2_com_atom.pos = group2->center_of_mass();
676
677   if (b_anisotropic) {
678     int const flags = coordnum::ef_gradients | coordnum::ef_anisotropic;
679     coordnum::switching_function<flags>(r0, r0_vec, en, ed,
680                                         group1_com_atom,
681                                         group2_com_atom,
682                                         NULL, 0.0);
683   } else {
684     int const flags = coordnum::ef_gradients;
685     coordnum::switching_function<flags>(r0, r0_vec, en, ed,
686                                         group1_com_atom,
687                                         group2_com_atom,
688                                         NULL, 0.0);
689   }
690
691   group1->set_weighted_gradient(group1_com_atom.grad);
692   group2->set_weighted_gradient(group2_com_atom.grad);
693 }
694
695
696 void colvar::groupcoordnum::apply_force(colvarvalue const &force)
697 {
698   if (!group1->noforce)
699     group1->apply_colvar_force(force.real_value);
700
701   if (!group2->noforce)
702     group2->apply_colvar_force(force.real_value);
703 }
704
705
706 simple_scalar_dist_functions(groupcoordnum)