Project

General

Profile

pme_pp.c

fixed src/mdlib/pme_pp.c for Gromacs 4.0 - Berk Hess, 08/23/2010 10:00 AM

 
1
/*
2
 * $Id$
3
 * 
4
 *                This source code is part of
5
 * 
6
 *                 G   R   O   M   A   C   S
7
 * 
8
 *          GROningen MAchine for Chemical Simulations
9
 * 
10
 *                        VERSION 3.2.0
11
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13
 * Copyright (c) 2001-2004, The GROMACS development team,
14
 * check out http://www.gromacs.org for more information.
15

16
 * This program is free software; you can redistribute it and/or
17
 * modify it under the terms of the GNU General Public License
18
 * as published by the Free Software Foundation; either version 2
19
 * of the License, or (at your option) any later version.
20
 * 
21
 * If you want to redistribute modifications, please consider that
22
 * scientific software is very special. Version control is crucial -
23
 * bugs must be traceable. We will be happy to consider code for
24
 * inclusion in the official distribution, but derived work must not
25
 * be called official GROMACS. Details are found in the README & COPYING
26
 * files - if they are missing, get the official version at www.gromacs.org.
27
 * 
28
 * To help us fund GROMACS development, we humbly ask that you cite
29
 * the papers on the package - you can find them in the top README file.
30
 * 
31
 * For more info, check our website at http://www.gromacs.org
32
 * 
33
 * And Hey:
34
 * GROwing Monsters And Cloning Shrimps
35
 */
36

    
37
#ifdef HAVE_CONFIG_H
38
#include <config.h>
39
#endif
40

    
41

    
42
#include <stdio.h>
43
#include <string.h>
44
#include <math.h>
45
#include "typedefs.h"
46
#include "smalloc.h"
47
#include "gmx_fatal.h"
48
#include "vec.h"
49
#include "pme.h"
50
#include "network.h"
51
#include "domdec.h"
52
#ifdef GMX_MPI
53
#include <mpi.h>
54
#endif
55
#include "mpelogging.h"
56

    
57
#define PP_PME_CHARGE   (1<<0)
58
#define PP_PME_CHARGEB  (1<<1)
59
#define PP_PME_COORD    (1<<2)
60
#define PP_PME_FEP      (1<<3)
61
#define PP_PME_FINISH   (1<<4)
62

    
63
#define PME_PP_TERM     (1<<0)
64
#define PME_PP_USR1     (1<<1)
65

    
66
typedef struct gmx_pme_pp {
67
#ifdef GMX_MPI
68
  MPI_Comm mpi_comm_mysim;
69
#endif
70
  int  nnode;        /* The number of PP node to communicate with  */
71
  int  *node;        /* The PP node ranks                          */
72
  int  node_peer;    /* The peer PP node rank                      */
73
  int  *nat;         /* The number of atom for each PP node        */
74
  int  flags_charge; /* The flags sent along with the last charges */
75
  real *chargeA;
76
  real *chargeB;
77
  rvec *x;
78
  rvec *f;
79
  int  nalloc;
80
#ifdef GMX_MPI
81
  MPI_Request *req;
82
  MPI_Status  *stat;
83
#endif
84
} t_gmx_pme_pp;
85

    
86
typedef struct gmx_pme_comm_n_box {
87
  int    natoms;
88
  matrix box;
89
  int    maxshift;
90
  real   lambda;
91
  int    flags;
92
} gmx_pme_comm_n_box_t;
93

    
94
typedef struct {
95
  matrix vir;
96
  real   energy;
97
  real   dvdlambda;
98
  float  cycles;
99
  int    flags;
100
} gmx_pme_comm_vir_ene_t;
101

    
102

    
103
/* The following stuff is needed for signal handling on the PME nodes.
104
 * The signal variables are defined in pme.c and also used in md.c.
105
 */ 
106
extern bool bGotTermSignal, bGotUsr1Signal; 
107

    
108

    
109
gmx_pme_pp_t gmx_pme_pp_init(t_commrec *cr)
110
{
111
  struct gmx_pme_pp *pme_pp;
112
  int rank;
113

    
114
  snew(pme_pp,1);
115

    
116
#ifdef GMX_MPI
117
  pme_pp->mpi_comm_mysim = cr->mpi_comm_mysim;
118
  MPI_Comm_rank(cr->mpi_comm_mygroup,&rank);
119
  get_pme_ddnodes(cr,rank,&pme_pp->nnode,&pme_pp->node,&pme_pp->node_peer);
120
  snew(pme_pp->nat,pme_pp->nnode);
121
  snew(pme_pp->req,2*pme_pp->nnode);
122
  snew(pme_pp->stat,2*pme_pp->nnode);
123
  pme_pp->nalloc = 0;
124
  pme_pp->flags_charge = 0;
125
#endif
126

    
127
  return pme_pp;
128
}
129

    
130
/* This should be faster with a real non-blocking MPI implementation */
131
/* #define GMX_PME_DELAYED_WAIT */
132

    
133
static void gmx_pme_send_q_x_wait(gmx_domdec_t *dd)
134
{
135
#ifdef GMX_MPI
136
  if (dd->nreq_pme) {
137
    MPI_Waitall(dd->nreq_pme,dd->req_pme,MPI_STATUSES_IGNORE);
138
    dd->nreq_pme = 0;
139
  }
140
#endif
141
}
142

    
143
static void gmx_pme_send_q_x(t_commrec *cr, int flags,
144
                             real *chargeA, real *chargeB,
145
                             matrix box, rvec *x,
146
                             real lambda, int maxshift)
147
{
148
  gmx_domdec_t *dd;
149
  gmx_pme_comm_n_box_t *cnb;
150
  int  n;
151

    
152
  dd = cr->dd;
153
  n = dd->nat_home;
154

    
155
  if (debug)
156
    fprintf(debug,"PP node %d sending to PME node %d: %d%s%s\n",
157
            cr->sim_nodeid,dd->pme_nodeid,n,
158
            flags & PP_PME_CHARGE ? " charges" : "",
159
            flags & PP_PME_COORD  ? " coordinates" : "");
160

    
161
#ifdef GMX_PME_DELAYED_WAIT
162
  /* When can not use cnb until pending communication has finished */
163
  gmx_pme_send_x_q_wait(dd);
164
#endif
165

    
166
  if (dd->pme_receive_vir_ener) {
167
    /* Peer PP node: communicate all data */
168
    if (dd->cnb == NULL)
169
      snew(dd->cnb,1);
170
    cnb = dd->cnb;
171

    
172
    cnb->flags    = flags;
173
    cnb->natoms   = n;
174
    cnb->maxshift = maxshift;
175
    cnb->lambda   = lambda;
176
    if (flags & PP_PME_COORD)
177
      copy_mat(box,cnb->box);
178
#ifdef GMX_MPI
179
    MPI_Isend(cnb,sizeof(*cnb),MPI_BYTE,
180
              dd->pme_nodeid,0,cr->mpi_comm_mysim,
181
              &dd->req_pme[dd->nreq_pme++]);
182
#endif
183
  } else if (flags & PP_PME_CHARGE) {
184
#ifdef GMX_MPI
185
    /* Communicate only the number of atoms */
186
    MPI_Isend(&n,sizeof(n),MPI_BYTE,
187
              dd->pme_nodeid,0,cr->mpi_comm_mysim,
188
              &dd->req_pme[dd->nreq_pme++]);
189
#endif
190
  }
191

    
192
#ifdef GMX_MPI
193
  if (n > 0) {
194
    if (flags & PP_PME_CHARGE) {
195
      MPI_Isend(chargeA,n*sizeof(real),MPI_BYTE,
196
                dd->pme_nodeid,1,cr->mpi_comm_mysim,
197
                &dd->req_pme[dd->nreq_pme++]);
198
    }
199
    if (flags & PP_PME_CHARGEB) {
200
      MPI_Isend(chargeB,n*sizeof(real),MPI_BYTE,
201
                dd->pme_nodeid,2,cr->mpi_comm_mysim,
202
                &dd->req_pme[dd->nreq_pme++]);
203
    }
204
    if (flags & PP_PME_COORD) {
205
      MPI_Isend(x[0],n*sizeof(rvec),MPI_BYTE,
206
                dd->pme_nodeid,3,cr->mpi_comm_mysim,
207
                &dd->req_pme[dd->nreq_pme++]);
208
    }
209
  }
210

    
211
#ifndef GMX_PME_DELAYED_WAIT
212
  /* Wait for the data to arrive */
213
  /* We can skip this wait as we are sure x and q will not be modified
214
   * before the next call to gmx_pme_send_x_q or gmx_pme_receive_f.
215
   */
216
  gmx_pme_send_q_x_wait(dd);
217
#endif
218
#endif
219
}
220

    
221
void gmx_pme_send_q(t_commrec *cr,
222
                    bool bFreeEnergy, real *chargeA, real *chargeB,
223
                    int maxshift)
224
{
225
  int flags;
226

    
227
  flags = PP_PME_CHARGE;
228
  if (bFreeEnergy)
229
    flags |= PP_PME_CHARGEB;
230

    
231
  gmx_pme_send_q_x(cr,flags,chargeA,chargeB,NULL,NULL,0,maxshift);
232
}
233

    
234
void gmx_pme_send_x(t_commrec *cr, matrix box, rvec *x,
235
                    bool bFreeEnergy, real lambda)
236
{
237
  int flags;
238

    
239
  flags = PP_PME_COORD;
240
  if (bFreeEnergy)
241
    flags |= PP_PME_FEP;
242

    
243
  gmx_pme_send_q_x(cr,flags,NULL,NULL,box,x,lambda,0);
244
}
245

    
246
void gmx_pme_finish(t_commrec *cr)
247
{
248
  int flags;
249

    
250
  flags = PP_PME_FINISH;
251

    
252
  gmx_pme_send_q_x(cr,flags,NULL,NULL,NULL,NULL,0,0);
253
}
254

    
255
int gmx_pme_recv_q_x(struct gmx_pme_pp *pme_pp,
256
                     real **chargeA, real **chargeB,
257
                     matrix box, rvec **x,rvec **f,
258
                     int *maxshift,
259
                     bool *bFreeEnergy,real *lambda)
260
{
261
  gmx_pme_comm_n_box_t cnb;
262
  int  nat=0,q,messages,sender;
263
  real *charge_pp;
264

    
265
  messages = 0;
266

    
267
  /* avoid compiler warning about unused variable without MPI support */
268
  cnb.flags = 0;        
269
#ifdef GMX_MPI
270
  do {
271
    /* Receive the send count and box from the peer PP node */
272
    MPI_Recv(&cnb,sizeof(cnb),MPI_BYTE,
273
              pme_pp->node_peer,0,
274
              pme_pp->mpi_comm_mysim,MPI_STATUS_IGNORE);
275

    
276
    if (debug)
277
      fprintf(debug,"PME only node receiving:%s%s%s\n",
278
              (cnb.flags & PP_PME_CHARGE) ? " charges" : "",
279
              (cnb.flags & PP_PME_COORD ) ? " coordinates" : "",
280
              (cnb.flags & PP_PME_FINISH) ? " finish" : "");
281

    
282
    if (cnb.flags & PP_PME_CHARGE) {
283
      /* Receive the send counts from the other PP nodes */
284
      for(sender=0; sender<pme_pp->nnode; sender++) {
285
        if (pme_pp->node[sender] == pme_pp->node_peer) {
286
          pme_pp->nat[sender] = cnb.natoms;
287
        } else {
288
          MPI_Irecv(&(pme_pp->nat[sender]),sizeof(pme_pp->nat[0]),MPI_BYTE,
289
                    pme_pp->node[sender],0,
290
                    pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
291
        }
292
      }
293
      MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
294
      messages = 0;
295

    
296
      nat = 0;
297
      for(sender=0; sender<pme_pp->nnode; sender++)
298
        nat += pme_pp->nat[sender];
299

    
300
      if (nat > pme_pp->nalloc) {
301
        pme_pp->nalloc = over_alloc_dd(nat);
302
        srenew(pme_pp->chargeA,pme_pp->nalloc);
303
        if (cnb.flags & PP_PME_CHARGEB)
304
          srenew(pme_pp->chargeB,pme_pp->nalloc);
305
        srenew(pme_pp->x,pme_pp->nalloc);
306
        srenew(pme_pp->f,pme_pp->nalloc);
307
      }
308

    
309
      /* maxshift is sent when the charges are sent */
310
      *maxshift = cnb.maxshift;
311

    
312
      /* Receive the charges in place */
313
      for(q=0; q<((cnb.flags & PP_PME_CHARGEB) ? 2 : 1); q++) {
314
        if (q == 0)
315
          charge_pp = pme_pp->chargeA;
316
        else
317
          charge_pp = pme_pp->chargeB;
318
        nat = 0;
319
        for(sender=0; sender<pme_pp->nnode; sender++) {
320
          if (pme_pp->nat[sender] > 0) {
321
            MPI_Irecv(charge_pp+nat,pme_pp->nat[sender]*sizeof(real),MPI_BYTE,
322
                      pme_pp->node[sender],1+q,
323
                      pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
324
            nat += pme_pp->nat[sender];
325
            if (debug)
326
              fprintf(debug,"Received from PP node %d: %d charges\n",
327
                      pme_pp->node[sender],pme_pp->nat[sender]);
328
          }
329
        }
330
      }
331
      
332
      pme_pp->flags_charge = cnb.flags;
333
    }
334

    
335
    if (cnb.flags & PP_PME_COORD) {
336
      if (!(pme_pp->flags_charge & PP_PME_CHARGE))
337
        gmx_incons("PME-only node received coordinates before charges");
338

    
339
      /* The box, FE flag and lambda are sent along with the coordinates */
340
      copy_mat(cnb.box,box);
341
      *bFreeEnergy = (cnb.flags & PP_PME_FEP);
342
      *lambda      = cnb.lambda;
343

    
344
      if (*bFreeEnergy && !(pme_pp->flags_charge & PP_PME_CHARGEB))
345
        gmx_incons("PME-only node received free energy request, but did not receive B-state charges");
346
      
347
      /* Receive the coordinates in place */
348
      nat = 0;
349
      for(sender=0; sender<pme_pp->nnode; sender++) {
350
        if (pme_pp->nat[sender] > 0) {
351
          MPI_Irecv(pme_pp->x[nat],pme_pp->nat[sender]*sizeof(rvec),MPI_BYTE,
352
                    pme_pp->node[sender],3,
353
                    pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
354
          nat += pme_pp->nat[sender];
355
          if (debug)
356
              fprintf(debug,"Received from PP node %d: %d coordinates\n",
357
                      pme_pp->node[sender],pme_pp->nat[sender]);
358
        }
359
      }
360
    }
361

    
362
    /* Wait for the coordinates and/or charges to arrive */
363
    MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
364
    messages = 0;
365
  } while (!(cnb.flags & (PP_PME_COORD | PP_PME_FINISH)));
366
#endif
367

    
368
  *chargeA = pme_pp->chargeA;
369
  *chargeB = pme_pp->chargeB;
370
  *x       = pme_pp->x;
371
  *f       = pme_pp->f;
372

    
373
  return ((cnb.flags & PP_PME_FINISH) ? -1 : nat);
374
}
375

    
376
static void receive_virial_energy(t_commrec *cr,
377
                                  matrix vir,real *energy,real *dvdlambda,
378
                                  float *pme_cycles) 
379
{
380
  gmx_pme_comm_vir_ene_t cve;
381

    
382
  if (cr->dd->pme_receive_vir_ener) {
383
    if (debug)
384
      fprintf(debug,
385
              "PP node %d receiving from PME node %d: virial and energy\n",
386
              cr->sim_nodeid,cr->dd->pme_nodeid);
387
#ifdef GMX_MPI
388
    MPI_Recv(&cve,sizeof(cve),MPI_BYTE,cr->dd->pme_nodeid,1,cr->mpi_comm_mysim,
389
             MPI_STATUS_IGNORE);
390
#else
391
    memset(&cve,0,sizeof(cve));
392
#endif
393
        
394
    m_add(vir,cve.vir,vir);
395
    *energy = cve.energy;
396
    *dvdlambda += cve.dvdlambda;
397
    *pme_cycles = cve.cycles;
398

    
399
    bGotTermSignal = (cve.flags & PME_PP_TERM);
400
    bGotUsr1Signal = (cve.flags & PME_PP_USR1);
401
  } else {
402
    *energy = 0;
403
    *pme_cycles = 0;
404
  }
405
}
406

    
407
void gmx_pme_receive_f(t_commrec *cr,
408
                       rvec f[], matrix vir, 
409
                       real *energy, real *dvdlambda,
410
                       float *pme_cycles)
411
{
412
  static rvec *f_rec=NULL;
413
  static int  nalloc=0;
414
  int natoms,i;
415

    
416
#ifdef GMX_PME_DELAYED_WAIT
417
  /* Wait for the x request to finish */
418
  gmx_pme_send_q_x_wait(cr->dd);
419
#endif
420

    
421
  natoms = cr->dd->nat_home;
422

    
423
  if (natoms > nalloc) {
424
    nalloc = over_alloc_dd(natoms);
425
    srenew(f_rec,nalloc);
426
  }
427

    
428
#ifdef GMX_MPI  
429
  MPI_Recv(f_rec[0],natoms*sizeof(rvec),MPI_BYTE,
430
           cr->dd->pme_nodeid,0,cr->mpi_comm_mysim,
431
           MPI_STATUS_IGNORE);
432
#endif
433

    
434
  for(i=0; i<natoms; i++)
435
    rvec_inc(f[i],f_rec[i]);
436
  
437
  receive_virial_energy(cr,vir,energy,dvdlambda,pme_cycles);
438
}
439

    
440
void gmx_pme_send_force_vir_ener(struct gmx_pme_pp *pme_pp,
441
                                 rvec *f, matrix vir,
442
                                 real energy, real dvdlambda,
443
                                 float cycles,
444
                                 bool bGotTermSignal,
445
                                 bool bGotUsr1Signal)
446
{
447
  gmx_pme_comm_vir_ene_t cve; 
448
  int messages,ind_start,ind_end,receiver;
449

    
450
  cve.cycles = cycles;
451

    
452
  /* Now the evaluated forces have to be transferred to the PP nodes */
453
  messages = 0;
454
  ind_end = 0;
455
  for (receiver=0; receiver<pme_pp->nnode; receiver++) {
456
    ind_start = ind_end;
457
    ind_end   = ind_start + pme_pp->nat[receiver];
458
#ifdef GMX_MPI
459
    if (MPI_Isend(f[ind_start],(ind_end-ind_start)*sizeof(rvec),MPI_BYTE,
460
                  pme_pp->node[receiver],0,
461
                  pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]) != 0)
462
      gmx_comm("MPI_Isend failed in do_pmeonly");
463
#endif
464
    }
465
  
466
  /* send virial and energy to our last PP node */
467
  copy_mat(vir,cve.vir);
468
  cve.energy    = energy;
469
  cve.dvdlambda = dvdlambda;
470
  cve.flags     = 0;
471
  if (bGotTermSignal)
472
    cve.flags |= PME_PP_TERM;
473
  if (bGotUsr1Signal)
474
    cve.flags |= PME_PP_USR1;
475
  
476
  cve.cycles = cycles;
477
  
478
  if (debug)
479
    fprintf(debug,"PME node sending to PP node %d: virial and energy\n",
480
            pme_pp->node_peer);
481
#ifdef GMX_MPI
482
  MPI_Isend(&cve,sizeof(cve),MPI_BYTE,
483
            pme_pp->node_peer,1,
484
            pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
485
  
486
  /* Wait for the forces to arrive */
487
  MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
488
#endif
489
}