From 000c6b9e5d94d9e2c5efb592cfb038d9775e0c21 Mon Sep 17 00:00:00 2001 From: Erik Lindahl Date: Fri, 6 Jun 2014 09:48:45 +0200 Subject: [PATCH] Process negative sigma correctly with combrule 2 or 3 The initial combination rule code (for sigma/epsilon) did not take negative sigma rules into account, which caused segfaults instead of the values reaching the code in convparm.c that handles sigma<0 during conversion to c6/c12 (where it is used to signal c6=0). Fixes #1391. Change-Id: I437f06d67c5ecfb58d236590288ad122bcdf2df0 --- src/kernel/toppush.c | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/src/kernel/toppush.c b/src/kernel/toppush.c index d2f142a025..c5c98ed589 100644 --- a/src/kernel/toppush.c +++ b/src/kernel/toppush.c @@ -95,7 +95,7 @@ void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atyp break; case eCOMB_ARITHMETIC: - /* c0 and c1 are epsilon and sigma */ + /* c0 and c1 are sigma and epsilon */ for (i = k = 0; (i < nr); i++) { for (j = 0; (j < nr); j++, k++) @@ -104,14 +104,21 @@ void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atyp cj0 = get_atomtype_nbparam(j, 0, atype); ci1 = get_atomtype_nbparam(i, 1, atype); cj1 = get_atomtype_nbparam(j, 1, atype); - plist->param[k].c[0] = (ci0+cj0)*0.5; + plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5; + /* Negative sigma signals that c6 should be set to zero later, + * so we need to propagate that through the combination rules. + */ + if (ci0 < 0 || cj0 < 0) + { + plist->param[k].c[0] *= -1; + } plist->param[k].c[1] = sqrt(ci1*cj1); - } + } } break; case eCOMB_GEOM_SIG_EPS: - /* c0 and c1 are epsilon and sigma */ + /* c0 and c1 are sigma and epsilon */ for (i = k = 0; (i < nr); i++) { for (j = 0; (j < nr); j++, k++) @@ -120,7 +127,14 @@ void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atyp cj0 = get_atomtype_nbparam(j, 0, atype); ci1 = get_atomtype_nbparam(i, 1, atype); cj1 = get_atomtype_nbparam(j, 1, atype); - plist->param[k].c[0] = sqrt(ci0*cj0); + plist->param[k].c[0] = sqrt(fabs(ci0*cj0)); + /* Negative sigma signals that c6 should be set to zero later, + * so we need to propagate that through the combination rules. + */ + if (ci0 < 0 || cj0 < 0) + { + plist->param[k].c[0] *= -1; + } plist->param[k].c[1] = sqrt(ci1*cj1); } } -- 2.22.0