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