Fix dEkin/dl handling with multiple ranks
[alexxy/gromacs.git] / src / gromacs / mdlib / stat.c
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5  * Copyright (c) 2001-2004, The GROMACS development team.
6  * Copyright (c) 2013,2014,2016, by the GROMACS development team, led by
7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8  * and including many others, as listed in the AUTHORS file in the
9  * top-level source directory and at http://www.gromacs.org.
10  *
11  * GROMACS is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU Lesser General Public License
13  * as published by the Free Software Foundation; either version 2.1
14  * of the License, or (at your option) any later version.
15  *
16  * GROMACS is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19  * Lesser General Public License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with GROMACS; if not, see
23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
25  *
26  * If you want to redistribute modifications to GROMACS, please
27  * consider that scientific software is very special. Version
28  * control is crucial - bugs must be traceable. We will be happy to
29  * consider code for inclusion in the official distribution, but
30  * derived work must not be called official GROMACS. Details are found
31  * in the README & COPYING files - if they are missing, get the
32  * official version at http://www.gromacs.org.
33  *
34  * To help us fund GROMACS development, we humbly ask that you cite
35  * the research papers on the package. Check out http://www.gromacs.org.
36  */
37 #ifdef HAVE_CONFIG_H
38 #include <config.h>
39 #endif
40
41 #include <string.h>
42 #include <stdio.h>
43 #include "typedefs.h"
44 #include "types/commrec.h"
45 #include "sysstuff.h"
46 #include "gmx_fatal.h"
47 #include "network.h"
48 #include "txtdump.h"
49 #include "names.h"
50 #include "physics.h"
51 #include "vec.h"
52 #include "gromacs/math/utilities.h"
53 #include "main.h"
54 #include "force.h"
55 #include "vcm.h"
56 #include "gromacs/utility/smalloc.h"
57 #include "gromacs/fileio/futil.h"
58 #include "network.h"
59 #include "rbin.h"
60 #include "tgroup.h"
61 #include "gromacs/fileio/xtcio.h"
62 #include "gromacs/fileio/gmxfio.h"
63 #include "gromacs/fileio/trnio.h"
64 #include "domdec.h"
65 #include "constr.h"
66 #include "checkpoint.h"
67 #include "xvgr.h"
68 #include "md_support.h"
69 #include "mdrun.h"
70 #include "sim_util.h"
71
72 typedef struct gmx_global_stat
73 {
74     t_bin *rb;
75     int   *itc0;
76     int   *itc1;
77 } t_gmx_global_stat;
78
79 gmx_global_stat_t global_stat_init(t_inputrec *ir)
80 {
81     gmx_global_stat_t gs;
82
83     snew(gs, 1);
84
85     gs->rb = mk_bin();
86     snew(gs->itc0, ir->opts.ngtc);
87     snew(gs->itc1, ir->opts.ngtc);
88
89     return gs;
90 }
91
92 void global_stat_destroy(gmx_global_stat_t gs)
93 {
94     destroy_bin(gs->rb);
95     sfree(gs->itc0);
96     sfree(gs->itc1);
97     sfree(gs);
98 }
99
100 static int filter_enerdterm(real *afrom, gmx_bool bToBuffer, real *ato,
101                             gmx_bool bTemp, gmx_bool bPres, gmx_bool bEner)
102 {
103     int i, to, from;
104
105     from = 0;
106     to   = 0;
107     for (i = 0; i < F_NRE; i++)
108     {
109         if (bToBuffer)
110         {
111             from = i;
112         }
113         else
114         {
115             to = i;
116         }
117         switch (i)
118         {
119             case F_EKIN:
120             case F_TEMP:
121             case F_DKDL:
122                 if (bTemp)
123                 {
124                     ato[to++] = afrom[from++];
125                 }
126                 break;
127             case F_PRES:
128             case F_PDISPCORR:
129                 if (bPres)
130                 {
131                     ato[to++] = afrom[from++];
132                 }
133                 break;
134             default:
135                 if (bEner)
136                 {
137                     ato[to++] = afrom[from++];
138                 }
139                 break;
140         }
141     }
142
143     return to;
144 }
145
146 void global_stat(FILE *fplog, gmx_global_stat_t gs,
147                  t_commrec *cr, gmx_enerdata_t *enerd,
148                  tensor fvir, tensor svir, rvec mu_tot,
149                  t_inputrec *inputrec,
150                  gmx_ekindata_t *ekind, gmx_constr_t constr,
151                  t_vcm *vcm,
152                  int nsig, real *sig,
153                  gmx_mtop_t *top_global, t_state *state_local,
154                  gmx_bool bSumEkinhOld, int flags)
155 /* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
156 {
157     t_bin     *rb;
158     int       *itc0, *itc1;
159     int        ie    = 0, ifv = 0, isv = 0, irmsd = 0, imu = 0;
160     int        idedl = 0, idedlo = 0, idvdll = 0, idvdlnl = 0, iepl = 0, icm = 0, imass = 0, ica = 0, inb = 0;
161     int        isig  = -1;
162     int        icj   = -1, ici = -1, icx = -1;
163     int        inn[egNR];
164     real       copyenerd[F_NRE];
165     int        nener, j;
166     real      *rmsd_data = NULL;
167     double     nb;
168     gmx_bool   bVV, bTemp, bEner, bPres, bConstrVir, bEkinAveVel, bFirstIterate, bReadEkin;
169
170     bVV           = EI_VV(inputrec->eI);
171     bTemp         = flags & CGLO_TEMPERATURE;
172     bEner         = flags & CGLO_ENERGY;
173     bPres         = (flags & CGLO_PRESSURE);
174     bConstrVir    = (flags & CGLO_CONSTRAINT);
175     bFirstIterate = (flags & CGLO_FIRSTITERATE);
176     bEkinAveVel   = (inputrec->eI == eiVV || (inputrec->eI == eiVVAK && bPres));
177     bReadEkin     = (flags & CGLO_READEKIN);
178
179     rb   = gs->rb;
180     itc0 = gs->itc0;
181     itc1 = gs->itc1;
182
183
184     reset_bin(rb);
185     /* This routine copies all the data to be summed to one big buffer
186      * using the t_bin struct.
187      */
188
189     /* First, we neeed to identify which enerd->term should be
190        communicated.  Temperature and pressure terms should only be
191        communicated and summed when they need to be, to avoid repeating
192        the sums and overcounting. */
193
194     nener = filter_enerdterm(enerd->term, TRUE, copyenerd, bTemp, bPres, bEner);
195
196     /* First, the data that needs to be communicated with velocity verlet every time
197        This is just the constraint virial.*/
198     if (bConstrVir)
199     {
200         isv = add_binr(rb, DIM*DIM, svir[0]);
201         where();
202     }
203
204 /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
205     if (bTemp || !bVV)
206     {
207         if (ekind)
208         {
209             for (j = 0; (j < inputrec->opts.ngtc); j++)
210             {
211                 if (bSumEkinhOld)
212                 {
213                     itc0[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinh_old[0]);
214                 }
215                 if (bEkinAveVel && !bReadEkin)
216                 {
217                     itc1[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinf[0]);
218                 }
219                 else if (!bReadEkin)
220                 {
221                     itc1[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinh[0]);
222                 }
223             }
224             /* these probably need to be put into one of these categories */
225             where();
226             idedl = add_binr(rb, 1, &(ekind->dekindl));
227             if (bSumEkinhOld)
228             {
229                 idedlo = add_binr(rb, 1, &(ekind->dekindl_old));
230             }
231             where();
232             ica   = add_binr(rb, 1, &(ekind->cosacc.mvcos));
233             where();
234         }
235     }
236     where();
237
238     if ((bPres || !bVV) && bFirstIterate)
239     {
240         ifv = add_binr(rb, DIM*DIM, fvir[0]);
241     }
242
243
244     if (bEner)
245     {
246         where();
247         if (bFirstIterate)
248         {
249             ie  = add_binr(rb, nener, copyenerd);
250         }
251         where();
252         if (constr)
253         {
254             rmsd_data = constr_rmsd_data(constr);
255             if (rmsd_data)
256             {
257                 irmsd = add_binr(rb, inputrec->eI == eiSD2 ? 3 : 2, rmsd_data);
258             }
259         }
260         if (!NEED_MUTOT(*inputrec))
261         {
262             imu = add_binr(rb, DIM, mu_tot);
263             where();
264         }
265
266         if (bFirstIterate)
267         {
268             for (j = 0; (j < egNR); j++)
269             {
270                 inn[j] = add_binr(rb, enerd->grpp.nener, enerd->grpp.ener[j]);
271             }
272             where();
273             if (inputrec->efep != efepNO)
274             {
275                 idvdll  = add_bind(rb, efptNR, enerd->dvdl_lin);
276                 idvdlnl = add_bind(rb, efptNR, enerd->dvdl_nonlin);
277                 if (enerd->n_lambda > 0)
278                 {
279                     iepl = add_bind(rb, enerd->n_lambda, enerd->enerpart_lambda);
280                 }
281             }
282         }
283     }
284
285     if (vcm)
286     {
287         icm   = add_binr(rb, DIM*vcm->nr, vcm->group_p[0]);
288         where();
289         imass = add_binr(rb, vcm->nr, vcm->group_mass);
290         where();
291         if (vcm->mode == ecmANGULAR)
292         {
293             icj   = add_binr(rb, DIM*vcm->nr, vcm->group_j[0]);
294             where();
295             icx   = add_binr(rb, DIM*vcm->nr, vcm->group_x[0]);
296             where();
297             ici   = add_binr(rb, DIM*DIM*vcm->nr, vcm->group_i[0][0]);
298             where();
299         }
300     }
301
302     if (DOMAINDECOMP(cr))
303     {
304         nb  = cr->dd->nbonded_local;
305         inb = add_bind(rb, 1, &nb);
306     }
307     where();
308     if (nsig > 0)
309     {
310         isig = add_binr(rb, nsig, sig);
311     }
312
313     /* Global sum it all */
314     if (debug)
315     {
316         fprintf(debug, "Summing %d energies\n", rb->maxreal);
317     }
318     sum_bin(rb, cr);
319     where();
320
321     /* Extract all the data locally */
322
323     if (bConstrVir)
324     {
325         extract_binr(rb, isv, DIM*DIM, svir[0]);
326     }
327
328     /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
329     if (bTemp || !bVV)
330     {
331         if (ekind)
332         {
333             for (j = 0; (j < inputrec->opts.ngtc); j++)
334             {
335                 if (bSumEkinhOld)
336                 {
337                     extract_binr(rb, itc0[j], DIM*DIM, ekind->tcstat[j].ekinh_old[0]);
338                 }
339                 if (bEkinAveVel && !bReadEkin)
340                 {
341                     extract_binr(rb, itc1[j], DIM*DIM, ekind->tcstat[j].ekinf[0]);
342                 }
343                 else if (!bReadEkin)
344                 {
345                     extract_binr(rb, itc1[j], DIM*DIM, ekind->tcstat[j].ekinh[0]);
346                 }
347             }
348             extract_binr(rb, idedl, 1, &(ekind->dekindl));
349             if (bSumEkinhOld)
350             {
351                 extract_binr(rb, idedlo, 1, &(ekind->dekindl_old));
352             }
353             extract_binr(rb, ica, 1, &(ekind->cosacc.mvcos));
354             where();
355         }
356     }
357     if ((bPres || !bVV) && bFirstIterate)
358     {
359         extract_binr(rb, ifv, DIM*DIM, fvir[0]);
360     }
361
362     if (bEner)
363     {
364         if (bFirstIterate)
365         {
366             extract_binr(rb, ie, nener, copyenerd);
367             if (rmsd_data)
368             {
369                 extract_binr(rb, irmsd, inputrec->eI == eiSD2 ? 3 : 2, rmsd_data);
370             }
371             if (!NEED_MUTOT(*inputrec))
372             {
373                 extract_binr(rb, imu, DIM, mu_tot);
374             }
375
376             for (j = 0; (j < egNR); j++)
377             {
378                 extract_binr(rb, inn[j], enerd->grpp.nener, enerd->grpp.ener[j]);
379             }
380             if (inputrec->efep != efepNO)
381             {
382                 extract_bind(rb, idvdll, efptNR, enerd->dvdl_lin);
383                 extract_bind(rb, idvdlnl, efptNR, enerd->dvdl_nonlin);
384                 if (enerd->n_lambda > 0)
385                 {
386                     extract_bind(rb, iepl, enerd->n_lambda, enerd->enerpart_lambda);
387                 }
388             }
389             if (DOMAINDECOMP(cr))
390             {
391                 extract_bind(rb, inb, 1, &nb);
392                 if ((int)(nb + 0.5) != cr->dd->nbonded_global)
393                 {
394                     dd_print_missing_interactions(fplog, cr, (int)(nb + 0.5), top_global, state_local);
395                 }
396             }
397             where();
398
399             filter_enerdterm(copyenerd, FALSE, enerd->term, bTemp, bPres, bEner);
400         }
401     }
402
403     if (vcm)
404     {
405         extract_binr(rb, icm, DIM*vcm->nr, vcm->group_p[0]);
406         where();
407         extract_binr(rb, imass, vcm->nr, vcm->group_mass);
408         where();
409         if (vcm->mode == ecmANGULAR)
410         {
411             extract_binr(rb, icj, DIM*vcm->nr, vcm->group_j[0]);
412             where();
413             extract_binr(rb, icx, DIM*vcm->nr, vcm->group_x[0]);
414             where();
415             extract_binr(rb, ici, DIM*DIM*vcm->nr, vcm->group_i[0][0]);
416             where();
417         }
418     }
419
420     if (nsig > 0)
421     {
422         extract_binr(rb, isig, nsig, sig);
423     }
424     where();
425 }
426
427 int do_per_step(gmx_int64_t step, gmx_int64_t nstep)
428 {
429     if (nstep != 0)
430     {
431         return ((step % nstep) == 0);
432     }
433     else
434     {
435         return 0;
436     }
437 }