src/kernel/runner.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  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <signal.h>
  41 #include <stdlib.h>
  42
  43 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
  44 /* _isnan() */
  45 #include <float.h>
  46 #endif
  47
  48 #include "typedefs.h"
  49 #include "smalloc.h"
  50 #include "sysstuff.h"
  51 #include "statutil.h"
  52 #include "mdrun.h"
  53 #include "network.h"
  54 #include "pull.h"
  55 #include "names.h"
  56 #include "disre.h"
  57 #include "orires.h"
  58 #include "dihre.h"
  59 #include "pppm.h"
  60 #include "pme.h"
  61 #include "mdatoms.h"
  62 #include "repl_ex.h"
  63 #include "qmmm.h"
  64 #include "mpelogging.h"
  65 #include "domdec.h"
  66 #include "partdec.h"
  67 #include "coulomb.h"
  68 #include "constr.h"
  69 #include "mvdata.h"
  70 #include "checkpoint.h"
  71 #include "mtop_util.h"
  72 #include "sighandler.h"
  73
  74 #include "md_openmm.h"
  75
  76 #ifdef GMX_LIB_MPI
  77 #include <mpi.h>
  78 #endif
  79 #ifdef GMX_THREADS
  80 #include "tmpi.h"
  81 #endif
  82
  83 #ifdef GMX_FAHCORE
  84 #include "corewrap.h"
  85 #endif
  86
  87 #ifdef GMX_OPENMM
  88 #include "md_openmm.h"
  89 #endif
  90
  91 /* The following two variables and the signal_handler function
  92  * are used from md.c and pme.c as well
  93  */
  94 extern bool bGotTermSignal, bGotUsr1Signal;
  95
  96 static RETSIGTYPE signal_handler(int n)
  97 {
  98     switch (n) {
  99         case SIGTERM:
 100             bGotTermSignal = TRUE;
 101             break;
 102 #ifdef HAVE_SIGUSR1
 103         case SIGUSR1:
 104             bGotUsr1Signal = TRUE;
 105             break;
 106 #endif
 107     }
 108 }
 109
 110 typedef struct {
 111     gmx_integrator_t *func;
 112 } gmx_intp_t;
 113
 114 /* The array should match the eI array in include/types/enums.h */
 115 #ifdef GMX_OPENMM  /* FIXME do_md_openmm needs fixing */
 116 const gmx_intp_t integrator[eiNR] = { {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm},{do_md_openmm}};
 117 #else
 118 const gmx_intp_t integrator[eiNR] = { {do_md}, {do_steep}, {do_cg}, {do_md}, {do_md}, {do_nm}, {do_lbfgs}, {do_tpi}, {do_tpi}, {do_md}, {do_md},{do_md}};
 119 #endif
 120
 121 gmx_large_int_t     deform_init_init_step_tpx;
 122 matrix              deform_init_box_tpx;
 123 #ifdef GMX_THREADS
 124 tMPI_Thread_mutex_t deform_init_box_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
 125 #endif
 126
 127
 128 #ifdef GMX_THREADS
 129 struct mdrunner_arglist
 130 {
 131     FILE *fplog;
 132     t_commrec *cr;
 133     int nfile;
 134     const t_filenm *fnm;
 135     output_env_t oenv;
 136     bool bVerbose;
 137     bool bCompact;
 138     int nstglobalcomm;
 139     ivec ddxyz;
 140     int dd_node_order;
 141     real rdd;
 142     real rconstr;
 143     const char *dddlb_opt;
 144     real dlb_scale;
 145     const char *ddcsx;
 146     const char *ddcsy;
 147     const char *ddcsz;
 148     int nstepout;
 149     int resetstep;
 150     int nmultisim;
 151     int repl_ex_nst;
 152     int repl_ex_seed;
 153     real pforce;
 154     real cpt_period;
 155     real max_hours;
 156     const char *deviceOptions;
 157     unsigned long Flags;
 158     int ret; /* return value */
 159 };
 160
 161
 162 static void mdrunner_start_fn(void *arg)
 163 {
 164     struct mdrunner_arglist *mda=(struct mdrunner_arglist*)arg;
 165     struct mdrunner_arglist mc=*mda; /* copy the arg list to make sure
 166                                         that it's thread-local. This doesn't
 167                                         copy pointed-to items, of course,
 168                                         but those are all const. */
 169     t_commrec *cr;  /* we need a local version of this */
 170     FILE *fplog=NULL;
 171     t_filenm *fnm=dup_tfn(mc.nfile, mc.fnm);
 172
 173     cr=init_par_threads(mc.cr);
 174     if (MASTER(cr))
 175     {
 176         fplog=mc.fplog;
 177     }
 178
 179
 180     mda->ret=mdrunner(fplog, cr, mc.nfile, mc.fnm, mc.oenv, mc.bVerbose,
 181                       mc.bCompact, mc.nstglobalcomm,
 182                       mc.ddxyz, mc.dd_node_order, mc.rdd,
 183                       mc.rconstr, mc.dddlb_opt, mc.dlb_scale,
 184                       mc.ddcsx, mc.ddcsy, mc.ddcsz, mc.nstepout, mc.resetstep, mc.nmultisim,
 185                       mc.repl_ex_nst, mc.repl_ex_seed, mc.pforce,
 186                       mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.Flags);
 187 }
 188
 189 #endif
 190
 191 int mdrunner_threads(int nthreads,
 192                      FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 193                      const output_env_t oenv, bool bVerbose,bool bCompact,
 194                      int nstglobalcomm,
 195                      ivec ddxyz,int dd_node_order,real rdd,real rconstr,
 196                      const char *dddlb_opt,real dlb_scale,
 197                      const char *ddcsx,const char *ddcsy,const char *ddcsz,
 198                      int nstepout,int resetstep,int nmultisim,int repl_ex_nst,
 199                      int repl_ex_seed, real pforce,real cpt_period,
 200                      real max_hours, const char *deviceOptions, unsigned long Flags)
 201 {
 202     int ret;
 203     /* first check whether we even need to start tMPI */
 204     if (nthreads < 2)
 205     {
 206         ret=mdrunner(fplog, cr, nfile, fnm, oenv, bVerbose, bCompact,
 207                      nstglobalcomm,
 208                      ddxyz, dd_node_order, rdd, rconstr, dddlb_opt, dlb_scale,
 209                      ddcsx, ddcsy, ddcsz, nstepout, resetstep, nmultisim, repl_ex_nst,
 210                      repl_ex_seed, pforce, cpt_period, max_hours, deviceOptions, Flags);
 211     }
 212     else
 213     {
 214 #ifdef GMX_THREADS
 215         struct mdrunner_arglist mda;
 216         /* fill the data structure to pass as void pointer to thread start fn */
 217         mda.fplog=fplog;
 218         mda.cr=cr;
 219         mda.nfile=nfile;
 220         mda.fnm=fnm;
 221         mda.oenv=oenv;
 222         mda.bVerbose=bVerbose;
 223         mda.bCompact=bCompact;
 224         mda.nstglobalcomm=nstglobalcomm;
 225         mda.ddxyz[XX]=ddxyz[XX];
 226         mda.ddxyz[YY]=ddxyz[YY];
 227         mda.ddxyz[ZZ]=ddxyz[ZZ];
 228         mda.dd_node_order=dd_node_order;
 229         mda.rdd=rdd;
 230         mda.rconstr=rconstr;
 231         mda.dddlb_opt=dddlb_opt;
 232         mda.dlb_scale=dlb_scale;
 233         mda.ddcsx=ddcsx;
 234         mda.ddcsy=ddcsy;
 235         mda.ddcsz=ddcsz;
 236         mda.nstepout=nstepout;
 237         mda.resetstep=resetstep;
 238         mda.nmultisim=nmultisim;
 239         mda.repl_ex_nst=repl_ex_nst;
 240         mda.repl_ex_seed=repl_ex_seed;
 241         mda.pforce=pforce;
 242         mda.cpt_period=cpt_period;
 243         mda.max_hours=max_hours;
 244         mda.deviceOptions=deviceOptions;
 245         mda.Flags=Flags;
 246
 247         fprintf(stderr, "Starting %d threads\n",nthreads);
 248         fflush(stderr);
 249         tMPI_Init_fn(nthreads, mdrunner_start_fn, (void*)(&mda) );
 250         ret=mda.ret;
 251 #else
 252         ret=-1;
 253         gmx_comm("Multiple threads requested but not compiled with threads");
 254 #endif
 255     }
 256     return ret;
 257 }
 258
 259
 260 int mdrunner(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
 261              const output_env_t oenv, bool bVerbose,bool bCompact,
 262              int nstglobalcomm,
 263              ivec ddxyz,int dd_node_order,real rdd,real rconstr,
 264              const char *dddlb_opt,real dlb_scale,
 265              const char *ddcsx,const char *ddcsy,const char *ddcsz,
 266              int nstepout,int resetstep,int nmultisim,int repl_ex_nst,int repl_ex_seed,
 267              real pforce,real cpt_period,real max_hours,
 268              const char *deviceOptions,
 269              unsigned long Flags)
 270 {
 271     double     nodetime=0,realtime;
 272     t_inputrec *inputrec;
 273     t_state    *state=NULL;
 274     matrix     box;
 275     gmx_ddbox_t ddbox;
 276     int        npme_major;
 277     real       tmpr1,tmpr2;
 278     t_nrnb     *nrnb;
 279     gmx_mtop_t *mtop=NULL;
 280     t_mdatoms  *mdatoms=NULL;
 281     t_forcerec *fr=NULL;
 282     t_fcdata   *fcd=NULL;
 283     real       ewaldcoeff=0;
 284     gmx_pme_t  *pmedata=NULL;
 285     gmx_vsite_t *vsite=NULL;
 286     gmx_constr_t constr;
 287     int        i,m,nChargePerturbed=-1,status,nalloc;
 288     char       *gro;
 289     gmx_wallcycle_t wcycle;
 290     bool       bReadRNG,bReadEkin;
 291     int        list;
 292     gmx_runtime_t runtime;
 293     int        rc;
 294     gmx_large_int_t reset_counters;
 295     gmx_edsam_t ed;
 296
 297     /* A parallel command line option consistency check */
 298     if (!PAR(cr) &&
 299         (ddxyz[XX] > 1 || ddxyz[YY] > 1 || ddxyz[ZZ] > 1 || cr->npmenodes > 0))
 300     {
 301         gmx_fatal(FARGS,
 302                   "The -dd or -npme option request a parallel simulation, "
 303 #ifndef GMX_MPI
 304                   "but mdrun was compiled without threads or MPI enabled"
 305 #else
 306 #ifdef GMX_THREADS
 307                   "but the number of threads (option -nt) is 1"
 308 #else
 309                   "but mdrun was not started through mpirun/mpiexec or only one process was requested through mpirun/mpiexec"
 310 #endif
 311 #endif
 312             );
 313     }
 314
 315     /* Essential dynamics */
 316     if (opt2bSet("-ei",nfile,fnm))
 317     {
 318         /* Open input and output files, allocate space for ED data structure */
 319         ed = ed_open(nfile,fnm,cr);
 320     }
 321     else
 322         ed=NULL;
 323
 324     snew(inputrec,1);
 325     snew(mtop,1);
 326
 327     if (bVerbose && SIMMASTER(cr))
 328     {
 329         fprintf(stderr,"Getting Loaded...\n");
 330     }
 331
 332     if (Flags & MD_APPENDFILES)
 333     {
 334         fplog = NULL;
 335     }
 336
 337     if (PAR(cr))
 338     {
 339         /* The master thread on the master node reads from disk,
 340          * then distributes everything to the other processors.
 341          */
 342
 343         list = (SIMMASTER(cr) && !(Flags & MD_APPENDFILES)) ?  (LIST_SCALARS | LIST_INPUTREC) : 0;
 344
 345         snew(state,1);
 346         init_parallel(fplog, opt2fn_master("-s",nfile,fnm,cr),cr,
 347                       inputrec,mtop,state,list);
 348
 349     }
 350     else
 351     {
 352         /* Read a file for a single processor */
 353         snew(state,1);
 354         init_single(fplog,inputrec,ftp2fn(efTPX,nfile,fnm),mtop,state);
 355     }
 356     if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
 357     {
 358         cr->npmenodes = 0;
 359     }
 360
 361 #ifdef GMX_FAHCORE
 362     fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
 363 #endif
 364
 365     /* NMR restraints must be initialized before load_checkpoint,
 366      * since with time averaging the history is added to t_state.
 367      * For proper consistency check we therefore need to extend
 368      * t_state here.
 369      * So the PME-only nodes (if present) will also initialize
 370      * the distance restraints.
 371      */
 372     snew(fcd,1);
 373
 374     /* This needs to be called before read_checkpoint to extend the state */
 375     init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
 376
 377     if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
 378     {
 379         if (PAR(cr) && !(Flags & MD_PARTDEC))
 380         {
 381             gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
 382         }
 383         /* Orientation restraints */
 384         if (MASTER(cr))
 385         {
 386             init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
 387                         state);
 388         }
 389     }
 390
 391     if (DEFORM(*inputrec))
 392     {
 393         /* Store the deform reference box before reading the checkpoint */
 394         if (SIMMASTER(cr))
 395         {
 396             copy_mat(state->box,box);
 397         }
 398         if (PAR(cr))
 399         {
 400             gmx_bcast(sizeof(box),box,cr);
 401         }
 402         /* Because we do not have the update struct available yet
 403          * in which the reference values should be stored,
 404          * we store them temporarily in static variables.
 405          * This should be thread safe, since they are only written once
 406          * and with identical values.
 407          */
 408 #ifdef GMX_THREADS
 409         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
 410 #endif
 411         deform_init_init_step_tpx = inputrec->init_step;
 412         copy_mat(box,deform_init_box_tpx);
 413 #ifdef GMX_THREADS
 414         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
 415 #endif
 416     }
 417
 418     if (opt2bSet("-cpi",nfile,fnm))
 419     {
 420         /* Check if checkpoint file exists before doing continuation.
 421          * This way we can use identical input options for the first and subsequent runs...
 422          */
 423         if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
 424         {
 425             load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
 426                             cr,Flags & MD_PARTDEC,ddxyz,
 427                             inputrec,state,&bReadRNG,&bReadEkin,
 428                             (Flags & MD_APPENDFILES));
 429
 430             if (bReadRNG)
 431             {
 432                 Flags |= MD_READ_RNG;
 433             }
 434             if (bReadEkin)
 435             {
 436                 Flags |= MD_READ_EKIN;
 437             }
 438         }
 439     }
 440
 441     if ((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
 442     {
 443         gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
 444                              Flags,&fplog);
 445     }
 446
 447     if (SIMMASTER(cr))
 448     {
 449         copy_mat(state->box,box);
 450     }
 451
 452     if (PAR(cr))
 453     {
 454         gmx_bcast(sizeof(box),box,cr);
 455     }
 456
 457     if (bVerbose && SIMMASTER(cr))
 458     {
 459         fprintf(stderr,"Loaded with Money\n\n");
 460     }
 461
 462     if (PAR(cr) && !((Flags & MD_PARTDEC) || EI_TPI(inputrec->eI)))
 463     {
 464         cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
 465                                            dddlb_opt,dlb_scale,
 466                                            ddcsx,ddcsy,ddcsz,
 467                                            mtop,inputrec,
 468                                            box,state->x,
 469                                            &ddbox,&npme_major);
 470
 471         make_dd_communicators(fplog,cr,dd_node_order);
 472
 473         /* Set overallocation to avoid frequent reallocation of arrays */
 474         set_over_alloc_dd(TRUE);
 475     }
 476     else
 477     {
 478         cr->duty = (DUTY_PP | DUTY_PME);
 479         npme_major = cr->nnodes;
 480
 481         if (inputrec->ePBC == epbcSCREW)
 482         {
 483             gmx_fatal(FARGS,
 484                       "pbc=%s is only implemented with domain decomposition",
 485                       epbc_names[inputrec->ePBC]);
 486         }
 487     }
 488
 489     if (PAR(cr))
 490     {
 491         /* After possible communicator splitting in make_dd_communicators.
 492          * we can set up the intra/inter node communication.
 493          */
 494         gmx_setup_nodecomm(fplog,cr);
 495     }
 496
 497     wcycle = wallcycle_init(fplog,resetstep,cr);
 498     if (PAR(cr))
 499     {
 500         /* Master synchronizes its value of reset_counters with all nodes
 501          * including PME only nodes */
 502         reset_counters = wcycle_get_reset_counters(wcycle);
 503         gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
 504         wcycle_set_reset_counters(wcycle, reset_counters);
 505     }
 506
 507
 508     snew(nrnb,1);
 509     if (cr->duty & DUTY_PP)
 510     {
 511         /* For domain decomposition we allocate dynamically
 512          * in dd_partition_system.
 513          */
 514         if (DOMAINDECOMP(cr))
 515         {
 516             bcast_state_setup(cr,state);
 517         }
 518         else
 519         {
 520             if (PAR(cr))
 521             {
 522                 if (!MASTER(cr))
 523                 {
 524                     snew(state,1);
 525                 }
 526                 bcast_state(cr,state,TRUE);
 527             }
 528         }
 529
 530         /* Dihedral Restraints */
 531         if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
 532         {
 533             init_dihres(fplog,mtop,inputrec,fcd);
 534         }
 535
 536         /* Initiate forcerecord */
 537         fr = mk_forcerec();
 538         init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
 539                       opt2fn("-table",nfile,fnm),
 540                       opt2fn("-tablep",nfile,fnm),
 541                       opt2fn("-tableb",nfile,fnm),FALSE,pforce);
 542
 543         /* version for PCA_NOT_READ_NODE (see md.c) */
 544         /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
 545           "nofile","nofile","nofile",FALSE,pforce);
 546           */
 547         fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
 548
 549         /* Initialize QM-MM */
 550         if(fr->bQMMM)
 551         {
 552             init_QMMMrec(cr,box,mtop,inputrec,fr);
 553         }
 554
 555         /* Initialize the mdatoms structure.
 556          * mdatoms is not filled with atom data,
 557          * as this can not be done now with domain decomposition.
 558          */
 559         mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
 560
 561         /* Initialize the virtual site communication */
 562         vsite = init_vsite(mtop,cr);
 563
 564         calc_shifts(box,fr->shift_vec);
 565
 566         /* With periodic molecules the charge groups should be whole at start up
 567          * and the virtual sites should not be far from their proper positions.
 568          */
 569         if (!inputrec->bContinuation && MASTER(cr) &&
 570             !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
 571         {
 572             /* Make molecules whole at start of run */
 573             if (fr->ePBC != epbcNONE)
 574             {
 575                 do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
 576             }
 577             if (vsite)
 578             {
 579                 /* Correct initial vsite positions are required
 580                  * for the initial distribution in the domain decomposition
 581                  * and for the initial shell prediction.
 582                  */
 583                 construct_vsites_mtop(fplog,vsite,mtop,state->x);
 584             }
 585         }
 586
 587         /* Initiate PPPM if necessary */
 588         if (fr->eeltype == eelPPPM)
 589         {
 590             if (mdatoms->nChargePerturbed)
 591             {
 592                 gmx_fatal(FARGS,"Free energy with %s is not implemented",
 593                           eel_names[fr->eeltype]);
 594             }
 595             status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
 596                                    getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
 597             if (status != 0)
 598             {
 599                 gmx_fatal(FARGS,"Error %d initializing PPPM",status);
 600             }
 601         }
 602
 603         if (EEL_PME(fr->eeltype))
 604         {
 605             ewaldcoeff = fr->ewaldcoeff;
 606             pmedata = &fr->pmedata;
 607         }
 608         else
 609         {
 610             pmedata = NULL;
 611         }
 612     }
 613     else
 614     {
 615         /* This is a PME only node */
 616
 617         /* We don't need the state */
 618         done_state(state);
 619
 620         ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
 621         snew(pmedata,1);
 622     }
 623
 624     /* Initiate PME if necessary,
 625      * either on all nodes or on dedicated PME nodes only. */
 626     if (EEL_PME(inputrec->coulombtype))
 627     {
 628         if (mdatoms)
 629         {
 630             nChargePerturbed = mdatoms->nChargePerturbed;
 631         }
 632         if (cr->npmenodes > 0)
 633         {
 634             /* The PME only nodes need to know nChargePerturbed */
 635             gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
 636         }
 637         if (cr->duty & DUTY_PME)
 638         {
 639             status = gmx_pme_init(pmedata,cr,npme_major,inputrec,
 640                                   mtop ? mtop->natoms : 0,nChargePerturbed,
 641                                   (Flags & MD_REPRODUCIBLE));
 642             if (status != 0)
 643             {
 644                 gmx_fatal(FARGS,"Error %d initializing PME",status);
 645             }
 646         }
 647     }
 648
 649
 650     if (integrator[inputrec->eI].func == do_md)
 651     {
 652         /* Turn on signal handling on all nodes */
 653         /*
 654          * (A user signal from the PME nodes (if any)
 655          * is communicated to the PP nodes.
 656          */
 657         signal_handler_install();
 658     }
 659
 660     if (cr->duty & DUTY_PP)
 661     {
 662         if (inputrec->ePull != epullNO)
 663         {
 664             /* Initialize pull code */
 665             init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
 666                       EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
 667         }
 668
 669         constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
 670
 671         if (DOMAINDECOMP(cr))
 672         {
 673             dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
 674                             Flags & MD_DDBONDCHECK,fr->cginfo_mb);
 675
 676             set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
 677
 678             setup_dd_grid(fplog,cr->dd);
 679         }
 680
 681         /* Now do whatever the user wants us to do (how flexible...) */
 682         integrator[inputrec->eI].func(fplog,cr,nfile,fnm,
 683                                       oenv,bVerbose,bCompact,
 684                                       nstglobalcomm,
 685                                       vsite,constr,
 686                                       nstepout,inputrec,mtop,
 687                                       fcd,state,
 688                                       mdatoms,nrnb,wcycle,ed,fr,
 689                                       repl_ex_nst,repl_ex_seed,
 690                                       cpt_period,max_hours,
 691                                       deviceOptions,
 692                                       Flags,
 693                                       &runtime);
 694
 695         if (inputrec->ePull != epullNO)
 696         {
 697             finish_pull(fplog,inputrec->pull);
 698         }
 699     }
 700     else
 701     {
 702         /* do PME only */
 703         gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
 704     }
 705
 706     if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
 707     {
 708         /* Some timing stats */
 709         if (MASTER(cr))
 710         {
 711             if (runtime.proc == 0)
 712             {
 713                 runtime.proc = runtime.real;
 714             }
 715         }
 716         else
 717         {
 718             runtime.real = 0;
 719         }
 720     }
 721
 722     wallcycle_stop(wcycle,ewcRUN);
 723
 724     /* Finish up, write some stuff
 725      * if rerunMD, don't write last frame again
 726      */
 727     finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
 728                inputrec,nrnb,wcycle,&runtime,
 729                EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
 730
 731     /* Does what it says */
 732     print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
 733
 734     /* Close logfile already here if we were appending to it */
 735     if (MASTER(cr) && (Flags & MD_APPENDFILES))
 736     {
 737         gmx_log_close(fplog);
 738     }
 739
 740     if(bGotStopNextStepSignal)
 741     {
 742         rc = 1;
 743     }
 744     else if(bGotStopNextNSStepSignal)
 745     {
 746         rc = 2;
 747     }
 748     else
 749     {
 750         rc = 0;
 751     }
 752
 753     return rc;
 754 }
 755
 756 void md_print_warning(const t_commrec *cr,FILE *fplog,const char *buf)
 757 {
 758     if (MASTER(cr))
 759     {
 760         fprintf(stderr,"\n%s\n",buf);
 761     }
 762     if (fplog)
 763     {
 764         fprintf(fplog,"\n%s\n",buf);
 765     }
 766 }