src/gromacs/mdlib/pme_pp.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   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 #include "sighandler.h"
  53
  54 #include "gromacs/utility/gmxmpi.h"
  55
  56 #define PP_PME_CHARGE         (1<<0)
  57 #define PP_PME_CHARGEB        (1<<1)
  58 #define PP_PME_COORD          (1<<2)
  59 #define PP_PME_FEP            (1<<3)
  60 #define PP_PME_ENER_VIR       (1<<4)
  61 #define PP_PME_FINISH         (1<<5)
  62 #define PP_PME_SWITCHGRID     (1<<6)
  63 #define PP_PME_RESETCOUNTERS  (1<<7)
  64
  65
  66 #define PME_PP_SIGSTOP     (1<<0)
  67 #define PME_PP_SIGSTOPNSS     (1<<1)
  68
  69 typedef struct gmx_pme_pp {
  70 #ifdef GMX_MPI
  71     MPI_Comm     mpi_comm_mysim;
  72 #endif
  73     int          nnode;        /* The number of PP node to communicate with  */
  74     int         *node;         /* The PP node ranks                          */
  75     int          node_peer;    /* The peer PP node rank                      */
  76     int         *nat;          /* The number of atom for each PP node        */
  77     int          flags_charge; /* The flags sent along with the last charges */
  78     real        *chargeA;
  79     real        *chargeB;
  80     rvec        *x;
  81     rvec        *f;
  82     int          nalloc;
  83 #ifdef GMX_MPI
  84     MPI_Request *req;
  85     MPI_Status  *stat;
  86 #endif
  87 } t_gmx_pme_pp;
  88
  89 typedef struct gmx_pme_comm_n_box {
  90     int             natoms;
  91     matrix          box;
  92     int             maxshift_x;
  93     int             maxshift_y;
  94     real            lambda;
  95     int             flags;
  96     gmx_large_int_t step;
  97     ivec            grid_size;  /* For PME grid tuning */
  98     real            ewaldcoeff; /* For PME grid tuning */
  99 } gmx_pme_comm_n_box_t;
 100
 101 typedef struct {
 102     matrix          vir;
 103     real            energy;
 104     real            dvdlambda;
 105     float           cycles;
 106     gmx_stop_cond_t stop_cond;
 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     {
 141         MPI_Waitall(dd->nreq_pme, dd->req_pme, MPI_STATUSES_IGNORE);
 142         dd->nreq_pme = 0;
 143     }
 144 #endif
 145 }
 146
 147 static void gmx_pme_send_q_x(t_commrec *cr, int flags,
 148                              real *chargeA, real *chargeB,
 149                              matrix box, rvec *x,
 150                              real lambda,
 151                              int maxshift_x, int maxshift_y,
 152                              gmx_large_int_t step)
 153 {
 154     gmx_domdec_t         *dd;
 155     gmx_pme_comm_n_box_t *cnb;
 156     int                   n;
 157
 158     dd = cr->dd;
 159     n  = dd->nat_home;
 160
 161     if (debug)
 162     {
 163         fprintf(debug, "PP node %d sending to PME node %d: %d%s%s\n",
 164                 cr->sim_nodeid, dd->pme_nodeid, n,
 165                 flags & PP_PME_CHARGE ? " charges" : "",
 166                 flags & PP_PME_COORD  ? " coordinates" : "");
 167     }
 168
 169 #ifdef GMX_PME_DELAYED_WAIT
 170     /* When can not use cnb until pending communication has finished */
 171     gmx_pme_send_x_q_wait(dd);
 172 #endif
 173
 174     if (dd->pme_receive_vir_ener)
 175     {
 176         /* Peer PP node: communicate all data */
 177         if (dd->cnb == NULL)
 178         {
 179             snew(dd->cnb, 1);
 180         }
 181         cnb = dd->cnb;
 182
 183         cnb->flags      = flags;
 184         cnb->natoms     = n;
 185         cnb->maxshift_x = maxshift_x;
 186         cnb->maxshift_y = maxshift_y;
 187         cnb->lambda     = lambda;
 188         cnb->step       = step;
 189         if (flags & PP_PME_COORD)
 190         {
 191             copy_mat(box, cnb->box);
 192         }
 193 #ifdef GMX_MPI
 194         MPI_Isend(cnb, sizeof(*cnb), MPI_BYTE,
 195                   dd->pme_nodeid, 0, cr->mpi_comm_mysim,
 196                   &dd->req_pme[dd->nreq_pme++]);
 197 #endif
 198     }
 199     else if (flags & PP_PME_CHARGE)
 200     {
 201 #ifdef GMX_MPI
 202         /* Communicate only the number of atoms */
 203         MPI_Isend(&n, sizeof(n), MPI_BYTE,
 204                   dd->pme_nodeid, 0, cr->mpi_comm_mysim,
 205                   &dd->req_pme[dd->nreq_pme++]);
 206 #endif
 207     }
 208
 209 #ifdef GMX_MPI
 210     if (n > 0)
 211     {
 212         if (flags & PP_PME_CHARGE)
 213         {
 214             MPI_Isend(chargeA, n*sizeof(real), MPI_BYTE,
 215                       dd->pme_nodeid, 1, cr->mpi_comm_mysim,
 216                       &dd->req_pme[dd->nreq_pme++]);
 217         }
 218         if (flags & PP_PME_CHARGEB)
 219         {
 220             MPI_Isend(chargeB, n*sizeof(real), MPI_BYTE,
 221                       dd->pme_nodeid, 2, cr->mpi_comm_mysim,
 222                       &dd->req_pme[dd->nreq_pme++]);
 223         }
 224         if (flags & PP_PME_COORD)
 225         {
 226             MPI_Isend(x[0], n*sizeof(rvec), MPI_BYTE,
 227                       dd->pme_nodeid, 3, cr->mpi_comm_mysim,
 228                       &dd->req_pme[dd->nreq_pme++]);
 229         }
 230     }
 231
 232 #ifndef GMX_PME_DELAYED_WAIT
 233     /* Wait for the data to arrive */
 234     /* We can skip this wait as we are sure x and q will not be modified
 235      * before the next call to gmx_pme_send_x_q or gmx_pme_receive_f.
 236      */
 237     gmx_pme_send_q_x_wait(dd);
 238 #endif
 239 #endif
 240 }
 241
 242 void gmx_pme_send_q(t_commrec *cr,
 243                     gmx_bool bFreeEnergy, real *chargeA, real *chargeB,
 244                     int maxshift_x, int maxshift_y)
 245 {
 246     int flags;
 247
 248     flags = PP_PME_CHARGE;
 249     if (bFreeEnergy)
 250     {
 251         flags |= PP_PME_CHARGEB;
 252     }
 253
 254     gmx_pme_send_q_x(cr, flags,
 255                      chargeA, chargeB, NULL, NULL, 0, maxshift_x, maxshift_y, -1);
 256 }
 257
 258 void gmx_pme_send_x(t_commrec *cr, matrix box, rvec *x,
 259                     gmx_bool bFreeEnergy, real lambda,
 260                     gmx_bool bEnerVir,
 261                     gmx_large_int_t step)
 262 {
 263     int flags;
 264
 265     flags = PP_PME_COORD;
 266     if (bFreeEnergy)
 267     {
 268         flags |= PP_PME_FEP;
 269     }
 270     if (bEnerVir)
 271     {
 272         flags |= PP_PME_ENER_VIR;
 273     }
 274
 275     gmx_pme_send_q_x(cr, flags, NULL, NULL, box, x, lambda, 0, 0, step);
 276 }
 277
 278 void gmx_pme_send_finish(t_commrec *cr)
 279 {
 280     int flags;
 281
 282     flags = PP_PME_FINISH;
 283
 284     gmx_pme_send_q_x(cr, flags, NULL, NULL, NULL, NULL, 0, 0, 0, -1);
 285 }
 286
 287 void gmx_pme_send_switchgrid(t_commrec *cr, ivec grid_size, real ewaldcoeff)
 288 {
 289 #ifdef GMX_MPI
 290     gmx_pme_comm_n_box_t cnb;
 291
 292     /* Only let one PP node signal each PME node */
 293     if (cr->dd->pme_receive_vir_ener)
 294     {
 295         cnb.flags = PP_PME_SWITCHGRID;
 296         copy_ivec(grid_size, cnb.grid_size);
 297         cnb.ewaldcoeff = ewaldcoeff;
 298
 299         /* We send this, uncommon, message blocking to simplify the code */
 300         MPI_Send(&cnb, sizeof(cnb), MPI_BYTE,
 301                  cr->dd->pme_nodeid, 0, cr->mpi_comm_mysim);
 302     }
 303 #endif
 304 }
 305
 306 void gmx_pme_send_resetcounters(t_commrec *cr, gmx_large_int_t step)
 307 {
 308 #ifdef GMX_MPI
 309     gmx_pme_comm_n_box_t cnb;
 310
 311     /* Only let one PP node signal each PME node */
 312     if (cr->dd->pme_receive_vir_ener)
 313     {
 314         cnb.flags = PP_PME_RESETCOUNTERS;
 315         cnb.step  = step;
 316
 317         /* We send this, uncommon, message blocking to simplify the code */
 318         MPI_Send(&cnb, sizeof(cnb), MPI_BYTE,
 319                  cr->dd->pme_nodeid, 0, cr->mpi_comm_mysim);
 320     }
 321 #endif
 322 }
 323
 324 int gmx_pme_recv_q_x(struct gmx_pme_pp *pme_pp,
 325                      int *natoms,
 326                      real **chargeA, real **chargeB,
 327                      matrix box, rvec **x, rvec **f,
 328                      int *maxshift_x, int *maxshift_y,
 329                      gmx_bool *bFreeEnergy, real *lambda,
 330                      gmx_bool *bEnerVir,
 331                      gmx_large_int_t *step,
 332                      ivec grid_size, real *ewaldcoeff)
 333 {
 334     gmx_pme_comm_n_box_t cnb;
 335     int                  nat = 0, q, messages, sender;
 336     real                *charge_pp;
 337
 338     messages = 0;
 339
 340     /* avoid compiler warning about unused variable without MPI support */
 341     cnb.flags = 0;
 342 #ifdef GMX_MPI
 343     do
 344     {
 345         /* Receive the send count, box and time step from the peer PP node */
 346         MPI_Recv(&cnb, sizeof(cnb), MPI_BYTE,
 347                  pme_pp->node_peer, 0,
 348                  pme_pp->mpi_comm_mysim, MPI_STATUS_IGNORE);
 349
 350         if (debug)
 351         {
 352             fprintf(debug, "PME only node receiving:%s%s%s%s%s\n",
 353                     (cnb.flags & PP_PME_CHARGE)        ? " charges" : "",
 354                     (cnb.flags & PP_PME_COORD )        ? " coordinates" : "",
 355                     (cnb.flags & PP_PME_FINISH)        ? " finish" : "",
 356                     (cnb.flags & PP_PME_SWITCHGRID)    ? " switch grid" : "",
 357                     (cnb.flags & PP_PME_RESETCOUNTERS) ? " reset counters" : "");
 358         }
 359
 360         if (cnb.flags & PP_PME_SWITCHGRID)
 361         {
 362             /* Special case, receive the new parameters and return */
 363             copy_ivec(cnb.grid_size, grid_size);
 364             *ewaldcoeff = cnb.ewaldcoeff;
 365
 366             return pmerecvqxSWITCHGRID;
 367         }
 368
 369         if (cnb.flags & PP_PME_RESETCOUNTERS)
 370         {
 371             /* Special case, receive the step and return */
 372             *step = cnb.step;
 373
 374             return pmerecvqxRESETCOUNTERS;
 375         }
 376
 377         if (cnb.flags & PP_PME_CHARGE)
 378         {
 379             /* Receive the send counts from the other PP nodes */
 380             for (sender = 0; sender < pme_pp->nnode; sender++)
 381             {
 382                 if (pme_pp->node[sender] == pme_pp->node_peer)
 383                 {
 384                     pme_pp->nat[sender] = cnb.natoms;
 385                 }
 386                 else
 387                 {
 388                     MPI_Irecv(&(pme_pp->nat[sender]), sizeof(pme_pp->nat[0]),
 389                               MPI_BYTE,
 390                               pme_pp->node[sender], 0,
 391                               pme_pp->mpi_comm_mysim, &pme_pp->req[messages++]);
 392                 }
 393             }
 394             MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 395             messages = 0;
 396
 397             nat = 0;
 398             for (sender = 0; sender < pme_pp->nnode; sender++)
 399             {
 400                 nat += pme_pp->nat[sender];
 401             }
 402
 403             if (nat > pme_pp->nalloc)
 404             {
 405                 pme_pp->nalloc = over_alloc_dd(nat);
 406                 srenew(pme_pp->chargeA, pme_pp->nalloc);
 407                 if (cnb.flags & PP_PME_CHARGEB)
 408                 {
 409                     srenew(pme_pp->chargeB, pme_pp->nalloc);
 410                 }
 411                 srenew(pme_pp->x, pme_pp->nalloc);
 412                 srenew(pme_pp->f, pme_pp->nalloc);
 413             }
 414
 415             /* maxshift is sent when the charges are sent */
 416             *maxshift_x = cnb.maxshift_x;
 417             *maxshift_y = cnb.maxshift_y;
 418
 419             /* Receive the charges in place */
 420             for (q = 0; q < ((cnb.flags & PP_PME_CHARGEB) ? 2 : 1); q++)
 421             {
 422                 if (q == 0)
 423                 {
 424                     charge_pp = pme_pp->chargeA;
 425                 }
 426                 else
 427                 {
 428                     charge_pp = pme_pp->chargeB;
 429                 }
 430                 nat = 0;
 431                 for (sender = 0; sender < pme_pp->nnode; sender++)
 432                 {
 433                     if (pme_pp->nat[sender] > 0)
 434                     {
 435                         MPI_Irecv(charge_pp+nat,
 436                                   pme_pp->nat[sender]*sizeof(real),
 437                                   MPI_BYTE,
 438                                   pme_pp->node[sender], 1+q,
 439                                   pme_pp->mpi_comm_mysim,
 440                                   &pme_pp->req[messages++]);
 441                         nat += pme_pp->nat[sender];
 442                         if (debug)
 443                         {
 444                             fprintf(debug, "Received from PP node %d: %d "
 445                                     "charges\n",
 446                                     pme_pp->node[sender], pme_pp->nat[sender]);
 447                         }
 448                     }
 449                 }
 450             }
 451
 452             pme_pp->flags_charge = cnb.flags;
 453         }
 454
 455         if (cnb.flags & PP_PME_COORD)
 456         {
 457             if (!(pme_pp->flags_charge & PP_PME_CHARGE))
 458             {
 459                 gmx_incons("PME-only node received coordinates before charges"
 460                            );
 461             }
 462
 463             /* The box, FE flag and lambda are sent along with the coordinates
 464              *  */
 465             copy_mat(cnb.box, box);
 466             *bFreeEnergy = (cnb.flags & PP_PME_FEP);
 467             *lambda      = cnb.lambda;
 468             *bEnerVir    = (cnb.flags & PP_PME_ENER_VIR);
 469
 470             if (*bFreeEnergy && !(pme_pp->flags_charge & PP_PME_CHARGEB))
 471             {
 472                 gmx_incons("PME-only node received free energy request, but "
 473                            "did not receive B-state charges");
 474             }
 475
 476             /* Receive the coordinates in place */
 477             nat = 0;
 478             for (sender = 0; sender < pme_pp->nnode; sender++)
 479             {
 480                 if (pme_pp->nat[sender] > 0)
 481                 {
 482                     MPI_Irecv(pme_pp->x[nat], pme_pp->nat[sender]*sizeof(rvec),
 483                               MPI_BYTE,
 484                               pme_pp->node[sender], 3,
 485                               pme_pp->mpi_comm_mysim, &pme_pp->req[messages++]);
 486                     nat += pme_pp->nat[sender];
 487                     if (debug)
 488                     {
 489                         fprintf(debug, "Received from PP node %d: %d "
 490                                 "coordinates\n",
 491                                 pme_pp->node[sender], pme_pp->nat[sender]);
 492                     }
 493                 }
 494             }
 495         }
 496
 497         /* Wait for the coordinates and/or charges to arrive */
 498         MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 499         messages = 0;
 500     }
 501     while (!(cnb.flags & (PP_PME_COORD | PP_PME_FINISH)));
 502
 503     *step = cnb.step;
 504 #endif
 505
 506     *natoms  = nat;
 507     *chargeA = pme_pp->chargeA;
 508     *chargeB = pme_pp->chargeB;
 509     *x       = pme_pp->x;
 510     *f       = pme_pp->f;
 511
 512     return ((cnb.flags & PP_PME_FINISH) ? pmerecvqxFINISH : pmerecvqxX);
 513 }
 514
 515 static void receive_virial_energy(t_commrec *cr,
 516                                   matrix vir, real *energy, real *dvdlambda,
 517                                   float *pme_cycles)
 518 {
 519     gmx_pme_comm_vir_ene_t cve;
 520
 521     if (cr->dd->pme_receive_vir_ener)
 522     {
 523         if (debug)
 524         {
 525             fprintf(debug,
 526                     "PP node %d receiving from PME node %d: virial and energy\n",
 527                     cr->sim_nodeid, cr->dd->pme_nodeid);
 528         }
 529 #ifdef GMX_MPI
 530         MPI_Recv(&cve, sizeof(cve), MPI_BYTE, cr->dd->pme_nodeid, 1, cr->mpi_comm_mysim,
 531                  MPI_STATUS_IGNORE);
 532 #else
 533         memset(&cve, 0, sizeof(cve));
 534 #endif
 535
 536         m_add(vir, cve.vir, vir);
 537         *energy     = cve.energy;
 538         *dvdlambda += cve.dvdlambda;
 539         *pme_cycles = cve.cycles;
 540
 541         if (cve.stop_cond != gmx_stop_cond_none)
 542         {
 543             gmx_set_stop_condition(cve.stop_cond);
 544         }
 545     }
 546     else
 547     {
 548         *energy     = 0;
 549         *pme_cycles = 0;
 550     }
 551 }
 552
 553 void gmx_pme_receive_f(t_commrec *cr,
 554                        rvec f[], matrix vir,
 555                        real *energy, real *dvdlambda,
 556                        float *pme_cycles)
 557 {
 558     int natoms, i;
 559
 560 #ifdef GMX_PME_DELAYED_WAIT
 561     /* Wait for the x request to finish */
 562     gmx_pme_send_q_x_wait(cr->dd);
 563 #endif
 564
 565     natoms = cr->dd->nat_home;
 566
 567     if (natoms > cr->dd->pme_recv_f_alloc)
 568     {
 569         cr->dd->pme_recv_f_alloc = over_alloc_dd(natoms);
 570         srenew(cr->dd->pme_recv_f_buf, cr->dd->pme_recv_f_alloc);
 571     }
 572
 573 #ifdef GMX_MPI
 574     MPI_Recv(cr->dd->pme_recv_f_buf[0],
 575              natoms*sizeof(rvec), MPI_BYTE,
 576              cr->dd->pme_nodeid, 0, cr->mpi_comm_mysim,
 577              MPI_STATUS_IGNORE);
 578 #endif
 579
 580     for (i = 0; i < natoms; i++)
 581     {
 582         rvec_inc(f[i], cr->dd->pme_recv_f_buf[i]);
 583     }
 584
 585
 586     receive_virial_energy(cr, vir, energy, dvdlambda, pme_cycles);
 587 }
 588
 589 void gmx_pme_send_force_vir_ener(struct gmx_pme_pp *pme_pp,
 590                                  rvec *f, matrix vir,
 591                                  real energy, real dvdlambda,
 592                                  float cycles)
 593 {
 594     gmx_pme_comm_vir_ene_t cve;
 595     int                    messages, ind_start, ind_end, receiver;
 596
 597     cve.cycles = cycles;
 598
 599     /* Now the evaluated forces have to be transferred to the PP nodes */
 600     messages = 0;
 601     ind_end  = 0;
 602     for (receiver = 0; receiver < pme_pp->nnode; receiver++)
 603     {
 604         ind_start = ind_end;
 605         ind_end   = ind_start + pme_pp->nat[receiver];
 606 #ifdef GMX_MPI
 607         if (MPI_Isend(f[ind_start], (ind_end-ind_start)*sizeof(rvec), MPI_BYTE,
 608                       pme_pp->node[receiver], 0,
 609                       pme_pp->mpi_comm_mysim, &pme_pp->req[messages++]) != 0)
 610         {
 611             gmx_comm("MPI_Isend failed in do_pmeonly");
 612         }
 613 #endif
 614     }
 615
 616     /* send virial and energy to our last PP node */
 617     copy_mat(vir, cve.vir);
 618     cve.energy    = energy;
 619     cve.dvdlambda = dvdlambda;
 620     /* check for the signals to send back to a PP node */
 621     cve.stop_cond = gmx_get_stop_condition();
 622
 623     cve.cycles = cycles;
 624
 625     if (debug)
 626     {
 627         fprintf(debug, "PME node sending to PP node %d: virial and energy\n",
 628                 pme_pp->node_peer);
 629     }
 630 #ifdef GMX_MPI
 631     MPI_Isend(&cve, sizeof(cve), MPI_BYTE,
 632               pme_pp->node_peer, 1,
 633               pme_pp->mpi_comm_mysim, &pme_pp->req[messages++]);
 634
 635     /* Wait for the forces to arrive */
 636     MPI_Waitall(messages, pme_pp->req, pme_pp->stat);
 637 #endif
 638 }