src/mdlib/pme_pp.c

   1 /*
   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
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40
  41 #include <stdio.h>
  42 #include <string.h>
  43 #include <math.h>
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "gmx_fatal.h"
  47 #include "vec.h"
  48 #include "pme.h"
  49 #include "network.h"
  50 #include "domdec.h"
  51 #include "sighandler.h"
  52
  53 #ifdef GMX_LIB_MPI
  54 #include <mpi.h>
  55 #endif
  56 #ifdef GMX_THREADS
  57 #include "tmpi.h"
  58 #endif
  59
  60 #include "mpelogging.h"
  61
  62 #define PP_PME_CHARGE   (1<<0)
  63 #define PP_PME_CHARGEB  (1<<1)
  64 #define PP_PME_COORD    (1<<2)
  65 #define PP_PME_FEP      (1<<3)
  66 #define PP_PME_FINISH   (1<<4)
  67
  68 #define PME_PP_SIGSTOP     (1<<0)
  69 #define PME_PP_SIGSTOPNSS     (1<<1)
  70
  71 typedef struct gmx_pme_pp {
  72 #ifdef GMX_MPI
  73   MPI_Comm mpi_comm_mysim;
  74 #endif
  75   int  nnode;        /* The number of PP node to communicate with  */
  76   int  *node;        /* The PP node ranks                          */
  77   int  node_peer;    /* The peer PP node rank                      */
  78   int  *nat;         /* The number of atom for each PP node        */
  79   int  flags_charge; /* The flags sent along with the last charges */
  80   real *chargeA;
  81   real *chargeB;
  82   rvec *x;
  83   rvec *f;
  84   int  nalloc;
  85 #ifdef GMX_MPI
  86   MPI_Request *req;
  87   MPI_Status  *stat;
  88 #endif
  89 } t_gmx_pme_pp;
  90
  91 typedef struct gmx_pme_comm_n_box {
  92   int    natoms;
  93   matrix box;
  94   int    maxshift0;
  95   int    maxshift1;
  96   real   lambda;
  97   int    flags;
  98   gmx_large_int_t step;
  99 } gmx_pme_comm_n_box_t;
 100
 101 typedef struct {
 102   matrix vir;
 103   real   energy;
 104   real   dvdlambda;
 105   float  cycles;
 106   int    flags;
 107 } gmx_pme_comm_vir_ene_t;
 108
 109
 110
 111
 112 gmx_pme_pp_t gmx_pme_pp_init(t_commrec *cr)
 113 {
 114   struct gmx_pme_pp *pme_pp;
 115   int rank;
 116
 117   snew(pme_pp,1);
 118
 119 #ifdef GMX_MPI
 120   pme_pp->mpi_comm_mysim = cr->mpi_comm_mysim;
 121   MPI_Comm_rank(cr->mpi_comm_mygroup,&rank);
 122   get_pme_ddnodes(cr,rank,&pme_pp->nnode,&pme_pp->node,&pme_pp->node_peer);
 123   snew(pme_pp->nat,pme_pp->nnode);
 124   snew(pme_pp->req,2*pme_pp->nnode);
 125   snew(pme_pp->stat,2*pme_pp->nnode);
 126   pme_pp->nalloc = 0;
 127   pme_pp->flags_charge = 0;
 128 #endif
 129
 130   return pme_pp;
 131 }
 132
 133 /* This should be faster with a real non-blocking MPI implementation */
 134 /* #define GMX_PME_DELAYED_WAIT */
 135
 136 static void gmx_pme_send_q_x_wait(gmx_domdec_t *dd)
 137 {
 138 #ifdef GMX_MPI
 139   if (dd->nreq_pme) {
 140     MPI_Waitall(dd->nreq_pme,dd->req_pme,MPI_STATUSES_IGNORE);
 141     dd->nreq_pme = 0;
 142   }
 143 #endif
 144 }
 145
 146 static void gmx_pme_send_q_x(t_commrec *cr, int flags,
 147                              real *chargeA, real *chargeB,
 148                              matrix box, rvec *x,
 149                              real lambda,
 150                              int maxshift0, int maxshift1,
 151                              gmx_large_int_t step)
 152 {
 153   gmx_domdec_t *dd;
 154   gmx_pme_comm_n_box_t *cnb;
 155   int  n;
 156
 157   dd = cr->dd;
 158   n = dd->nat_home;
 159
 160   if (debug)
 161     fprintf(debug,"PP node %d sending to PME node %d: %d%s%s\n",
 162             cr->sim_nodeid,dd->pme_nodeid,n,
 163             flags & PP_PME_CHARGE ? " charges" : "",
 164             flags & PP_PME_COORD  ? " coordinates" : "");
 165
 166 #ifdef GMX_PME_DELAYED_WAIT
 167   /* When can not use cnb until pending communication has finished */
 168   gmx_pme_send_x_q_wait(dd);
 169 #endif
 170
 171   if (dd->pme_receive_vir_ener) {
 172     /* Peer PP node: communicate all data */
 173     if (dd->cnb == NULL)
 174       snew(dd->cnb,1);
 175     cnb = dd->cnb;
 176
 177     cnb->flags     = flags;
 178     cnb->natoms    = n;
 179     cnb->maxshift0 = maxshift0;
 180     cnb->maxshift1 = maxshift1;
 181     cnb->lambda    = lambda;
 182     cnb->step      = step;
 183     if (flags & PP_PME_COORD)
 184       copy_mat(box,cnb->box);
 185 #ifdef GMX_MPI
 186     MPI_Isend(cnb,sizeof(*cnb),MPI_BYTE,
 187               dd->pme_nodeid,0,cr->mpi_comm_mysim,
 188               &dd->req_pme[dd->nreq_pme++]);
 189 #endif
 190   } else if (flags & PP_PME_CHARGE) {
 191 #ifdef GMX_MPI
 192     /* Communicate only the number of atoms */
 193     MPI_Isend(&n,sizeof(n),MPI_BYTE,
 194               dd->pme_nodeid,0,cr->mpi_comm_mysim,
 195               &dd->req_pme[dd->nreq_pme++]);
 196 #endif
 197   }
 198
 199 #ifdef GMX_MPI
 200   if (flags & PP_PME_CHARGE) {
 201     MPI_Isend(chargeA,n*sizeof(real),MPI_BYTE,
 202               dd->pme_nodeid,1,cr->mpi_comm_mysim,
 203               &dd->req_pme[dd->nreq_pme++]);
 204   }
 205   if (flags & PP_PME_CHARGEB) {
 206     MPI_Isend(chargeB,n*sizeof(real),MPI_BYTE,
 207               dd->pme_nodeid,2,cr->mpi_comm_mysim,
 208               &dd->req_pme[dd->nreq_pme++]);
 209   }
 210   if (flags & PP_PME_COORD) {
 211     MPI_Isend(x[0],n*sizeof(rvec),MPI_BYTE,
 212               dd->pme_nodeid,3,cr->mpi_comm_mysim,
 213               &dd->req_pme[dd->nreq_pme++]);
 214   }
 215
 216 #ifndef GMX_PME_DELAYED_WAIT
 217   /* Wait for the data to arrive */
 218   /* We can skip this wait as we are sure x and q will not be modified
 219    * before the next call to gmx_pme_send_x_q or gmx_pme_receive_f.
 220    */
 221   gmx_pme_send_q_x_wait(dd);
 222 #endif
 223 #endif
 224 }
 225
 226 void gmx_pme_send_q(t_commrec *cr,
 227                     bool bFreeEnergy, real *chargeA, real *chargeB,
 228                     int maxshift0, int maxshift1)
 229 {
 230   int flags;
 231
 232   flags = PP_PME_CHARGE;
 233   if (bFreeEnergy)
 234     flags |= PP_PME_CHARGEB;
 235
 236   gmx_pme_send_q_x(cr,flags,chargeA,chargeB,NULL,NULL,0,maxshift0,maxshift1,-1);
 237 }
 238
 239 void gmx_pme_send_x(t_commrec *cr, matrix box, rvec *x,
 240                     bool bFreeEnergy, real lambda,gmx_large_int_t step)
 241 {
 242   int flags;
 243
 244   flags = PP_PME_COORD;
 245   if (bFreeEnergy)
 246     flags |= PP_PME_FEP;
 247
 248   gmx_pme_send_q_x(cr,flags,NULL,NULL,box,x,lambda,0,0,step);
 249 }
 250
 251 void gmx_pme_finish(t_commrec *cr)
 252 {
 253   int flags;
 254
 255   flags = PP_PME_FINISH;
 256
 257   gmx_pme_send_q_x(cr,flags,NULL,NULL,NULL,NULL,0,0,0,-1);
 258 }
 259
 260 int gmx_pme_recv_q_x(struct gmx_pme_pp *pme_pp,
 261                      real **chargeA, real **chargeB,
 262                      matrix box, rvec **x,rvec **f,
 263                      int *maxshift0, int *maxshift1,
 264                      bool *bFreeEnergy,real *lambda,
 265                      gmx_large_int_t *step)
 266 {
 267     gmx_pme_comm_n_box_t cnb;
 268     int  nat=0,q,messages,sender;
 269     real *charge_pp;
 270
 271     messages = 0;
 272
 273     /* avoid compiler warning about unused variable without MPI support */
 274     cnb.flags = 0;
 275 #ifdef GMX_MPI
 276     do {
 277         /* Receive the send count, box and time step from the peer PP node */
 278         MPI_Recv(&cnb,sizeof(cnb),MPI_BYTE,
 279                  pme_pp->node_peer,0,
 280                  pme_pp->mpi_comm_mysim,MPI_STATUS_IGNORE);
 281
 282         if (debug)
 283             fprintf(debug,"PME only node receiving:%s%s%s\n",
 284                     (cnb.flags & PP_PME_CHARGE) ? " charges" : "",
 285                         (cnb.flags & PP_PME_COORD ) ? " coordinates" : "",
 286                             (cnb.flags & PP_PME_FINISH) ? " finish" : "");
 287
 288         if (cnb.flags & PP_PME_CHARGE) {
 289             /* Receive the send counts from the other PP nodes */
 290             for(sender=0; sender<pme_pp->nnode; sender++) {
 291                 if (pme_pp->node[sender] == pme_pp->node_peer) {
 292                     pme_pp->nat[sender] = cnb.natoms;
 293                 } else {
 294                     MPI_Irecv(&(pme_pp->nat[sender]),sizeof(pme_pp->nat[0]),
 295                               MPI_BYTE,
 296                               pme_pp->node[sender],0,
 297                               pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
 298                 }
 299             }
 300             MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 301             messages = 0;
 302
 303             nat = 0;
 304             for(sender=0; sender<pme_pp->nnode; sender++)
 305                 nat += pme_pp->nat[sender];
 306
 307             if (nat > pme_pp->nalloc) {
 308                 pme_pp->nalloc = over_alloc_dd(nat);
 309                 srenew(pme_pp->chargeA,pme_pp->nalloc);
 310                 if (cnb.flags & PP_PME_CHARGEB)
 311                     srenew(pme_pp->chargeB,pme_pp->nalloc);
 312                 srenew(pme_pp->x,pme_pp->nalloc);
 313                 srenew(pme_pp->f,pme_pp->nalloc);
 314             }
 315
 316             /* maxshift is sent when the charges are sent */
 317             *maxshift0 = cnb.maxshift0;
 318             *maxshift1 = cnb.maxshift1;
 319
 320             /* Receive the charges in place */
 321             for(q=0; q<((cnb.flags & PP_PME_CHARGEB) ? 2 : 1); q++) {
 322                 if (q == 0)
 323                     charge_pp = pme_pp->chargeA;
 324                 else
 325                     charge_pp = pme_pp->chargeB;
 326                 nat = 0;
 327                 for(sender=0; sender<pme_pp->nnode; sender++) {
 328                     if (pme_pp->nat[sender] > 0) {
 329                         MPI_Irecv(charge_pp+nat,
 330                                   pme_pp->nat[sender]*sizeof(real),
 331                                   MPI_BYTE,
 332                                   pme_pp->node[sender],1+q,
 333                                   pme_pp->mpi_comm_mysim,
 334                                   &pme_pp->req[messages++]);
 335                         nat += pme_pp->nat[sender];
 336                         if (debug)
 337                             fprintf(debug,"Received from PP node %d: %d "
 338                                 "charges\n",
 339                                     pme_pp->node[sender],pme_pp->nat[sender]);
 340                     }
 341                 }
 342             }
 343
 344             pme_pp->flags_charge = cnb.flags;
 345         }
 346
 347         if (cnb.flags & PP_PME_COORD) {
 348             if (!(pme_pp->flags_charge & PP_PME_CHARGE))
 349                 gmx_incons("PME-only node received coordinates before charges"
 350                     );
 351
 352             /* The box, FE flag and lambda are sent along with the coordinates
 353              *  */
 354             copy_mat(cnb.box,box);
 355             *bFreeEnergy = (cnb.flags & PP_PME_FEP);
 356             *lambda      = cnb.lambda;
 357
 358             if (*bFreeEnergy && !(pme_pp->flags_charge & PP_PME_CHARGEB))
 359                 gmx_incons("PME-only node received free energy request, but "
 360                     "did not receive B-state charges");
 361
 362             /* Receive the coordinates in place */
 363             nat = 0;
 364             for(sender=0; sender<pme_pp->nnode; sender++) {
 365                 if (pme_pp->nat[sender] > 0) {
 366                     MPI_Irecv(pme_pp->x[nat],pme_pp->nat[sender]*sizeof(rvec),
 367                               MPI_BYTE,
 368                               pme_pp->node[sender],3,
 369                               pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
 370                     nat += pme_pp->nat[sender];
 371                     if (debug)
 372                         fprintf(debug,"Received from PP node %d: %d "
 373                             "coordinates\n",
 374                                 pme_pp->node[sender],pme_pp->nat[sender]);
 375                 }
 376             }
 377         }
 378
 379         /* Wait for the coordinates and/or charges to arrive */
 380         MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 381         messages = 0;
 382     } while (!(cnb.flags & (PP_PME_COORD | PP_PME_FINISH)));
 383
 384     *step = cnb.step;
 385 #endif
 386
 387     *chargeA = pme_pp->chargeA;
 388     *chargeB = pme_pp->chargeB;
 389     *x       = pme_pp->x;
 390     *f       = pme_pp->f;
 391
 392
 393     return ((cnb.flags & PP_PME_FINISH) ? -1 : nat);
 394 }
 395
 396 static void receive_virial_energy(t_commrec *cr,
 397                                   matrix vir,real *energy,real *dvdlambda,
 398                                   float *pme_cycles)
 399 {
 400   gmx_pme_comm_vir_ene_t cve;
 401
 402   if (cr->dd->pme_receive_vir_ener) {
 403     if (debug)
 404       fprintf(debug,
 405               "PP node %d receiving from PME node %d: virial and energy\n",
 406               cr->sim_nodeid,cr->dd->pme_nodeid);
 407 #ifdef GMX_MPI
 408     MPI_Recv(&cve,sizeof(cve),MPI_BYTE,cr->dd->pme_nodeid,1,cr->mpi_comm_mysim,
 409              MPI_STATUS_IGNORE);
 410 #else
 411     memset(&cve,0,sizeof(cve));
 412 #endif
 413
 414     m_add(vir,cve.vir,vir);
 415     *energy = cve.energy;
 416     *dvdlambda += cve.dvdlambda;
 417     *pme_cycles = cve.cycles;
 418
 419     bGotStopNextStepSignal = (cve.flags & PME_PP_SIGSTOP);
 420     bGotStopNextNSStepSignal = (cve.flags & PME_PP_SIGSTOPNSS);
 421   } else {
 422     *energy = 0;
 423     *pme_cycles = 0;
 424   }
 425 }
 426
 427 void gmx_pme_receive_f(t_commrec *cr,
 428                        rvec f[], matrix vir,
 429                        real *energy, real *dvdlambda,
 430                        float *pme_cycles)
 431 {
 432   int natoms,i;
 433
 434 #ifdef GMX_PME_DELAYED_WAIT
 435   /* Wait for the x request to finish */
 436   gmx_pme_send_q_x_wait(cr->dd);
 437 #endif
 438
 439   natoms = cr->dd->nat_home;
 440
 441   if (natoms > cr->dd->pme_recv_f_alloc)
 442   {
 443       cr->dd->pme_recv_f_alloc = over_alloc_dd(natoms);
 444       srenew(cr->dd->pme_recv_f_buf, cr->dd->pme_recv_f_alloc);
 445   }
 446
 447 #ifdef GMX_MPI
 448   MPI_Recv(cr->dd->pme_recv_f_buf[0],
 449            natoms*sizeof(rvec),MPI_BYTE,
 450            cr->dd->pme_nodeid,0,cr->mpi_comm_mysim,
 451            MPI_STATUS_IGNORE);
 452 #endif
 453
 454   for(i=0; i<natoms; i++)
 455       rvec_inc(f[i],cr->dd->pme_recv_f_buf[i]);
 456
 457
 458   receive_virial_energy(cr,vir,energy,dvdlambda,pme_cycles);
 459 }
 460
 461 void gmx_pme_send_force_vir_ener(struct gmx_pme_pp *pme_pp,
 462                                  rvec *f, matrix vir,
 463                                  real energy, real dvdlambda,
 464                                  float cycles,
 465                                  bool bGotStopNextStepSignal,
 466                                  bool bGotStopNextNSStepSignal)
 467 {
 468   gmx_pme_comm_vir_ene_t cve;
 469   int messages,ind_start,ind_end,receiver;
 470
 471   cve.cycles = cycles;
 472
 473   /* Now the evaluated forces have to be transferred to the PP nodes */
 474   messages = 0;
 475   ind_end = 0;
 476   for (receiver=0; receiver<pme_pp->nnode; receiver++) {
 477     ind_start = ind_end;
 478     ind_end   = ind_start + pme_pp->nat[receiver];
 479 #ifdef GMX_MPI
 480     if (MPI_Isend(f[ind_start],(ind_end-ind_start)*sizeof(rvec),MPI_BYTE,
 481                   pme_pp->node[receiver],0,
 482                   pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]) != 0)
 483       gmx_comm("MPI_Isend failed in do_pmeonly");
 484 #endif
 485     }
 486
 487   /* send virial and energy to our last PP node */
 488   copy_mat(vir,cve.vir);
 489   cve.energy    = energy;
 490   cve.dvdlambda = dvdlambda;
 491   cve.flags     = 0;
 492   if (bGotStopNextStepSignal)
 493     cve.flags |= PME_PP_SIGSTOP;
 494   if (bGotStopNextNSStepSignal)
 495     cve.flags |= PME_PP_SIGSTOPNSS;
 496
 497   cve.cycles = cycles;
 498
 499   if (debug)
 500     fprintf(debug,"PME node sending to PP node %d: virial and energy\n",
 501             pme_pp->node_peer);
 502 #ifdef GMX_MPI
 503   MPI_Isend(&cve,sizeof(cve),MPI_BYTE,
 504             pme_pp->node_peer,1,
 505             pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
 506
 507   /* Wait for the forces to arrive */
 508   MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 509 #endif
 510 }