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