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 #ifdef __linux
  40 #define _GNU_SOURCE
  41 #include <sched.h>
  42 #include <sys/syscall.h>
  43 #endif
  44 #include <signal.h>
  45 #include <stdlib.h>
  46
  47 #include "typedefs.h"
  48 #include "smalloc.h"
  49 #include "sysstuff.h"
  50 #include "statutil.h"
  51 #include "mdrun.h"
  52 #include "network.h"
  53 #include "pull.h"
  54 #include "names.h"
  55 #include "disre.h"
  56 #include "orires.h"
  57 #include "dihre.h"
  58 #include "pme.h"
  59 #include "mdatoms.h"
  60 #include "repl_ex.h"
  61 #include "qmmm.h"
  62 #include "mpelogging.h"
  63 #include "domdec.h"
  64 #include "partdec.h"
  65 #include "coulomb.h"
  66 #include "constr.h"
  67 #include "mvdata.h"
  68 #include "checkpoint.h"
  69 #include "mtop_util.h"
  70 #include "sighandler.h"
  71 #include "tpxio.h"
  72 #include "txtdump.h"
  73 #include "pull_rotation.h"
  74 #include "md_openmm.h"
  75
  76 #ifdef GMX_LIB_MPI
  77 #include <mpi.h>
  78 #endif
  79 #ifdef GMX_THREAD_MPI
  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 #ifdef GMX_OPENMP
  92 #include <omp.h>
  93 #endif
  94
  95
  96 typedef struct {
  97     gmx_integrator_t *func;
  98 } gmx_intp_t;
  99
 100 /* The array should match the eI array in include/types/enums.h */
 101 #ifdef GMX_OPENMM  /* FIXME do_md_openmm needs fixing */
 102 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}};
 103 #else
 104 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}};
 105 #endif
 106
 107 gmx_large_int_t     deform_init_init_step_tpx;
 108 matrix              deform_init_box_tpx;
 109 #ifdef GMX_THREAD_MPI
 110 tMPI_Thread_mutex_t deform_init_box_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
 111 #endif
 112
 113
 114 #ifdef GMX_THREAD_MPI
 115 struct mdrunner_arglist
 116 {
 117     FILE *fplog;
 118     t_commrec *cr;
 119     int nfile;
 120     const t_filenm *fnm;
 121     output_env_t oenv;
 122     gmx_bool bVerbose;
 123     gmx_bool bCompact;
 124     int nstglobalcomm;
 125     ivec ddxyz;
 126     int dd_node_order;
 127     real rdd;
 128     real rconstr;
 129     const char *dddlb_opt;
 130     real dlb_scale;
 131     const char *ddcsx;
 132     const char *ddcsy;
 133     const char *ddcsz;
 134     int nstepout;
 135     int resetstep;
 136     int nmultisim;
 137     int repl_ex_nst;
 138     int repl_ex_seed;
 139     real pforce;
 140     real cpt_period;
 141     real max_hours;
 142     const char *deviceOptions;
 143     unsigned long Flags;
 144     int ret; /* return value */
 145 };
 146
 147
 148 /* The function used for spawning threads. Extracts the mdrunner()
 149    arguments from its one argument and calls mdrunner(), after making
 150    a commrec. */
 151 static void mdrunner_start_fn(void *arg)
 152 {
 153     struct mdrunner_arglist *mda=(struct mdrunner_arglist*)arg;
 154     struct mdrunner_arglist mc=*mda; /* copy the arg list to make sure
 155                                         that it's thread-local. This doesn't
 156                                         copy pointed-to items, of course,
 157                                         but those are all const. */
 158     t_commrec *cr;  /* we need a local version of this */
 159     FILE *fplog=NULL;
 160     t_filenm *fnm;
 161
 162     fnm = dup_tfn(mc.nfile, mc.fnm);
 163
 164     cr = init_par_threads(mc.cr);
 165
 166     if (MASTER(cr))
 167     {
 168         fplog=mc.fplog;
 169     }
 170
 171     mda->ret=mdrunner(cr->nnodes, fplog, cr, mc.nfile, fnm, mc.oenv,
 172                       mc.bVerbose, mc.bCompact, mc.nstglobalcomm,
 173                       mc.ddxyz, mc.dd_node_order, mc.rdd,
 174                       mc.rconstr, mc.dddlb_opt, mc.dlb_scale,
 175                       mc.ddcsx, mc.ddcsy, mc.ddcsz, mc.nstepout, mc.resetstep,
 176                       mc.nmultisim, mc.repl_ex_nst, mc.repl_ex_seed, mc.pforce,
 177                       mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.Flags);
 178 }
 179
 180 /* called by mdrunner() to start a specific number of threads (including
 181    the main thread) for thread-parallel runs. This in turn calls mdrunner()
 182    for each thread.
 183    All options besides nthreads are the same as for mdrunner(). */
 184 static t_commrec *mdrunner_start_threads(int nthreads,
 185               FILE *fplog,t_commrec *cr,int nfile,
 186               const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
 187               gmx_bool bCompact, int nstglobalcomm,
 188               ivec ddxyz,int dd_node_order,real rdd,real rconstr,
 189               const char *dddlb_opt,real dlb_scale,
 190               const char *ddcsx,const char *ddcsy,const char *ddcsz,
 191               int nstepout,int resetstep,int nmultisim,int repl_ex_nst,
 192               int repl_ex_seed, real pforce,real cpt_period, real max_hours,
 193               const char *deviceOptions, unsigned long Flags)
 194 {
 195     int ret;
 196     struct mdrunner_arglist *mda;
 197     t_commrec *crn; /* the new commrec */
 198     t_filenm *fnmn;
 199
 200     /* first check whether we even need to start tMPI */
 201     if (nthreads<2)
 202         return cr;
 203
 204     /* a few small, one-time, almost unavoidable memory leaks: */
 205     snew(mda,1);
 206     fnmn=dup_tfn(nfile, fnm);
 207
 208     /* fill the data structure to pass as void pointer to thread start fn */
 209     mda->fplog=fplog;
 210     mda->cr=cr;
 211     mda->nfile=nfile;
 212     mda->fnm=fnmn;
 213     mda->oenv=oenv;
 214     mda->bVerbose=bVerbose;
 215     mda->bCompact=bCompact;
 216     mda->nstglobalcomm=nstglobalcomm;
 217     mda->ddxyz[XX]=ddxyz[XX];
 218     mda->ddxyz[YY]=ddxyz[YY];
 219     mda->ddxyz[ZZ]=ddxyz[ZZ];
 220     mda->dd_node_order=dd_node_order;
 221     mda->rdd=rdd;
 222     mda->rconstr=rconstr;
 223     mda->dddlb_opt=dddlb_opt;
 224     mda->dlb_scale=dlb_scale;
 225     mda->ddcsx=ddcsx;
 226     mda->ddcsy=ddcsy;
 227     mda->ddcsz=ddcsz;
 228     mda->nstepout=nstepout;
 229     mda->resetstep=resetstep;
 230     mda->nmultisim=nmultisim;
 231     mda->repl_ex_nst=repl_ex_nst;
 232     mda->repl_ex_seed=repl_ex_seed;
 233     mda->pforce=pforce;
 234     mda->cpt_period=cpt_period;
 235     mda->max_hours=max_hours;
 236     mda->deviceOptions=deviceOptions;
 237     mda->Flags=Flags;
 238
 239     fprintf(stderr, "Starting %d threads\n",nthreads);
 240     fflush(stderr);
 241     /* now spawn new threads that start mdrunner_start_fn(), while
 242        the main thread returns */
 243     ret=tMPI_Init_fn(TRUE, nthreads, mdrunner_start_fn, (void*)(mda) );
 244     if (ret!=TMPI_SUCCESS)
 245         return NULL;
 246
 247     /* make a new comm_rec to reflect the new situation */
 248     crn=init_par_threads(cr);
 249     return crn;
 250 }
 251
 252
 253 /* Get the number of threads to use for thread-MPI based on how many
 254  * were requested, which algorithms we're using,
 255  * and how many particles there are.
 256  */
 257 static int get_nthreads_mpi(int nthreads_requested, t_inputrec *inputrec,
 258                             gmx_mtop_t *mtop)
 259 {
 260     int nthreads,nthreads_new;
 261     int min_atoms_per_thread;
 262     char *env;
 263
 264     nthreads = nthreads_requested;
 265
 266     /* determine # of hardware threads. */
 267     if (nthreads_requested < 1)
 268     {
 269         if ((env = getenv("GMX_MAX_THREADS")) != NULL)
 270         {
 271             nthreads = 0;
 272             sscanf(env,"%d",&nthreads);
 273             if (nthreads < 1)
 274             {
 275                 gmx_fatal(FARGS,"GMX_MAX_THREADS (%d) should be larger than 0",
 276                           nthreads);
 277             }
 278         }
 279         else
 280         {
 281             nthreads = tMPI_Thread_get_hw_number();
 282         }
 283     }
 284
 285     if (inputrec->eI == eiNM || EI_TPI(inputrec->eI))
 286     {
 287         /* Steps are divided over the nodes iso splitting the atoms */
 288         min_atoms_per_thread = 0;
 289     }
 290     else
 291     {
 292         min_atoms_per_thread = MIN_ATOMS_PER_THREAD;
 293     }
 294
 295     /* Check if an algorithm does not support parallel simulation.  */
 296     if (nthreads != 1 &&
 297         ( inputrec->eI == eiLBFGS ||
 298           inputrec->coulombtype == eelEWALD ) )
 299     {
 300         fprintf(stderr,"\nThe integration or electrostatics algorithm doesn't support parallel runs. Not starting any threads.\n");
 301         nthreads = 1;
 302     }
 303     else if (nthreads_requested < 1 &&
 304              mtop->natoms/nthreads < min_atoms_per_thread)
 305     {
 306         /* the thread number was chosen automatically, but there are too many
 307            threads (too few atoms per thread) */
 308         nthreads_new = max(1,mtop->natoms/min_atoms_per_thread);
 309
 310         if (nthreads_new > 8 || (nthreads == 8 && nthreads_new > 4))
 311         {
 312             /* Use only multiples of 4 above 8 threads
 313              * or with an 8-core processor
 314              * (to avoid 6 threads on 8 core processors with 4 real cores).
 315              */
 316             nthreads_new = (nthreads_new/4)*4;
 317         }
 318         else if (nthreads_new > 4)
 319         {
 320             /* Avoid 5 or 7 threads */
 321             nthreads_new = (nthreads_new/2)*2;
 322         }
 323
 324         nthreads = nthreads_new;
 325
 326         fprintf(stderr,"\n");
 327         fprintf(stderr,"NOTE: Parallelization is limited by the small number of atoms,\n");
 328         fprintf(stderr,"      only starting %d threads.\n",nthreads);
 329         fprintf(stderr,"      You can use the -nt option to optimize the number of threads.\n\n");
 330     }
 331     return nthreads;
 332 }
 333 #endif
 334
 335
 336 int mdrunner(int nthreads_requested, FILE *fplog,t_commrec *cr,int nfile,
 337              const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
 338              gmx_bool bCompact, int nstglobalcomm,
 339              ivec ddxyz,int dd_node_order,real rdd,real rconstr,
 340              const char *dddlb_opt,real dlb_scale,
 341              const char *ddcsx,const char *ddcsy,const char *ddcsz,
 342              int nstepout,int resetstep,int nmultisim,int repl_ex_nst,
 343              int repl_ex_seed, real pforce,real cpt_period,real max_hours,
 344              const char *deviceOptions, unsigned long Flags)
 345 {
 346     double     nodetime=0,realtime;
 347     t_inputrec *inputrec;
 348     t_state    *state=NULL;
 349     matrix     box;
 350     gmx_ddbox_t ddbox={0};
 351     int        npme_major,npme_minor;
 352     real       tmpr1,tmpr2;
 353     t_nrnb     *nrnb;
 354     gmx_mtop_t *mtop=NULL;
 355     t_mdatoms  *mdatoms=NULL;
 356     t_forcerec *fr=NULL;
 357     t_fcdata   *fcd=NULL;
 358     real       ewaldcoeff=0;
 359     gmx_pme_t  *pmedata=NULL;
 360     gmx_vsite_t *vsite=NULL;
 361     gmx_constr_t constr;
 362     int        i,m,nChargePerturbed=-1,status,nalloc;
 363     char       *gro;
 364     gmx_wallcycle_t wcycle;
 365     gmx_bool       bReadRNG,bReadEkin;
 366     int        list;
 367     gmx_runtime_t runtime;
 368     int        rc;
 369     gmx_large_int_t reset_counters;
 370     gmx_edsam_t ed=NULL;
 371     t_commrec   *cr_old=cr;
 372     int         nthreads_mpi=1;
 373     int         nthreads_pme=1;
 374
 375     /* CAUTION: threads may be started later on in this function, so
 376        cr doesn't reflect the final parallel state right now */
 377     snew(inputrec,1);
 378     snew(mtop,1);
 379
 380     if (bVerbose && SIMMASTER(cr))
 381     {
 382         fprintf(stderr,"Getting Loaded...\n");
 383     }
 384
 385     if (Flags & MD_APPENDFILES)
 386     {
 387         fplog = NULL;
 388     }
 389
 390     snew(state,1);
 391     if (MASTER(cr))
 392     {
 393         /* Read (nearly) all data required for the simulation */
 394         read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
 395
 396         /* NOW the threads will be started: */
 397 #ifdef GMX_THREAD_MPI
 398         nthreads_mpi = get_nthreads_mpi(nthreads_requested, inputrec, mtop);
 399
 400         if (nthreads_mpi > 1)
 401         {
 402             /* now start the threads. */
 403             cr=mdrunner_start_threads(nthreads_mpi, fplog, cr_old, nfile, fnm,
 404                                       oenv, bVerbose, bCompact, nstglobalcomm,
 405                                       ddxyz, dd_node_order, rdd, rconstr,
 406                                       dddlb_opt, dlb_scale, ddcsx, ddcsy, ddcsz,
 407                                       nstepout, resetstep, nmultisim,
 408                                       repl_ex_nst, repl_ex_seed, pforce,
 409                                       cpt_period, max_hours, deviceOptions,
 410                                       Flags);
 411             /* the main thread continues here with a new cr. We don't deallocate
 412                the old cr because other threads may still be reading it. */
 413             if (cr == NULL)
 414             {
 415                 gmx_comm("Failed to spawn threads");
 416             }
 417         }
 418 #endif
 419     }
 420     /* END OF CAUTION: cr is now reliable */
 421
 422     if (PAR(cr))
 423     {
 424         /* now broadcast everything to the non-master nodes/threads: */
 425         init_parallel(fplog, cr, inputrec, mtop);
 426     }
 427     if (fplog != NULL)
 428     {
 429         pr_inputrec(fplog,0,"Input Parameters",inputrec,FALSE);
 430     }
 431
 432     /* now make sure the state is initialized and propagated */
 433     set_state_entries(state,inputrec,cr->nnodes);
 434
 435     /* A parallel command line option consistency check that we can
 436        only do after any threads have started. */
 437     if (!PAR(cr) &&
 438         (ddxyz[XX] > 1 || ddxyz[YY] > 1 || ddxyz[ZZ] > 1 || cr->npmenodes > 0))
 439     {
 440         gmx_fatal(FARGS,
 441                   "The -dd or -npme option request a parallel simulation, "
 442 #ifndef GMX_MPI
 443                   "but mdrun was compiled without threads or MPI enabled"
 444 #else
 445 #ifdef GMX_THREAD_MPI
 446                   "but the number of threads (option -nt) is 1"
 447 #else
 448                   "but mdrun was not started through mpirun/mpiexec or only one process was requested through mpirun/mpiexec"
 449 #endif
 450 #endif
 451             );
 452     }
 453
 454     if ((Flags & MD_RERUN) &&
 455         (EI_ENERGY_MINIMIZATION(inputrec->eI) || eiNM == inputrec->eI))
 456     {
 457         gmx_fatal(FARGS, "The .mdp file specified an energy mininization or normal mode algorithm, and these are not compatible with mdrun -rerun");
 458     }
 459
 460     if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
 461     {
 462         /* All-vs-all loops do not work with domain decomposition */
 463         Flags |= MD_PARTDEC;
 464     }
 465
 466     if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
 467     {
 468         if (cr->npmenodes > 0)
 469         {
 470             if (!EEL_PME(inputrec->coulombtype))
 471             {
 472                 gmx_fatal_collective(FARGS,cr,NULL,
 473                                      "PME nodes are requested, but the system does not use PME electrostatics");
 474             }
 475             if (Flags & MD_PARTDEC)
 476             {
 477                 gmx_fatal_collective(FARGS,cr,NULL,
 478                                      "PME nodes are requested, but particle decomposition does not support separate PME nodes");
 479             }
 480         }
 481
 482         cr->npmenodes = 0;
 483     }
 484
 485 #ifdef GMX_FAHCORE
 486     fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
 487 #endif
 488
 489     /* NMR restraints must be initialized before load_checkpoint,
 490      * since with time averaging the history is added to t_state.
 491      * For proper consistency check we therefore need to extend
 492      * t_state here.
 493      * So the PME-only nodes (if present) will also initialize
 494      * the distance restraints.
 495      */
 496     snew(fcd,1);
 497
 498     /* This needs to be called before read_checkpoint to extend the state */
 499     init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
 500
 501     if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
 502     {
 503         if (PAR(cr) && !(Flags & MD_PARTDEC))
 504         {
 505             gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
 506         }
 507         /* Orientation restraints */
 508         if (MASTER(cr))
 509         {
 510             init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
 511                         state);
 512         }
 513     }
 514
 515     if (DEFORM(*inputrec))
 516     {
 517         /* Store the deform reference box before reading the checkpoint */
 518         if (SIMMASTER(cr))
 519         {
 520             copy_mat(state->box,box);
 521         }
 522         if (PAR(cr))
 523         {
 524             gmx_bcast(sizeof(box),box,cr);
 525         }
 526         /* Because we do not have the update struct available yet
 527          * in which the reference values should be stored,
 528          * we store them temporarily in static variables.
 529          * This should be thread safe, since they are only written once
 530          * and with identical values.
 531          */
 532 #ifdef GMX_THREAD_MPI
 533         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
 534 #endif
 535         deform_init_init_step_tpx = inputrec->init_step;
 536         copy_mat(box,deform_init_box_tpx);
 537 #ifdef GMX_THREAD_MPI
 538         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
 539 #endif
 540     }
 541
 542     if (opt2bSet("-cpi",nfile,fnm))
 543     {
 544         /* Check if checkpoint file exists before doing continuation.
 545          * This way we can use identical input options for the first and subsequent runs...
 546          */
 547         if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
 548         {
 549             load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
 550                             cr,Flags & MD_PARTDEC,ddxyz,
 551                             inputrec,state,&bReadRNG,&bReadEkin,
 552                             (Flags & MD_APPENDFILES),
 553                             (Flags & MD_APPENDFILESSET));
 554
 555             if (bReadRNG)
 556             {
 557                 Flags |= MD_READ_RNG;
 558             }
 559             if (bReadEkin)
 560             {
 561                 Flags |= MD_READ_EKIN;
 562             }
 563         }
 564     }
 565
 566     if (((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
 567 #ifdef GMX_THREAD_MPI
 568         /* With thread MPI only the master node/thread exists in mdrun.c,
 569          * therefore non-master nodes need to open the "seppot" log file here.
 570          */
 571         || (!MASTER(cr) && (Flags & MD_SEPPOT))
 572 #endif
 573         )
 574     {
 575         gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
 576                              Flags,&fplog);
 577     }
 578
 579     if (SIMMASTER(cr))
 580     {
 581         copy_mat(state->box,box);
 582     }
 583
 584     if (PAR(cr))
 585     {
 586         gmx_bcast(sizeof(box),box,cr);
 587     }
 588
 589     /* Essential dynamics */
 590     if (opt2bSet("-ei",nfile,fnm))
 591     {
 592         /* Open input and output files, allocate space for ED data structure */
 593         ed = ed_open(nfile,fnm,Flags,cr);
 594     }
 595
 596     if (bVerbose && SIMMASTER(cr))
 597     {
 598         fprintf(stderr,"Loaded with Money\n\n");
 599     }
 600
 601     if (PAR(cr) && !((Flags & MD_PARTDEC) ||
 602                      EI_TPI(inputrec->eI) ||
 603                      inputrec->eI == eiNM))
 604     {
 605         cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
 606                                            dddlb_opt,dlb_scale,
 607                                            ddcsx,ddcsy,ddcsz,
 608                                            mtop,inputrec,
 609                                            box,state->x,
 610                                            &ddbox,&npme_major,&npme_minor);
 611
 612         make_dd_communicators(fplog,cr,dd_node_order);
 613
 614         /* Set overallocation to avoid frequent reallocation of arrays */
 615         set_over_alloc_dd(TRUE);
 616     }
 617     else
 618     {
 619         /* PME, if used, is done on all nodes with 1D decomposition */
 620         cr->npmenodes = 0;
 621         cr->duty = (DUTY_PP | DUTY_PME);
 622         npme_major = 1;
 623         npme_minor = 1;
 624         if (!EI_TPI(inputrec->eI))
 625         {
 626             npme_major = cr->nnodes;
 627         }
 628
 629         if (inputrec->ePBC == epbcSCREW)
 630         {
 631             gmx_fatal(FARGS,
 632                       "pbc=%s is only implemented with domain decomposition",
 633                       epbc_names[inputrec->ePBC]);
 634         }
 635     }
 636
 637     if (PAR(cr))
 638     {
 639         /* After possible communicator splitting in make_dd_communicators.
 640          * we can set up the intra/inter node communication.
 641          */
 642         gmx_setup_nodecomm(fplog,cr);
 643     }
 644
 645     /* get number of OpenMP/PME threads
 646      * env variable should be read only on one node to make sure it is identical everywhere */
 647 #ifdef GMX_OPENMP
 648     if (EEL_PME(inputrec->coulombtype))
 649     {
 650         if (MASTER(cr))
 651         {
 652             char *ptr;
 653             if ((ptr=getenv("GMX_PME_NTHREADS")) != NULL)
 654             {
 655                 sscanf(ptr,"%d",&nthreads_pme);
 656             }
 657             if (fplog != NULL && nthreads_pme > 1)
 658             {
 659                 fprintf(fplog,"Using %d threads for PME\n",nthreads_pme);
 660             }
 661         }
 662         if (PAR(cr))
 663         {
 664             gmx_bcast_sim(sizeof(nthreads_pme),&nthreads_pme,cr);
 665         }
 666     }
 667 #endif
 668
 669     wcycle = wallcycle_init(fplog,resetstep,cr,nthreads_pme);
 670     if (PAR(cr))
 671     {
 672         /* Master synchronizes its value of reset_counters with all nodes
 673          * including PME only nodes */
 674         reset_counters = wcycle_get_reset_counters(wcycle);
 675         gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
 676         wcycle_set_reset_counters(wcycle, reset_counters);
 677     }
 678
 679
 680     snew(nrnb,1);
 681     if (cr->duty & DUTY_PP)
 682     {
 683         /* For domain decomposition we allocate dynamically
 684          * in dd_partition_system.
 685          */
 686         if (DOMAINDECOMP(cr))
 687         {
 688             bcast_state_setup(cr,state);
 689         }
 690         else
 691         {
 692             if (PAR(cr))
 693             {
 694                 bcast_state(cr,state,TRUE);
 695             }
 696         }
 697
 698         /* Dihedral Restraints */
 699         if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
 700         {
 701             init_dihres(fplog,mtop,inputrec,fcd);
 702         }
 703
 704         /* Initiate forcerecord */
 705         fr = mk_forcerec();
 706         init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
 707                       opt2fn("-table",nfile,fnm),
 708                       opt2fn("-tabletf",nfile,fnm),
 709                       opt2fn("-tablep",nfile,fnm),
 710                       opt2fn("-tableb",nfile,fnm),FALSE,pforce);
 711
 712         /* version for PCA_NOT_READ_NODE (see md.c) */
 713         /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
 714           "nofile","nofile","nofile","nofile",FALSE,pforce);
 715           */
 716         fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
 717
 718         /* Initialize QM-MM */
 719         if(fr->bQMMM)
 720         {
 721             init_QMMMrec(cr,box,mtop,inputrec,fr);
 722         }
 723
 724         /* Initialize the mdatoms structure.
 725          * mdatoms is not filled with atom data,
 726          * as this can not be done now with domain decomposition.
 727          */
 728         mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
 729
 730         /* Initialize the virtual site communication */
 731         vsite = init_vsite(mtop,cr);
 732
 733         calc_shifts(box,fr->shift_vec);
 734
 735         /* With periodic molecules the charge groups should be whole at start up
 736          * and the virtual sites should not be far from their proper positions.
 737          */
 738         if (!inputrec->bContinuation && MASTER(cr) &&
 739             !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
 740         {
 741             /* Make molecules whole at start of run */
 742             if (fr->ePBC != epbcNONE)
 743             {
 744                 do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
 745             }
 746             if (vsite)
 747             {
 748                 /* Correct initial vsite positions are required
 749                  * for the initial distribution in the domain decomposition
 750                  * and for the initial shell prediction.
 751                  */
 752                 construct_vsites_mtop(fplog,vsite,mtop,state->x);
 753             }
 754         }
 755
 756         if (EEL_PME(fr->eeltype))
 757         {
 758             ewaldcoeff = fr->ewaldcoeff;
 759             pmedata = &fr->pmedata;
 760         }
 761         else
 762         {
 763             pmedata = NULL;
 764         }
 765     }
 766     else
 767     {
 768         /* This is a PME only node */
 769
 770         /* We don't need the state */
 771         done_state(state);
 772
 773         ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
 774         snew(pmedata,1);
 775     }
 776
 777     /* Initiate PME if necessary,
 778      * either on all nodes or on dedicated PME nodes only. */
 779     if (EEL_PME(inputrec->coulombtype))
 780     {
 781         if (mdatoms)
 782         {
 783             nChargePerturbed = mdatoms->nChargePerturbed;
 784         }
 785         if (cr->npmenodes > 0)
 786         {
 787             /* The PME only nodes need to know nChargePerturbed */
 788             gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
 789         }
 790
 791
 792         /* Set CPU affinity. Can be important for performance.
 793            On some systems (e.g. Cray) CPU Affinity is set by default.
 794            But default assigning doesn't work (well) with only some ranks
 795            having threads. This causes very low performance.
 796            External tools have cumbersome syntax for setting affinity
 797            in the case that only some ranks have threads.
 798            Thus it is important that GROMACS sets the affinity internally at
 799            if only PME is using threads.
 800         */
 801
 802 #ifdef GMX_OPENMP
 803 #ifdef __linux
 804 #ifdef GMX_LIB_MPI
 805         {
 806             int core;
 807             MPI_Comm comm_intra; /* intra communicator (but different to nc.comm_intra includes PME nodes) */
 808             MPI_Comm_split(MPI_COMM_WORLD,gmx_hostname_num(),gmx_node_rank(),&comm_intra);
 809             int local_omp_nthreads = (cr->duty & DUTY_PME) ? nthreads_pme : 1; /* threads on this node */
 810             MPI_Scan(&local_omp_nthreads,&core, 1, MPI_INT, MPI_SUM, comm_intra);
 811             core-=local_omp_nthreads; /* make exclusive scan */
 812 #pragma omp parallel firstprivate(core) num_threads(local_omp_nthreads)
 813             {
 814                 cpu_set_t mask;
 815                 CPU_ZERO(&mask);
 816                 core+=omp_get_thread_num();
 817                 CPU_SET(core,&mask);
 818                 sched_setaffinity((pid_t) syscall (SYS_gettid),sizeof(cpu_set_t),&mask);
 819             }
 820         }
 821 #endif /*GMX_MPI*/
 822 #endif /*__linux*/
 823 #endif /*GMX_OPENMP*/
 824
 825         if (cr->duty & DUTY_PME)
 826         {
 827             status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
 828                                   mtop ? mtop->natoms : 0,nChargePerturbed,
 829                                   (Flags & MD_REPRODUCIBLE),nthreads_pme);
 830             if (status != 0)
 831             {
 832                 gmx_fatal(FARGS,"Error %d initializing PME",status);
 833             }
 834         }
 835     }
 836
 837
 838     if (integrator[inputrec->eI].func == do_md
 839 #ifdef GMX_OPENMM
 840         ||
 841         integrator[inputrec->eI].func == do_md_openmm
 842 #endif
 843         )
 844     {
 845         /* Turn on signal handling on all nodes */
 846         /*
 847          * (A user signal from the PME nodes (if any)
 848          * is communicated to the PP nodes.
 849          */
 850         signal_handler_install();
 851     }
 852
 853     if (cr->duty & DUTY_PP)
 854     {
 855         if (inputrec->ePull != epullNO)
 856         {
 857             /* Initialize pull code */
 858             init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
 859                       EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
 860         }
 861
 862         if (inputrec->bRot)
 863         {
 864            /* Initialize enforced rotation code */
 865            init_rot(fplog,inputrec,nfile,fnm,cr,state->x,state->box,mtop,oenv,
 866                     bVerbose,Flags);
 867         }
 868
 869         constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
 870
 871         if (DOMAINDECOMP(cr))
 872         {
 873             dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
 874                             Flags & MD_DDBONDCHECK,fr->cginfo_mb);
 875
 876             set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
 877
 878             setup_dd_grid(fplog,cr->dd);
 879         }
 880
 881         /* Now do whatever the user wants us to do (how flexible...) */
 882         integrator[inputrec->eI].func(fplog,cr,nfile,fnm,
 883                                       oenv,bVerbose,bCompact,
 884                                       nstglobalcomm,
 885                                       vsite,constr,
 886                                       nstepout,inputrec,mtop,
 887                                       fcd,state,
 888                                       mdatoms,nrnb,wcycle,ed,fr,
 889                                       repl_ex_nst,repl_ex_seed,
 890                                       cpt_period,max_hours,
 891                                       deviceOptions,
 892                                       Flags,
 893                                       &runtime);
 894
 895         if (inputrec->ePull != epullNO)
 896         {
 897             finish_pull(fplog,inputrec->pull);
 898         }
 899
 900         if (inputrec->bRot)
 901         {
 902             finish_rot(fplog,inputrec->rot);
 903         }
 904
 905     }
 906     else
 907     {
 908         /* do PME only */
 909         gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
 910     }
 911
 912     if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
 913     {
 914         /* Some timing stats */
 915         if (SIMMASTER(cr))
 916         {
 917             if (runtime.proc == 0)
 918             {
 919                 runtime.proc = runtime.real;
 920             }
 921         }
 922         else
 923         {
 924             runtime.real = 0;
 925         }
 926     }
 927
 928     wallcycle_stop(wcycle,ewcRUN);
 929
 930     /* Finish up, write some stuff
 931      * if rerunMD, don't write last frame again
 932      */
 933     finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
 934                inputrec,nrnb,wcycle,&runtime,
 935                EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
 936
 937     /* Does what it says */
 938     print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
 939
 940     /* Close logfile already here if we were appending to it */
 941     if (MASTER(cr) && (Flags & MD_APPENDFILES))
 942     {
 943         gmx_log_close(fplog);
 944     }
 945
 946     rc=(int)gmx_get_stop_condition();
 947
 948 #ifdef GMX_THREAD_MPI
 949     /* we need to join all threads. The sub-threads join when they
 950        exit this function, but the master thread needs to be told to
 951        wait for that. */
 952     if (PAR(cr) && MASTER(cr))
 953     {
 954         tMPI_Finalize();
 955     }
 956 #endif
 957
 958     return rc;
 959 }