src/mdlib/pme_pp.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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43
  44 #include <stdio.h>
  45 #include <string.h>
  46 #include <math.h>
  47 #include "typedefs.h"
  48 #include "smalloc.h"
  49 #include "gmx_fatal.h"
  50 #include "vec.h"
  51 #include "pme.h"
  52 #include "network.h"
  53 #include "domdec.h"
  54 #include "sighandler.h"
  55
  56 #ifdef GMX_LIB_MPI
  57 #include <mpi.h>
  58 #endif
  59 #ifdef GMX_THREAD_MPI
  60 #include "tmpi.h"
  61 #endif
  62
  63 #include "mpelogging.h"
  64
  65 #define PP_PME_CHARGE   (1<<0)
  66 #define PP_PME_CHARGEB  (1<<1)
  67 #define PP_PME_COORD    (1<<2)
  68 #define PP_PME_FEP      (1<<3)
  69 #define PP_PME_ENER_VIR (1<<4)
  70 #define PP_PME_FINISH   (1<<5)
  71 #define PP_PME_SWITCH   (1<<6)
  72
  73 #define PME_PP_SIGSTOP     (1<<0)
  74 #define PME_PP_SIGSTOPNSS     (1<<1)
  75
  76 typedef struct gmx_pme_pp {
  77 #ifdef GMX_MPI
  78   MPI_Comm mpi_comm_mysim;
  79 #endif
  80   int  nnode;        /* The number of PP node to communicate with  */
  81   int  *node;        /* The PP node ranks                          */
  82   int  node_peer;    /* The peer PP node rank                      */
  83   int  *nat;         /* The number of atom for each PP node        */
  84   int  flags_charge; /* The flags sent along with the last charges */
  85   real *chargeA;
  86   real *chargeB;
  87   rvec *x;
  88   rvec *f;
  89   int  nalloc;
  90 #ifdef GMX_MPI
  91   MPI_Request *req;
  92   MPI_Status  *stat;
  93 #endif
  94 } t_gmx_pme_pp;
  95
  96 typedef struct gmx_pme_comm_n_box {
  97     int    natoms;
  98     matrix box;
  99     int    maxshift_x;
 100     int    maxshift_y;
 101     real   lambda;
 102     int    flags;
 103     gmx_large_int_t step;
 104     ivec   grid_size;     /* For PME grid tuning */
 105     real   ewaldcoeff;    /* For PME grid tuning */
 106 } gmx_pme_comm_n_box_t;
 107
 108 typedef struct {
 109   matrix vir;
 110   real   energy;
 111   real   dvdlambda;
 112   float  cycles;
 113   gmx_stop_cond_t stop_cond;
 114 } gmx_pme_comm_vir_ene_t;
 115
 116
 117
 118
 119 gmx_pme_pp_t gmx_pme_pp_init(t_commrec *cr)
 120 {
 121   struct gmx_pme_pp *pme_pp;
 122   int rank;
 123
 124   snew(pme_pp,1);
 125
 126 #ifdef GMX_MPI
 127   pme_pp->mpi_comm_mysim = cr->mpi_comm_mysim;
 128   MPI_Comm_rank(cr->mpi_comm_mygroup,&rank);
 129   get_pme_ddnodes(cr,rank,&pme_pp->nnode,&pme_pp->node,&pme_pp->node_peer);
 130   snew(pme_pp->nat,pme_pp->nnode);
 131   snew(pme_pp->req,2*pme_pp->nnode);
 132   snew(pme_pp->stat,2*pme_pp->nnode);
 133   pme_pp->nalloc = 0;
 134   pme_pp->flags_charge = 0;
 135 #endif
 136
 137   return pme_pp;
 138 }
 139
 140 /* This should be faster with a real non-blocking MPI implementation */
 141 /* #define GMX_PME_DELAYED_WAIT */
 142
 143 static void gmx_pme_send_q_x_wait(gmx_domdec_t *dd)
 144 {
 145 #ifdef GMX_MPI
 146   if (dd->nreq_pme) {
 147     MPI_Waitall(dd->nreq_pme,dd->req_pme,MPI_STATUSES_IGNORE);
 148     dd->nreq_pme = 0;
 149   }
 150 #endif
 151 }
 152
 153 static void gmx_pme_send_q_x(t_commrec *cr, int flags,
 154                              real *chargeA, real *chargeB,
 155                              matrix box, rvec *x,
 156                              real lambda,
 157                              int maxshift_x, int maxshift_y,
 158                              gmx_large_int_t step)
 159 {
 160   gmx_domdec_t *dd;
 161   gmx_pme_comm_n_box_t *cnb;
 162   int  n;
 163
 164   dd = cr->dd;
 165   n = dd->nat_home;
 166
 167   if (debug)
 168     fprintf(debug,"PP node %d sending to PME node %d: %d%s%s\n",
 169             cr->sim_nodeid,dd->pme_nodeid,n,
 170             flags & PP_PME_CHARGE ? " charges" : "",
 171             flags & PP_PME_COORD  ? " coordinates" : "");
 172
 173 #ifdef GMX_PME_DELAYED_WAIT
 174   /* When can not use cnb until pending communication has finished */
 175   gmx_pme_send_x_q_wait(dd);
 176 #endif
 177
 178   if (dd->pme_receive_vir_ener) {
 179     /* Peer PP node: communicate all data */
 180     if (dd->cnb == NULL)
 181       snew(dd->cnb,1);
 182     cnb = dd->cnb;
 183
 184     cnb->flags      = flags;
 185     cnb->natoms     = n;
 186     cnb->maxshift_x = maxshift_x;
 187     cnb->maxshift_y = maxshift_y;
 188     cnb->lambda     = lambda;
 189     cnb->step       = step;
 190     if (flags & PP_PME_COORD)
 191       copy_mat(box,cnb->box);
 192 #ifdef GMX_MPI
 193     MPI_Isend(cnb,sizeof(*cnb),MPI_BYTE,
 194               dd->pme_nodeid,0,cr->mpi_comm_mysim,
 195               &dd->req_pme[dd->nreq_pme++]);
 196 #endif
 197   } else if (flags & PP_PME_CHARGE) {
 198 #ifdef GMX_MPI
 199     /* Communicate only the number of atoms */
 200     MPI_Isend(&n,sizeof(n),MPI_BYTE,
 201               dd->pme_nodeid,0,cr->mpi_comm_mysim,
 202               &dd->req_pme[dd->nreq_pme++]);
 203 #endif
 204   }
 205
 206 #ifdef GMX_MPI
 207   if (n > 0) {
 208     if (flags & PP_PME_CHARGE) {
 209       MPI_Isend(chargeA,n*sizeof(real),MPI_BYTE,
 210                 dd->pme_nodeid,1,cr->mpi_comm_mysim,
 211                 &dd->req_pme[dd->nreq_pme++]);
 212     }
 213     if (flags & PP_PME_CHARGEB) {
 214       MPI_Isend(chargeB,n*sizeof(real),MPI_BYTE,
 215                 dd->pme_nodeid,2,cr->mpi_comm_mysim,
 216                 &dd->req_pme[dd->nreq_pme++]);
 217     }
 218     if (flags & PP_PME_COORD) {
 219       MPI_Isend(x[0],n*sizeof(rvec),MPI_BYTE,
 220                 dd->pme_nodeid,3,cr->mpi_comm_mysim,
 221                 &dd->req_pme[dd->nreq_pme++]);
 222     }
 223   }
 224
 225 #ifndef GMX_PME_DELAYED_WAIT
 226   /* Wait for the data to arrive */
 227   /* We can skip this wait as we are sure x and q will not be modified
 228    * before the next call to gmx_pme_send_x_q or gmx_pme_receive_f.
 229    */
 230   gmx_pme_send_q_x_wait(dd);
 231 #endif
 232 #endif
 233 }
 234
 235 void gmx_pme_send_q(t_commrec *cr,
 236                     gmx_bool bFreeEnergy, real *chargeA, real *chargeB,
 237                     int maxshift_x, int maxshift_y)
 238 {
 239   int flags;
 240
 241   flags = PP_PME_CHARGE;
 242   if (bFreeEnergy)
 243     flags |= PP_PME_CHARGEB;
 244
 245   gmx_pme_send_q_x(cr,flags,
 246                    chargeA,chargeB,NULL,NULL,0,maxshift_x,maxshift_y,-1);
 247 }
 248
 249 void gmx_pme_send_x(t_commrec *cr, matrix box, rvec *x,
 250                     gmx_bool bFreeEnergy, real lambda,
 251                     gmx_bool bEnerVir,
 252                     gmx_large_int_t step)
 253 {
 254   int flags;
 255
 256   flags = PP_PME_COORD;
 257   if (bFreeEnergy)
 258     flags |= PP_PME_FEP;
 259   if (bEnerVir)
 260     flags |= PP_PME_ENER_VIR;
 261
 262   gmx_pme_send_q_x(cr,flags,NULL,NULL,box,x,lambda,0,0,step);
 263 }
 264
 265 void gmx_pme_send_finish(t_commrec *cr)
 266 {
 267   int flags;
 268
 269   flags = PP_PME_FINISH;
 270
 271   gmx_pme_send_q_x(cr,flags,NULL,NULL,NULL,NULL,0,0,0,-1);
 272 }
 273
 274 void gmx_pme_send_switch(t_commrec *cr, ivec grid_size, real ewaldcoeff)
 275 {
 276 #ifdef GMX_MPI
 277     gmx_pme_comm_n_box_t cnb;
 278
 279     if (cr->dd->pme_receive_vir_ener)
 280     {
 281         cnb.flags = PP_PME_SWITCH;
 282         copy_ivec(grid_size,cnb.grid_size);
 283         cnb.ewaldcoeff = ewaldcoeff;
 284
 285         /* We send this, uncommon, message blocking to simplify the code */
 286         MPI_Send(&cnb,sizeof(cnb),MPI_BYTE,
 287                  cr->dd->pme_nodeid,0,cr->mpi_comm_mysim);
 288     }
 289 #endif
 290 }
 291
 292 int gmx_pme_recv_q_x(struct gmx_pme_pp *pme_pp,
 293                      real **chargeA, real **chargeB,
 294                      matrix box, rvec **x,rvec **f,
 295                      int *maxshift_x, int *maxshift_y,
 296                      gmx_bool *bFreeEnergy,real *lambda,
 297                      gmx_bool *bEnerVir,
 298                      gmx_large_int_t *step,
 299                      ivec grid_size, real *ewaldcoeff)
 300 {
 301     gmx_pme_comm_n_box_t cnb;
 302     int  nat=0,q,messages,sender;
 303     real *charge_pp;
 304
 305     messages = 0;
 306
 307     /* avoid compiler warning about unused variable without MPI support */
 308     cnb.flags = 0;
 309 #ifdef GMX_MPI
 310     do {
 311         /* Receive the send count, box and time step from the peer PP node */
 312         MPI_Recv(&cnb,sizeof(cnb),MPI_BYTE,
 313                  pme_pp->node_peer,0,
 314                  pme_pp->mpi_comm_mysim,MPI_STATUS_IGNORE);
 315
 316         if (debug)
 317         {
 318             fprintf(debug,"PME only node receiving:%s%s%s%s\n",
 319                     (cnb.flags & PP_PME_CHARGE) ? " charges" : "",
 320                     (cnb.flags & PP_PME_COORD ) ? " coordinates" : "",
 321                     (cnb.flags & PP_PME_FINISH) ? " finish" : "",
 322                     (cnb.flags & PP_PME_SWITCH) ? " switch" : "");
 323         }
 324
 325         if (cnb.flags & PP_PME_SWITCH)
 326         {
 327             /* Special case, receive the new parameters and return */
 328             copy_ivec(cnb.grid_size,grid_size);
 329             *ewaldcoeff = cnb.ewaldcoeff;
 330
 331             return -2;
 332         }
 333
 334         if (cnb.flags & PP_PME_CHARGE) {
 335             /* Receive the send counts from the other PP nodes */
 336             for(sender=0; sender<pme_pp->nnode; sender++) {
 337                 if (pme_pp->node[sender] == pme_pp->node_peer) {
 338                     pme_pp->nat[sender] = cnb.natoms;
 339                 } else {
 340                     MPI_Irecv(&(pme_pp->nat[sender]),sizeof(pme_pp->nat[0]),
 341                               MPI_BYTE,
 342                               pme_pp->node[sender],0,
 343                               pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
 344                 }
 345             }
 346             MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 347             messages = 0;
 348
 349             nat = 0;
 350             for(sender=0; sender<pme_pp->nnode; sender++)
 351                 nat += pme_pp->nat[sender];
 352
 353             if (nat > pme_pp->nalloc) {
 354                 pme_pp->nalloc = over_alloc_dd(nat);
 355                 srenew(pme_pp->chargeA,pme_pp->nalloc);
 356                 if (cnb.flags & PP_PME_CHARGEB)
 357                     srenew(pme_pp->chargeB,pme_pp->nalloc);
 358                 srenew(pme_pp->x,pme_pp->nalloc);
 359                 srenew(pme_pp->f,pme_pp->nalloc);
 360             }
 361
 362             /* maxshift is sent when the charges are sent */
 363             *maxshift_x = cnb.maxshift_x;
 364             *maxshift_y = cnb.maxshift_y;
 365
 366             /* Receive the charges in place */
 367             for(q=0; q<((cnb.flags & PP_PME_CHARGEB) ? 2 : 1); q++) {
 368                 if (q == 0)
 369                     charge_pp = pme_pp->chargeA;
 370                 else
 371                     charge_pp = pme_pp->chargeB;
 372                 nat = 0;
 373                 for(sender=0; sender<pme_pp->nnode; sender++) {
 374                     if (pme_pp->nat[sender] > 0) {
 375                         MPI_Irecv(charge_pp+nat,
 376                                   pme_pp->nat[sender]*sizeof(real),
 377                                   MPI_BYTE,
 378                                   pme_pp->node[sender],1+q,
 379                                   pme_pp->mpi_comm_mysim,
 380                                   &pme_pp->req[messages++]);
 381                         nat += pme_pp->nat[sender];
 382                         if (debug)
 383                             fprintf(debug,"Received from PP node %d: %d "
 384                                 "charges\n",
 385                                     pme_pp->node[sender],pme_pp->nat[sender]);
 386                     }
 387                 }
 388             }
 389
 390             pme_pp->flags_charge = cnb.flags;
 391         }
 392
 393         if (cnb.flags & PP_PME_COORD) {
 394             if (!(pme_pp->flags_charge & PP_PME_CHARGE))
 395                 gmx_incons("PME-only node received coordinates before charges"
 396                     );
 397
 398             /* The box, FE flag and lambda are sent along with the coordinates
 399              *  */
 400             copy_mat(cnb.box,box);
 401             *bFreeEnergy = (cnb.flags & PP_PME_FEP);
 402             *lambda      = cnb.lambda;
 403             *bEnerVir    = (cnb.flags & PP_PME_ENER_VIR);
 404
 405             if (*bFreeEnergy && !(pme_pp->flags_charge & PP_PME_CHARGEB))
 406                 gmx_incons("PME-only node received free energy request, but "
 407                     "did not receive B-state charges");
 408
 409             /* Receive the coordinates in place */
 410             nat = 0;
 411             for(sender=0; sender<pme_pp->nnode; sender++) {
 412                 if (pme_pp->nat[sender] > 0) {
 413                     MPI_Irecv(pme_pp->x[nat],pme_pp->nat[sender]*sizeof(rvec),
 414                               MPI_BYTE,
 415                               pme_pp->node[sender],3,
 416                               pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
 417                     nat += pme_pp->nat[sender];
 418                     if (debug)
 419                         fprintf(debug,"Received from PP node %d: %d "
 420                             "coordinates\n",
 421                                 pme_pp->node[sender],pme_pp->nat[sender]);
 422                 }
 423             }
 424         }
 425
 426         /* Wait for the coordinates and/or charges to arrive */
 427         MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 428         messages = 0;
 429     } while (!(cnb.flags & (PP_PME_COORD | PP_PME_FINISH)));
 430
 431     *step = cnb.step;
 432 #endif
 433
 434     *chargeA = pme_pp->chargeA;
 435     *chargeB = pme_pp->chargeB;
 436     *x       = pme_pp->x;
 437     *f       = pme_pp->f;
 438
 439
 440     return ((cnb.flags & PP_PME_FINISH) ? -1 : nat);
 441 }
 442
 443 static void receive_virial_energy(t_commrec *cr,
 444                                   matrix vir,real *energy,real *dvdlambda,
 445                                   float *pme_cycles)
 446 {
 447   gmx_pme_comm_vir_ene_t cve;
 448
 449   if (cr->dd->pme_receive_vir_ener) {
 450     if (debug)
 451       fprintf(debug,
 452               "PP node %d receiving from PME node %d: virial and energy\n",
 453               cr->sim_nodeid,cr->dd->pme_nodeid);
 454 #ifdef GMX_MPI
 455     MPI_Recv(&cve,sizeof(cve),MPI_BYTE,cr->dd->pme_nodeid,1,cr->mpi_comm_mysim,
 456              MPI_STATUS_IGNORE);
 457 #else
 458     memset(&cve,0,sizeof(cve));
 459 #endif
 460
 461     m_add(vir,cve.vir,vir);
 462     *energy = cve.energy;
 463     *dvdlambda += cve.dvdlambda;
 464     *pme_cycles = cve.cycles;
 465
 466     if ( cve.stop_cond != gmx_stop_cond_none )
 467     {
 468         gmx_set_stop_condition(cve.stop_cond);
 469     }
 470   } else {
 471     *energy = 0;
 472     *pme_cycles = 0;
 473   }
 474 }
 475
 476 void gmx_pme_receive_f(t_commrec *cr,
 477                        rvec f[], matrix vir,
 478                        real *energy, real *dvdlambda,
 479                        float *pme_cycles)
 480 {
 481   int natoms,i;
 482
 483 #ifdef GMX_PME_DELAYED_WAIT
 484   /* Wait for the x request to finish */
 485   gmx_pme_send_q_x_wait(cr->dd);
 486 #endif
 487
 488   natoms = cr->dd->nat_home;
 489
 490   if (natoms > cr->dd->pme_recv_f_alloc)
 491   {
 492       cr->dd->pme_recv_f_alloc = over_alloc_dd(natoms);
 493       srenew(cr->dd->pme_recv_f_buf, cr->dd->pme_recv_f_alloc);
 494   }
 495
 496 #ifdef GMX_MPI
 497   MPI_Recv(cr->dd->pme_recv_f_buf[0],
 498            natoms*sizeof(rvec),MPI_BYTE,
 499            cr->dd->pme_nodeid,0,cr->mpi_comm_mysim,
 500            MPI_STATUS_IGNORE);
 501 #endif
 502
 503   for(i=0; i<natoms; i++)
 504       rvec_inc(f[i],cr->dd->pme_recv_f_buf[i]);
 505
 506
 507   receive_virial_energy(cr,vir,energy,dvdlambda,pme_cycles);
 508 }
 509
 510 void gmx_pme_send_force_vir_ener(struct gmx_pme_pp *pme_pp,
 511                                  rvec *f, matrix vir,
 512                                  real energy, real dvdlambda,
 513                                  float cycles)
 514 {
 515   gmx_pme_comm_vir_ene_t cve;
 516   int messages,ind_start,ind_end,receiver;
 517
 518   cve.cycles = cycles;
 519
 520   /* Now the evaluated forces have to be transferred to the PP nodes */
 521   messages = 0;
 522   ind_end = 0;
 523   for (receiver=0; receiver<pme_pp->nnode; receiver++) {
 524     ind_start = ind_end;
 525     ind_end   = ind_start + pme_pp->nat[receiver];
 526 #ifdef GMX_MPI
 527     if (MPI_Isend(f[ind_start],(ind_end-ind_start)*sizeof(rvec),MPI_BYTE,
 528                   pme_pp->node[receiver],0,
 529                   pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]) != 0)
 530       gmx_comm("MPI_Isend failed in do_pmeonly");
 531 #endif
 532     }
 533
 534   /* send virial and energy to our last PP node */
 535   copy_mat(vir,cve.vir);
 536   cve.energy    = energy;
 537   cve.dvdlambda = dvdlambda;
 538   /* check for the signals to send back to a PP node */
 539   cve.stop_cond = gmx_get_stop_condition();
 540
 541   cve.cycles = cycles;
 542
 543   if (debug)
 544     fprintf(debug,"PME node sending to PP node %d: virial and energy\n",
 545             pme_pp->node_peer);
 546 #ifdef GMX_MPI
 547   MPI_Isend(&cve,sizeof(cve),MPI_BYTE,
 548             pme_pp->node_peer,1,
 549             pme_pp->mpi_comm_mysim,&pme_pp->req[messages++]);
 550
 551   /* Wait for the forces to arrive */
 552   MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 553 #endif
 554 }