src/gromacs/mdlib/tgroup.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #include "gmxpre.h"
  40
  41 #include "tgroup.h"
  42
  43 #include "gromacs/math/vec.h"
  44 #include "gromacs/mdlib/coupling.h"
  45 #include "gromacs/mdtypes/group.h"
  46 #include "gromacs/mdtypes/inputrec.h"
  47
  48
  49 real sum_ekin(const t_grpopts* opts, gmx_ekindata_t* ekind, real* dekindlambda, bool bEkinAveVel, bool bScaleEkin)
  50 {
  51     int           i, j, m, ngtc;
  52     real          T;
  53     t_grp_tcstat* tcstat;
  54     real          nrdf, nd, *ndf;
  55
  56     ngtc = opts->ngtc;
  57     ndf  = opts->nrdf;
  58
  59     T    = 0;
  60     nrdf = 0;
  61
  62     clear_mat(ekind->ekin);
  63
  64     for (i = 0; (i < ngtc); i++)
  65     {
  66
  67         nd     = ndf[i];
  68         tcstat = &ekind->tcstat[i];
  69         /* Sometimes a group does not have degrees of freedom, e.g.
  70          * when it consists of shells and virtual sites, then we just
  71          * set the temperatue to 0 and also neglect the kinetic
  72          * energy, which should be  zero anyway.
  73          */
  74
  75         if (nd > 0)
  76         {
  77             if (bEkinAveVel)
  78             {
  79                 if (!bScaleEkin)
  80                 {
  81                     /* in this case, kinetic energy is from the current velocities already */
  82                     msmul(tcstat->ekinf, tcstat->ekinscalef_nhc, tcstat->ekinf);
  83                 }
  84             }
  85             else
  86             {
  87                 /* Calculate the full step Ekin as the average of the half steps */
  88                 for (j = 0; (j < DIM); j++)
  89                 {
  90                     for (m = 0; (m < DIM); m++)
  91                     {
  92                         tcstat->ekinf[j][m] = 0.5
  93                                               * (tcstat->ekinh[j][m] * tcstat->ekinscaleh_nhc
  94                                                  + tcstat->ekinh_old[j][m]);
  95                     }
  96                 }
  97             }
  98             m_add(tcstat->ekinf, ekind->ekin, ekind->ekin);
  99
 100             tcstat->Th = calc_temp(trace(tcstat->ekinh), nd);
 101             tcstat->T  = calc_temp(trace(tcstat->ekinf), nd);
 102
 103             /* after the scaling factors have been multiplied in, we can remove them */
 104             if (bEkinAveVel)
 105             {
 106                 tcstat->ekinscalef_nhc = 1.0;
 107             }
 108             else
 109             {
 110                 tcstat->ekinscaleh_nhc = 1.0;
 111             }
 112         }
 113         else
 114         {
 115             tcstat->T  = 0;
 116             tcstat->Th = 0;
 117         }
 118         T += nd * tcstat->T;
 119         nrdf += nd;
 120     }
 121     if (nrdf > 0)
 122     {
 123         T /= nrdf;
 124     }
 125     if (dekindlambda)
 126     {
 127         if (bEkinAveVel)
 128         {
 129             *dekindlambda = ekind->dekindl;
 130         }
 131         else
 132         {
 133             *dekindlambda = 0.5 * (ekind->dekindl + ekind->dekindl_old);
 134         }
 135     }
 136     return T;
 137 }