Skip to content

Commit b3df77f

Browse files
giacomofioringiacomofiorin
committed
cleanup: Use vector cutoff for coordNum even in the isotropic case
1 parent e89d39e commit b3df77f

File tree

2 files changed

+113
-141
lines changed

2 files changed

+113
-141
lines changed

‎src/colvarcomp_coordnums.cpp‎

Lines changed: 107 additions & 125 deletions
Original file line numberDiff line numberDiff line change
@@ -19,11 +19,9 @@ colvar::coordnum::coordnum()
1919
{
2020
set_function_type("coordNum");
2121
x.type(colvarvalue::type_scalar);
22-
colvarproxy *proxy = cvm::main()->proxy;
23-
r0 = proxy->angstrom_to_internal(4.0);
24-
r0_vec = cvm::rvector(proxy->angstrom_to_internal(4.0),
25-
proxy->angstrom_to_internal(4.0),
26-
proxy->angstrom_to_internal(4.0));
22+
cvm::real const r0 = cvm::main()->proxy->angstrom_to_internal(4.0);
23+
r0_vec = {r0, r0, r0};
24+
// Default upper boundary not yet known
2725
}
2826

2927

@@ -49,20 +47,25 @@ int colvar::coordnum::init(std::string const &conf)
4947
cvm::error("Error: only group2 is allowed to be a dummy atom\n", COLVARS_INPUT_ERROR);
5048
}
5149

52-
bool const b_isotropic = get_keyval(conf, "cutoff", r0, r0);
50+
// Get the default value from r0_vec to report it
51+
cvm::real r0 = r0_vec[0];
52+
bool const b_redefined_cutoff = get_keyval(conf, "cutoff", r0, r0);
5353

5454
if (get_keyval(conf, "cutoff3", r0_vec, r0_vec)) {
55-
if (b_isotropic) {
56-
error_code |= cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" "
57-
"at the same time.\n",
58-
COLVARS_INPUT_ERROR);
55+
if (b_redefined_cutoff) {
56+
error_code |=
57+
cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
58+
COLVARS_INPUT_ERROR);
5959
}
6060

61-
b_anisotropic = true;
6261
// remove meaningless negative signs
6362
if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
6463
if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
6564
if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
65+
} else {
66+
if (b_redefined_cutoff) {
67+
r0_vec = {r0, r0, r0};
68+
}
6669
}
6770

6871
get_keyval(conf, "expNumer", en, en);
@@ -122,14 +125,14 @@ colvar::coordnum::~coordnum()
122125

123126
template<int flags> void colvar::coordnum::main_loop(bool **pairlist_elem)
124127
{
125-
const cvm::rvector inv_r0_vec(
126-
1.0 / ((flags & ef_anisotropic) ? r0_vec.x : r0),
127-
1.0 / ((flags & ef_anisotropic) ? r0_vec.y : r0),
128-
1.0 / ((flags & ef_anisotropic) ? r0_vec.z : r0));
129-
cvm::rvector const inv_r0sq_vec(
128+
const cvm::rvector inv_r0_vec{
129+
1.0 / r0_vec.x,
130+
1.0 / r0_vec.y,
131+
1.0 / r0_vec.z};
132+
cvm::rvector const inv_r0sq_vec{
130133
inv_r0_vec.x*inv_r0_vec.x,
131134
inv_r0_vec.y*inv_r0_vec.y,
132-
inv_r0_vec.z*inv_r0_vec.z);
135+
inv_r0_vec.z*inv_r0_vec.z};
133136
if (b_group2_center_only) {
134137
const cvm::atom_pos group2_com = group2->center_of_mass();
135138
cvm::rvector group2_com_grad(0, 0, 0);
@@ -187,41 +190,17 @@ template<int compute_flags> int colvar::coordnum::compute_coordnum()
187190

188191
bool *pairlist_elem = use_pairlist ? pairlist : NULL;
189192

190-
if (b_anisotropic) {
191-
192-
if (use_pairlist) {
193-
if (rebuild_pairlist) {
194-
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist |
195-
ef_rebuild_pairlist;
196-
main_loop<flags>(&pairlist_elem);
197-
} else {
198-
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist;
199-
main_loop<flags>(&pairlist_elem);
200-
}
201-
193+
if (use_pairlist) {
194+
if (rebuild_pairlist) {
195+
int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
196+
main_loop<flags>(&pairlist_elem);
202197
} else {
203-
204-
int const flags = compute_flags | ef_anisotropic;
205-
main_loop<flags>(NULL);
198+
int const flags = compute_flags | ef_use_pairlist;
199+
main_loop<flags>(&pairlist_elem);
206200
}
207-
208201
} else {
209-
210-
if (use_pairlist) {
211-
212-
if (rebuild_pairlist) {
213-
int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
214-
main_loop<flags>(&pairlist_elem);
215-
} else {
216-
int const flags = compute_flags | ef_use_pairlist;
217-
main_loop<flags>(&pairlist_elem);
218-
}
219-
220-
} else {
221-
222-
int const flags = compute_flags;
223-
main_loop<flags>(NULL);
224-
}
202+
int const flags = compute_flags;
203+
main_loop<flags>(NULL);
225204
}
226205

227206
return COLVARS_OK;
@@ -250,8 +229,8 @@ void colvar::coordnum::calc_gradients()
250229

251230
colvar::h_bond::h_bond()
252231
{
253-
colvarproxy *proxy = cvm::main()->proxy;
254-
r0 = proxy->angstrom_to_internal(3.3);
232+
cvm::real const r0 = cvm::main()->proxy->angstrom_to_internal(3.3);
233+
r0_vec = {r0, r0, r0};
255234
set_function_type("hBond");
256235
x.type(colvarvalue::type_scalar);
257236
init_scalar_boundaries(0.0, 1.0);
@@ -285,7 +264,11 @@ int colvar::h_bond::init(std::string const &conf)
285264
modify_atom.add_atom(cvm::atom_group::init_atom_from_proxy(p, d_num));
286265
}
287266

288-
get_keyval(conf, "cutoff", r0, r0);
267+
cvm::real r0 = r0_vec[0];
268+
bool const b_redefined_cutoff = get_keyval(conf, "cutoff", r0, r0);
269+
if (b_redefined_cutoff) {
270+
r0_vec = {r0, r0, r0};
271+
}
289272
get_keyval(conf, "expNumer", en, en);
290273
get_keyval(conf, "expDenom", ed, ed);
291274

@@ -309,7 +292,7 @@ colvar::h_bond::h_bond(cvm::atom_group::simple_atom const &acceptor,
309292
cvm::real r0_i, int en_i, int ed_i)
310293
: h_bond()
311294
{
312-
r0 = r0_i;
295+
r0_vec = {r0_i, r0_i, r0_i};
313296
en = en_i;
314297
ed = ed_i;
315298
register_atom_group(new cvm::atom_group);
@@ -322,22 +305,25 @@ colvar::h_bond::h_bond(cvm::atom_group::simple_atom const &acceptor,
322305
void colvar::h_bond::calc_value()
323306
{
324307
int const flags = coordnum::ef_null;
325-
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
326308
cvm::rvector G1, G2;
327309
const cvm::atom_pos A1{atom_groups[0]->pos_x(0),
328310
atom_groups[0]->pos_y(0),
329311
atom_groups[0]->pos_z(0)};
330312
const cvm::atom_pos A2{atom_groups[0]->pos_x(1),
331313
atom_groups[0]->pos_y(1),
332314
atom_groups[0]->pos_z(1)};
333-
const cvm::rvector inv_r0_vec(
334-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0),
335-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0),
336-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0));
337-
cvm::rvector const inv_r0sq_vec(
338-
inv_r0_vec.x*inv_r0_vec.x,
339-
inv_r0_vec.y*inv_r0_vec.y,
340-
inv_r0_vec.z*inv_r0_vec.z);
315+
316+
const cvm::rvector inv_r0_vec{
317+
1.0 / r0_vec.x,
318+
1.0 / r0_vec.y,
319+
1.0 / r0_vec.z
320+
};
321+
cvm::rvector const inv_r0sq_vec{
322+
inv_r0_vec.x * inv_r0_vec.x,
323+
inv_r0_vec.y * inv_r0_vec.y,
324+
inv_r0_vec.z * inv_r0_vec.z
325+
};
326+
341327
x.real_value =
342328
coordnum::switching_function<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
343329
atom_groups[0]->pos_x(0),
@@ -360,15 +346,16 @@ void colvar::h_bond::calc_value()
360346
void colvar::h_bond::calc_gradients()
361347
{
362348
int const flags = coordnum::ef_gradients;
363-
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
364-
const cvm::rvector inv_r0_vec(
365-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0),
366-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0),
367-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0));
368-
cvm::rvector const inv_r0sq_vec(
349+
const cvm::rvector inv_r0_vec{
350+
1.0 / r0_vec.x,
351+
1.0 / r0_vec.y,
352+
1.0 / r0_vec.z
353+
};
354+
cvm::rvector const inv_r0sq_vec{
369355
inv_r0_vec.x*inv_r0_vec.x,
370356
inv_r0_vec.y*inv_r0_vec.y,
371-
inv_r0_vec.z*inv_r0_vec.z);
357+
inv_r0_vec.z*inv_r0_vec.z
358+
};
372359
coordnum::switching_function<flags>(inv_r0_vec, inv_r0sq_vec, en, ed,
373360
atom_groups[0]->pos_x(0),
374361
atom_groups[0]->pos_y(0),
@@ -386,12 +373,12 @@ void colvar::h_bond::calc_gradients()
386373
}
387374

388375

389-
390376
colvar::selfcoordnum::selfcoordnum()
391377
{
392378
set_function_type("selfCoordNum");
393379
x.type(colvarvalue::type_scalar);
394-
r0 = cvm::main()->proxy->angstrom_to_internal(4.0);
380+
cvm::real const r0 = cvm::main()->proxy->angstrom_to_internal(4.0);
381+
r0_vec = {r0, r0, r0};
395382
}
396383

397384

@@ -405,7 +392,11 @@ int colvar::selfcoordnum::init(std::string const &conf)
405392
return error_code | COLVARS_INPUT_ERROR;
406393
}
407394

408-
get_keyval(conf, "cutoff", r0, r0);
395+
cvm::real r0 = r0_vec[0];
396+
bool const b_redefined_cutoff = get_keyval(conf, "cutoff", r0, r0);
397+
if (b_redefined_cutoff) {
398+
r0_vec = {r0, r0, r0};
399+
}
409400
get_keyval(conf, "expNumer", en, en);
410401
get_keyval(conf, "expDenom", ed, ed);
411402

@@ -452,8 +443,6 @@ colvar::selfcoordnum::~selfcoordnum()
452443

453444
template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
454445
{
455-
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
456-
457446
bool const use_pairlist = (pairlist != NULL);
458447
bool const rebuild_pairlist = (pairlist != NULL) &&
459448
(cvm::step_relative() % pairlist_freq == 0);
@@ -464,14 +453,16 @@ template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
464453

465454
// Always isotropic (TODO: enable the ellipsoid?)
466455
#define CALL_KERNEL(flags) do { \
467-
const cvm::rvector inv_r0_vec( \
468-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), \
469-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), \
470-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); \
471-
cvm::rvector const inv_r0sq_vec( \
472-
inv_r0_vec.x*inv_r0_vec.x, \
473-
inv_r0_vec.y*inv_r0_vec.y, \
474-
inv_r0_vec.z*inv_r0_vec.z); \
456+
const cvm::rvector inv_r0_vec{ \
457+
1.0 / r0_vec.x, \
458+
1.0 / r0_vec.y, \
459+
1.0 / r0_vec.z \
460+
}; \
461+
cvm::rvector const inv_r0sq_vec{ \
462+
inv_r0_vec.x*inv_r0_vec.x, \
463+
inv_r0_vec.y*inv_r0_vec.y, \
464+
inv_r0_vec.z*inv_r0_vec.z \
465+
}; \
475466
for (i = 0; i < n - 1; i++) { \
476467
for (j = i + 1; j < n; j++) { \
477468
x.real_value += \
@@ -535,11 +526,8 @@ colvar::groupcoordnum::groupcoordnum()
535526
set_function_type("groupCoord");
536527
x.type(colvarvalue::type_scalar);
537528
init_scalar_boundaries(0.0, 1.0);
538-
colvarproxy *proxy = cvm::main()->proxy;
539-
r0 = proxy->angstrom_to_internal(4.0);
540-
r0_vec = cvm::rvector(proxy->angstrom_to_internal(4.0),
541-
proxy->angstrom_to_internal(4.0),
542-
proxy->angstrom_to_internal(4.0));
529+
cvm::real const r0 = cvm::main()->proxy->angstrom_to_internal(4.0);
530+
r0_vec = {r0, r0, r0};
543531
}
544532

545533

@@ -552,19 +540,23 @@ int colvar::groupcoordnum::init(std::string const &conf)
552540
return cvm::error("Error: neither group can be a dummy atom\n", COLVARS_INPUT_ERROR);
553541
}
554542

555-
bool const b_scale = get_keyval(conf, "cutoff", r0, r0);
543+
cvm::real r0 = r0_vec[0];
544+
bool const b_redefined_cutoff = get_keyval(conf, "cutoff", r0, r0);
556545

557546
if (get_keyval(conf, "cutoff3", r0_vec, r0_vec)) {
558-
if (b_scale) {
547+
if (b_redefined_cutoff) {
559548
error_code |=
560549
cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
561550
COLVARS_INPUT_ERROR);
562551
}
563-
b_anisotropic = true;
564552
// remove meaningless negative signs
565553
if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
566554
if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
567555
if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
556+
} else {
557+
if (b_redefined_cutoff) {
558+
r0_vec = {r0, r0, r0};
559+
}
568560
}
569561

570562
get_keyval(conf, "expNumer", en, en);
@@ -591,32 +583,25 @@ void colvar::groupcoordnum::calc_value()
591583
{
592584
const cvm::atom_pos A1 = group1->center_of_mass();
593585
const cvm::atom_pos A2 = group2->center_of_mass();
594-
#define CALL_KERNEL(flags) do { \
595-
const cvm::rvector inv_r0_vec( \
596-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), \
597-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), \
598-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); \
599-
cvm::rvector const inv_r0sq_vec( \
600-
inv_r0_vec.x*inv_r0_vec.x, \
601-
inv_r0_vec.y*inv_r0_vec.y, \
602-
inv_r0_vec.z*inv_r0_vec.z); \
603-
cvm::rvector G1, G2; \
604-
const cvm::rvector r0sq_vec(r0_vec.x*r0_vec.x, \
605-
r0_vec.y*r0_vec.y, \
606-
r0_vec.z*r0_vec.z); \
586+
#define CALL_KERNEL(flags) do { \
587+
const cvm::rvector inv_r0_vec{ \
588+
1.0 / r0_vec.x, \
589+
1.0 / r0_vec.y, \
590+
1.0 / r0_vec.z \
591+
}; \
592+
cvm::rvector const inv_r0sq_vec{ \
593+
inv_r0_vec.x*inv_r0_vec.x, \
594+
inv_r0_vec.y*inv_r0_vec.y, \
595+
inv_r0_vec.z*inv_r0_vec.z \
596+
}; \
597+
cvm::rvector G1, G2; \
607598
x.real_value = coordnum::switching_function<flags>(inv_r0_vec, inv_r0sq_vec, en, ed, \
608599
A1.x, A1.y, A1.z, \
609600
A2.x, A2.y, A2.z, \
610601
G1.x, G1.y, G1.z, \
611602
G2.x, G2.y, G2.z, NULL, 0.0); \
612603
} while (0);
613-
if (b_anisotropic) {
614-
int const flags = coordnum::ef_anisotropic;
615-
CALL_KERNEL(flags);
616-
} else {
617-
int const flags = coordnum::ef_null;
618-
CALL_KERNEL(flags);
619-
}
604+
CALL_KERNEL(coordnum::ef_null);
620605
#undef CALL_KERNEL
621606
}
622607

@@ -627,27 +612,24 @@ void colvar::groupcoordnum::calc_gradients()
627612
const cvm::atom_pos A2 = group2->center_of_mass();
628613
cvm::rvector G1(0, 0, 0), G2(0, 0, 0);
629614
#define CALL_KERNEL(flags) do { \
630-
const cvm::rvector inv_r0_vec( \
631-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), \
632-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), \
633-
1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); \
634-
cvm::rvector const inv_r0sq_vec( \
635-
inv_r0_vec.x*inv_r0_vec.x, \
636-
inv_r0_vec.y*inv_r0_vec.y, \
637-
inv_r0_vec.z*inv_r0_vec.z); \
615+
const cvm::rvector inv_r0_vec{ \
616+
1.0 / r0_vec.x, \
617+
1.0 / r0_vec.y, \
618+
1.0 / r0_vec.z \
619+
}; \
620+
cvm::rvector const inv_r0sq_vec{ \
621+
inv_r0_vec.x*inv_r0_vec.x, \
622+
inv_r0_vec.y*inv_r0_vec.y, \
623+
inv_r0_vec.z*inv_r0_vec.z \
624+
}; \
638625
coordnum::switching_function<flags>(inv_r0_vec, inv_r0sq_vec, en, ed, \
639626
A1.x, A1.y, A1.z, \
640627
A2.x, A2.y, A2.z, \
641628
G1.x, G1.y, G1.z, \
642629
G2.x, G2.y, G2.z, NULL, 0.0); \
643630
} while (0);
644-
if (b_anisotropic) {
645-
int const flags = coordnum::ef_gradients | coordnum::ef_anisotropic;
646-
CALL_KERNEL(flags);
647-
} else {
648-
int const flags = coordnum::ef_gradients;
649-
CALL_KERNEL(flags);
650-
}
631+
int const flags = coordnum::ef_gradients;
632+
CALL_KERNEL(flags);
651633
group1->set_weighted_gradient(G1);
652634
group2->set_weighted_gradient(G2);
653635
}

0 commit comments

Comments
 (0)