src/gromacs/mdlib/tgroup.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROwing Monsters And Cloning Shrimps
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40
  41 #include <math.h>
  42 #include "macros.h"
  43 #include "main.h"
  44 #include "smalloc.h"
  45 #include "futil.h"
  46 #include "tgroup.h"
  47 #include "vec.h"
  48 #include "network.h"
  49 #include "smalloc.h"
  50 #include "update.h"
  51 #include "rbin.h"
  52 #include "mtop_util.h"
  53 #include "gmx_omp_nthreads.h"
  54
  55 static void init_grptcstat(int ngtc, t_grp_tcstat tcstat[])
  56 {
  57     int i, j;
  58
  59     for (i = 0; (i < ngtc); i++)
  60     {
  61         tcstat[i].T = 0;
  62         clear_mat(tcstat[i].ekinh);
  63         clear_mat(tcstat[i].ekinh_old);
  64         clear_mat(tcstat[i].ekinf);
  65     }
  66 }
  67
  68 static void init_grpstat(FILE *log,
  69                          gmx_mtop_t *mtop, int ngacc, t_grp_acc gstat[])
  70 {
  71     gmx_groups_t           *groups;
  72     gmx_mtop_atomloop_all_t aloop;
  73     int                     i, grp;
  74     t_atom                 *atom;
  75
  76     if (ngacc > 0)
  77     {
  78         groups = &mtop->groups;
  79         aloop  = gmx_mtop_atomloop_all_init(mtop);
  80         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
  81         {
  82             grp = ggrpnr(groups, egcACC, i);
  83             if ((grp < 0) && (grp >= ngacc))
  84             {
  85                 gmx_incons("Input for acceleration groups wrong");
  86             }
  87             gstat[grp].nat++;
  88             /* This will not work for integrator BD */
  89             gstat[grp].mA += atom->m;
  90             gstat[grp].mB += atom->mB;
  91         }
  92     }
  93 }
  94
  95 void init_ekindata(FILE *log, gmx_mtop_t *mtop, t_grpopts *opts,
  96                    gmx_ekindata_t *ekind)
  97 {
  98     int i;
  99     int nthread, thread;
 100 #ifdef DEBUG
 101     fprintf(log, "ngtc: %d, ngacc: %d, ngener: %d\n", opts->ngtc, opts->ngacc,
 102             opts->ngener);
 103 #endif
 104
 105     /* bNEMD tells if we should remove remove the COM velocity
 106      * from the velocities during velocity scaling in T-coupling.
 107      * Turn this on when we have multiple acceleration groups
 108      * or one accelerated group.
 109      */
 110     ekind->bNEMD = (opts->ngacc > 1 || norm(opts->acc[0]) > 0);
 111
 112     ekind->ngtc = opts->ngtc;
 113     snew(ekind->tcstat, opts->ngtc);
 114     init_grptcstat(opts->ngtc, ekind->tcstat);
 115     /* Set Berendsen tcoupl lambda's to 1,
 116      * so runs without Berendsen coupling are not affected.
 117      */
 118     for (i = 0; i < opts->ngtc; i++)
 119     {
 120         ekind->tcstat[i].lambda         = 1.0;
 121         ekind->tcstat[i].vscale_nhc     = 1.0;
 122         ekind->tcstat[i].ekinscaleh_nhc = 1.0;
 123         ekind->tcstat[i].ekinscalef_nhc = 1.0;
 124     }
 125
 126     nthread = gmx_omp_nthreads_get(emntUpdate);
 127
 128     snew(ekind->ekin_work_alloc, nthread);
 129     snew(ekind->ekin_work, nthread);
 130 #pragma omp parallel for num_threads(nthread) schedule(static)
 131     for (thread = 0; thread < nthread; thread++)
 132     {
 133         /* Allocate 2 elements extra on both sides,
 134          * so in single precision we have 2*3*3*4=72 bytes buffer
 135          * on both sides to avoid cache pollution.
 136          */
 137         snew(ekind->ekin_work_alloc[thread], ekind->ngtc+4);
 138         ekind->ekin_work[thread] = ekind->ekin_work_alloc[thread] + 2;
 139     }
 140
 141     ekind->ngacc = opts->ngacc;
 142     snew(ekind->grpstat, opts->ngacc);
 143     init_grpstat(log, mtop, opts->ngacc, ekind->grpstat);
 144 }
 145
 146 void accumulate_u(t_commrec *cr, t_grpopts *opts, gmx_ekindata_t *ekind)
 147 {
 148     /* This routine will only be called when it's necessary */
 149     t_bin *rb;
 150     int    g;
 151
 152     rb = mk_bin();
 153
 154     for (g = 0; (g < opts->ngacc); g++)
 155     {
 156         add_binr(rb, DIM, ekind->grpstat[g].u);
 157     }
 158     sum_bin(rb, cr);
 159
 160     for (g = 0; (g < opts->ngacc); g++)
 161     {
 162         extract_binr(rb, DIM*g, DIM, ekind->grpstat[g].u);
 163     }
 164     destroy_bin(rb);
 165 }
 166
 167 /* I don't think accumulate_ekin is used anymore? */
 168
 169 #if 0
 170 static void accumulate_ekin(t_commrec *cr, t_grpopts *opts,
 171                             gmx_ekindata_t *ekind)
 172 {
 173     int g;
 174
 175     if (PAR(cr))
 176     {
 177         for (g = 0; (g < opts->ngtc); g++)
 178         {
 179             gmx_sum(DIM*DIM, ekind->tcstat[g].ekinf[0], cr);
 180         }
 181     }
 182 }
 183 #endif
 184
 185 void update_ekindata(int start, int homenr, gmx_ekindata_t *ekind,
 186                      t_grpopts *opts, rvec v[], t_mdatoms *md, real lambda)
 187 {
 188     int  d, g, n;
 189     real mv;
 190
 191     /* calculate mean velocities at whole timestep */
 192     for (g = 0; (g < opts->ngtc); g++)
 193     {
 194         ekind->tcstat[g].T = 0;
 195     }
 196
 197     if (ekind->bNEMD)
 198     {
 199         for (g = 0; (g < opts->ngacc); g++)
 200         {
 201             clear_rvec(ekind->grpstat[g].u);
 202         }
 203
 204         g = 0;
 205         for (n = start; (n < start+homenr); n++)
 206         {
 207             if (md->cACC)
 208             {
 209                 g = md->cACC[n];
 210             }
 211             for (d = 0; (d < DIM); d++)
 212             {
 213                 mv                      = md->massT[n]*v[n][d];
 214                 ekind->grpstat[g].u[d] += mv;
 215             }
 216         }
 217
 218         for (g = 0; (g < opts->ngacc); g++)
 219         {
 220             for (d = 0; (d < DIM); d++)
 221             {
 222                 ekind->grpstat[g].u[d] /=
 223                     (1-lambda)*ekind->grpstat[g].mA + lambda*ekind->grpstat[g].mB;
 224             }
 225         }
 226     }
 227 }
 228
 229 real sum_ekin(t_grpopts *opts, gmx_ekindata_t *ekind, real *dekindlambda,
 230               gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld, gmx_bool bScaleEkin)
 231 {
 232     int           i, j, m, ngtc;
 233     real          T, ek;
 234     t_grp_tcstat *tcstat;
 235     real          nrdf, nd, *ndf;
 236
 237     ngtc = opts->ngtc;
 238     ndf  = opts->nrdf;
 239
 240     T    = 0;
 241     nrdf = 0;
 242
 243     clear_mat(ekind->ekin);
 244
 245     for (i = 0; (i < ngtc); i++)
 246     {
 247
 248         nd     = ndf[i];
 249         tcstat = &ekind->tcstat[i];
 250         /* Sometimes a group does not have degrees of freedom, e.g.
 251          * when it consists of shells and virtual sites, then we just
 252          * set the temperatue to 0 and also neglect the kinetic
 253          * energy, which should be  zero anyway.
 254          */
 255
 256         if (nd > 0)
 257         {
 258             if (bEkinAveVel)
 259             {
 260                 if (!bScaleEkin)
 261                 {
 262                     /* in this case, kinetic energy is from the current velocities already */
 263                     msmul(tcstat->ekinf, tcstat->ekinscalef_nhc, tcstat->ekinf);
 264                 }
 265             }
 266             else
 267
 268             {
 269                 /* Calculate the full step Ekin as the average of the half steps */
 270                 for (j = 0; (j < DIM); j++)
 271                 {
 272                     for (m = 0; (m < DIM); m++)
 273                     {
 274                         tcstat->ekinf[j][m] =
 275                             0.5*(tcstat->ekinh[j][m]*tcstat->ekinscaleh_nhc + tcstat->ekinh_old[j][m]);
 276                     }
 277                 }
 278             }
 279             m_add(tcstat->ekinf, ekind->ekin, ekind->ekin);
 280
 281             tcstat->Th = calc_temp(trace(tcstat->ekinh), nd);
 282             tcstat->T  = calc_temp(trace(tcstat->ekinf), nd);
 283
 284             /* after the scaling factors have been multiplied in, we can remove them */
 285             if (bEkinAveVel)
 286             {
 287                 tcstat->ekinscalef_nhc = 1.0;
 288             }
 289             else
 290             {
 291                 tcstat->ekinscaleh_nhc = 1.0;
 292             }
 293         }
 294         else
 295         {
 296             tcstat->T  = 0;
 297             tcstat->Th = 0;
 298         }
 299         T    += nd*tcstat->T;
 300         nrdf += nd;
 301     }
 302     if (nrdf > 0)
 303     {
 304         T /= nrdf;
 305     }
 306     if (dekindlambda)
 307     {
 308         *dekindlambda = 0.5*(ekind->dekindl + ekind->dekindl_old);
 309     }
 310     return T;
 311 }