Add missing Ewald charge term with PME on GPU

With PME on GPU and a system with net charge, the net charge energy
contribution was missing.

This is a backport of 55231153a622f83f053da0b24cd0f2e8c9752962
(MR !3350).

Fixes #4668

docs/release-notes/2021/2021.7.rst

@@ -16,6 +16,18 @@ in the :ref:`release-notes`.
 Fixes where mdrun could behave incorrectly
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
+Add missing net charge term when running PME on a GPU
+"""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When PME was running on a GPU, the term due to a net charge of the system was missing.
+In normal runs this only changed the potential energy by a constant, which is usually
+not relevant. In free-energy calculations where the net charge of the systen changes,
+the would lead to incorrect dV/dlambda and Delta lambda values (but one should anyhow
+avoid changing the net charge of a system with free-energy calculations).
+
+:issue:`4668`
+
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

src/gromacs/mdlib/force.cpp

@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2022, 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.
@@ -122,7 +122,8 @@ void calculateLongRangeNonbondeds(t_forcerec*                   fr,
     /* Do long-range electrostatics and/or LJ-PME
      * and compute PME surface terms when necessary.
      */
-    if ((computePmeOnCpu || fr->ic->eeltype == eelEWALD || haveEwaldSurfaceTerm)
+    if ((computePmeOnCpu || fr->ic->eeltype == eelEWALD || haveEwaldSurfaceTerm || fr->qsum[0] != 0
+         || fr->qsum[1] != 0)
         && stepWork.computeNonbondedForces)
     {
         int  status = 0;