1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
35 /* This file is completely threadsafe - keep it that way! */
52 #include "mtop_util.h"
54 static void init_grptcstat(int ngtc,t_grp_tcstat tcstat[])
58 for(i=0; (i<ngtc); i++) {
60 clear_mat(tcstat[i].ekinh);
61 clear_mat(tcstat[i].ekinh_old);
62 clear_mat(tcstat[i].ekinf);
66 static void init_grpstat(FILE *log,
67 gmx_mtop_t *mtop,int ngacc,t_grp_acc gstat[])
70 gmx_mtop_atomloop_all_t aloop;
75 groups = &mtop->groups;
76 aloop = gmx_mtop_atomloop_all_init(mtop);
77 while (gmx_mtop_atomloop_all_next(aloop,&i,&atom))
79 grp = ggrpnr(groups,egcACC,i);
80 if ((grp < 0) && (grp >= ngacc))
82 gmx_incons("Input for acceleration groups wrong");
85 /* This will not work for integrator BD */
86 gstat[grp].mA += atom->m;
87 gstat[grp].mB += atom->mB;
92 void init_ekindata(FILE *log,gmx_mtop_t *mtop,t_grpopts *opts,
93 gmx_ekindata_t *ekind)
97 fprintf(log,"ngtc: %d, ngacc: %d, ngener: %d\n",opts->ngtc,opts->ngacc,
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.
106 ekind->bNEMD = (opts->ngacc > 1 || norm(opts->acc[0]) > 0);
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.
114 for(i=0; i<opts->ngtc; i++)
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;
122 ekind->ngacc = opts->ngacc;
123 snew(ekind->grpstat,opts->ngacc);
124 init_grpstat(log,mtop,opts->ngacc,ekind->grpstat);
127 void accumulate_u(t_commrec *cr,t_grpopts *opts,gmx_ekindata_t *ekind)
129 /* This routine will only be called when it's necessary */
135 for(g=0; (g<opts->ngacc); g++)
137 add_binr(rb,DIM,ekind->grpstat[g].u);
141 for(g=0; (g<opts->ngacc); g++)
143 extract_binr(rb,DIM*g,DIM,ekind->grpstat[g].u);
148 /* I don't think accumulate_ekin is used anymore? */
151 static void accumulate_ekin(t_commrec *cr,t_grpopts *opts,
152 gmx_ekindata_t *ekind)
158 for(g=0; (g<opts->ngtc); g++)
160 gmx_sum(DIM*DIM,ekind->tcstat[g].ekinf[0],cr);
166 void update_ekindata(int start,int homenr,gmx_ekindata_t *ekind,
167 t_grpopts *opts,rvec v[],t_mdatoms *md,real lambda)
172 /* calculate mean velocities at whole timestep */
173 for(g=0; (g<opts->ngtc); g++) {
174 ekind->tcstat[g].T = 0;
178 for (g=0; (g<opts->ngacc); g++)
179 clear_rvec(ekind->grpstat[g].u);
182 for(n=start; (n<start+homenr); n++) {
185 for(d=0; (d<DIM);d++) {
186 mv = md->massT[n]*v[n][d];
187 ekind->grpstat[g].u[d] += mv;
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;
200 real sum_ekin(t_grpopts *opts,gmx_ekindata_t *ekind,real *dekindlambda,
201 gmx_bool bEkinAveVel, gmx_bool bSaveEkinOld, gmx_bool bScaleEkin)
205 t_grp_tcstat *tcstat;
214 clear_mat(ekind->ekin);
216 for(i=0; (i<ngtc); 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.
232 /* in this case, kinetic energy is from the current velocities already */
233 msmul(tcstat->ekinf,tcstat->ekinscalef_nhc,tcstat->ekinf);
239 /* Calculate the full step Ekin as the average of the half steps */
240 for(j=0; (j<DIM); j++)
242 for(m=0; (m<DIM); m++)
244 tcstat->ekinf[j][m] =
245 0.5*(tcstat->ekinh[j][m]*tcstat->ekinscaleh_nhc + tcstat->ekinh_old[j][m]);
249 m_add(tcstat->ekinf,ekind->ekin,ekind->ekin);
251 tcstat->Th = calc_temp(trace(tcstat->ekinh),nd);
252 tcstat->T = calc_temp(trace(tcstat->ekinf),nd);
254 /* after the scaling factors have been multiplied in, we can remove them */
257 tcstat->ekinscalef_nhc = 1.0;
261 tcstat->ekinscaleh_nhc = 1.0;
278 *dekindlambda = 0.5*(ekind->dekindl + ekind->dekindl_old);