Skip to content

Commit 1cfa768

Browse files
HanatoKgiacomofiorin
HanatoK
authored and
giacomofiorin
committed
refactor: remove some redundant code
1 parent 69b0a4f commit 1cfa768

File tree

2 files changed

+19
-115
lines changed

2 files changed

+19
-115
lines changed

‎src/colvarcomp_coordnums.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -176,7 +176,7 @@ void colvar::coordnum::compute_tolerance_l2_max()
176176
size_t i;
177177
// Find the value of l2 such that F(l2) = 0 using the Newton method
178178
for (i = 0; i < num_iters_max; i++) {
179-
F = switching_function<ef_use_pairlist | ef_gradients>(l2, dFdl2, en, ed, tolerance);
179+
F = switching_function<ef_use_pairlist | ef_gradients, 0, 0>(l2, dFdl2, en, ed, tolerance);
180180
if ((std::fabs(F) < result_tol) || (std::fabs(dFdl2) < dF_tol)) {
181181
break;
182182
}
@@ -532,7 +532,7 @@ template <int flags, int n, int m> inline void colvar::selfcoordnum::selfcoordnu
532532
diff = cvm::main()->proxy->position_distance(cvm::rvector{x1, y1, z1}, cvm::rvector{x2, y2, z2});
533533
}
534534
partial = compute_pair_coordnum<flags, n, m>(
535-
inv_r0_vec, inv_r0sq_vec, diff,
535+
inv_r0_vec, inv_r0sq_vec, diff, n, m,
536536
g1.x, g1.y, g1.z, gx2, gy2, gz2,
537537
tolerance, tolerance_l2_max);
538538
}

‎src/colvarcomp_coordnums.h‎

Lines changed: 17 additions & 113 deletions
Original file line numberDiff line numberDiff line change
@@ -43,9 +43,6 @@ class colvar::coordnum : public colvar::cvc {
4343
/// @param en Numerator exponent
4444
/// @param ed Denominator exponent
4545
/// @param pairlist_tol Pairlist tolerance
46-
template <int flags>
47-
static cvm::real switching_function(cvm::real const &l2, cvm::real &dFdl2, int en, int ed,
48-
cvm::real pairlist_tol);
4946

5047
/// Main kernel for the coordination number
5148
template <int flags>
@@ -57,14 +54,16 @@ class colvar::coordnum : public colvar::cvc {
5754
cvm::real &g2x, cvm::real &g2y, cvm::real &g2z,
5855
cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max);
5956

60-
template <int flags, int en, int ed>
57+
template <int flags, int t_en, int t_ed>
6158
static cvm::real switching_function(cvm::real const &l2, cvm::real &dFdl2,
59+
int en, int ed,
6260
cvm::real pairlist_tol);
6361

64-
template <int flags, int en, int ed>
62+
template <int flags, int t_en, int t_ed>
6563
static cvm::real compute_pair_coordnum(cvm::rvector const &inv_r0_vec,
6664
cvm::rvector const &inv_r0sq_vec,
6765
const cvm::rvector& diff,
66+
int en, int ed,
6867
cvm::real &g1x, cvm::real &g1y, cvm::real &g1z,
6968
cvm::real &g2x, cvm::real &g2y, cvm::real &g2z,
7069
cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max);
@@ -182,74 +181,15 @@ class colvar::h_bond : public colvar::cvc {
182181
};
183182

184183

185-
template <int flags>
186-
inline cvm::real colvar::coordnum::switching_function(cvm::real const &l2, cvm::real &dFdl2,
187-
int en, int ed,
188-
cvm::real pairlist_tol)
189-
{
190-
// Assume en and ed are even integers, and avoid sqrt in the following
191-
int const en2 = en/2;
192-
int const ed2 = ed/2;
193-
194-
cvm::real const xn = cvm::integer_power(l2, en2);
195-
cvm::real const xd = cvm::integer_power(l2, ed2);
196-
cvm::real const eps_l2 = 1.0e-7;
197-
cvm::real const h = l2 - 1.0;
198-
cvm::real const en2_r = (cvm::real) en2;
199-
cvm::real const ed2_r = (cvm::real) ed2;
200-
cvm::real func_no_pairlist;
201-
202-
if (std::abs(h) < eps_l2) {
203-
// Order-2 Taylor expansion: c0 + c1*h + c2*h^2
204-
cvm::real const c0 = en2_r / ed2_r;
205-
cvm::real const c1 = (en2_r * (en2_r - ed2_r)) / (2.0 * ed2_r);
206-
cvm::real const c2 = (en2_r * (en2_r - ed2_r) * (2.0 * en2_r - ed2_r - 3.0)) / (12.0 * ed2_r);
207-
func_no_pairlist = c0 + h * (c1 + h * c2);
208-
} else {
209-
func_no_pairlist = (1.0 - xn) / (1.0 - xd);
210-
}
211-
212-
cvm::real func, inv_one_pairlist_tol;
213-
if (flags & ef_use_pairlist) {
214-
inv_one_pairlist_tol = 1 / (1.0-pairlist_tol);
215-
func = (func_no_pairlist - pairlist_tol) * inv_one_pairlist_tol;
216-
} else {
217-
func = func_no_pairlist;
218-
}
219-
220-
// If the value is too small and we are correcting for the tolerance, the result is negative
221-
// and we need to exclude it rather than let it contribute to the sum or the gradients.
222-
if (func < 0)
223-
return 0;
224-
225-
if (flags & ef_gradients) {
226-
// Logarithmic derivative: 1st-order Taylor expansion around l2 = 1
227-
cvm::real log_deriv;
228-
if (std::abs(h) < eps_l2) {
229-
cvm::real const g0 = 0.5 * (en2_r - ed2_r);
230-
cvm::real const g1 = ((en2_r - ed2_r) * (en2_r + ed2_r - 6.0)) / 12.0;
231-
log_deriv = g0 + h * g1;
232-
} else {
233-
log_deriv = (ed2_r * xd / ((1.0 - xd) * l2)) - (en2_r * xn / ((1.0 - xn) * l2));
234-
}
235-
dFdl2 = (flags & ef_use_pairlist) ?
236-
func_no_pairlist * inv_one_pairlist_tol * log_deriv :
237-
func * log_deriv;
238-
}
239-
240-
return func;
241-
}
242-
243-
template <int flags, int en, int ed>
184+
template <int flags, int t_en, int t_ed>
244185
inline cvm::real colvar::coordnum::switching_function(
245-
cvm::real const &l2, cvm::real &dFdl2, cvm::real pairlist_tol)
186+
cvm::real const &l2, cvm::real &dFdl2, int en, int ed, cvm::real pairlist_tol)
246187
{
247188
// Assume en and ed are even integers, and avoid sqrt in the following
248-
const int constexpr en2 = en/2;
249-
const int constexpr ed2 = ed/2;
250-
251-
cvm::real const xn = cvm::integer_power<en2>(l2);
252-
cvm::real const xd = cvm::integer_power<ed2>(l2);
189+
const int en2 = t_en != 0 ? t_en / 2 : en / 2;
190+
const int ed2 = t_ed != 0 ? t_ed / 2 : ed / 2;
191+
cvm::real const xn = t_en != 0 ? cvm::integer_power<t_en / 2>(l2) : cvm::integer_power(l2, en2);
192+
cvm::real const xd = t_ed != 0 ? cvm::integer_power<t_ed / 2>(l2) : cvm::integer_power(l2, ed2);
253193
cvm::real const eps_l2 = 1.0e-7;
254194
cvm::real const h = l2 - 1.0;
255195
cvm::real const en2_r = (cvm::real) en2;
@@ -323,47 +263,17 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
323263
? cvm::main()->proxy->position_distance_internal(pos1, pos2)
324264
: cvm::main()->proxy->position_distance(pos1, pos2);
325265

326-
cvm::rvector const scal_diff(diff.x * inv_r0_vec.x,
327-
diff.y * inv_r0_vec.y,
328-
diff.z * inv_r0_vec.z);
329-
cvm::real const l2 = scal_diff.norm2();
330-
if (flags & ef_use_pairlist) {
331-
if (l2 > pairlist_tol_l2_max) {
332-
// Exit if the distance is such that F(l2) < pairlist_tol
333-
return 0.0;
334-
}
335-
}
336-
337-
cvm::real dFdl2 = 0.0;
338-
cvm::real F = switching_function<flags>(l2, dFdl2, en, ed, pairlist_tol);
339-
340-
if ((flags & ef_gradients) && (F > 0.0)) {
341-
cvm::rvector const dl2dx((2.0 * inv_r0sq_vec.x) * diff.x,
342-
(2.0 * inv_r0sq_vec.y) * diff.y,
343-
(2.0 * inv_r0sq_vec.z) * diff.z);
344-
345-
const cvm::rvector G = dFdl2*dl2dx;
346-
g1x += -1.0*G.x;
347-
g1y += -1.0*G.y;
348-
g1z += -1.0*G.z;
349-
g2x += G.x;
350-
g2y += G.y;
351-
g2z += G.z;
352-
}
353-
354-
return F;
266+
return compute_pair_coordnum<flags, 0, 0>(
267+
inv_r0_vec, inv_r0sq_vec, diff, en, ed,
268+
g1x, g1y, g1z, g2x, g2y, g2z,
269+
pairlist_tol, pairlist_tol_l2_max);
355270
}
356271

357-
template<int flags, int en, int ed>
272+
template<int flags, int t_en, int t_ed>
358273
inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv_r0_vec,
359274
cvm::rvector const &inv_r0sq_vec,
360275
const cvm::rvector& diff,
361-
// const cvm::real a1x,
362-
// const cvm::real a1y,
363-
// const cvm::real a1z,
364-
// const cvm::real a2x,
365-
// const cvm::real a2y,
366-
// const cvm::real a2z,
276+
int en, int ed,
367277
cvm::real& g1x,
368278
cvm::real& g1y,
369279
cvm::real& g1z,
@@ -373,12 +283,6 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
373283
cvm::real pairlist_tol,
374284
cvm::real pairlist_tol_l2_max)
375285
{
376-
// const cvm::atom_pos pos1{a1x, a1y, a1z};
377-
// const cvm::atom_pos pos2{a2x, a2y, a2z};
378-
// cvm::rvector const diff = (flags & ef_use_internal_pbc)
379-
// ? cvm::main()->proxy->position_distance_internal(pos1, pos2)
380-
// : cvm::main()->proxy->position_distance(pos1, pos2);
381-
382286
cvm::rvector const scal_diff(diff.x * inv_r0_vec.x,
383287
diff.y * inv_r0_vec.y,
384288
diff.z * inv_r0_vec.z);
@@ -391,7 +295,7 @@ inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv
391295
}
392296

393297
cvm::real dFdl2 = 0.0;
394-
cvm::real F = switching_function<flags, en, ed>(l2, dFdl2, pairlist_tol);
298+
const cvm::real F = switching_function<flags, t_en, t_ed>(l2, dFdl2, en, ed, pairlist_tol);
395299

396300
if ((flags & ef_gradients) && (F > 0.0)) {
397301
cvm::rvector const dl2dx((2.0 * inv_r0sq_vec.x) * diff.x,

0 commit comments

Comments
 (0)