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