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
  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
 101   /* bNEMD tells if we should remove remove the COM velocity
 102    * from the velocities during velocity scaling in T-coupling.
 103    * Turn this on when we have multiple acceleration groups
 104    * or one accelerated group.
 105    */
 106   ekind->bNEMD = (opts->ngacc > 1 || norm(opts->acc[0]) > 0);
 107
 108   ekind->ngtc = opts->ngtc;
 109   snew(ekind->tcstat,opts->ngtc);
 110   init_grptcstat(opts->ngtc,ekind->tcstat);
 111   /* Set Berendsen tcoupl lambda's to 1,
 112    * so runs without Berendsen coupling are not affected.
 113    */
 114   for(i=0; i<opts->ngtc; i++)
 115   {
 116       ekind->tcstat[i].lambda = 1.0;
 117       ekind->tcstat[i].vscale_nhc = 1.0;
 118       ekind->tcstat[i].ekinscaleh_nhc = 1.0;
 119       ekind->tcstat[i].ekinscalef_nhc = 1.0;
 120   }
 121
 122   ekind->ngacc = opts->ngacc;
 123   snew(ekind->grpstat,opts->ngacc);
 124   init_grpstat(log,mtop,opts->ngacc,ekind->grpstat);
 125 }
 126
 127 void accumulate_u(t_commrec *cr,t_grpopts *opts,gmx_ekindata_t *ekind)
 128 {
 129     /* This routine will only be called when it's necessary */
 130     t_bin *rb;
 131     int   g;
 132
 133     rb = mk_bin();
 134
 135     for(g=0; (g<opts->ngacc); g++)
 136     {
 137         add_binr(rb,DIM,ekind->grpstat[g].u);
 138     }
 139     sum_bin(rb,cr);
 140
 141     for(g=0; (g<opts->ngacc); g++)
 142     {
 143         extract_binr(rb,DIM*g,DIM,ekind->grpstat[g].u);
 144     }
 145     destroy_bin(rb);
 146 }
 147
 148 /* I don't think accumulate_ekin is used anymore? */
 149
 150 #if 0
 151 static void accumulate_ekin(t_commrec *cr,t_grpopts *opts,
 152                             gmx_ekindata_t *ekind)
 153 {
 154     int g;
 155
 156     if(PAR(cr))
 157     {
 158         for(g=0; (g<opts->ngtc); g++)
 159         {
 160             gmx_sum(DIM*DIM,ekind->tcstat[g].ekinf[0],cr);
 161         }
 162     }
 163 }
 164 #endif
 165
 166 void update_ekindata(int start,int homenr,gmx_ekindata_t *ekind,
 167                      t_grpopts *opts,rvec v[],t_mdatoms *md,real lambda)
 168 {
 169   int  d,g,n;
 170   real mv;
 171
 172   /* calculate mean velocities at whole timestep */
 173   for(g=0; (g<opts->ngtc); g++) {
 174     ekind->tcstat[g].T = 0;
 175   }
 176
 177   if (ekind->bNEMD) {
 178     for (g=0; (g<opts->ngacc); g++)
 179       clear_rvec(ekind->grpstat[g].u);
 180
 181     g = 0;
 182     for(n=start; (n<start+homenr); n++) {
 183       if (md->cACC)
 184         g = md->cACC[n];
 185       for(d=0; (d<DIM);d++) {
 186           mv = md->massT[n]*v[n][d];
 187           ekind->grpstat[g].u[d] += mv;
 188       }
 189     }
 190
 191     for (g=0; (g < opts->ngacc); g++) {
 192       for(d=0; (d<DIM);d++) {
 193           ekind->grpstat[g].u[d] /=
 194               (1-lambda)*ekind->grpstat[g].mA + lambda*ekind->grpstat[g].mB;
 195       }
 196     }
 197   }
 198 }
 199
 200 real sum_ekin(t_grpopts *opts,gmx_ekindata_t *ekind,real *dekindlambda,
 201               gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld, gmx_bool bScaleEkin)
 202 {
 203     int          i,j,m,ngtc;
 204     real         T,ek;
 205     t_grp_tcstat *tcstat;
 206     real         nrdf,nd,*ndf;
 207
 208     ngtc = opts->ngtc;
 209     ndf  = opts->nrdf;
 210
 211     T = 0;
 212     nrdf = 0;
 213
 214     clear_mat(ekind->ekin);
 215
 216     for(i=0; (i<ngtc); i++)
 217     {
 218
 219         nd = ndf[i];
 220         tcstat = &ekind->tcstat[i];
 221         /* Sometimes a group does not have degrees of freedom, e.g.
 222          * when it consists of shells and virtual sites, then we just
 223          * set the temperatue to 0 and also neglect the kinetic
 224          * energy, which should be  zero anyway.
 225          */
 226
 227         if (nd > 0) {
 228             if (bEkinAveVel)
 229             {
 230                 if (!bScaleEkin)
 231                 {
 232                     /* in this case, kinetic energy is from the current velocities already */
 233                     msmul(tcstat->ekinf,tcstat->ekinscalef_nhc,tcstat->ekinf);
 234                 }
 235             }
 236             else
 237
 238             {
 239                 /* Calculate the full step Ekin as the average of the half steps */
 240                 for(j=0; (j<DIM); j++)
 241                 {
 242                     for(m=0; (m<DIM); m++)
 243                     {
 244                         tcstat->ekinf[j][m] =
 245                             0.5*(tcstat->ekinh[j][m]*tcstat->ekinscaleh_nhc + tcstat->ekinh_old[j][m]);
 246                     }
 247                 }
 248             }
 249             m_add(tcstat->ekinf,ekind->ekin,ekind->ekin);
 250
 251             tcstat->Th = calc_temp(trace(tcstat->ekinh),nd);
 252             tcstat->T  = calc_temp(trace(tcstat->ekinf),nd);
 253
 254             /* after the scaling factors have been multiplied in, we can remove them */
 255             if (bEkinAveVel)
 256             {
 257                 tcstat->ekinscalef_nhc = 1.0;
 258             }
 259             else
 260             {
 261                 tcstat->ekinscaleh_nhc = 1.0;
 262             }
 263         }
 264         else
 265         {
 266             tcstat->T  = 0;
 267             tcstat->Th = 0;
 268         }
 269         T    += nd*tcstat->T;
 270         nrdf += nd;
 271     }
 272     if (nrdf > 0)
 273     {
 274         T/=nrdf;
 275     }
 276     if (dekindlambda)
 277     {
 278         *dekindlambda = 0.5*(ekind->dekindl + ekind->dekindl_old);
 279     }
 280     return T;
 281 }