src/gromacs/mdlib/tpi.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  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <string.h>
  41 #include <time.h>
  42 #include <math.h>
  43 #include "sysstuff.h"
  44 #include "string2.h"
  45 #include "network.h"
  46 #include "confio.h"
  47 #include "smalloc.h"
  48 #include "nrnb.h"
  49 #include "main.h"
  50 #include "chargegroup.h"
  51 #include "force.h"
  52 #include "macros.h"
  53 #include "random.h"
  54 #include "names.h"
  55 #include "gmx_fatal.h"
  56 #include "txtdump.h"
  57 #include "typedefs.h"
  58 #include "update.h"
  59 #include "random.h"
  60 #include "constr.h"
  61 #include "vec.h"
  62 #include "statutil.h"
  63 #include "tgroup.h"
  64 #include "mdebin.h"
  65 #include "vsite.h"
  66 #include "force.h"
  67 #include "mdrun.h"
  68 #include "domdec.h"
  69 #include "partdec.h"
  70 #include "gmx_random.h"
  71 #include "physics.h"
  72 #include "xvgr.h"
  73 #include "mdatoms.h"
  74 #include "ns.h"
  75 #include "gmx_wallcycle.h"
  76 #include "mtop_util.h"
  77 #include "gmxfio.h"
  78 #include "pme.h"
  79 #include "gbutil.h"
  80
  81 #ifdef GMX_X86_SSE2
  82 #include "gmx_x86_sse2.h"
  83 #endif
  84
  85
  86 static void global_max(t_commrec *cr, int *n)
  87 {
  88     int *sum, i;
  89
  90     snew(sum, cr->nnodes);
  91     sum[cr->nodeid] = *n;
  92     gmx_sumi(cr->nnodes, sum, cr);
  93     for (i = 0; i < cr->nnodes; i++)
  94     {
  95         *n = max(*n, sum[i]);
  96     }
  97
  98     sfree(sum);
  99 }
 100
 101 static void realloc_bins(double **bin, int *nbin, int nbin_new)
 102 {
 103     int i;
 104
 105     if (nbin_new != *nbin)
 106     {
 107         srenew(*bin, nbin_new);
 108         for (i = *nbin; i < nbin_new; i++)
 109         {
 110             (*bin)[i] = 0;
 111         }
 112         *nbin = nbin_new;
 113     }
 114 }
 115
 116 double do_tpi(FILE *fplog, t_commrec *cr,
 117               int nfile, const t_filenm fnm[],
 118               const output_env_t oenv, gmx_bool bVerbose, gmx_bool gmx_unused bCompact,
 119               int gmx_unused nstglobalcomm,
 120               gmx_vsite_t gmx_unused *vsite, gmx_constr_t gmx_unused constr,
 121               int gmx_unused stepout,
 122               t_inputrec *inputrec,
 123               gmx_mtop_t *top_global, t_fcdata *fcd,
 124               t_state *state,
 125               t_mdatoms *mdatoms,
 126               t_nrnb *nrnb, gmx_wallcycle_t wcycle,
 127               gmx_edsam_t gmx_unused ed,
 128               t_forcerec *fr,
 129               int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
 130               gmx_membed_t gmx_unused membed,
 131               real gmx_unused cpt_period, real gmx_unused max_hours,
 132               const char gmx_unused *deviceOptions,
 133               unsigned long gmx_unused Flags,
 134               gmx_runtime_t *runtime)
 135 {
 136     const char     *TPI = "Test Particle Insertion";
 137     gmx_localtop_t *top;
 138     gmx_groups_t   *groups;
 139     gmx_enerdata_t *enerd;
 140     rvec           *f;
 141     real            lambda, t, temp, beta, drmax, epot;
 142     double          embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
 143     t_trxstatus    *status;
 144     t_trxframe      rerun_fr;
 145     gmx_bool        bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS, bOurStep;
 146     tensor          force_vir, shake_vir, vir, pres;
 147     int             cg_tp, a_tp0, a_tp1, ngid, gid_tp, nener, e;
 148     rvec           *x_mol;
 149     rvec            mu_tot, x_init, dx, x_tp;
 150     int             nnodes, frame, nsteps, step;
 151     int             i, start, end;
 152     gmx_rng_t       tpi_rand;
 153     FILE           *fp_tpi = NULL;
 154     char           *ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN];
 155     double          dbl, dump_ener;
 156     gmx_bool        bCavity;
 157     int             nat_cavity  = 0, d;
 158     real           *mass_cavity = NULL, mass_tot;
 159     int             nbin;
 160     double          invbinw, *bin, refvolshift, logV, bUlogV;
 161     real            dvdl, prescorr, enercorr, dvdlcorr;
 162     gmx_bool        bEnergyOutOfBounds;
 163     const char     *tpid_leg[2] = {"direct", "reweighted"};
 164
 165     /* Since there is no upper limit to the insertion energies,
 166      * we need to set an upper limit for the distribution output.
 167      */
 168     real bU_bin_limit      = 50;
 169     real bU_logV_bin_limit = bU_bin_limit + 10;
 170
 171     nnodes = cr->nnodes;
 172
 173     top = gmx_mtop_generate_local_top(top_global, inputrec);
 174
 175     groups = &top_global->groups;
 176
 177     bCavity = (inputrec->eI == eiTPIC);
 178     if (bCavity)
 179     {
 180         ptr = getenv("GMX_TPIC_MASSES");
 181         if (ptr == NULL)
 182         {
 183             nat_cavity = 1;
 184         }
 185         else
 186         {
 187             /* Read (multiple) masses from env var GMX_TPIC_MASSES,
 188              * The center of mass of the last atoms is then used for TPIC.
 189              */
 190             nat_cavity = 0;
 191             while (sscanf(ptr, "%lf%n", &dbl, &i) > 0)
 192             {
 193                 srenew(mass_cavity, nat_cavity+1);
 194                 mass_cavity[nat_cavity] = dbl;
 195                 fprintf(fplog, "mass[%d] = %f\n",
 196                         nat_cavity+1, mass_cavity[nat_cavity]);
 197                 nat_cavity++;
 198                 ptr += i;
 199             }
 200             if (nat_cavity == 0)
 201             {
 202                 gmx_fatal(FARGS, "Found %d masses in GMX_TPIC_MASSES", nat_cavity);
 203             }
 204         }
 205     }
 206
 207     /*
 208        init_em(fplog,TPI,inputrec,&lambda,nrnb,mu_tot,
 209        state->box,fr,mdatoms,top,cr,nfile,fnm,NULL,NULL);*/
 210     /* We never need full pbc for TPI */
 211     fr->ePBC = epbcXYZ;
 212     /* Determine the temperature for the Boltzmann weighting */
 213     temp = inputrec->opts.ref_t[0];
 214     if (fplog)
 215     {
 216         for (i = 1; (i < inputrec->opts.ngtc); i++)
 217         {
 218             if (inputrec->opts.ref_t[i] != temp)
 219             {
 220                 fprintf(fplog, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 221                 fprintf(stderr, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 222             }
 223         }
 224         fprintf(fplog,
 225                 "\n  The temperature for test particle insertion is %.3f K\n\n",
 226                 temp);
 227     }
 228     beta = 1.0/(BOLTZ*temp);
 229
 230     /* Number of insertions per frame */
 231     nsteps = inputrec->nsteps;
 232
 233     /* Use the same neighborlist with more insertions points
 234      * in a sphere of radius drmax around the initial point
 235      */
 236     /* This should be a proper mdp parameter */
 237     drmax = inputrec->rtpi;
 238
 239     /* An environment variable can be set to dump all configurations
 240      * to pdb with an insertion energy <= this value.
 241      */
 242     dump_pdb  = getenv("GMX_TPI_DUMP");
 243     dump_ener = 0;
 244     if (dump_pdb)
 245     {
 246         sscanf(dump_pdb, "%lf", &dump_ener);
 247     }
 248
 249     atoms2md(top_global, inputrec, 0, NULL, 0, top_global->natoms, mdatoms);
 250     update_mdatoms(mdatoms, inputrec->fepvals->init_lambda);
 251
 252     snew(enerd, 1);
 253     init_enerdata(groups->grps[egcENER].nr, inputrec->fepvals->n_lambda, enerd);
 254     snew(f, top_global->natoms);
 255
 256     /* Print to log file  */
 257     runtime_start(runtime);
 258     print_date_and_time(fplog, cr->nodeid,
 259                         "Started Test Particle Insertion", runtime);
 260     wallcycle_start(wcycle, ewcRUN);
 261
 262     /* The last charge group is the group to be inserted */
 263     cg_tp = top->cgs.nr - 1;
 264     a_tp0 = top->cgs.index[cg_tp];
 265     a_tp1 = top->cgs.index[cg_tp+1];
 266     if (debug)
 267     {
 268         fprintf(debug, "TPI cg %d, atoms %d-%d\n", cg_tp, a_tp0, a_tp1);
 269     }
 270     if (a_tp1 - a_tp0 > 1 &&
 271         (inputrec->rlist < inputrec->rcoulomb ||
 272          inputrec->rlist < inputrec->rvdw))
 273     {
 274         gmx_fatal(FARGS, "Can not do TPI for multi-atom molecule with a twin-range cut-off");
 275     }
 276     snew(x_mol, a_tp1-a_tp0);
 277
 278     bDispCorr = (inputrec->eDispCorr != edispcNO);
 279     bCharge   = FALSE;
 280     for (i = a_tp0; i < a_tp1; i++)
 281     {
 282         /* Copy the coordinates of the molecule to be insterted */
 283         copy_rvec(state->x[i], x_mol[i-a_tp0]);
 284         /* Check if we need to print electrostatic energies */
 285         bCharge |= (mdatoms->chargeA[i] != 0 ||
 286                     (mdatoms->chargeB && mdatoms->chargeB[i] != 0));
 287     }
 288     bRFExcl = (bCharge && EEL_RF(fr->eeltype) && fr->eeltype != eelRF_NEC);
 289
 290     calc_cgcm(fplog, cg_tp, cg_tp+1, &(top->cgs), state->x, fr->cg_cm);
 291     if (bCavity)
 292     {
 293         if (norm(fr->cg_cm[cg_tp]) > 0.5*inputrec->rlist && fplog)
 294         {
 295             fprintf(fplog, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
 296             fprintf(stderr, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
 297         }
 298     }
 299     else
 300     {
 301         /* Center the molecule to be inserted at zero */
 302         for (i = 0; i < a_tp1-a_tp0; i++)
 303         {
 304             rvec_dec(x_mol[i], fr->cg_cm[cg_tp]);
 305         }
 306     }
 307
 308     if (fplog)
 309     {
 310         fprintf(fplog, "\nWill insert %d atoms %s partial charges\n",
 311                 a_tp1-a_tp0, bCharge ? "with" : "without");
 312
 313         fprintf(fplog, "\nWill insert %d times in each frame of %s\n",
 314                 nsteps, opt2fn("-rerun", nfile, fnm));
 315     }
 316
 317     if (!bCavity)
 318     {
 319         if (inputrec->nstlist > 1)
 320         {
 321             if (drmax == 0 && a_tp1-a_tp0 == 1)
 322             {
 323                 gmx_fatal(FARGS, "Re-using the neighborlist %d times for insertions of a single atom in a sphere of radius %f does not make sense", inputrec->nstlist, drmax);
 324             }
 325             if (fplog)
 326             {
 327                 fprintf(fplog, "Will use the same neighborlist for %d insertions in a sphere of radius %f\n", inputrec->nstlist, drmax);
 328             }
 329         }
 330     }
 331     else
 332     {
 333         if (fplog)
 334         {
 335             fprintf(fplog, "Will insert randomly in a sphere of radius %f around the center of the cavity\n", drmax);
 336         }
 337     }
 338
 339     ngid   = groups->grps[egcENER].nr;
 340     gid_tp = GET_CGINFO_GID(fr->cginfo[cg_tp]);
 341     nener  = 1 + ngid;
 342     if (bDispCorr)
 343     {
 344         nener += 1;
 345     }
 346     if (bCharge)
 347     {
 348         nener += ngid;
 349         if (bRFExcl)
 350         {
 351             nener += 1;
 352         }
 353         if (EEL_FULL(fr->eeltype))
 354         {
 355             nener += 1;
 356         }
 357     }
 358     snew(sum_UgembU, nener);
 359
 360     /* Initialize random generator */
 361     tpi_rand = gmx_rng_init(inputrec->ld_seed);
 362
 363     if (MASTER(cr))
 364     {
 365         fp_tpi = xvgropen(opt2fn("-tpi", nfile, fnm),
 366                           "TPI energies", "Time (ps)",
 367                           "(kJ mol\\S-1\\N) / (nm\\S3\\N)", oenv);
 368         xvgr_subtitle(fp_tpi, "f. are averages over one frame", oenv);
 369         snew(leg, 4+nener);
 370         e = 0;
 371         sprintf(str, "-kT log(<Ve\\S-\\betaU\\N>/<V>)");
 372         leg[e++] = strdup(str);
 373         sprintf(str, "f. -kT log<e\\S-\\betaU\\N>");
 374         leg[e++] = strdup(str);
 375         sprintf(str, "f. <e\\S-\\betaU\\N>");
 376         leg[e++] = strdup(str);
 377         sprintf(str, "f. V");
 378         leg[e++] = strdup(str);
 379         sprintf(str, "f. <Ue\\S-\\betaU\\N>");
 380         leg[e++] = strdup(str);
 381         for (i = 0; i < ngid; i++)
 382         {
 383             sprintf(str, "f. <U\\sVdW %s\\Ne\\S-\\betaU\\N>",
 384                     *(groups->grpname[groups->grps[egcENER].nm_ind[i]]));
 385             leg[e++] = strdup(str);
 386         }
 387         if (bDispCorr)
 388         {
 389             sprintf(str, "f. <U\\sdisp c\\Ne\\S-\\betaU\\N>");
 390             leg[e++] = strdup(str);
 391         }
 392         if (bCharge)
 393         {
 394             for (i = 0; i < ngid; i++)
 395             {
 396                 sprintf(str, "f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
 397                         *(groups->grpname[groups->grps[egcENER].nm_ind[i]]));
 398                 leg[e++] = strdup(str);
 399             }
 400             if (bRFExcl)
 401             {
 402                 sprintf(str, "f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
 403                 leg[e++] = strdup(str);
 404             }
 405             if (EEL_FULL(fr->eeltype))
 406             {
 407                 sprintf(str, "f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
 408                 leg[e++] = strdup(str);
 409             }
 410         }
 411         xvgr_legend(fp_tpi, 4+nener, (const char**)leg, oenv);
 412         for (i = 0; i < 4+nener; i++)
 413         {
 414             sfree(leg[i]);
 415         }
 416         sfree(leg);
 417     }
 418     clear_rvec(x_init);
 419     V_all     = 0;
 420     VembU_all = 0;
 421
 422     invbinw = 10;
 423     nbin    = 10;
 424     snew(bin, nbin);
 425
 426     bNotLastFrame = read_first_frame(oenv, &status, opt2fn("-rerun", nfile, fnm),
 427                                      &rerun_fr, TRX_NEED_X);
 428     frame = 0;
 429
 430     if (rerun_fr.natoms - (bCavity ? nat_cavity : 0) !=
 431         mdatoms->nr - (a_tp1 - a_tp0))
 432     {
 433         gmx_fatal(FARGS, "Number of atoms in trajectory (%d)%s "
 434                   "is not equal the number in the run input file (%d) "
 435                   "minus the number of atoms to insert (%d)\n",
 436                   rerun_fr.natoms, bCavity ? " minus one" : "",
 437                   mdatoms->nr, a_tp1-a_tp0);
 438     }
 439
 440     refvolshift = log(det(rerun_fr.box));
 441
 442 #ifdef GMX_X86_SSE2
 443     /* Make sure we don't detect SSE overflow generated before this point */
 444     gmx_mm_check_and_reset_overflow();
 445 #endif
 446
 447     while (bNotLastFrame)
 448     {
 449         lambda = rerun_fr.lambda;
 450         t      = rerun_fr.time;
 451
 452         sum_embU = 0;
 453         for (e = 0; e < nener; e++)
 454         {
 455             sum_UgembU[e] = 0;
 456         }
 457
 458         /* Copy the coordinates from the input trajectory */
 459         for (i = 0; i < rerun_fr.natoms; i++)
 460         {
 461             copy_rvec(rerun_fr.x[i], state->x[i]);
 462         }
 463         copy_mat(rerun_fr.box, state->box);
 464
 465         V    = det(state->box);
 466         logV = log(V);
 467
 468         bStateChanged = TRUE;
 469         bNS           = TRUE;
 470         for (step = 0; step < nsteps; step++)
 471         {
 472             /* In parallel all nodes generate all random configurations.
 473              * In that way the result is identical to a single cpu tpi run.
 474              */
 475             if (!bCavity)
 476             {
 477                 /* Random insertion in the whole volume */
 478                 bNS = (step % inputrec->nstlist == 0);
 479                 if (bNS)
 480                 {
 481                     /* Generate a random position in the box */
 482                     x_init[XX] = gmx_rng_uniform_real(tpi_rand)*state->box[XX][XX];
 483                     x_init[YY] = gmx_rng_uniform_real(tpi_rand)*state->box[YY][YY];
 484                     x_init[ZZ] = gmx_rng_uniform_real(tpi_rand)*state->box[ZZ][ZZ];
 485                 }
 486                 if (inputrec->nstlist == 1)
 487                 {
 488                     copy_rvec(x_init, x_tp);
 489                 }
 490                 else
 491                 {
 492                     /* Generate coordinates within |dx|=drmax of x_init */
 493                     do
 494                     {
 495                         dx[XX] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
 496                         dx[YY] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
 497                         dx[ZZ] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
 498                     }
 499                     while (norm2(dx) > drmax*drmax);
 500                     rvec_add(x_init, dx, x_tp);
 501                 }
 502             }
 503             else
 504             {
 505                 /* Random insertion around a cavity location
 506                  * given by the last coordinate of the trajectory.
 507                  */
 508                 if (step == 0)
 509                 {
 510                     if (nat_cavity == 1)
 511                     {
 512                         /* Copy the location of the cavity */
 513                         copy_rvec(rerun_fr.x[rerun_fr.natoms-1], x_init);
 514                     }
 515                     else
 516                     {
 517                         /* Determine the center of mass of the last molecule */
 518                         clear_rvec(x_init);
 519                         mass_tot = 0;
 520                         for (i = 0; i < nat_cavity; i++)
 521                         {
 522                             for (d = 0; d < DIM; d++)
 523                             {
 524                                 x_init[d] +=
 525                                     mass_cavity[i]*rerun_fr.x[rerun_fr.natoms-nat_cavity+i][d];
 526                             }
 527                             mass_tot += mass_cavity[i];
 528                         }
 529                         for (d = 0; d < DIM; d++)
 530                         {
 531                             x_init[d] /= mass_tot;
 532                         }
 533                     }
 534                 }
 535                 /* Generate coordinates within |dx|=drmax of x_init */
 536                 do
 537                 {
 538                     dx[XX] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
 539                     dx[YY] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
 540                     dx[ZZ] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
 541                 }
 542                 while (norm2(dx) > drmax*drmax);
 543                 rvec_add(x_init, dx, x_tp);
 544             }
 545
 546             if (a_tp1 - a_tp0 == 1)
 547             {
 548                 /* Insert a single atom, just copy the insertion location */
 549                 copy_rvec(x_tp, state->x[a_tp0]);
 550             }
 551             else
 552             {
 553                 /* Copy the coordinates from the top file */
 554                 for (i = a_tp0; i < a_tp1; i++)
 555                 {
 556                     copy_rvec(x_mol[i-a_tp0], state->x[i]);
 557                 }
 558                 /* Rotate the molecule randomly */
 559                 rotate_conf(a_tp1-a_tp0, state->x+a_tp0, NULL,
 560                             2*M_PI*gmx_rng_uniform_real(tpi_rand),
 561                             2*M_PI*gmx_rng_uniform_real(tpi_rand),
 562                             2*M_PI*gmx_rng_uniform_real(tpi_rand));
 563                 /* Shift to the insertion location */
 564                 for (i = a_tp0; i < a_tp1; i++)
 565                 {
 566                     rvec_inc(state->x[i], x_tp);
 567                 }
 568             }
 569
 570             /* Check if this insertion belongs to this node */
 571             bOurStep = TRUE;
 572             if (PAR(cr))
 573             {
 574                 switch (inputrec->eI)
 575                 {
 576                     case eiTPI:
 577                         bOurStep = ((step / inputrec->nstlist) % nnodes == cr->nodeid);
 578                         break;
 579                     case eiTPIC:
 580                         bOurStep = (step % nnodes == cr->nodeid);
 581                         break;
 582                     default:
 583                         gmx_fatal(FARGS, "Unknown integrator %s", ei_names[inputrec->eI]);
 584                 }
 585             }
 586             if (bOurStep)
 587             {
 588                 /* Clear some matrix variables  */
 589                 clear_mat(force_vir);
 590                 clear_mat(shake_vir);
 591                 clear_mat(vir);
 592                 clear_mat(pres);
 593
 594                 /* Set the charge group center of mass of the test particle */
 595                 copy_rvec(x_init, fr->cg_cm[top->cgs.nr-1]);
 596
 597                 /* Calc energy (no forces) on new positions.
 598                  * Since we only need the intermolecular energy
 599                  * and the RF exclusion terms of the inserted molecule occur
 600                  * within a single charge group we can pass NULL for the graph.
 601                  * This also avoids shifts that would move charge groups
 602                  * out of the box.
 603                  *
 604                  * Some checks above ensure than we can not have
 605                  * twin-range interactions together with nstlist > 1,
 606                  * therefore we do not need to remember the LR energies.
 607                  */
 608                 /* Make do_force do a single node force calculation */
 609                 cr->nnodes = 1;
 610                 do_force(fplog, cr, inputrec,
 611                          step, nrnb, wcycle, top, &top_global->groups,
 612                          state->box, state->x, &state->hist,
 613                          f, force_vir, mdatoms, enerd, fcd,
 614                          state->lambda,
 615                          NULL, fr, NULL, mu_tot, t, NULL, NULL, FALSE,
 616                          GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
 617                          (bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS | GMX_FORCE_DO_LR : 0) |
 618                          (bStateChanged ? GMX_FORCE_STATECHANGED : 0));
 619                 cr->nnodes    = nnodes;
 620                 bStateChanged = FALSE;
 621                 bNS           = FALSE;
 622
 623                 /* Calculate long range corrections to pressure and energy */
 624                 calc_dispcorr(fplog, inputrec, fr, step, top_global->natoms, state->box,
 625                               lambda, pres, vir, &prescorr, &enercorr, &dvdlcorr);
 626                 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
 627                 enerd->term[F_DISPCORR]  = enercorr;
 628                 enerd->term[F_EPOT]     += enercorr;
 629                 enerd->term[F_PRES]     += prescorr;
 630                 enerd->term[F_DVDL_VDW] += dvdlcorr;
 631
 632                 epot               = enerd->term[F_EPOT];
 633                 bEnergyOutOfBounds = FALSE;
 634 #ifdef GMX_X86_SSE2
 635                 /* With SSE the energy can overflow, check for this */
 636                 if (gmx_mm_check_and_reset_overflow())
 637                 {
 638                     if (debug)
 639                     {
 640                         fprintf(debug, "Found an SSE overflow, assuming the energy is out of bounds\n");
 641                     }
 642                     bEnergyOutOfBounds = TRUE;
 643                 }
 644 #endif
 645                 /* If the compiler doesn't optimize this check away
 646                  * we catch the NAN energies.
 647                  * The epot>GMX_REAL_MAX check catches inf values,
 648                  * which should nicely result in embU=0 through the exp below,
 649                  * but it does not hurt to check anyhow.
 650                  */
 651                 /* Non-bonded Interaction usually diverge at r=0.
 652                  * With tabulated interaction functions the first few entries
 653                  * should be capped in a consistent fashion between
 654                  * repulsion, dispersion and Coulomb to avoid accidental
 655                  * negative values in the total energy.
 656                  * The table generation code in tables.c does this.
 657                  * With user tbales the user should take care of this.
 658                  */
 659                 if (epot != epot || epot > GMX_REAL_MAX)
 660                 {
 661                     bEnergyOutOfBounds = TRUE;
 662                 }
 663                 if (bEnergyOutOfBounds)
 664                 {
 665                     if (debug)
 666                     {
 667                         fprintf(debug, "\n  time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n", t, step, epot);
 668                     }
 669                     embU = 0;
 670                 }
 671                 else
 672                 {
 673                     embU      = exp(-beta*epot);
 674                     sum_embU += embU;
 675                     /* Determine the weighted energy contributions of each energy group */
 676                     e                = 0;
 677                     sum_UgembU[e++] += epot*embU;
 678                     if (fr->bBHAM)
 679                     {
 680                         for (i = 0; i < ngid; i++)
 681                         {
 682                             sum_UgembU[e++] +=
 683                                 (enerd->grpp.ener[egBHAMSR][GID(i, gid_tp, ngid)] +
 684                                  enerd->grpp.ener[egBHAMLR][GID(i, gid_tp, ngid)])*embU;
 685                         }
 686                     }
 687                     else
 688                     {
 689                         for (i = 0; i < ngid; i++)
 690                         {
 691                             sum_UgembU[e++] +=
 692                                 (enerd->grpp.ener[egLJSR][GID(i, gid_tp, ngid)] +
 693                                  enerd->grpp.ener[egLJLR][GID(i, gid_tp, ngid)])*embU;
 694                         }
 695                     }
 696                     if (bDispCorr)
 697                     {
 698                         sum_UgembU[e++] += enerd->term[F_DISPCORR]*embU;
 699                     }
 700                     if (bCharge)
 701                     {
 702                         for (i = 0; i < ngid; i++)
 703                         {
 704                             sum_UgembU[e++] +=
 705                                 (enerd->grpp.ener[egCOULSR][GID(i, gid_tp, ngid)] +
 706                                  enerd->grpp.ener[egCOULLR][GID(i, gid_tp, ngid)])*embU;
 707                         }
 708                         if (bRFExcl)
 709                         {
 710                             sum_UgembU[e++] += enerd->term[F_RF_EXCL]*embU;
 711                         }
 712                         if (EEL_FULL(fr->eeltype))
 713                         {
 714                             sum_UgembU[e++] += enerd->term[F_COUL_RECIP]*embU;
 715                         }
 716                     }
 717                 }
 718
 719                 if (embU == 0 || beta*epot > bU_bin_limit)
 720                 {
 721                     bin[0]++;
 722                 }
 723                 else
 724                 {
 725                     i = (int)((bU_logV_bin_limit
 726                                - (beta*epot - logV + refvolshift))*invbinw
 727                               + 0.5);
 728                     if (i < 0)
 729                     {
 730                         i = 0;
 731                     }
 732                     if (i >= nbin)
 733                     {
 734                         realloc_bins(&bin, &nbin, i+10);
 735                     }
 736                     bin[i]++;
 737                 }
 738
 739                 if (debug)
 740                 {
 741                     fprintf(debug, "TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
 742                             step, epot, x_tp[XX], x_tp[YY], x_tp[ZZ]);
 743                 }
 744
 745                 if (dump_pdb && epot <= dump_ener)
 746                 {
 747                     sprintf(str, "t%g_step%d.pdb", t, step);
 748                     sprintf(str2, "t: %f step %d ener: %f", t, step, epot);
 749                     write_sto_conf_mtop(str, str2, top_global, state->x, state->v,
 750                                         inputrec->ePBC, state->box);
 751                 }
 752             }
 753         }
 754
 755         if (PAR(cr))
 756         {
 757             /* When running in parallel sum the energies over the processes */
 758             gmx_sumd(1,    &sum_embU, cr);
 759             gmx_sumd(nener, sum_UgembU, cr);
 760         }
 761
 762         frame++;
 763         V_all     += V;
 764         VembU_all += V*sum_embU/nsteps;
 765
 766         if (fp_tpi)
 767         {
 768             if (bVerbose || frame%10 == 0 || frame < 10)
 769             {
 770                 fprintf(stderr, "mu %10.3e <mu> %10.3e\n",
 771                         -log(sum_embU/nsteps)/beta, -log(VembU_all/V_all)/beta);
 772             }
 773
 774             fprintf(fp_tpi, "%10.3f %12.5e %12.5e %12.5e %12.5e",
 775                     t,
 776                     VembU_all == 0 ? 20/beta : -log(VembU_all/V_all)/beta,
 777                     sum_embU == 0  ? 20/beta : -log(sum_embU/nsteps)/beta,
 778                     sum_embU/nsteps, V);
 779             for (e = 0; e < nener; e++)
 780             {
 781                 fprintf(fp_tpi, " %12.5e", sum_UgembU[e]/nsteps);
 782             }
 783             fprintf(fp_tpi, "\n");
 784             fflush(fp_tpi);
 785         }
 786
 787         bNotLastFrame = read_next_frame(oenv, status, &rerun_fr);
 788     } /* End of the loop  */
 789     runtime_end(runtime);
 790
 791     close_trj(status);
 792
 793     if (fp_tpi != NULL)
 794     {
 795         gmx_fio_fclose(fp_tpi);
 796     }
 797
 798     if (fplog != NULL)
 799     {
 800         fprintf(fplog, "\n");
 801         fprintf(fplog, "  <V>  = %12.5e nm^3\n", V_all/frame);
 802         fprintf(fplog, "  <mu> = %12.5e kJ/mol\n", -log(VembU_all/V_all)/beta);
 803     }
 804
 805     /* Write the Boltzmann factor histogram */
 806     if (PAR(cr))
 807     {
 808         /* When running in parallel sum the bins over the processes */
 809         i = nbin;
 810         global_max(cr, &i);
 811         realloc_bins(&bin, &nbin, i);
 812         gmx_sumd(nbin, bin, cr);
 813     }
 814     if (MASTER(cr))
 815     {
 816         fp_tpi = xvgropen(opt2fn("-tpid", nfile, fnm),
 817                           "TPI energy distribution",
 818                           "\\betaU - log(V/<V>)", "count", oenv);
 819         sprintf(str, "number \\betaU > %g: %9.3e", bU_bin_limit, bin[0]);
 820         xvgr_subtitle(fp_tpi, str, oenv);
 821         xvgr_legend(fp_tpi, 2, (const char **)tpid_leg, oenv);
 822         for (i = nbin-1; i > 0; i--)
 823         {
 824             bUlogV = -i/invbinw + bU_logV_bin_limit - refvolshift + log(V_all/frame);
 825             fprintf(fp_tpi, "%6.2f %10d %12.5e\n",
 826                     bUlogV,
 827                     (int)(bin[i]+0.5),
 828                     bin[i]*exp(-bUlogV)*V_all/VembU_all);
 829         }
 830         gmx_fio_fclose(fp_tpi);
 831     }
 832     sfree(bin);
 833
 834     sfree(sum_UgembU);
 835
 836     runtime->nsteps_done = frame*inputrec->nsteps;
 837
 838     return 0;
 839 }