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