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

diff --git a/src/kernel/toppush.c b/src/kernel/toppush.c
index d2f142a02511edab1dba6a179f3f44d29c925614..c5c98ed58943ccb5f882f7c58d2f74403fcf6af3 100644 (file)
--- 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);
                         }
                     }