src/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         tcstat[i].T = 0;
  61         clear_mat(tcstat[i].ekinh);
  62         clear_mat(tcstat[i].ekinh_old);
  63         clear_mat(tcstat[i].ekinf);
  64     }
  65 }
  66
  67 static void init_grpstat(FILE *log,
  68                          gmx_mtop_t *mtop,int ngacc,t_grp_acc gstat[])
  69 {
  70     gmx_groups_t *groups;
  71     gmx_mtop_atomloop_all_t aloop;
  72     int    i,grp;
  73     t_atom *atom;
  74
  75     if (ngacc > 0) {
  76         groups = &mtop->groups;
  77         aloop = gmx_mtop_atomloop_all_init(mtop);
  78         while (gmx_mtop_atomloop_all_next(aloop,&i,&atom))
  79         {
  80             grp = ggrpnr(groups,egcACC,i);
  81             if ((grp < 0) && (grp >= ngacc))
  82             {
  83                 gmx_incons("Input for acceleration groups wrong");
  84             }
  85             gstat[grp].nat++;
  86             /* This will not work for integrator BD */
  87             gstat[grp].mA += atom->m;
  88             gstat[grp].mB += atom->mB;
  89         }
  90     }
  91 }
  92
  93 void init_ekindata(FILE *log,gmx_mtop_t *mtop,t_grpopts *opts,
  94                    gmx_ekindata_t *ekind)
  95 {
  96   int i;
  97   int nthread,thread;
  98 #ifdef DEBUG
  99   fprintf(log,"ngtc: %d, ngacc: %d, ngener: %d\n",opts->ngtc,opts->ngacc,
 100           opts->ngener);
 101 #endif
 102
 103   /* bNEMD tells if we should remove remove the COM velocity
 104    * from the velocities during velocity scaling in T-coupling.
 105    * Turn this on when we have multiple acceleration groups
 106    * or one accelerated group.
 107    */
 108   ekind->bNEMD = (opts->ngacc > 1 || norm(opts->acc[0]) > 0);
 109
 110   ekind->ngtc = opts->ngtc;
 111   snew(ekind->tcstat,opts->ngtc);
 112   init_grptcstat(opts->ngtc,ekind->tcstat);
 113   /* Set Berendsen tcoupl lambda's to 1,
 114    * so runs without Berendsen coupling are not affected.
 115    */
 116   for(i=0; i<opts->ngtc; i++)
 117   {
 118       ekind->tcstat[i].lambda = 1.0;
 119       ekind->tcstat[i].vscale_nhc = 1.0;
 120       ekind->tcstat[i].ekinscaleh_nhc = 1.0;
 121       ekind->tcstat[i].ekinscalef_nhc = 1.0;
 122   }
 123
 124     nthread = gmx_omp_nthreads_get(emntUpdate);
 125
 126     snew(ekind->ekin_work_alloc,nthread);
 127     snew(ekind->ekin_work,nthread);
 128 #pragma omp parallel for num_threads(nthread) schedule(static)
 129     for(thread=0; thread<nthread; thread++)
 130     {
 131         /* Allocate 2 elements extra on both sides,
 132          * so in single precision we have 2*3*3*4=72 bytes buffer
 133          * on both sides to avoid cache pollution.
 134          */
 135         snew(ekind->ekin_work_alloc[thread],ekind->ngtc+4);
 136         ekind->ekin_work[thread] = ekind->ekin_work_alloc[thread] + 2;
 137     }
 138
 139   ekind->ngacc = opts->ngacc;
 140   snew(ekind->grpstat,opts->ngacc);
 141   init_grpstat(log,mtop,opts->ngacc,ekind->grpstat);
 142 }
 143
 144 void accumulate_u(t_commrec *cr,t_grpopts *opts,gmx_ekindata_t *ekind)
 145 {
 146     /* This routine will only be called when it's necessary */
 147     t_bin *rb;
 148     int   g;
 149
 150     rb = mk_bin();
 151
 152     for(g=0; (g<opts->ngacc); g++)
 153     {
 154         add_binr(rb,DIM,ekind->grpstat[g].u);
 155     }
 156     sum_bin(rb,cr);
 157
 158     for(g=0; (g<opts->ngacc); g++)
 159     {
 160         extract_binr(rb,DIM*g,DIM,ekind->grpstat[g].u);
 161     }
 162     destroy_bin(rb);
 163 }
 164
 165 /* I don't think accumulate_ekin is used anymore? */
 166
 167 #if 0
 168 static void accumulate_ekin(t_commrec *cr,t_grpopts *opts,
 169                             gmx_ekindata_t *ekind)
 170 {
 171     int g;
 172
 173     if(PAR(cr))
 174     {
 175         for(g=0; (g<opts->ngtc); g++)
 176         {
 177             gmx_sum(DIM*DIM,ekind->tcstat[g].ekinf[0],cr);
 178         }
 179     }
 180 }
 181 #endif
 182
 183 void update_ekindata(int start,int homenr,gmx_ekindata_t *ekind,
 184                      t_grpopts *opts,rvec v[],t_mdatoms *md,real lambda)
 185 {
 186   int  d,g,n;
 187   real mv;
 188
 189   /* calculate mean velocities at whole timestep */
 190   for(g=0; (g<opts->ngtc); g++) {
 191     ekind->tcstat[g].T = 0;
 192   }
 193
 194   if (ekind->bNEMD) {
 195     for (g=0; (g<opts->ngacc); g++)
 196       clear_rvec(ekind->grpstat[g].u);
 197
 198     g = 0;
 199     for(n=start; (n<start+homenr); n++) {
 200       if (md->cACC)
 201         g = md->cACC[n];
 202       for(d=0; (d<DIM);d++) {
 203           mv = md->massT[n]*v[n][d];
 204           ekind->grpstat[g].u[d] += mv;
 205       }
 206     }
 207
 208     for (g=0; (g < opts->ngacc); g++) {
 209       for(d=0; (d<DIM);d++) {
 210           ekind->grpstat[g].u[d] /=
 211               (1-lambda)*ekind->grpstat[g].mA + lambda*ekind->grpstat[g].mB;
 212       }
 213     }
 214   }
 215 }
 216
 217 real sum_ekin(t_grpopts *opts,gmx_ekindata_t *ekind,real *dekindlambda,
 218               gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld, gmx_bool bScaleEkin)
 219 {
 220     int          i,j,m,ngtc;
 221     real         T,ek;
 222     t_grp_tcstat *tcstat;
 223     real         nrdf,nd,*ndf;
 224
 225     ngtc = opts->ngtc;
 226     ndf  = opts->nrdf;
 227
 228     T = 0;
 229     nrdf = 0;
 230
 231     clear_mat(ekind->ekin);
 232
 233     for(i=0; (i<ngtc); i++)
 234     {
 235
 236         nd = ndf[i];
 237         tcstat = &ekind->tcstat[i];
 238         /* Sometimes a group does not have degrees of freedom, e.g.
 239          * when it consists of shells and virtual sites, then we just
 240          * set the temperatue to 0 and also neglect the kinetic
 241          * energy, which should be  zero anyway.
 242          */
 243
 244         if (nd > 0) {
 245             if (bEkinAveVel)
 246             {
 247                 if (!bScaleEkin)
 248                 {
 249                     /* in this case, kinetic energy is from the current velocities already */
 250                     msmul(tcstat->ekinf,tcstat->ekinscalef_nhc,tcstat->ekinf);
 251                 }
 252             }
 253             else
 254
 255             {
 256                 /* Calculate the full step Ekin as the average of the half steps */
 257                 for(j=0; (j<DIM); j++)
 258                 {
 259                     for(m=0; (m<DIM); m++)
 260                     {
 261                         tcstat->ekinf[j][m] =
 262                             0.5*(tcstat->ekinh[j][m]*tcstat->ekinscaleh_nhc + tcstat->ekinh_old[j][m]);
 263                     }
 264                 }
 265             }
 266             m_add(tcstat->ekinf,ekind->ekin,ekind->ekin);
 267
 268             tcstat->Th = calc_temp(trace(tcstat->ekinh),nd);
 269             tcstat->T  = calc_temp(trace(tcstat->ekinf),nd);
 270
 271             /* after the scaling factors have been multiplied in, we can remove them */
 272             if (bEkinAveVel)
 273             {
 274                 tcstat->ekinscalef_nhc = 1.0;
 275             }
 276             else
 277             {
 278                 tcstat->ekinscaleh_nhc = 1.0;
 279             }
 280         }
 281         else
 282         {
 283             tcstat->T  = 0;
 284             tcstat->Th = 0;
 285         }
 286         T    += nd*tcstat->T;
 287         nrdf += nd;
 288     }
 289     if (nrdf > 0)
 290     {
 291         T/=nrdf;
 292     }
 293     if (dekindlambda)
 294     {
 295         *dekindlambda = 0.5*(ekind->dekindl + ekind->dekindl_old);
 296     }
 297     return T;
 298 }