Move QMMM source to C++
authorErik Lindahl <erik@kth.se>
Tue, 7 Jul 2015 20:50:41 +0000 (22:50 +0200)
committerErik Lindahl <erik@kth.se>
Tue, 7 Jul 2015 22:45:55 +0000 (00:45 +0200)
Tested by enabling the four different QM/MM interfaces,
but not linking to the actual programs since those calls
are still unaltered.

Change-Id: I86163a0db56003687e6226a78169141af98940b7

src/gromacs/mdlib/qm_gamess.cpp [moved from src/gromacs/mdlib/qm_gamess.c with 95% similarity]
src/gromacs/mdlib/qm_gaussian.cpp [moved from src/gromacs/mdlib/qm_gaussian.c with 97% similarity]
src/gromacs/mdlib/qm_mopac.cpp [moved from src/gromacs/mdlib/qm_mopac.c with 91% similarity]
src/gromacs/mdlib/qm_orca.cpp [moved from src/gromacs/mdlib/qm_orca.c with 98% similarity]
src/gromacs/mdlib/qmmm.cpp [moved from src/gromacs/mdlib/qmmm.c with 95% similarity]

similarity index 95%
rename from src/gromacs/mdlib/qm_gamess.c
rename to src/gromacs/mdlib/qm_gamess.cpp
index d913f7f48925abd33818aa8e90317dac5d590b2a..b7c729b29d34b8b7c9c896eda8e5ca5f9970e1ac 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -55,6 +55,7 @@
 #include "gromacs/legacyheaders/qmmm.h"
 #include "gromacs/legacyheaders/txtdump.h"
 #include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/utility/fatalerror.h"
 
 
 void
-F77_FUNC(inigms, IMIGMS) (void);
+    F77_FUNC(inigms, IMIGMS) (void);
 
 void
-F77_FUNC(endgms, ENDGMS) (void);
+    F77_FUNC(endgms, ENDGMS) (void);
 
 void
-F77_FUNC(grads, GRADS) (int *nrqmat, real *qmcrd, int *nrmmat, real *mmchrg,
-                        real *mmcrd, real *qmgrad, real *mmgrad, real *energy);
+    F77_FUNC(grads, GRADS) (int *nrqmat, real *qmcrd, int *nrmmat, real *mmchrg,
+                            real *mmcrd, real *qmgrad, real *mmgrad, real *energy);
 
 
 
@@ -88,7 +89,7 @@ void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
      * dynamics simulations. 7-6-2002 (London)
      */
     int
-        i, j, rank;
+        i, j;
     FILE
        *out;
     char
@@ -228,7 +229,7 @@ void init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
     }
 }
 
-real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
+real call_gamess(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
                  rvec f[], rvec fshift[])
 {
     /* do the actual QMMM calculation using GAMESS-UK. In this
@@ -237,7 +238,7 @@ real call_gamess(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
      * gradient routines linked directly
      */
     int
-        i, j, rank;
+        i, j;
     real
         QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy;
     t_QMMMrec
similarity index 97%
rename from src/gromacs/mdlib/qm_gaussian.c
rename to src/gromacs/mdlib/qm_gaussian.cpp
index d90232c7e86a55fae13a1b0ca5d9a31f4c1fb51f..6f519c26b3f64f322f1d876215fa04b548e8f9ea 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #include "gromacs/legacyheaders/typedefs.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
+#include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
-
 /* TODO: this should be made thread-safe */
 
 /* Gaussian interface routines */
 
-void init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
+void init_gaussian(t_QMrec *qm)
 {
-    FILE
-       *rffile = NULL, *out = NULL;
+    FILE *out = NULL;
     ivec
-        basissets[eQMbasisNR] = {{0, 3, 0},
+          basissets[eQMbasisNR] = {{0, 3, 0},
                                  {0, 3, 0}, /*added for double sto-3g entry in names.c*/
                                  {5, 0, 0},
                                  {5, 0, 1},
@@ -81,9 +80,9 @@ void init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
                                  {1, 6, 11},
                                  {4, 6, 0}};
     char
-       *buf = NULL;
+         *buf = NULL;
     int
-        i;
+          i;
 
     /* using the ivec above to convert the basis read form the mdp file
      * in a human readable format into some numbers for the gaussian
@@ -695,8 +694,7 @@ void write_gaussian_input(int step, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
 
 }  /* write_gaussian_input */
 
-real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], int step,
-                          t_QMrec *qm, t_MMrec *mm)
+real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], t_QMrec *qm, t_MMrec *mm)
 {
     int
         i, j, atnum;
@@ -801,8 +799,7 @@ real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], int step,
     return(QMener);
 }
 
-real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step,
-                             gmx_bool swapped, t_QMrec *qm, t_MMrec *mm)
+real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step, t_QMrec *qm, t_MMrec *mm)
 {
     int
         i;
@@ -1023,8 +1020,7 @@ void do_gaussian(int step, char *exe)
     }
 }
 
-real call_gaussian(t_commrec *cr,  t_forcerec *fr,
-                   t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
+real call_gaussian(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 {
     /* normal gaussian jobs */
     static int
@@ -1045,7 +1041,7 @@ real call_gaussian(t_commrec *cr,  t_forcerec *fr,
 
     write_gaussian_input(step, fr, qm, mm);
     do_gaussian(step, exe);
-    QMener = read_gaussian_output(QMgrad, MMgrad, step, qm, mm);
+    QMener = read_gaussian_output(QMgrad, MMgrad, qm, mm);
     /* put the QMMM forces in the force array and to the fshift
      */
     for (i = 0; i < qm->nrQMatoms; i++)
@@ -1071,8 +1067,7 @@ real call_gaussian(t_commrec *cr,  t_forcerec *fr,
 
 } /* call_gaussian */
 
-real call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
-                      rvec f[], rvec fshift[])
+real call_gaussian_SH(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 {
     /* a gaussian call routine intended for doing diabatic surface
      * "sliding". See the manual for the theoretical background of this
@@ -1128,7 +1123,7 @@ real call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
     write_gaussian_SH_input(step, swapped, fr, qm, mm);
 
     do_gaussian(step, exe);
-    QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, swapped, qm, mm);
+    QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
 
     /* check for a surface hop. Only possible if we were already state
      * averaging.
@@ -1149,7 +1144,7 @@ real call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
         {
             write_gaussian_SH_input(step, swapped, fr, qm, mm);
             do_gaussian(step, exe);
-            QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, swapped, qm, mm);
+            QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
         }
     }
     /* add the QMMM forces to the gmx force array and fshift
similarity index 91%
rename from src/gromacs/mdlib/qm_mopac.c
rename to src/gromacs/mdlib/qm_mopac.cpp
index 6eb689cd3ea5d5f4002da26408af314e9154625c..7496965bf09be4db5fcce49db1edc2266f4ff30f 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 
 /* mopac interface routines */
 void
-F77_FUNC(domldt, DOMLDT) (int *nrqmat, int labels[], char keywords[]);
+    F77_FUNC(domldt, DOMLDT) (int *nrqmat, int labels[], char keywords[]);
 
 void
-F77_FUNC(domop, DOMOP) (int *nrqmat, double qmcrd[], int *nrmmat,
-                        double mmchrg[], double mmcrd[], double qmgrad[],
-                        double mmgrad[], double *energy, double qmcharges[]);
+    F77_FUNC(domop, DOMOP) (int *nrqmat, double qmcrd[], int *nrmmat,
+                            double mmchrg[], double mmcrd[], double qmgrad[],
+                            double mmgrad[], double *energy, double qmcharges[]);
 
 
 
-void init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
+void init_mopac(t_QMrec *qm)
 {
     /* initializes the mopac routines ans sets up the semiempirical
      * computation by calling moldat(). The inline mopac routines can
@@ -104,8 +104,7 @@ void init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm)
 
 } /* init_mopac */
 
-real call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
-                rvec f[], rvec fshift[])
+real call_mopac(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 {
     /* do the actual QMMM calculation using directly linked mopac subroutines
      */
@@ -167,8 +166,7 @@ real call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
     return (QMener);
 }
 
-real call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm, t_MMrec *mm,
-                   rvec f[], rvec fshift[])
+real call_mopac_SH(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[])
 {
     /* do the actual SH QMMM calculation using directly linked mopac
        subroutines */
similarity index 98%
rename from src/gromacs/mdlib/qm_orca.c
rename to src/gromacs/mdlib/qm_orca.cpp
index 11ddfd325b901c82b2b25aa211aa228f59c9546b..2f121806fd3ef98fec128f7ad797d0208067ca73 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -284,7 +284,7 @@ real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr,
     int
         i, j, atnum;
     char
-        buf[300], tmp[300], orca_xyzFilename[300], orca_pcgradFilename[300], orca_engradFilename[300];
+        buf[300], orca_xyzFilename[300], orca_pcgradFilename[300], orca_engradFilename[300];
     real
         QMener;
     FILE
@@ -320,8 +320,8 @@ real read_orca_output(rvec QMgrad[], rvec MMgrad[], t_forcerec *fr,
                 gmx_fatal(FARGS, "Unexpected end of ORCA output");
             }
 #ifdef GMX_DOUBLE
-            sscanf(buf, "%s%lf%lf%lf\n",
-                   tmp,
+            sscanf(buf, "%d%lf%lf%lf\n",
+                   &atnum,
                    &qm->xQM[i][XX],
                    &qm->xQM[i][YY],
                    &qm->xQM[i][ZZ]);
similarity index 95%
rename from src/gromacs/mdlib/qmmm.c
rename to src/gromacs/mdlib/qmmm.cpp
index bb83405a5548e083d481bcf88616aed3c422a2b5..cf7f7b344082ac850676c7b20f2412245ac43d3d 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #include <stdlib.h>
 #include <string.h>
 
+#include <cmath>
+
+#include <algorithm>
+
 #include "gromacs/fileio/confio.h"
 #include "gromacs/legacyheaders/force.h"
 #include "gromacs/legacyheaders/macros.h"
@@ -73,36 +77,32 @@ void
 init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
 
 real
-call_gamess(t_commrec *cr, t_forcerec *fr,
+call_gamess(t_forcerec *fr,
             t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
 
 #elif defined GMX_QMMM_MOPAC
 /* MOPAC interface */
 
 void
-init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
+init_mopac(t_QMrec *qm);
 
 real
-call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
-           t_MMrec *mm, rvec f[], rvec fshift[]);
+call_mopac(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
 
 real
-call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
-              t_MMrec *mm, rvec f[], rvec fshift[]);
+call_mopac_SH(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
 
 #elif defined GMX_QMMM_GAUSSIAN
 /* GAUSSIAN interface */
 
 void
-init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
+init_gaussian(t_QMrec *qm);
 
 real
-call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
-                 t_MMrec *mm, rvec f[], rvec fshift[]);
+call_gaussian_SH(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
 
 real
-call_gaussian(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
-              t_MMrec *mm, rvec f[], rvec fshift[]);
+call_gaussian(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
 
 #elif defined GMX_QMMM_ORCA
 /* ORCA interface */
@@ -165,11 +165,11 @@ real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec
 #ifdef GMX_QMMM_MOPAC
         if (qm->bSH)
         {
-            QMener = call_mopac_SH(cr, fr, qm, mm, f, fshift);
+            QMener = call_mopac_SH(qm, mm, f, fshift);
         }
         else
         {
-            QMener = call_mopac(cr, fr, qm, mm, f, fshift);
+            QMener = call_mopac(qm, mm, f, fshift);
         }
 #else
         gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
@@ -181,7 +181,7 @@ real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec
         if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
         {
 #ifdef GMX_QMMM_GAUSSIAN
-            QMener = call_gaussian_SH(cr, fr, qm, mm, f, fshift);
+            QMener = call_gaussian_SH(fr, qm, mm, f, fshift);
 #else
             gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian.");
 #endif
@@ -189,9 +189,9 @@ real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec
         else
         {
 #ifdef GMX_QMMM_GAMESS
-            QMener = call_gamess(cr, fr, qm, mm, f, fshift);
+            QMener = call_gamess(fr, qm, mm, f, fshift);
 #elif defined GMX_QMMM_GAUSSIAN
-            QMener = call_gaussian(cr, fr, qm, mm, f, fshift);
+            QMener = call_gaussian(fr, qm, mm, f, fshift);
 #elif defined GMX_QMMM_ORCA
             QMener = call_orca(fr, qm, mm, f, fshift);
 #else
@@ -210,7 +210,7 @@ void init_QMroutine(t_commrec gmx_unused *cr, t_QMrec gmx_unused *qm, t_MMrec gm
     {
 #ifdef GMX_QMMM_MOPAC
         /* do a semi-empiprical calculation */
-        init_mopac(cr, qm, mm);
+        init_mopac(qm);
 #else
         gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 #endif
@@ -221,7 +221,7 @@ void init_QMroutine(t_commrec gmx_unused *cr, t_QMrec gmx_unused *qm, t_MMrec gm
 #ifdef GMX_QMMM_GAMESS
         init_gamess(cr, qm, mm);
 #elif defined GMX_QMMM_GAUSSIAN
-        init_gaussian(cr, qm, mm);
+        init_gaussian(qm);
 #elif defined GMX_QMMM_ORCA
         init_orca(qm);
 #else
@@ -262,7 +262,7 @@ static void punch_QMMM_excl(t_QMrec *qm, t_MMrec *mm, t_blocka *excls)
     FILE
        *out = NULL;
     int
-        i, j, k, nrexcl = 0, *excluded = NULL, max = 0;
+        i, j, k, nrexcl = 0, *excluded = NULL, max_excl = 0;
 
 
     out = fopen("QMMMexcl.dat", "w");
@@ -280,10 +280,10 @@ static void punch_QMMM_excl(t_QMrec *qm, t_MMrec *mm, t_blocka *excls)
             {
                 if (mm->indexMM[k] == excls->a[j]) /* the excluded MM atom */
                 {
-                    if (nrexcl >= max)
+                    if (nrexcl >= max_excl)
                     {
-                        max += 1000;
-                        srenew(excluded, max);
+                        max_excl += 1000;
+                        srenew(excluded, max_excl);
                     }
                     excluded[nrexcl++] = k;
                     continue;
@@ -485,8 +485,8 @@ void init_QMMMrec(t_commrec  *cr,
     t_ilist                 *ilist_mol;
     gmx_mtop_atomlookup_t    alook;
 
-    c6au  = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 6));
-    c12au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 12));
+    c6au  = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 6));
+    c12au = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 12));
     /* issue a fatal if the user wants to run with more than one node */
     if (PAR(cr))
     {
@@ -753,7 +753,7 @@ void init_QMMMrec(t_commrec  *cr,
         {
 #ifdef GMX_QMMM_MOPAC
             /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
-            init_mopac(cr, qr->qm[0], qr->mm);
+            init_mopac(qr->qm[0]);
 #else
             gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 #endif
@@ -764,7 +764,7 @@ void init_QMMMrec(t_commrec  *cr,
 #ifdef GMX_QMMM_GAMESS
             init_gamess(cr, qr->qm[0], qr->mm);
 #elif defined GMX_QMMM_GAUSSIAN
-            init_gaussian(cr, qr->qm[0], qr->mm);
+            init_gaussian(qr->qm[0]);
 #elif defined GMX_QMMM_ORCA
             init_orca(qr->qm[0]);
 #else
@@ -797,8 +797,6 @@ void update_QMMMrec(t_commrec      *cr,
         QMMMlist;
     rvec
         dx, crd;
-    int
-       *MMatoms;
     t_QMrec
        *qm;
     t_MMrec
@@ -810,8 +808,8 @@ void update_QMMMrec(t_commrec      *cr,
     real
         c12au, c6au;
 
-    c6au  = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 6));
-    c12au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 12));
+    c6au  = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 6));
+    c12au = (HARTREE2KJ*AVOGADRO*std::pow(BOHR2NM, 12));
 
     /* every cpu has this array. On every processor we fill this array
      * with 1's and 0's. 1's indicate the atoms is a QM atom on the
@@ -873,7 +871,7 @@ void update_QMMMrec(t_commrec      *cr,
                 crd[0] = IS2X(QMMMlist.shift[i]) + IS2X(qm_i_particles[i].shift);
                 crd[1] = IS2Y(QMMMlist.shift[i]) + IS2Y(qm_i_particles[i].shift);
                 crd[2] = IS2Z(QMMMlist.shift[i]) + IS2Z(qm_i_particles[i].shift);
-                is     = XYZ2IS(crd[0], crd[1], crd[2]);
+                is     = static_cast<int>(XYZ2IS(crd[0], crd[1], crd[2]));
                 for (j = QMMMlist.jindex[i];
                      j < QMMMlist.jindex[i+1];
                      j++)
@@ -897,9 +895,13 @@ void update_QMMMrec(t_commrec      *cr,
             qsort(qm_i_particles, QMMMlist.nri,
                   (size_t)sizeof(qm_i_particles[0]),
                   struct_comp);
-            qsort(mm_j_particles, mm_nr,
-                  (size_t)sizeof(mm_j_particles[0]),
-                  struct_comp);
+            /* The mm_j_particles argument to qsort is not allowed to be NULL */
+            if (mm_nr > 0)
+            {
+                qsort(mm_j_particles, mm_nr,
+                      (size_t)sizeof(mm_j_particles[0]),
+                      struct_comp);
+            }
             /* remove multiples in the QM shift array, since in init_QMMM() we
              * went through the atom numbers from 0 to md.nr, the order sorted
              * here matches the one of QMindex already.