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