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