src/gromacs/imd/imd.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  *
  38  * \brief
  39  * Implements functions of imd.h.
  40  *
  41  * Re-implementation of basic IMD functions from NAMD/VMD from scratch,
  42  * see imdsocket.h for references to the IMD API.
  43  *
  44  * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
  45  *
  46  * \ingroup module_imd
  47  */
  48
  49 #ifdef HAVE_CONFIG_H
  50 #include <config.h>
  51 #endif
  52
  53
  54 #include <string.h>
  55
  56 #ifdef GMX_NATIVE_WINDOWS
  57 #include <windows.h>
  58 #else
  59 #include <unistd.h>
  60 #endif
  61
  62 #include "imd.h"
  63 #include "imdsocket.h"
  64 #include "gromacs/utility/smalloc.h"
  65 #include "network.h"
  66 #include "mdrun.h"
  67 #include "sighandler.h"
  68 #include "gmx_ga2la.h"
  69 #include "xvgr.h"
  70 #include "gromacs/mdlib/groupcoord.h"
  71 #include "gromacs/fileio/confio.h"
  72 #include "mtop_util.h"
  73 #include "names.h"
  74 #include "gromacs/timing/wallcycle.h"
  75
  76 #include "gmx_fatal.h"
  77
  78 /*! \brief How long shall we wait in seconds until we check for a connection again? */
  79 #define IMDLOOPWAIT 1
  80
  81 /*! \brief How long shall we check for the IMD_GO? */
  82 #define IMDCONNECTWAIT 2
  83
  84 /*! \brief IMD Header Size. */
  85 #define HEADERSIZE 8
  86 /*! \brief IMD Protocol Version. */
  87 #define IMDVERSION 2
  88
  89 /*! \brief Broadcast d to all nodes */
  90 #define  block_bc(cr, d) gmx_bcast(sizeof(d), &(d), (cr))
  91
  92 /*! \brief Broadcast nr elements of d to all nodes */
  93 #define  nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
  94
  95
  96 /*! \internal
  97  * \brief
  98  * IMD (interactive molecular dynamics) energy record.
  99  *
 100  * As in the original IMD implementation. Energies in kcal/mol.
 101  * NOTE: We return the energies in GROMACS / SI units,
 102  * so they also show up as SI in VMD.
 103  *
 104  */
 105 typedef struct
 106 {
 107     gmx_int32_t tstep;     /**< time step                                     */
 108     float       T_abs;     /**< absolute temperature                          */
 109     float       E_tot;     /**< total energy                                  */
 110     float       E_pot;     /**< potential energy                              */
 111     float       E_vdw;     /**< van der Waals energy                          */
 112     float       E_coul;    /**< Coulomb interaction energy                    */
 113     float       E_bond;    /**< bonds energy                                  */
 114     float       E_angle;   /**< angles energy                                 */
 115     float       E_dihe;    /**< dihedrals energy                              */
 116     float       E_impr;    /**< improper dihedrals energy                     */
 117 } IMDEnergyBlock;
 118
 119
 120 /*! \internal
 121  * \brief IMD (interactive molecular dynamics) communication structure.
 122  *
 123  * This structure defines the IMD communication message header & protocol version.
 124  */
 125 typedef struct
 126 {
 127     gmx_int32_t type;      /**< Type of IMD message, see IMDType_t above      */
 128     gmx_int32_t length;    /**< Length                                        */
 129 } IMDHeader;
 130
 131
 132 /*! \internal
 133  * \brief IMD (interactive molecular dynamics) main data structure.
 134  *
 135  * Contains private IMD data
 136  */
 137 typedef struct gmx_IMD
 138 {
 139     FILE      *outf;                 /**< Output file for IMD data, mainly forces.    */
 140
 141     int        nat;                  /**< Number of atoms that can be pulled via IMD. */
 142     int        nat_loc;              /**< Part of the atoms that are local.           */
 143     atom_id   *ind;                  /**< Global indices of the IMD atoms.            */
 144     atom_id   *ind_loc;              /**< Local indices of the IMD atoms.             */
 145     int        nalloc_loc;           /**< Allocation size for ind_loc.                */
 146     rvec      *xa;                   /**< Positions for all IMD atoms assembled on
 147                                           the master node.                            */
 148     ivec      *xa_shifts;            /**< Shifts for all IMD atoms, to make
 149                                           molecule(s) whole.                          */
 150     ivec      *xa_eshifts;           /**< Extra shifts since last DD step.            */
 151     rvec      *xa_old;               /**< Old positions for all IMD atoms on master.  */
 152     int       *xa_ind;               /**< Position of each local atom in the
 153                                           collective array.                           */
 154
 155     int             nstimd;          /**< Global IMD frequency, known to all nodes.   */
 156     int             nstimd_new;      /**< New frequency from IMD client, master only. */
 157     int             nstimd_def;      /**< Default IMD frequency when disconnected.    */
 158
 159     int             port;            /**< Port to use for network socket.             */
 160     IMDSocket      *socket;          /**< The IMD socket on the master node.          */
 161     IMDSocket      *clientsocket;    /**< The IMD socket on the client.               */
 162     int             length;          /**< Length we got with last header.             */
 163
 164     gmx_bool        bWConnect;       /**< Shall we block and wait for connection?     */
 165     gmx_bool        bTerminated;     /**< Set if MD is terminated.                    */
 166     gmx_bool        bTerminatable;   /**< Set if MD can be terminated.                */
 167     gmx_bool        bConnected;      /**< Set if connection is present.               */
 168     gmx_bool        bNewForces;      /**< Set if we received new forces.              */
 169     gmx_bool        bForceActivated; /**< Set if pulling from VMD is allowed.         */
 170
 171     IMDEnergyBlock *energies;        /**< Pointer to energies we send back.           */
 172
 173     gmx_int32_t     vmd_nforces;     /**< Number of VMD forces.                       */
 174     gmx_int32_t    *vmd_f_ind;       /**< VMD forces indices.                         */
 175     float          *vmd_forces;      /**< The VMD forces flat in memory.              */
 176     int             nforces;         /**< Number of actual MD forces;
 177                                           this gets communicated to the clients.      */
 178     atom_id        *f_ind;           /**< Force indices.                              */
 179     rvec           *f;               /**< The IMD pulling forces.                     */
 180
 181     char           *forcesendbuf;    /**< Buffer for force sending.                   */
 182     char           *coordsendbuf;    /**< Buffer for coordinate sending.              */
 183     char           *energysendbuf;   /**< Send buffer for energies.                   */
 184     rvec           *sendxbuf;        /**< Buffer to make molecules whole before
 185                                           sending.                                    */
 186
 187     t_block         mols;            /**< Molecules block in IMD group.               */
 188
 189     /* The next block is used on the master node only to reduce the output
 190      * without sacrificing information. If any of these values changes,
 191      * we need to write output */
 192     int       old_nforces;           /**< Old value for nforces.                      */
 193     atom_id  *old_f_ind;             /**< Old values for force indices.               */
 194     rvec     *old_forces;            /**< Old values for IMD pulling forces.          */
 195
 196 } t_gmx_IMD_setup;
 197
 198
 199 /*! \internal
 200  * \brief Enum for types of IMD messages.
 201  *
 202  * We use the same records as the NAMD/VMD IMD implementation.
 203  */
 204 typedef enum IMDType_t
 205 {
 206     IMD_DISCONNECT,      /**< client disconnect                               */
 207     IMD_ENERGIES,        /**< energy data                                     */
 208     IMD_FCOORDS,         /**< atomic coordinates                              */
 209     IMD_GO,              /**< start command for the simulation                */
 210     IMD_HANDSHAKE,       /**< handshake to determine little/big endianness    */
 211     IMD_KILL,            /**< terminates the simulation                       */
 212     IMD_MDCOMM,          /**< force data                                      */
 213     IMD_PAUSE,           /**< pauses the simulation                           */
 214     IMD_TRATE,           /**< sets the IMD transmission and processing rate   */
 215     IMD_IOERROR,         /**< I/O error                                       */
 216     IMD_NR               /**< number of entries                               */
 217 } IMDMessageType;
 218
 219
 220 /*! \internal
 221  * \brief Names of the IMDType for error messages.
 222  */
 223 const char *eIMDType_names[IMD_NR + 1] = {
 224     "IMD_DISCONNECT",
 225     "IMD_ENERGIES",
 226     "IMD_FCOORDS",
 227     "IMD_GO",
 228     "IMD_HANDSHAKE",
 229     "IMD_KILL",
 230     "IMD_MDCOMM",
 231     "IMD_PAUSE",
 232     "IMD_TRATE",
 233     "IMD_IOERROR",
 234     NULL
 235 };
 236
 237
 238 #ifdef GMX_IMD
 239
 240 /*! \brief Fills the header with message and the length argument. */
 241 static void fill_header(IMDHeader *header, IMDMessageType type, gmx_int32_t length)
 242 {
 243     /* We (ab-)use htonl network function for the correct endianness */
 244     header->type   = htonl((gmx_int32_t) type);
 245     header->length = htonl(length);
 246 }
 247
 248
 249 /*! \brief Swaps the endianess of the header. */
 250 static void swap_header(IMDHeader *header)
 251 {
 252     /* and vice versa... */
 253     header->type   = ntohl(header->type);
 254     header->length = ntohl(header->length);
 255 }
 256
 257
 258 /*! \brief Reads multiple bytes from socket. */
 259 static gmx_int32_t imd_read_multiple(IMDSocket *socket, char *datptr, gmx_int32_t toread)
 260 {
 261     gmx_int32_t leftcount, countread;
 262
 263
 264     leftcount = toread;
 265     /* Try to read while we haven't reached toread */
 266     while (leftcount != 0)
 267     {
 268         if ((countread = imdsock_read(socket, datptr, leftcount)) < 0)
 269         {
 270             /* interrupted function call, try again... */
 271             if (errno == EINTR)
 272             {
 273                 countread = 0;
 274             }
 275             /* this is an unexpected error, return what we got */
 276             else
 277             {
 278                 return toread - leftcount;
 279             }
 280
 281             /* nothing read, finished */
 282         }
 283         else if (countread == 0)
 284         {
 285             break;
 286         }
 287         leftcount -= countread;
 288         datptr    += countread;
 289     } /* end while */
 290
 291     /* return nr of bytes read */
 292     return toread - leftcount;
 293 }
 294
 295
 296 /*! \brief Writes multiple bytes to socket in analogy to imd_read_multiple. */
 297 static gmx_int32_t imd_write_multiple(IMDSocket *socket, const char *datptr, gmx_int32_t towrite)
 298 {
 299     gmx_int32_t leftcount, countwritten;
 300
 301
 302     leftcount = towrite;
 303     while (leftcount != 0)
 304     {
 305         if ((countwritten = imdsock_write(socket, datptr, leftcount)) <= 0)
 306         {
 307             if (errno == EINTR)
 308             {
 309                 countwritten = 0;
 310             }
 311             else
 312             {
 313                 return towrite - leftcount;
 314             }
 315         }
 316         leftcount -= countwritten;
 317         datptr    += countwritten;
 318     } /* end while */
 319
 320     return towrite - leftcount;
 321 }
 322
 323
 324 /*! \brief Handshake with IMD client. */
 325 static int imd_handshake(IMDSocket *socket)
 326 {
 327     IMDHeader header;
 328
 329
 330     fill_header(&header, IMD_HANDSHAKE, 1);
 331     header.length = IMDVERSION; /* client wants unswapped version */
 332
 333     return (imd_write_multiple(socket, (char *) &header, HEADERSIZE) != HEADERSIZE);
 334 }
 335
 336
 337 /*! \brief Send energies using the energy block and the send buffer. */
 338 static int imd_send_energies(IMDSocket *socket, const IMDEnergyBlock *energies, char *buffer)
 339 {
 340     gmx_int32_t recsize;
 341
 342
 343     recsize = HEADERSIZE + sizeof(IMDEnergyBlock);
 344     fill_header((IMDHeader *) buffer, IMD_ENERGIES, 1);
 345     memcpy(buffer + HEADERSIZE, energies, sizeof(IMDEnergyBlock));
 346
 347     return (imd_write_multiple(socket, buffer, recsize) != recsize);
 348 }
 349
 350
 351 /*! \brief Receive IMD header from socket, sets the length and returns the IMD message. */
 352 static IMDMessageType imd_recv_header(IMDSocket *socket, gmx_int32_t *length)
 353 {
 354     IMDHeader header;
 355
 356
 357     if (imd_read_multiple(socket, (char *) &header, HEADERSIZE) != HEADERSIZE)
 358     {
 359         return IMD_IOERROR;
 360     }
 361     swap_header(&header);
 362     *length = header.length;
 363
 364     return (IMDMessageType) header.type;
 365 }
 366
 367
 368 /*! \brief Receive force indices and forces.
 369  *
 370  * The number of forces was previously communicated via the header.
 371  */
 372 static int imd_recv_mdcomm(IMDSocket *socket, gmx_int32_t nforces, gmx_int32_t *forcendx, float *forces)
 373 {
 374     int retsize, retbytes;
 375
 376
 377     /* reading indices */
 378     retsize  = sizeof(gmx_int32_t) * nforces;
 379     retbytes = imd_read_multiple(socket, (char *) forcendx, retsize);
 380     if (retbytes != retsize)
 381     {
 382         return FALSE;
 383     }
 384
 385     /* reading forces as float array */
 386     retsize  = 3 * sizeof(float) * nforces;
 387     retbytes = imd_read_multiple(socket, (char *) forces, retsize);
 388     if (retbytes != retsize)
 389     {
 390         return FALSE;
 391     }
 392
 393     return TRUE;
 394 }
 395
 396 #endif
 397
 398 /* GROMACS specific functions for the IMD implementation */
 399
 400 extern void write_IMDgroup_to_file(gmx_bool bIMD, t_inputrec *ir, t_state *state,
 401                                    gmx_mtop_t *sys, int nfile, const t_filenm fnm[])
 402 {
 403     t_atoms IMDatoms;
 404
 405
 406     if (bIMD)
 407     {
 408         IMDatoms = gmx_mtop_global_atoms(sys);
 409         write_sto_conf_indexed(opt2fn("-imd", nfile, fnm), "IMDgroup", &IMDatoms,
 410                                state->x, state->v, ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
 411     }
 412 }
 413
 414
 415 extern void dd_make_local_IMD_atoms(gmx_bool bIMD, gmx_domdec_t *dd, t_IMD *imd)
 416 {
 417     gmx_ga2la_t         ga2la;
 418     t_gmx_IMD_setup    *IMDsetup;
 419
 420
 421     if (bIMD)
 422     {
 423         IMDsetup = imd->setup;
 424         ga2la    = dd->ga2la;
 425
 426         dd_make_local_group_indices(
 427                 ga2la, IMDsetup->nat, IMDsetup->ind, &IMDsetup->nat_loc,
 428                 &IMDsetup->ind_loc, &IMDsetup->nalloc_loc, IMDsetup->xa_ind);
 429     }
 430 }
 431
 432
 433 #ifdef GMX_IMD
 434 /*! \brief Send positions from rvec.
 435  *
 436  * We need a separate send buffer and conversion to Angstrom.
 437  */
 438 static int imd_send_rvecs(IMDSocket *socket, int nat, rvec *x, char *buffer)
 439 {
 440     gmx_int32_t size;
 441     int         i;
 442     float       sendx[3];
 443     int         tuplesize = 3 * sizeof(float);
 444
 445
 446     /* Required size for the send buffer */
 447     size = HEADERSIZE + 3 * sizeof(float) * nat;
 448
 449     /* Prepare header */
 450     fill_header((IMDHeader *) buffer, IMD_FCOORDS, (gmx_int32_t) nat);
 451     for (i = 0; i < nat; i++)
 452     {
 453         sendx[0] = (float) x[i][0] * NM2A;
 454         sendx[1] = (float) x[i][1] * NM2A;
 455         sendx[2] = (float) x[i][2] * NM2A;
 456         memcpy(buffer + HEADERSIZE + i * tuplesize, sendx, tuplesize);
 457     }
 458
 459     return (imd_write_multiple(socket, buffer, size) != size);
 460 }
 461
 462
 463 /*! \brief Initializes the IMD private data. */
 464 static t_gmx_IMD_setup* imd_create(int imdatoms, int nstimddef, int imdport)
 465 {
 466     t_gmx_IMD_setup *IMDsetup = NULL;
 467
 468
 469     snew(IMDsetup, 1);
 470     IMDsetup->nat             = imdatoms;
 471     IMDsetup->bTerminated     = FALSE;
 472     IMDsetup->bTerminatable   = FALSE;
 473     IMDsetup->bWConnect       = FALSE;
 474     IMDsetup->bConnected      = FALSE;
 475     IMDsetup->bForceActivated = FALSE;
 476     IMDsetup->bNewForces      = FALSE;
 477     IMDsetup->bForceActivated = FALSE;
 478     IMDsetup->nstimd          = 1;
 479     IMDsetup->nstimd_new      = 1;
 480     IMDsetup->nstimd_def      = nstimddef;
 481     if (imdport < 1)
 482     {
 483         IMDsetup->port        = 0;
 484         fprintf(stderr, "%s You chose a port number < 1. Will automatically assign a free port.\n", IMDstr);
 485     }
 486     else
 487     {
 488         IMDsetup->port        = imdport;
 489     }
 490
 491     return IMDsetup;
 492 }
 493
 494
 495 /*! \brief Prepare the socket on the MASTER. */
 496 static void imd_prepare_master_socket(t_gmx_IMD_setup *IMDsetup)
 497 {
 498     int ret;
 499
 500
 501 #ifdef GMX_NATIVE_WINDOWS
 502     /* Winsock requires separate initialization */
 503     fprintf(stderr, "%s Initializing winsock.\n", IMDstr);
 504 #ifdef GMX_HAVE_WINSOCK
 505     if (imdsock_winsockinit())
 506     {
 507         gmx_fatal(FARGS, "%s Failed to initialize winsock.\n", IMDstr);
 508     }
 509 #endif
 510 #endif
 511
 512     /* The rest is identical, first create and bind a socket and set to listen then. */
 513     fprintf(stderr, "%s Setting up incoming socket.\n", IMDstr);
 514     IMDsetup->socket = imdsock_create();
 515     if (!IMDsetup->socket)
 516     {
 517         gmx_fatal(FARGS, "%s Failed to create socket.", IMDstr);
 518     }
 519
 520     /* Bind to port */
 521     ret = imdsock_bind(IMDsetup->socket, IMDsetup->port);
 522     if (ret)
 523     {
 524         gmx_fatal(FARGS, "%s binding socket to port %d failed with error %d.\n", IMDstr, IMDsetup->port, ret);
 525     }
 526
 527     if (imd_sock_listen(IMDsetup->socket))
 528     {
 529         gmx_fatal(FARGS, "%s socket listen failed with error %d.\n", IMDstr, ret);
 530     }
 531
 532     if (imdsock_getport(IMDsetup->socket, &IMDsetup->port))
 533     {
 534         gmx_fatal(FARGS, "%s Could not determine port number.\n", IMDstr, ret);
 535     }
 536
 537     fprintf(stderr, "%s Listening for IMD connection on port %d.\n", IMDstr, IMDsetup->port);
 538 }
 539
 540
 541 /*! \brief Disconnect the client. */
 542 static void imd_disconnect(t_gmx_IMD_setup *IMDsetup)
 543 {
 544     /* Write out any buffered pulling data */
 545     fflush(IMDsetup->outf);
 546
 547     /* we first try to shut down the clientsocket */
 548     imdsock_shutdown(IMDsetup->clientsocket);
 549     if (!imdsock_destroy(IMDsetup->clientsocket))
 550     {
 551         fprintf(stderr, "%s Failed to destroy socket.\n", IMDstr);
 552     }
 553
 554     /* then we reset the IMD step to its default, and reset the connection boolean */
 555     IMDsetup->nstimd_new   = IMDsetup->nstimd_def;
 556     IMDsetup->clientsocket = NULL;
 557     IMDsetup->bConnected   = FALSE;
 558 }
 559
 560
 561 /*! \brief Prints an error message and disconnects the client.
 562  *
 563  *  Does not terminate mdrun!
 564  */
 565 static void imd_fatal(t_gmx_IMD_setup *IMDsetup, const char *msg)
 566 {
 567     fprintf(stderr, "%s %s", IMDstr, msg);
 568     imd_disconnect(IMDsetup);
 569     fprintf(stderr, "%s disconnected.\n", IMDstr);
 570 }
 571
 572
 573 /*! \brief Check whether we got an incoming connection. */
 574 static gmx_bool imd_tryconnect(t_gmx_IMD_setup *IMDsetup)
 575 {
 576     if (imdsock_tryread(IMDsetup->socket, 0, 0) > 0)
 577     {
 578         /* yes, we got something, accept on clientsocket */
 579         IMDsetup->clientsocket = imdsock_accept(IMDsetup->socket);
 580         if (!IMDsetup->clientsocket)
 581         {
 582             fprintf(stderr, "%s Accepting the connection on the socket failed.\n", IMDstr);
 583             return FALSE;
 584         }
 585
 586         /* handshake with client */
 587         if (imd_handshake(IMDsetup->clientsocket))
 588         {
 589             imd_fatal(IMDsetup, "Connection failed.\n");
 590             return FALSE;
 591         }
 592
 593         fprintf(stderr, "%s Connection established, checking if I got IMD_GO orders.\n", IMDstr);
 594
 595         /* Check if we get the proper "GO" command from client. */
 596         if (imdsock_tryread(IMDsetup->clientsocket, IMDCONNECTWAIT, 0) != 1 || imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length)) != IMD_GO)
 597         {
 598             imd_fatal(IMDsetup, "No IMD_GO order received. IMD connection failed.\n");
 599         }
 600
 601         /* IMD connected */
 602         IMDsetup->bConnected = TRUE;
 603
 604         return TRUE;
 605     }
 606
 607     return FALSE;
 608 }
 609
 610
 611 /*! \brief Wrap imd_tryconnect in order to make it blocking.
 612  *
 613  * Used when the simulation should wait for an incoming connection.
 614  */
 615 static void imd_blockconnect(t_gmx_IMD_setup *IMDsetup)
 616 {
 617     /* do not wait for connection, when e.g. ctrl+c is pressed and we will terminate anyways. */
 618     if (!(int) gmx_get_stop_condition() == gmx_stop_cond_none)
 619     {
 620         return;
 621     }
 622
 623     fprintf(stderr, "%s Will wait until I have a connection and IMD_GO orders.\n", IMDstr);
 624
 625     /* while we have no clientsocket... 2nd part: we should still react on ctrl+c */
 626     while ((!IMDsetup->clientsocket) && ((int) gmx_get_stop_condition() == gmx_stop_cond_none))
 627     {
 628         imd_tryconnect(IMDsetup);
 629 #ifdef GMX_NATIVE_WINDOWS
 630         /* for whatever reason, it is called Sleep on windows */
 631         Sleep(IMDLOOPWAIT);
 632 #else
 633         sleep(IMDLOOPWAIT);
 634 #endif
 635     }
 636 }
 637
 638
 639 /*! \brief Make sure that our array holding the forces received via IMD is large enough. */
 640 static void imd_prepare_vmd_Forces(t_gmx_IMD_setup *IMDsetup)
 641 {
 642     srenew((IMDsetup->vmd_f_ind), IMDsetup->vmd_nforces);
 643     srenew((IMDsetup->vmd_forces), 3*IMDsetup->vmd_nforces);
 644 }
 645
 646
 647 /*! \brief Reads forces received via IMD. */
 648 static void imd_read_vmd_Forces(t_gmx_IMD_setup *IMDsetup)
 649 {
 650     /* the length of the previously received header tells us the nr of forces we will receive */
 651     IMDsetup->vmd_nforces = IMDsetup->length;
 652     /* prepare the arrays */
 653     imd_prepare_vmd_Forces(IMDsetup);
 654     /* Now we read the forces... */
 655     if (!(imd_recv_mdcomm(IMDsetup->clientsocket, IMDsetup->vmd_nforces, IMDsetup->vmd_f_ind, IMDsetup->vmd_forces)))
 656     {
 657         imd_fatal(IMDsetup, "Error while reading forces from remote. Disconnecting\n");
 658     }
 659 }
 660
 661
 662 /*! \brief Prepares the MD force arrays. */
 663 static void imd_prepare_MD_Forces(t_gmx_IMD_setup *IMDsetup)
 664 {
 665     srenew((IMDsetup->f_ind), IMDsetup->nforces);
 666     srenew((IMDsetup->f    ), IMDsetup->nforces);
 667 }
 668
 669
 670 /*! \brief Copy IMD forces to MD forces.
 671  *
 672  * Do conversion from Cal->Joule and from
 673  * Angstrom -> nm and from a pointer array to arrays to 3*N array.
 674  */
 675 static void imd_copyto_MD_Forces(t_gmx_IMD_setup *IMDsetup)
 676 {
 677     int  i;
 678     real conversion = CAL2JOULE * NM2A;
 679
 680
 681     for (i = 0; i < IMDsetup->nforces; i++)
 682     {
 683         /* Copy the indices, a copy is important because we may update the incoming forces
 684          * whenever we receive new forces while the MD forces are only communicated upon
 685          * IMD communication.*/
 686         IMDsetup->f_ind[i] = IMDsetup->vmd_f_ind[i];
 687
 688         /* Convert to rvecs and do a proper unit conversion */
 689         IMDsetup->f[i][0] = IMDsetup->vmd_forces[3*i    ] * conversion;
 690         IMDsetup->f[i][1] = IMDsetup->vmd_forces[3*i + 1] * conversion;
 691         IMDsetup->f[i][2] = IMDsetup->vmd_forces[3*i + 2] * conversion;
 692     }
 693 }
 694
 695
 696 /*! \brief Return TRUE if any of the forces or indices changed. */
 697 static gmx_bool bForcesChanged(t_gmx_IMD_setup *IMDsetup)
 698 {
 699     int i;
 700
 701
 702     /* First, check whether the number of pulled atoms changed */
 703     if (IMDsetup->nforces != IMDsetup->old_nforces)
 704     {
 705         return TRUE;
 706     }
 707
 708     /* Second, check whether any of the involved atoms changed */
 709     for (i = 0; i < IMDsetup->nforces; i++)
 710     {
 711         if (IMDsetup->f_ind[i] != IMDsetup->old_f_ind[i])
 712         {
 713             return TRUE;
 714         }
 715     }
 716
 717     /* Third, check whether all forces are the same */
 718     for (i = 0; i < IMDsetup->nforces; i++)
 719     {
 720         if (IMDsetup->f[i][XX] != IMDsetup->old_forces[i][XX])
 721         {
 722             return TRUE;
 723         }
 724         if (IMDsetup->f[i][YY] != IMDsetup->old_forces[i][YY])
 725         {
 726             return TRUE;
 727         }
 728         if (IMDsetup->f[i][ZZ] != IMDsetup->old_forces[i][ZZ])
 729         {
 730             return TRUE;
 731         }
 732     }
 733
 734     /* All old and new forces are identical! */
 735     return FALSE;
 736 }
 737
 738
 739 /*! \brief Fill the old_f_ind and old_forces arrays with the new, old values. */
 740 static void keep_old_values(t_gmx_IMD_setup *IMDsetup)
 741 {
 742     int i;
 743
 744
 745     IMDsetup->old_nforces = IMDsetup->nforces;
 746
 747     for (i = 0; i < IMDsetup->nforces; i++)
 748     {
 749         IMDsetup->old_f_ind[i] = IMDsetup->f_ind[i];
 750         copy_rvec(IMDsetup->f[i], IMDsetup->old_forces[i]);
 751     }
 752 }
 753
 754
 755 /*! \brief Returns TRUE if any component of the two rvecs differs. */
 756 static gmx_inline gmx_bool rvecs_differ(const rvec v1, const rvec v2)
 757 {
 758     int i;
 759
 760
 761     for (i = 0; i < DIM; i++)
 762     {
 763         if (v1[i] != v2[i])
 764         {
 765             return TRUE;
 766         }
 767     }
 768
 769     return FALSE;
 770 }
 771
 772
 773 /*! \brief Write the applied pull forces to logfile.
 774  *
 775  * Call on master only!
 776  */
 777 static void output_imd_forces(t_inputrec *ir, double time)
 778 {
 779     t_gmx_IMD_setup *IMDsetup;
 780     int              i;
 781
 782
 783     IMDsetup = ir->imd->setup;
 784
 785     if (bForcesChanged(IMDsetup))
 786     {
 787         /* Write time and total number of applied IMD forces */
 788         fprintf(IMDsetup->outf, "%14.6e%6d", time, IMDsetup->nforces);
 789
 790         /* Write out the global atom indices of the pulled atoms and the forces itself,
 791          * write out a force only if it has changed since the last output */
 792         for (i = 0; i < IMDsetup->nforces; i++)
 793         {
 794             if (rvecs_differ(IMDsetup->f[i], IMDsetup->old_forces[i]))
 795             {
 796                 fprintf(IMDsetup->outf, "%9d", IMDsetup->ind[IMDsetup->f_ind[i]] + 1);
 797                 fprintf(IMDsetup->outf, "%12.4e%12.4e%12.4e", IMDsetup->f[i][0], IMDsetup->f[i][1], IMDsetup->f[i][2]);
 798             }
 799         }
 800         fprintf(IMDsetup->outf, "\n");
 801
 802         keep_old_values(IMDsetup);
 803     }
 804 }
 805
 806
 807 /*! \brief Synchronize the nodes. */
 808 static void imd_sync_nodes(t_inputrec *ir, t_commrec *cr, double t)
 809 {
 810     int              new_nforces = 0;
 811     t_gmx_IMD_setup *IMDsetup;
 812     int              start, end, i;
 813
 814
 815     IMDsetup = ir->imd->setup;
 816
 817     /* Notify the other nodes whether we are still connected. */
 818     if (PAR(cr))
 819     {
 820         block_bc(cr, IMDsetup->bConnected);
 821     }
 822
 823     /* ...if not connected, the job is done here. */
 824     if (!IMDsetup->bConnected)
 825     {
 826         return;
 827     }
 828
 829     /* Let the other nodes know whether we got a new IMD synchronization frequency. */
 830     if (PAR(cr))
 831     {
 832         block_bc(cr, IMDsetup->nstimd_new);
 833     }
 834
 835     /* Now we all set the (new) nstimd communication time step */
 836     IMDsetup->nstimd = IMDsetup->nstimd_new;
 837
 838     /* We're done if we don't allow pulling at all */
 839     if (!(IMDsetup->bForceActivated))
 840     {
 841         return;
 842     }
 843
 844     /* OK, let's check if we have received forces which we need to communicate
 845      * to the other nodes */
 846     if (MASTER(cr))
 847     {
 848         if (IMDsetup->bNewForces)
 849         {
 850             new_nforces = IMDsetup->vmd_nforces;
 851         }
 852         /* make the "new_forces" negative, when we did not receive new ones */
 853         else
 854         {
 855             new_nforces = IMDsetup->vmd_nforces * -1;
 856         }
 857     }
 858
 859     /* make new_forces known to the clients */
 860     if (PAR(cr))
 861     {
 862         block_bc(cr, new_nforces);
 863     }
 864
 865     /* When new_natoms < 0 then we know that these are still the same forces
 866      * so we don't communicate them, otherwise... */
 867     if (new_nforces >= 0)
 868     {
 869         /* set local VMD and nforces */
 870         IMDsetup->vmd_nforces = new_nforces;
 871         IMDsetup->nforces     = IMDsetup->vmd_nforces;
 872
 873         /* now everybody knows the number of forces in f_ind, so we can prepare
 874          * the target arrays for indices and forces */
 875         imd_prepare_MD_Forces(IMDsetup);
 876
 877         /* we first update the MD forces on the master by converting the VMD forces */
 878         if (MASTER(cr))
 879         {
 880             imd_copyto_MD_Forces(IMDsetup);
 881             /* We also write out forces on every update, so that we know which
 882              * forces are applied for every step */
 883             if (IMDsetup->outf)
 884             {
 885                 output_imd_forces(ir, t);
 886             }
 887         }
 888
 889         /* In parallel mode we communicate the to-be-applied forces to the other nodes */
 890         if (PAR(cr))
 891         {
 892             nblock_bc(cr, IMDsetup->nforces, IMDsetup->f_ind);
 893             nblock_bc(cr, IMDsetup->nforces, IMDsetup->f    );
 894         }
 895
 896         /* done communicating the forces, reset bNewForces */
 897         IMDsetup->bNewForces = FALSE;
 898     }
 899 }
 900
 901
 902 /*! \brief Reads header from the client and decides what to do. */
 903 static void imd_readcommand(t_gmx_IMD_setup *IMDsetup)
 904 {
 905     gmx_bool       IMDpaused = FALSE;
 906     IMDMessageType itype;
 907
 908
 909     while (IMDsetup->clientsocket && (imdsock_tryread(IMDsetup->clientsocket, 0, 0) > 0 || IMDpaused))
 910     {
 911         itype = imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length));
 912         /* let's see what we got: */
 913         switch (itype)
 914         {
 915             /* IMD asks us to terminate the simulation, check if the user allowed this */
 916             case IMD_KILL:
 917                 if (IMDsetup->bTerminatable)
 918                 {
 919                     fprintf(stderr, " %s Terminating connection and running simulation (if supported by integrator).\n", IMDstr);
 920                     IMDsetup->bTerminated = TRUE;
 921                     IMDsetup->bWConnect   = FALSE;
 922                     gmx_set_stop_condition(gmx_stop_cond_next);
 923                 }
 924                 else
 925                 {
 926                     fprintf(stderr, " %s Set -imdterm command line switch to allow mdrun termination from within IMD.\n", IMDstr);
 927                 }
 928
 929                 break;
 930
 931             /* the client doen't want to talk to us anymore */
 932             case IMD_DISCONNECT:
 933                 fprintf(stderr, " %s Disconnecting client.\n", IMDstr);
 934                 imd_disconnect(IMDsetup);
 935                 break;
 936
 937             /* we got new forces, read them and set bNewForces flag */
 938             case IMD_MDCOMM:
 939                 imd_read_vmd_Forces(IMDsetup);
 940                 IMDsetup->bNewForces = TRUE;
 941                 break;
 942
 943             /* the client asks us to (un)pause the simulation. So we toggle the IMDpaused state */
 944             case IMD_PAUSE:
 945                 if (IMDpaused)
 946                 {
 947                     fprintf(stderr, " %s Un-pause command received.\n", IMDstr);
 948                     IMDpaused = FALSE;
 949                 }
 950                 else
 951                 {
 952                     fprintf(stderr, " %s Pause command received.\n", IMDstr);
 953                     IMDpaused = TRUE;
 954                 }
 955
 956                 break;
 957
 958             /* the client sets a new transfer rate, if we get 0, we reset the rate
 959              * to the default. VMD filters 0 however */
 960             case IMD_TRATE:
 961                 IMDsetup->nstimd_new = (IMDsetup->length > 0) ? IMDsetup->length : IMDsetup->nstimd_def;
 962                 fprintf(stderr, " %s Update frequency will be set to %d.\n", IMDstr, IMDsetup->nstimd_new);
 963                 break;
 964
 965             /* Catch all rule for the remaining IMD types which we don't expect */
 966             default:
 967                 fprintf(stderr, " %s Received unexpected %s.\n", IMDstr, ENUM_NAME((int)itype, IMD_NR, eIMDType_names));
 968                 imd_fatal(IMDsetup, "Terminating connection\n");
 969                 break;
 970         } /* end switch */
 971     }     /* end while  */
 972 }
 973
 974
 975 /*! \brief Open IMD output file and write header information.
 976  *
 977  * Call on master only.
 978  */
 979 static FILE *open_imd_out(
 980         const char           *fn,
 981         t_gmx_IMD_setup      *IMDsetup,
 982         int                   nat_total,
 983         output_env_t          oenv,
 984         unsigned long         Flags)
 985 {
 986     FILE       *fp;
 987
 988
 989     /* Open log file of applied IMD forces if requested */
 990     if (fn && oenv)
 991     {
 992         /* If we append to an existing file, all the header information is already there */
 993         if (Flags & MD_APPENDFILES)
 994         {
 995             fp = gmx_fio_fopen(fn, "a+");
 996         }
 997         else
 998         {
 999             fp = gmx_fio_fopen(fn, "w+");
1000             if (IMDsetup->nat == nat_total)
1001             {
1002                 fprintf(fp, "# Note that you can select an IMD index group in the .mdp file if a subset of the atoms suffices.\n");
1003             }
1004
1005             xvgr_header(fp, "IMD Pull Forces", "Time (ps)", "# of Forces / Atom IDs / Forces (kJ/mol)", exvggtNONE, oenv);
1006
1007             fprintf(fp, "# Can display and manipulate %d (of a total of %d) atoms via IMD.\n", IMDsetup->nat, nat_total);
1008             fprintf(fp, "# column 1    : time (ps)\n");
1009             fprintf(fp, "# column 2    : total number of atoms feeling an IMD pulling force at that time\n");
1010             fprintf(fp, "# cols. 3.-6  : global atom number of pulled atom, x-force, y-force, z-force (kJ/mol)\n");
1011             fprintf(fp, "# then follow : atom-ID, f[x], f[y], f[z] for more atoms in case the force on multiple atoms is changed simultaneously.\n");
1012             fprintf(fp, "# Note that the force on any atom is always equal to the last value for that atom-ID found in the data.\n");
1013             fflush(fp);
1014         }
1015
1016         /* To reduce the output file size we remember the old values and output only
1017          * when something changed */
1018         snew(IMDsetup->old_f_ind, IMDsetup->nat);  /* One can never pull on more atoms */
1019         snew(IMDsetup->old_forces, IMDsetup->nat);
1020
1021         return fp;
1022     }
1023
1024     fprintf(stdout, "%s For a log of the IMD pull forces explicitly specify '-if' on the command line.\n"
1025             "%s (Not possible with energy minimization.)\n", IMDstr, IMDstr);
1026
1027     return NULL;
1028 }
1029 #endif
1030
1031
1032 extern void IMD_finalize(gmx_bool bIMD, t_IMD *imd)
1033 {
1034     if (bIMD)
1035     {
1036         if (imd->setup->outf)
1037         {
1038             gmx_fio_fclose(imd->setup->outf);
1039         }
1040     }
1041 }
1042
1043
1044 #ifdef GMX_IMD
1045 /*! \brief Creates the molecule start-end position array of molecules in the IMD group. */
1046 static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, gmx_mtop_t *top_global)
1047 {
1048     int      i, ii;
1049     int      gstart, gend, count;
1050     t_block  gmols, lmols;
1051     int      nat;
1052     atom_id *ind;
1053
1054     gmols = top_global->mols;
1055     nat   = IMDsetup->nat;
1056     ind   = IMDsetup->ind;
1057
1058     lmols.nr = 0;
1059
1060     /* check whether index is sorted */
1061     for (i = 0; i < nat-1; i++)
1062     {
1063         if (ind[i] > ind[i+1])
1064         {
1065             gmx_fatal(FARGS, "%s IMD index is not sorted. This is currently not supported.\n", IMDstr);
1066         }
1067     }
1068
1069     snew(lmols.index, gmols.nr+1);
1070     lmols.index[0] = 0;
1071
1072     for (i = 0; i < gmols.nr; i++)
1073     {
1074         gstart = gmols.index[i];
1075         gend   = gmols.index[i+1];
1076         count  = 0;
1077         for (ii = 0; ii < nat; ii++)
1078         {
1079             if ((ind[ii] >= gstart) && (ind[ii] < gend))
1080             {
1081                 count += 1;
1082             }
1083         }
1084         if (count > 0)
1085         {
1086             lmols.index[lmols.nr+1] = lmols.index[lmols.nr]+count;
1087             lmols.nr               += 1;
1088         }
1089     }
1090     srenew(lmols.index, lmols.nr+1);
1091     lmols.nalloc_index = lmols.nr+1;
1092     IMDsetup->mols     = lmols;
1093 }
1094
1095
1096 /*! \brief Copied and modified from groupcoord.c shift_positions_group(). */
1097 static void shift_positions(
1098         matrix box,
1099         rvec   x[],      /* The positions [0..nr] */
1100         ivec   is,       /* The shift [0..nr] */
1101         int    nr)       /* The number of positions */
1102 {
1103     int      i, tx, ty, tz;
1104
1105     /* Loop over the group's atoms */
1106     if (TRICLINIC(box))
1107     {
1108         for (i = 0; i < nr; i++)
1109         {
1110             tx = is[XX];
1111             ty = is[YY];
1112             tz = is[ZZ];
1113
1114             x[i][XX] = x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1115             x[i][YY] = x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1116             x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
1117         }
1118     }
1119     else
1120     {
1121         for (i = 0; i < nr; i++)
1122         {
1123             tx = is[XX];
1124             ty = is[YY];
1125             tz = is[ZZ];
1126
1127             x[i][XX] = x[i][XX]-tx*box[XX][XX];
1128             x[i][YY] = x[i][YY]-ty*box[YY][YY];
1129             x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
1130         }
1131     }
1132 }
1133
1134
1135 /*! \brief Removes shifts of molecules diffused outside of the box. */
1136 static void imd_remove_molshifts(t_gmx_IMD_setup *IMDsetup, matrix box)
1137 {
1138     int     i, ii, molsize;
1139     ivec    largest, smallest, shift;
1140     t_block mols;
1141
1142
1143     mols = IMDsetup->mols;
1144
1145     /* for each molecule also present in IMD group */
1146     for (i = 0; i < mols.nr; i++)
1147     {
1148         /* first we determine the minimum and maximum shifts for each molecule */
1149
1150         clear_ivec(largest);
1151         clear_ivec(smallest);
1152         clear_ivec(shift);
1153
1154         copy_ivec(IMDsetup->xa_shifts[mols.index[i]], largest);
1155         copy_ivec(IMDsetup->xa_shifts[mols.index[i]], smallest);
1156
1157         for (ii = mols.index[i]+1; ii < mols.index[i+1]; ii++)
1158         {
1159             if (IMDsetup->xa_shifts[ii][XX] > largest[XX])
1160             {
1161                 largest[XX]  = IMDsetup->xa_shifts[ii][XX];
1162             }
1163             if (IMDsetup->xa_shifts[ii][XX] < smallest[XX])
1164             {
1165                 smallest[XX] = IMDsetup->xa_shifts[ii][XX];
1166             }
1167
1168             if (IMDsetup->xa_shifts[ii][YY] > largest[YY])
1169             {
1170                 largest[YY]  = IMDsetup->xa_shifts[ii][YY];
1171             }
1172             if (IMDsetup->xa_shifts[ii][YY] < smallest[YY])
1173             {
1174                 smallest[YY] = IMDsetup->xa_shifts[ii][YY];
1175             }
1176
1177             if (IMDsetup->xa_shifts[ii][ZZ] > largest[ZZ])
1178             {
1179                 largest[ZZ]  = IMDsetup->xa_shifts[ii][ZZ];
1180             }
1181             if (IMDsetup->xa_shifts[ii][ZZ] < smallest[ZZ])
1182             {
1183                 smallest[ZZ] = IMDsetup->xa_shifts[ii][ZZ];
1184             }
1185
1186         }
1187
1188         /* check if we what we can subtract/add to the positions
1189          * to put them back into the central box */
1190         if (smallest[XX] > 0)
1191         {
1192             shift[XX] = smallest[XX];
1193         }
1194         if (smallest[YY] > 0)
1195         {
1196             shift[YY] = smallest[YY];
1197         }
1198         if (smallest[ZZ] > 0)
1199         {
1200             shift[ZZ] = smallest[ZZ];
1201         }
1202
1203         if (largest[XX] < 0)
1204         {
1205             shift[XX] = largest[XX];
1206         }
1207         if (largest[YY] < 0)
1208         {
1209             shift[YY] = largest[YY];
1210         }
1211         if (largest[ZZ] < 0)
1212         {
1213             shift[ZZ] = largest[ZZ];
1214         }
1215
1216         /* is there a shift at all? */
1217         if ((shift[XX]) || (shift[YY]) || (shift[ZZ]))
1218         {
1219             molsize = mols.index[i+1]-mols.index[i];
1220             /* shift the positions */
1221             shift_positions(box, &(IMDsetup->xa[mols.index[i]]), shift, molsize);
1222         }
1223
1224     }
1225 }
1226
1227
1228 /*! \brief Initialize arrays used to assemble the positions from the other nodes. */
1229 static void init_imd_prepare_for_x_assembly(t_commrec *cr, rvec x[], t_gmx_IMD_setup *IMDsetup)
1230 {
1231     int i, ii;
1232
1233
1234     snew(IMDsetup->xa,         IMDsetup->nat);
1235     snew(IMDsetup->xa_ind,     IMDsetup->nat);
1236     snew(IMDsetup->xa_shifts,  IMDsetup->nat);
1237     snew(IMDsetup->xa_eshifts, IMDsetup->nat);
1238     snew(IMDsetup->xa_old,     IMDsetup->nat);
1239
1240     /* Save the original (whole) set of positions such that later the
1241      * molecule can always be made whole again */
1242     if (MASTER(cr))
1243     {
1244         for (i = 0; i < IMDsetup->nat; i++)
1245         {
1246             ii = IMDsetup->ind[i];
1247             copy_rvec(x[ii], IMDsetup->xa_old[i]);
1248         }
1249     }
1250
1251     if (!PAR(cr))
1252     {
1253         IMDsetup->nat_loc = IMDsetup->nat;
1254         IMDsetup->ind_loc = IMDsetup->ind;
1255
1256         /* xa_ind[i] needs to be set to i for serial runs */
1257         for (i = 0; i < IMDsetup->nat; i++)
1258         {
1259             IMDsetup->xa_ind[i] = i;
1260         }
1261     }
1262
1263     /* Communicate initial coordinates xa_old to all processes */
1264 #ifdef GMX_MPI
1265     if (PAR(cr))
1266     {
1267         gmx_bcast(IMDsetup->nat * sizeof(IMDsetup->xa_old[0]), IMDsetup->xa_old, cr);
1268     }
1269 #endif
1270 }
1271 #endif
1272
1273
1274 /*! \brief Check for non-working integrator / parallel options. */
1275 static void imd_check_integrator_parallel(t_inputrec *ir, t_commrec *cr)
1276 {
1277     if (PAR(cr))
1278     {
1279         if (((ir->eI) == eiSteep) || ((ir->eI) == eiCG) || ((ir->eI) == eiLBFGS) || ((ir->eI) == eiNM))
1280         {
1281             gmx_fatal(FARGS, "%s Energy minimization via steep, CG, lbfgs and nm in parallel is currently not supported by IMD.\n", IMDstr);
1282             return;
1283         }
1284     }
1285 }
1286
1287
1288 extern void init_IMD(
1289         t_inputrec    *ir,
1290         t_commrec     *cr,
1291         gmx_mtop_t    *top_global,
1292         FILE          *fplog,
1293         int            defnstimd,
1294         rvec           x[],
1295         int            nfile,
1296         const t_filenm fnm[],
1297         output_env_t   oenv,
1298         int            imdport,
1299         unsigned long  Flags)
1300 {
1301     int              i;
1302     int              nat_total;
1303     t_gmx_IMD_setup *IMDsetup;
1304     gmx_int32_t      bufxsize;
1305     gmx_bool         bIMD = FALSE;
1306
1307
1308     /* We will allow IMD sessions only if explicitly enabled in the .tpr file */
1309     if (FALSE == ir->bIMD)
1310     {
1311         return;
1312     }
1313
1314     /* It seems we have a .tpr file that defines an IMD group and thus allows IMD sessions.
1315      * Check whether we can actually provide the IMD functionality for this setting: */
1316     if (MASTER(cr))
1317     {
1318         /* Check whether IMD was enabled by one of the command line switches: */
1319         if ((Flags & MD_IMDWAIT) || (Flags & MD_IMDTERM) || (Flags & MD_IMDPULL))
1320         {
1321             /* Multiple simulations or replica exchange */
1322             if (MULTISIM(cr))
1323             {
1324                 fprintf(stderr, "%s Cannot use IMD for multiple simulations or replica exchange.\n", IMDstr);
1325             }
1326             /* OK, IMD seems to be allowed and turned on... */
1327             else
1328             {
1329                 fprintf(stderr, "%s Enabled. This simulation will accept incoming IMD connections.\n", IMDstr);
1330                 bIMD = TRUE;
1331             }
1332         }
1333         else
1334         {
1335             fprintf(stderr, "%s None of the -imd switches was used.\n"
1336                     "%s This run will not accept incoming IMD connections\n", IMDstr, IMDstr);
1337         }
1338     } /* end master only */
1339
1340     /* Disable IMD if not all the needed functionality is there! */
1341 #if defined(GMX_NATIVE_WINDOWS) && !defined(GMX_HAVE_WINSOCK)
1342     bIMD = FALSE;
1343     fprintf(stderr, "Disabling IMD because the winsock library was not found at compile time.\n");
1344 #endif
1345
1346     /* Let the other nodes know whether we want IMD */
1347     if (PAR(cr))
1348     {
1349         block_bc(cr, bIMD);
1350     }
1351     /* ... and update our local inputrec accordingly. */
1352     ir->bIMD = bIMD;
1353
1354     /*... if not we are done.*/
1355     if (!ir->bIMD)
1356     {
1357         return;
1358     }
1359
1360
1361     /* check if we're using a sane integrator / parallel combination */
1362     imd_check_integrator_parallel(ir, cr);
1363
1364
1365     /*
1366      *****************************************************
1367      * From here on we assume that IMD is turned on      *
1368      *****************************************************
1369      */
1370
1371 #ifdef GMX_IMD
1372     nat_total = top_global->natoms;
1373
1374     /* Initialize IMD setup structure. If we read in a pre-IMD .tpr file, imd->nat
1375      * will be zero. For those cases we transfer _all_ atomic positions */
1376     ir->imd->setup = imd_create(ir->imd->nat > 0 ? ir->imd->nat : nat_total,
1377                                 defnstimd, imdport);
1378     IMDsetup       = ir->imd->setup;
1379
1380     /* We might need to open an output file for IMD forces data */
1381     if (MASTER(cr))
1382     {
1383         IMDsetup->outf = open_imd_out(opt2fn("-if", nfile, fnm), ir->imd->setup, nat_total, oenv, Flags);
1384     }
1385
1386     /* Make sure that we operate with a valid atom index array for the IMD atoms */
1387     if (ir->imd->nat > 0)
1388     {
1389         /* Point to the user-supplied array of atom numbers */
1390         IMDsetup->ind = ir->imd->ind;
1391     }
1392     else
1393     {
1394         /* Make a dummy (ind[i] = i) array of all atoms */
1395         snew(IMDsetup->ind, nat_total);
1396         for (i = 0; i < nat_total; i++)
1397         {
1398             IMDsetup->ind[i] = i;
1399         }
1400     }
1401
1402     /* read environment on master and prepare socket for incoming connections */
1403     if (MASTER(cr))
1404     {
1405         /* we allocate memory for our IMD energy structure */
1406         gmx_int32_t recsize = HEADERSIZE + sizeof(IMDEnergyBlock);
1407         snew(IMDsetup->energysendbuf, recsize);
1408
1409         /* Shall we wait for a connection? */
1410         if (Flags & MD_IMDWAIT)
1411         {
1412             IMDsetup->bWConnect = TRUE;
1413             fprintf(stderr, "%s Pausing simulation while no IMD connection present (-imdwait).\n", IMDstr);
1414         }
1415
1416         /* Will the IMD clients be able to terminate the simulation? */
1417         if (Flags & MD_IMDTERM)
1418         {
1419             IMDsetup->bTerminatable = TRUE;
1420             fprintf(stderr, "%s Allow termination of the simulation from IMD client (-imdterm).\n", IMDstr);
1421         }
1422
1423         /* Is pulling from IMD client allowed? */
1424         if (Flags & MD_IMDPULL)
1425         {
1426             IMDsetup->bForceActivated = TRUE;
1427             fprintf(stderr, "%s Pulling from IMD remote is enabled (-imdpull).\n", IMDstr);
1428         }
1429
1430         /* Initialize send buffers with constant size */
1431         snew(IMDsetup->sendxbuf, IMDsetup->nat);
1432         snew(IMDsetup->energies, 1);
1433         bufxsize = HEADERSIZE + 3 * sizeof(float) * IMDsetup->nat;
1434         snew(IMDsetup->coordsendbuf, bufxsize);
1435     }
1436
1437     /* do we allow interactive pulling? If so let the other nodes know. */
1438     if (PAR(cr))
1439     {
1440         block_bc(cr, IMDsetup->bForceActivated);
1441     }
1442
1443     /* setup the listening socket on master process */
1444     if (MASTER(cr))
1445     {
1446         fprintf(fplog, "%s Setting port for connection requests to %d.\n", IMDstr, IMDsetup->port);
1447         fprintf(stderr, "%s Turning on IMD - port for incoming requests is %d.\n", IMDstr, IMDsetup->port);
1448         imd_prepare_master_socket(IMDsetup);
1449         /* Wait until we have a connection if specified before */
1450         if (IMDsetup->bWConnect)
1451         {
1452             imd_blockconnect(IMDsetup);
1453         }
1454         else
1455         {
1456             fprintf(stderr, "%s -imdwait not set, starting simulation.\n", IMDstr);
1457         }
1458     }
1459     /* Let the other nodes know whether we are connected */
1460     imd_sync_nodes(ir, cr, 0);
1461
1462     /* Initialize arrays used to assemble the positions from the other nodes */
1463     init_imd_prepare_for_x_assembly(cr, x, IMDsetup);
1464
1465     /* Initialize molecule blocks to make them whole later...*/
1466     if (MASTER(cr))
1467     {
1468         init_imd_prepare_mols_in_imdgroup(IMDsetup, top_global);
1469     }
1470 #else
1471     gmx_incons("init_IMD: this GROMACS version was not compiled with IMD support!");
1472 #endif
1473 }
1474
1475
1476 extern gmx_bool do_IMD(
1477         gmx_bool        bIMD,
1478         gmx_int64_t     step,
1479         t_commrec      *cr,
1480         gmx_bool        bNS,
1481         matrix          box,
1482         rvec            x[],
1483         t_inputrec     *ir,
1484         double          t,
1485         gmx_wallcycle_t wcycle)
1486 {
1487     gmx_bool         imdstep = FALSE;
1488     t_gmx_IMD_setup *IMDsetup;
1489
1490
1491     /* IMD at all? */
1492     if (!bIMD)
1493     {
1494         return FALSE;
1495     }
1496
1497 #ifdef GMX_IMD
1498     wallcycle_start(wcycle, ewcIMD);
1499
1500     IMDsetup = ir->imd->setup;
1501
1502     /* read command from client and check if new incoming connection */
1503     if (MASTER(cr))
1504     {
1505         /* If not already connected, check for new connections */
1506         if (!IMDsetup->clientsocket)
1507         {
1508             if (IMDsetup->bWConnect)
1509             {
1510                 imd_blockconnect(IMDsetup);
1511             }
1512             else
1513             {
1514                 imd_tryconnect(IMDsetup);
1515             }
1516         }
1517
1518         /* Let's see if we have new IMD messages for us */
1519         if (IMDsetup->clientsocket)
1520         {
1521             imd_readcommand(IMDsetup);
1522         }
1523     }
1524
1525     /* is this an IMD communication step? */
1526     imdstep = do_per_step(step, IMDsetup->nstimd);
1527
1528     /* OK so this is an IMD step ... */
1529     if (imdstep)
1530     {
1531         /* First we sync all nodes to let everybody know whether we are connected to VMD */
1532         imd_sync_nodes(ir, cr, t);
1533     }
1534
1535     /* If a client is connected, we collect the positions
1536      * and put molecules back into the box before transfer */
1537     if ((imdstep && IMDsetup->bConnected)
1538         || bNS)            /* independent of imdstep, we communicate positions at each NS step */
1539     {
1540         /* Transfer the IMD positions to the master node. Every node contributes
1541          * its local positions x and stores them in the assembled xa array. */
1542         communicate_group_positions(cr, IMDsetup->xa, IMDsetup->xa_shifts, IMDsetup->xa_eshifts,
1543                                     TRUE, x, IMDsetup->nat, IMDsetup->nat_loc,
1544                                     IMDsetup->ind_loc, IMDsetup->xa_ind, IMDsetup->xa_old, box);
1545
1546         /* If connected and master -> remove shifts */
1547         if ((imdstep && IMDsetup->bConnected) && MASTER(cr))
1548         {
1549             imd_remove_molshifts(IMDsetup, box);
1550         }
1551     }
1552
1553     wallcycle_stop(wcycle, ewcIMD);
1554 #else
1555     gmx_incons("do_IMD called without IMD support!");
1556 #endif
1557
1558     return imdstep;
1559 }
1560
1561
1562 extern void IMD_fill_energy_record(gmx_bool bIMD, t_IMD *imd, gmx_enerdata_t *enerd,
1563                                    gmx_int64_t step, gmx_bool bHaveNewEnergies)
1564 {
1565     IMDEnergyBlock *ene;
1566     t_gmx_IMD       IMDsetup;
1567
1568
1569     if (bIMD)
1570     {
1571 #ifdef GMX_IMD
1572         IMDsetup = imd->setup;
1573
1574         if (IMDsetup->clientsocket)
1575         {
1576             ene = IMDsetup->energies;
1577
1578             ene->tstep = step;
1579
1580             /* In MPI-parallel simulations the energies are not accessible a at every time step.
1581              * We update them if we have new values, otherwise, the energy values from the
1582              * last global communication step are still on display in the viewer. */
1583             if (bHaveNewEnergies)
1584             {
1585                 ene->T_abs   = (float)  enerd->term[F_TEMP   ];
1586                 ene->E_pot   = (float)  enerd->term[F_EPOT   ];
1587                 ene->E_tot   = (float)  enerd->term[F_ETOT   ];
1588                 ene->E_bond  = (float)  enerd->term[F_BONDS  ];
1589                 ene->E_angle = (float)  enerd->term[F_ANGLES ];
1590                 ene->E_dihe  = (float)  enerd->term[F_PDIHS  ];
1591                 ene->E_impr  = (float)  enerd->term[F_IDIHS  ];
1592                 ene->E_vdw   = (float)  enerd->term[F_LJ     ];
1593                 ene->E_coul  = (float) (enerd->term[F_COUL_SR] + enerd->term[F_COUL_LR]);
1594             }
1595         }
1596 #else
1597         gmx_incons("IMD_fill_energy_record called without IMD support.");
1598 #endif
1599     }
1600 }
1601
1602
1603 extern void IMD_send_positions(t_IMD *imd)
1604 {
1605 #ifdef GMX_IMD
1606     t_gmx_IMD IMDsetup;
1607
1608
1609     IMDsetup = imd->setup;
1610
1611     if (IMDsetup->clientsocket)
1612     {
1613
1614         if (imd_send_energies(IMDsetup->clientsocket, IMDsetup->energies, IMDsetup->energysendbuf))
1615         {
1616             imd_fatal(IMDsetup, "Error sending updated energies. Disconnecting client.\n");
1617         }
1618
1619         if (imd_send_rvecs(IMDsetup->clientsocket, IMDsetup->nat, IMDsetup->xa, IMDsetup->coordsendbuf))
1620         {
1621             imd_fatal(IMDsetup, "Error sending updated positions. Disconnecting client.\n");
1622         }
1623     }
1624 #else
1625     gmx_incons("IMD_send_positions called without IMD support.");
1626 #endif
1627 }
1628
1629
1630 extern void IMD_prep_energies_send_positions(gmx_bool bIMD, gmx_bool bIMDstep,
1631                                              t_IMD *imd, gmx_enerdata_t *enerd,
1632                                              gmx_int64_t step, gmx_bool bHaveNewEnergies,
1633                                              gmx_wallcycle_t wcycle)
1634 {
1635     if (bIMD)
1636     {
1637 #ifdef GMX_IMD
1638         wallcycle_start(wcycle, ewcIMD);
1639
1640         /* Update time step for IMD and prepare IMD energy record if we have new energies. */
1641         IMD_fill_energy_record(TRUE, imd, enerd, step, bHaveNewEnergies);
1642
1643         if (bIMDstep)
1644         {
1645             /* Send positions and energies to VMD client via IMD */
1646             IMD_send_positions(imd);
1647         }
1648
1649         wallcycle_stop(wcycle, ewcIMD);
1650 #else
1651         gmx_incons("IMD_prep_energies_send_positions called without IMD support.");
1652 #endif
1653     }
1654 }
1655
1656
1657 extern int IMD_get_step(t_gmx_IMD IMDsetup)
1658 {
1659     return IMDsetup->nstimd;
1660 }
1661
1662
1663 extern void IMD_apply_forces(gmx_bool bIMD, t_IMD *imd, t_commrec *cr, rvec *f,
1664                              gmx_wallcycle_t wcycle)
1665 {
1666     int              i, j;
1667     int              locndx;
1668     t_gmx_IMD_setup *IMDsetup;
1669
1670
1671     if (bIMD)
1672     {
1673 #ifdef GMX_IMD
1674         wallcycle_start(wcycle, ewcIMD);
1675
1676         IMDsetup = imd->setup;
1677
1678         /* Are forces allowed at all? If not we're done */
1679         if (!IMDsetup->bForceActivated)
1680         {
1681             return;
1682         }
1683
1684         for (i = 0; i < IMDsetup->nforces; i++)
1685         {
1686             /* j are the indices in the "System group".*/
1687             j = IMDsetup->ind[IMDsetup->f_ind[i]];
1688
1689             /* check if this is a local atom and find out locndx */
1690             if (PAR(cr) && ga2la_get_home(cr->dd->ga2la, j, &locndx))
1691             {
1692                 j = locndx;
1693             }
1694
1695             rvec_inc(f[j], IMDsetup->f[i]);
1696         }
1697
1698         wallcycle_start(wcycle, ewcIMD);
1699 #else
1700         gmx_incons("IMD_apply_forces called without IMD support.");
1701 #endif
1702     }
1703 }