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 "pme.h"
  62 #include "mdatoms.h"
  63 #include "repl_ex.h"
  64 #include "qmmm.h"
  65 #include "domdec.h"
  66 #include "partdec.h"
  67 #include "coulomb.h"
  68 #include "constr.h"
  69 #include "mvdata.h"
  70 #include "checkpoint.h"
  71 #include "mtop_util.h"
  72 #include "sighandler.h"
  73 #include "tpxio.h"
  74 #include "txtdump.h"
  75 #include "pull_rotation.h"
  76 #include "membed.h"
  77 #include "macros.h"
  78
  79 #include "gmx_omp.h"
  80
  81 #ifdef GMX_LIB_MPI
  82 #include <mpi.h>
  83 #endif
  84 #ifdef GMX_THREAD_MPI
  85 #include "tmpi.h"
  86 #endif
  87
  88 #ifdef GMX_FAHCORE
  89 #include "corewrap.h"
  90 #endif
  91
  92 #ifdef GMX_OPENMM
  93 #include "md_openmm.h"
  94 #endif
  95
  96
  97 typedef struct {
  98     gmx_integrator_t *func;
  99 } gmx_intp_t;
 100
 101 /* The array should match the eI array in include/types/enums.h */
 102 #ifdef GMX_OPENMM  /* FIXME do_md_openmm needs fixing */
 103 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}};
 104 #else
 105 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}};
 106 #endif
 107
 108 gmx_large_int_t     deform_init_init_step_tpx;
 109 matrix              deform_init_box_tpx;
 110 #ifdef GMX_THREAD_MPI
 111 tMPI_Thread_mutex_t deform_init_box_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
 112 #endif
 113
 114
 115 #ifdef GMX_THREAD_MPI
 116 struct mdrunner_arglist
 117 {
 118     FILE *fplog;
 119     t_commrec *cr;
 120     int nfile;
 121     const t_filenm *fnm;
 122     output_env_t oenv;
 123     gmx_bool bVerbose;
 124     gmx_bool bCompact;
 125     int nstglobalcomm;
 126     ivec ddxyz;
 127     int dd_node_order;
 128     real rdd;
 129     real rconstr;
 130     const char *dddlb_opt;
 131     real dlb_scale;
 132     const char *ddcsx;
 133     const char *ddcsy;
 134     const char *ddcsz;
 135     int nstepout;
 136     int resetstep;
 137     int nmultisim;
 138     int repl_ex_nst;
 139     int repl_ex_nex;
 140     int repl_ex_seed;
 141     real pforce;
 142     real cpt_period;
 143     real max_hours;
 144     const char *deviceOptions;
 145     unsigned long Flags;
 146     int ret; /* return value */
 147 };
 148
 149
 150 /* The function used for spawning threads. Extracts the mdrunner()
 151    arguments from its one argument and calls mdrunner(), after making
 152    a commrec. */
 153 static void mdrunner_start_fn(void *arg)
 154 {
 155     struct mdrunner_arglist *mda=(struct mdrunner_arglist*)arg;
 156     struct mdrunner_arglist mc=*mda; /* copy the arg list to make sure
 157                                         that it's thread-local. This doesn't
 158                                         copy pointed-to items, of course,
 159                                         but those are all const. */
 160     t_commrec *cr;  /* we need a local version of this */
 161     FILE *fplog=NULL;
 162     t_filenm *fnm;
 163
 164     fnm = dup_tfn(mc.nfile, mc.fnm);
 165
 166     cr = init_par_threads(mc.cr);
 167
 168     if (MASTER(cr))
 169     {
 170         fplog=mc.fplog;
 171     }
 172
 173     mda->ret=mdrunner(cr->nnodes, fplog, cr, mc.nfile, fnm, mc.oenv,
 174                       mc.bVerbose, mc.bCompact, mc.nstglobalcomm,
 175                       mc.ddxyz, mc.dd_node_order, mc.rdd,
 176                       mc.rconstr, mc.dddlb_opt, mc.dlb_scale,
 177                       mc.ddcsx, mc.ddcsy, mc.ddcsz, mc.nstepout, mc.resetstep,
 178                       mc.nmultisim, mc.repl_ex_nst, mc.repl_ex_nex, mc.repl_ex_seed, mc.pforce,
 179                       mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.Flags);
 180 }
 181
 182 /* called by mdrunner() to start a specific number of threads (including
 183    the main thread) for thread-parallel runs. This in turn calls mdrunner()
 184    for each thread.
 185    All options besides nthreads are the same as for mdrunner(). */
 186 static t_commrec *mdrunner_start_threads(int nthreads,
 187               FILE *fplog,t_commrec *cr,int nfile,
 188               const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
 189               gmx_bool bCompact, int nstglobalcomm,
 190               ivec ddxyz,int dd_node_order,real rdd,real rconstr,
 191               const char *dddlb_opt,real dlb_scale,
 192               const char *ddcsx,const char *ddcsy,const char *ddcsz,
 193               int nstepout,int resetstep,int nmultisim,int repl_ex_nst,
 194               int repl_ex_nex, int repl_ex_seed, real pforce,real cpt_period, real max_hours,
 195               const char *deviceOptions, unsigned long Flags)
 196 {
 197     int ret;
 198     struct mdrunner_arglist *mda;
 199     t_commrec *crn; /* the new commrec */
 200     t_filenm *fnmn;
 201
 202     /* first check whether we even need to start tMPI */
 203     if (nthreads<2)
 204         return cr;
 205
 206     /* a few small, one-time, almost unavoidable memory leaks: */
 207     snew(mda,1);
 208     fnmn=dup_tfn(nfile, fnm);
 209
 210     /* fill the data structure to pass as void pointer to thread start fn */
 211     mda->fplog=fplog;
 212     mda->cr=cr;
 213     mda->nfile=nfile;
 214     mda->fnm=fnmn;
 215     mda->oenv=oenv;
 216     mda->bVerbose=bVerbose;
 217     mda->bCompact=bCompact;
 218     mda->nstglobalcomm=nstglobalcomm;
 219     mda->ddxyz[XX]=ddxyz[XX];
 220     mda->ddxyz[YY]=ddxyz[YY];
 221     mda->ddxyz[ZZ]=ddxyz[ZZ];
 222     mda->dd_node_order=dd_node_order;
 223     mda->rdd=rdd;
 224     mda->rconstr=rconstr;
 225     mda->dddlb_opt=dddlb_opt;
 226     mda->dlb_scale=dlb_scale;
 227     mda->ddcsx=ddcsx;
 228     mda->ddcsy=ddcsy;
 229     mda->ddcsz=ddcsz;
 230     mda->nstepout=nstepout;
 231     mda->resetstep=resetstep;
 232     mda->nmultisim=nmultisim;
 233     mda->repl_ex_nst=repl_ex_nst;
 234     mda->repl_ex_nex=repl_ex_nex;
 235     mda->repl_ex_seed=repl_ex_seed;
 236     mda->pforce=pforce;
 237     mda->cpt_period=cpt_period;
 238     mda->max_hours=max_hours;
 239     mda->deviceOptions=deviceOptions;
 240     mda->Flags=Flags;
 241
 242     fprintf(stderr, "Starting %d threads\n",nthreads);
 243     fflush(stderr);
 244     /* now spawn new threads that start mdrunner_start_fn(), while
 245        the main thread returns */
 246     ret=tMPI_Init_fn(TRUE, nthreads, mdrunner_start_fn, (void*)(mda) );
 247     if (ret!=TMPI_SUCCESS)
 248         return NULL;
 249
 250     /* make a new comm_rec to reflect the new situation */
 251     crn=init_par_threads(cr);
 252     return crn;
 253 }
 254
 255
 256 /* Get the number of threads to use for thread-MPI based on how many
 257  * were requested, which algorithms we're using,
 258  * and how many particles there are.
 259  */
 260 static int get_nthreads_mpi(int nthreads_requested, t_inputrec *inputrec,
 261                             gmx_mtop_t *mtop)
 262 {
 263     int nthreads,nthreads_new;
 264     int min_atoms_per_thread;
 265     char *env;
 266
 267     nthreads = nthreads_requested;
 268
 269     /* determine # of hardware threads. */
 270     if (nthreads_requested < 1)
 271     {
 272         if ((env = getenv("GMX_MAX_THREADS")) != NULL)
 273         {
 274             nthreads = 0;
 275             sscanf(env,"%d",&nthreads);
 276             if (nthreads < 1)
 277             {
 278                 gmx_fatal(FARGS,"GMX_MAX_THREADS (%d) should be larger than 0",
 279                           nthreads);
 280             }
 281         }
 282         else
 283         {
 284             nthreads = tMPI_Thread_get_hw_number();
 285         }
 286     }
 287
 288     if (inputrec->eI == eiNM || EI_TPI(inputrec->eI))
 289     {
 290         /* Steps are divided over the nodes iso splitting the atoms */
 291         min_atoms_per_thread = 0;
 292     }
 293     else
 294     {
 295         min_atoms_per_thread = MIN_ATOMS_PER_THREAD;
 296     }
 297
 298     /* Check if an algorithm does not support parallel simulation.  */
 299     if (nthreads != 1 &&
 300         ( inputrec->eI == eiLBFGS ||
 301           inputrec->coulombtype == eelEWALD ) )
 302     {
 303         fprintf(stderr,"\nThe integration or electrostatics algorithm doesn't support parallel runs. Not starting any threads.\n");
 304         nthreads = 1;
 305     }
 306     else if (nthreads_requested < 1 &&
 307              mtop->natoms/nthreads < min_atoms_per_thread)
 308     {
 309         /* the thread number was chosen automatically, but there are too many
 310            threads (too few atoms per thread) */
 311         nthreads_new = max(1,mtop->natoms/min_atoms_per_thread);
 312
 313         if (nthreads_new > 8 || (nthreads == 8 && nthreads_new > 4))
 314         {
 315             /* Use only multiples of 4 above 8 threads
 316              * or with an 8-core processor
 317              * (to avoid 6 threads on 8 core processors with 4 real cores).
 318              */
 319             nthreads_new = (nthreads_new/4)*4;
 320         }
 321         else if (nthreads_new > 4)
 322         {
 323             /* Avoid 5 or 7 threads */
 324             nthreads_new = (nthreads_new/2)*2;
 325         }
 326
 327         nthreads = nthreads_new;
 328
 329         fprintf(stderr,"\n");
 330         fprintf(stderr,"NOTE: Parallelization is limited by the small number of atoms,\n");
 331         fprintf(stderr,"      only starting %d threads.\n",nthreads);
 332         fprintf(stderr,"      You can use the -nt option to optimize the number of threads.\n\n");
 333     }
 334     return nthreads;
 335 }
 336 #endif
 337
 338
 339 int mdrunner(int nthreads_requested, FILE *fplog,t_commrec *cr,int nfile,
 340              const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
 341              gmx_bool bCompact, int nstglobalcomm,
 342              ivec ddxyz,int dd_node_order,real rdd,real rconstr,
 343              const char *dddlb_opt,real dlb_scale,
 344              const char *ddcsx,const char *ddcsy,const char *ddcsz,
 345              int nstepout,int resetstep,int nmultisim, int repl_ex_nst, int repl_ex_nex,
 346              int repl_ex_seed, real pforce,real cpt_period,real max_hours,
 347              const char *deviceOptions, unsigned long Flags)
 348 {
 349     double     nodetime=0,realtime;
 350     t_inputrec *inputrec;
 351     t_state    *state=NULL;
 352     matrix     box;
 353     gmx_ddbox_t ddbox={0};
 354     int        npme_major,npme_minor;
 355     real       tmpr1,tmpr2;
 356     t_nrnb     *nrnb;
 357     gmx_mtop_t *mtop=NULL;
 358     t_mdatoms  *mdatoms=NULL;
 359     t_forcerec *fr=NULL;
 360     t_fcdata   *fcd=NULL;
 361     real       ewaldcoeff=0;
 362     gmx_pme_t  *pmedata=NULL;
 363     gmx_vsite_t *vsite=NULL;
 364     gmx_constr_t constr;
 365     int        i,m,nChargePerturbed=-1,status,nalloc;
 366     char       *gro;
 367     gmx_wallcycle_t wcycle;
 368     gmx_bool       bReadRNG,bReadEkin;
 369     int        list;
 370     gmx_runtime_t runtime;
 371     int        rc;
 372     gmx_large_int_t reset_counters;
 373     gmx_edsam_t ed=NULL;
 374     t_commrec   *cr_old=cr;
 375     int         nthreads_mpi=1;
 376     int         nthreads_pme=1;
 377     gmx_membed_t membed=NULL;
 378
 379     /* CAUTION: threads may be started later on in this function, so
 380        cr doesn't reflect the final parallel state right now */
 381     snew(inputrec,1);
 382     snew(mtop,1);
 383
 384     if (bVerbose && SIMMASTER(cr))
 385     {
 386         fprintf(stderr,"Getting Loaded...\n");
 387     }
 388
 389     if (Flags & MD_APPENDFILES)
 390     {
 391         fplog = NULL;
 392     }
 393
 394     snew(state,1);
 395     if (MASTER(cr))
 396     {
 397         /* Read (nearly) all data required for the simulation */
 398         read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
 399
 400         /* NOW the threads will be started: */
 401 #ifdef GMX_THREAD_MPI
 402         nthreads_mpi = get_nthreads_mpi(nthreads_requested, inputrec, mtop);
 403
 404         if (nthreads_mpi > 1)
 405         {
 406             /* now start the threads. */
 407             cr=mdrunner_start_threads(nthreads_mpi, fplog, cr_old, nfile, fnm,
 408                                       oenv, bVerbose, bCompact, nstglobalcomm,
 409                                       ddxyz, dd_node_order, rdd, rconstr,
 410                                       dddlb_opt, dlb_scale, ddcsx, ddcsy, ddcsz,
 411                                       nstepout, resetstep, nmultisim,
 412                                       repl_ex_nst, repl_ex_nex, repl_ex_seed, pforce,
 413                                       cpt_period, max_hours, deviceOptions,
 414                                       Flags);
 415             /* the main thread continues here with a new cr. We don't deallocate
 416                the old cr because other threads may still be reading it. */
 417             if (cr == NULL)
 418             {
 419                 gmx_comm("Failed to spawn threads");
 420             }
 421         }
 422 #endif
 423     }
 424     /* END OF CAUTION: cr is now reliable */
 425
 426     /* g_membed initialisation *
 427      * Because we change the mtop, init_membed is called before the init_parallel *
 428      * (in case we ever want to make it run in parallel) */
 429     if (opt2bSet("-membed",nfile,fnm))
 430     {
 431         if (MASTER(cr))
 432         {
 433             fprintf(stderr,"Initializing membed");
 434         }
 435         membed = init_membed(fplog,nfile,fnm,mtop,inputrec,state,cr,&cpt_period);
 436     }
 437
 438     if (PAR(cr))
 439     {
 440         /* now broadcast everything to the non-master nodes/threads: */
 441         init_parallel(fplog, cr, inputrec, mtop);
 442     }
 443     if (fplog != NULL)
 444     {
 445         pr_inputrec(fplog,0,"Input Parameters",inputrec,FALSE);
 446     }
 447
 448     /* now make sure the state is initialized and propagated */
 449     set_state_entries(state,inputrec,cr->nnodes);
 450
 451     /* remove when vv and rerun works correctly! */
 452     if (PAR(cr) && EI_VV(inputrec->eI) && ((Flags & MD_RERUN) || (Flags & MD_RERUN_VSITE)))
 453     {
 454         gmx_fatal(FARGS,
 455                   "Currently can't do velocity verlet with rerun in parallel.");
 456     }
 457
 458     /* A parallel command line option consistency check that we can
 459        only do after any threads have started. */
 460     if (!PAR(cr) &&
 461         (ddxyz[XX] > 1 || ddxyz[YY] > 1 || ddxyz[ZZ] > 1 || cr->npmenodes > 0))
 462     {
 463         gmx_fatal(FARGS,
 464                   "The -dd or -npme option request a parallel simulation, "
 465 #ifndef GMX_MPI
 466                   "but mdrun was compiled without threads or MPI enabled"
 467 #else
 468 #ifdef GMX_THREAD_MPI
 469                   "but the number of threads (option -nt) is 1"
 470 #else
 471                   "but mdrun was not started through mpirun/mpiexec or only one process was requested through mpirun/mpiexec"
 472 #endif
 473 #endif
 474             );
 475     }
 476
 477     if ((Flags & MD_RERUN) &&
 478         (EI_ENERGY_MINIMIZATION(inputrec->eI) || eiNM == inputrec->eI))
 479     {
 480         gmx_fatal(FARGS, "The .mdp file specified an energy mininization or normal mode algorithm, and these are not compatible with mdrun -rerun");
 481     }
 482
 483     if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
 484     {
 485         /* All-vs-all loops do not work with domain decomposition */
 486         Flags |= MD_PARTDEC;
 487     }
 488
 489     if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
 490     {
 491         if (cr->npmenodes > 0)
 492         {
 493             if (!EEL_PME(inputrec->coulombtype))
 494             {
 495                 gmx_fatal_collective(FARGS,cr,NULL,
 496                                      "PME nodes are requested, but the system does not use PME electrostatics");
 497             }
 498             if (Flags & MD_PARTDEC)
 499             {
 500                 gmx_fatal_collective(FARGS,cr,NULL,
 501                                      "PME nodes are requested, but particle decomposition does not support separate PME nodes");
 502             }
 503         }
 504
 505         cr->npmenodes = 0;
 506     }
 507
 508 #ifdef GMX_FAHCORE
 509     fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
 510 #endif
 511
 512     /* NMR restraints must be initialized before load_checkpoint,
 513      * since with time averaging the history is added to t_state.
 514      * For proper consistency check we therefore need to extend
 515      * t_state here.
 516      * So the PME-only nodes (if present) will also initialize
 517      * the distance restraints.
 518      */
 519     snew(fcd,1);
 520
 521     /* This needs to be called before read_checkpoint to extend the state */
 522     init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
 523
 524     if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
 525     {
 526         if (PAR(cr) && !(Flags & MD_PARTDEC))
 527         {
 528             gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
 529         }
 530         /* Orientation restraints */
 531         if (MASTER(cr))
 532         {
 533             init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
 534                         state);
 535         }
 536     }
 537
 538     if (DEFORM(*inputrec))
 539     {
 540         /* Store the deform reference box before reading the checkpoint */
 541         if (SIMMASTER(cr))
 542         {
 543             copy_mat(state->box,box);
 544         }
 545         if (PAR(cr))
 546         {
 547             gmx_bcast(sizeof(box),box,cr);
 548         }
 549         /* Because we do not have the update struct available yet
 550          * in which the reference values should be stored,
 551          * we store them temporarily in static variables.
 552          * This should be thread safe, since they are only written once
 553          * and with identical values.
 554          */
 555 #ifdef GMX_THREAD_MPI
 556         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
 557 #endif
 558         deform_init_init_step_tpx = inputrec->init_step;
 559         copy_mat(box,deform_init_box_tpx);
 560 #ifdef GMX_THREAD_MPI
 561         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
 562 #endif
 563     }
 564
 565     if (opt2bSet("-cpi",nfile,fnm))
 566     {
 567         /* Check if checkpoint file exists before doing continuation.
 568          * This way we can use identical input options for the first and subsequent runs...
 569          */
 570         if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
 571         {
 572             load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
 573                             cr,Flags & MD_PARTDEC,ddxyz,
 574                             inputrec,state,&bReadRNG,&bReadEkin,
 575                             (Flags & MD_APPENDFILES),
 576                             (Flags & MD_APPENDFILESSET));
 577
 578             if (bReadRNG)
 579             {
 580                 Flags |= MD_READ_RNG;
 581             }
 582             if (bReadEkin)
 583             {
 584                 Flags |= MD_READ_EKIN;
 585             }
 586         }
 587     }
 588
 589     if (((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
 590 #ifdef GMX_THREAD_MPI
 591         /* With thread MPI only the master node/thread exists in mdrun.c,
 592          * therefore non-master nodes need to open the "seppot" log file here.
 593          */
 594         || (!MASTER(cr) && (Flags & MD_SEPPOT))
 595 #endif
 596         )
 597     {
 598         gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
 599                              Flags,&fplog);
 600     }
 601
 602     if (SIMMASTER(cr))
 603     {
 604         copy_mat(state->box,box);
 605     }
 606
 607     if (PAR(cr))
 608     {
 609         gmx_bcast(sizeof(box),box,cr);
 610     }
 611
 612     /* Essential dynamics */
 613     if (opt2bSet("-ei",nfile,fnm))
 614     {
 615         /* Open input and output files, allocate space for ED data structure */
 616         ed = ed_open(nfile,fnm,Flags,cr);
 617     }
 618
 619     if (bVerbose && SIMMASTER(cr))
 620     {
 621         fprintf(stderr,"Loaded with Money\n\n");
 622     }
 623
 624     if (PAR(cr) && !((Flags & MD_PARTDEC) ||
 625                      EI_TPI(inputrec->eI) ||
 626                      inputrec->eI == eiNM))
 627     {
 628         cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
 629                                            dddlb_opt,dlb_scale,
 630                                            ddcsx,ddcsy,ddcsz,
 631                                            mtop,inputrec,
 632                                            box,state->x,
 633                                            &ddbox,&npme_major,&npme_minor);
 634
 635         make_dd_communicators(fplog,cr,dd_node_order);
 636
 637         /* Set overallocation to avoid frequent reallocation of arrays */
 638         set_over_alloc_dd(TRUE);
 639     }
 640     else
 641     {
 642         /* PME, if used, is done on all nodes with 1D decomposition */
 643         cr->npmenodes = 0;
 644         cr->duty = (DUTY_PP | DUTY_PME);
 645         npme_major = 1;
 646         npme_minor = 1;
 647         if (!EI_TPI(inputrec->eI))
 648         {
 649             npme_major = cr->nnodes;
 650         }
 651
 652         if (inputrec->ePBC == epbcSCREW)
 653         {
 654             gmx_fatal(FARGS,
 655                       "pbc=%s is only implemented with domain decomposition",
 656                       epbc_names[inputrec->ePBC]);
 657         }
 658     }
 659
 660     if (PAR(cr))
 661     {
 662         /* After possible communicator splitting in make_dd_communicators.
 663          * we can set up the intra/inter node communication.
 664          */
 665         gmx_setup_nodecomm(fplog,cr);
 666     }
 667
 668     /* get number of OpenMP/PME threads
 669      * env variable should be read only on one node to make sure it is identical everywhere */
 670 #ifdef GMX_OPENMP
 671     if (EEL_PME(inputrec->coulombtype))
 672     {
 673         if (MASTER(cr))
 674         {
 675             char *ptr;
 676             if ((ptr=getenv("GMX_PME_NTHREADS")) != NULL)
 677             {
 678                 sscanf(ptr,"%d",&nthreads_pme);
 679             }
 680             if (fplog != NULL && nthreads_pme > 1)
 681             {
 682                 fprintf(fplog,"Using %d threads for PME\n",nthreads_pme);
 683             }
 684         }
 685         if (PAR(cr))
 686         {
 687             gmx_bcast_sim(sizeof(nthreads_pme),&nthreads_pme,cr);
 688         }
 689     }
 690 #endif
 691
 692     wcycle = wallcycle_init(fplog,resetstep,cr,nthreads_pme);
 693     if (PAR(cr))
 694     {
 695         /* Master synchronizes its value of reset_counters with all nodes
 696          * including PME only nodes */
 697         reset_counters = wcycle_get_reset_counters(wcycle);
 698         gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
 699         wcycle_set_reset_counters(wcycle, reset_counters);
 700     }
 701
 702
 703     snew(nrnb,1);
 704     if (cr->duty & DUTY_PP)
 705     {
 706         /* For domain decomposition we allocate dynamically
 707          * in dd_partition_system.
 708          */
 709         if (DOMAINDECOMP(cr))
 710         {
 711             bcast_state_setup(cr,state);
 712         }
 713         else
 714         {
 715             if (PAR(cr))
 716             {
 717                 bcast_state(cr,state,TRUE);
 718             }
 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+=gmx_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, inputrec->fepvals->init_lambda,
 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,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_nex,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 }