src/gromacs/mdrun/tpi.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  *
  39  * \brief This file defines the integrator for test particle insertion
  40  *
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_mdrun
  43  */
  44 #include "gmxpre.h"
  45
  46 #include <cmath>
  47 #include <cstdlib>
  48 #include <cstring>
  49 #include <ctime>
  50
  51 #include <algorithm>
  52
  53 #include <cfenv>
  54
  55 #include "gromacs/commandline/filenm.h"
  56 #include "gromacs/domdec/dlbtiming.h"
  57 #include "gromacs/domdec/domdec.h"
  58 #include "gromacs/ewald/pme.h"
  59 #include "gromacs/fileio/confio.h"
  60 #include "gromacs/fileio/trxio.h"
  61 #include "gromacs/fileio/xvgr.h"
  62 #include "gromacs/gmxlib/chargegroup.h"
  63 #include "gromacs/gmxlib/conformation_utilities.h"
  64 #include "gromacs/gmxlib/network.h"
  65 #include "gromacs/gmxlib/nrnb.h"
  66 #include "gromacs/math/units.h"
  67 #include "gromacs/math/vec.h"
  68 #include "gromacs/mdlib/constr.h"
  69 #include "gromacs/mdlib/dispersioncorrection.h"
  70 #include "gromacs/mdlib/energyoutput.h"
  71 #include "gromacs/mdlib/force.h"
  72 #include "gromacs/mdlib/force_flags.h"
  73 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  74 #include "gromacs/mdlib/mdatoms.h"
  75 #include "gromacs/mdlib/tgroup.h"
  76 #include "gromacs/mdlib/update.h"
  77 #include "gromacs/mdlib/vsite.h"
  78 #include "gromacs/mdrunutility/printtime.h"
  79 #include "gromacs/mdtypes/commrec.h"
  80 #include "gromacs/mdtypes/forcerec.h"
  81 #include "gromacs/mdtypes/group.h"
  82 #include "gromacs/mdtypes/inputrec.h"
  83 #include "gromacs/mdtypes/md_enums.h"
  84 #include "gromacs/mdtypes/mdrunoptions.h"
  85 #include "gromacs/mdtypes/state.h"
  86 #include "gromacs/nbnxm/nbnxm.h"
  87 #include "gromacs/pbcutil/pbc.h"
  88 #include "gromacs/random/threefry.h"
  89 #include "gromacs/random/uniformrealdistribution.h"
  90 #include "gromacs/timing/wallcycle.h"
  91 #include "gromacs/timing/walltime_accounting.h"
  92 #include "gromacs/topology/mtop_util.h"
  93 #include "gromacs/trajectory/trajectoryframe.h"
  94 #include "gromacs/utility/cstringutil.h"
  95 #include "gromacs/utility/fatalerror.h"
  96 #include "gromacs/utility/gmxassert.h"
  97 #include "gromacs/utility/logger.h"
  98 #include "gromacs/utility/smalloc.h"
  99
 100 #include "legacysimulator.h"
 101
 102 //! Global max algorithm
 103 static void global_max(t_commrec *cr, int *n)
 104 {
 105     int *sum, i;
 106
 107     snew(sum, cr->nnodes);
 108     sum[cr->nodeid] = *n;
 109     gmx_sumi(cr->nnodes, sum, cr);
 110     for (i = 0; i < cr->nnodes; i++)
 111     {
 112         *n = std::max(*n, sum[i]);
 113     }
 114
 115     sfree(sum);
 116 }
 117
 118 //! Reallocate arrays.
 119 static void realloc_bins(double **bin, int *nbin, int nbin_new)
 120 {
 121     int i;
 122
 123     if (nbin_new != *nbin)
 124     {
 125         srenew(*bin, nbin_new);
 126         for (i = *nbin; i < nbin_new; i++)
 127         {
 128             (*bin)[i] = 0;
 129         }
 130         *nbin = nbin_new;
 131     }
 132 }
 133
 134 //! Computes and returns the RF exclusion energy for the last molecule starting at \p beginAtom
 135 static real
 136 reactionFieldExclusionCorrection(gmx::ArrayRef<const gmx::RVec> x,
 137                                  const t_mdatoms               &mdatoms,
 138                                  const interaction_const_t     &ic,
 139                                  const int                      beginAtom)
 140 {
 141     real energy = 0;
 142
 143     for (int i = beginAtom; i < mdatoms.homenr; i++)
 144     {
 145         const real qi = mdatoms.chargeA[i];
 146         energy -= 0.5*qi*qi*ic.c_rf;
 147
 148         for (int j = i + 1; j < mdatoms.homenr; j++)
 149         {
 150             const real qj  = mdatoms.chargeA[j];
 151             const real rsq = distance2(x[i], x[j]);
 152             energy += qi*qj*(ic.k_rf*rsq - ic.c_rf);
 153         }
 154     }
 155
 156     return ic.epsfac*energy;
 157 }
 158
 159 namespace gmx
 160 {
 161
 162 // TODO: Convert to use the nbnxm kernels by putting the system and the teset molecule on two separate search grids
 163 void
 164 LegacySimulator::do_tpi()
 165 {
 166     GMX_RELEASE_ASSERT(gmx_omp_nthreads_get(emntDefault) == 1, "TPI does not support OpenMP");
 167
 168     gmx_localtop_t          top;
 169     PaddedVector<gmx::RVec> f {};
 170     real                    lambda, t, temp, beta, drmax, epot;
 171     double                  embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
 172     t_trxstatus            *status;
 173     t_trxframe              rerun_fr;
 174     gmx_bool                bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS;
 175     tensor                  force_vir, shake_vir, vir, pres;
 176     int                     a_tp0, a_tp1, ngid, gid_tp, nener, e;
 177     rvec                   *x_mol;
 178     rvec                    mu_tot, x_init, dx;
 179     int                     nnodes, frame;
 180     int64_t                 frame_step_prev, frame_step;
 181     int64_t                 nsteps, stepblocksize = 0, step;
 182     int64_t                 seed;
 183     int                     i;
 184     FILE                   *fp_tpi = nullptr;
 185     char                   *ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN];
 186     double                  dbl, dump_ener;
 187     gmx_bool                bCavity;
 188     int                     nat_cavity  = 0, d;
 189     real                   *mass_cavity = nullptr, mass_tot;
 190     int                     nbin;
 191     double                  invbinw, *bin, refvolshift, logV, bUlogV;
 192     gmx_bool                bEnergyOutOfBounds;
 193     const char             *tpid_leg[2] = {"direct", "reweighted"};
 194     auto                    mdatoms     = mdAtoms->mdatoms();
 195
 196     GMX_UNUSED_VALUE(outputProvider);
 197
 198     GMX_LOG(mdlog.info).asParagraph().
 199         appendText("Note that it is planned to change the command gmx mdrun -tpi "
 200                    "(and -tpic) to make the functionality available in a different "
 201                    "form in a future version of GROMACS, e.g. gmx test-particle-insertion.");
 202
 203     /* Since there is no upper limit to the insertion energies,
 204      * we need to set an upper limit for the distribution output.
 205      */
 206     real bU_bin_limit      = 50;
 207     real bU_logV_bin_limit = bU_bin_limit + 10;
 208
 209     nnodes = cr->nnodes;
 210
 211     gmx_mtop_generate_local_top(*top_global, &top, inputrec->efep != efepNO);
 212
 213     SimulationGroups *groups = &top_global->groups;
 214
 215     bCavity = (inputrec->eI == eiTPIC);
 216     if (bCavity)
 217     {
 218         ptr = getenv("GMX_TPIC_MASSES");
 219         if (ptr == nullptr)
 220         {
 221             nat_cavity = 1;
 222         }
 223         else
 224         {
 225             /* Read (multiple) masses from env var GMX_TPIC_MASSES,
 226              * The center of mass of the last atoms is then used for TPIC.
 227              */
 228             nat_cavity = 0;
 229             while (sscanf(ptr, "%20lf%n", &dbl, &i) > 0)
 230             {
 231                 srenew(mass_cavity, nat_cavity+1);
 232                 mass_cavity[nat_cavity] = dbl;
 233                 fprintf(fplog, "mass[%d] = %f\n",
 234                         nat_cavity+1, mass_cavity[nat_cavity]);
 235                 nat_cavity++;
 236                 ptr += i;
 237             }
 238             if (nat_cavity == 0)
 239             {
 240                 gmx_fatal(FARGS, "Found %d masses in GMX_TPIC_MASSES", nat_cavity);
 241             }
 242         }
 243     }
 244
 245     /*
 246        init_em(fplog,TPI,inputrec,&lambda,nrnb,mu_tot,
 247        state_global->box,fr,mdatoms,top,cr,nfile,fnm,NULL,NULL);*/
 248     /* We never need full pbc for TPI */
 249     fr->ePBC = epbcXYZ;
 250     /* Determine the temperature for the Boltzmann weighting */
 251     temp = inputrec->opts.ref_t[0];
 252     if (fplog)
 253     {
 254         for (i = 1; (i < inputrec->opts.ngtc); i++)
 255         {
 256             if (inputrec->opts.ref_t[i] != temp)
 257             {
 258                 fprintf(fplog, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 259                 fprintf(stderr, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 260             }
 261         }
 262         fprintf(fplog,
 263                 "\n  The temperature for test particle insertion is %.3f K\n\n",
 264                 temp);
 265     }
 266     beta = 1.0/(BOLTZ*temp);
 267
 268     /* Number of insertions per frame */
 269     nsteps = inputrec->nsteps;
 270
 271     /* Use the same neighborlist with more insertions points
 272      * in a sphere of radius drmax around the initial point
 273      */
 274     /* This should be a proper mdp parameter */
 275     drmax = inputrec->rtpi;
 276
 277     /* An environment variable can be set to dump all configurations
 278      * to pdb with an insertion energy <= this value.
 279      */
 280     dump_pdb  = getenv("GMX_TPI_DUMP");
 281     dump_ener = 0;
 282     if (dump_pdb)
 283     {
 284         sscanf(dump_pdb, "%20lf", &dump_ener);
 285     }
 286
 287     atoms2md(top_global, inputrec, -1, nullptr, top_global->natoms, mdAtoms);
 288     update_mdatoms(mdatoms, inputrec->fepvals->init_lambda);
 289
 290     f.resizeWithPadding(top_global->natoms);
 291
 292     /* Print to log file  */
 293     walltime_accounting_start_time(walltime_accounting);
 294     wallcycle_start(wcycle, ewcRUN);
 295     print_start(fplog, cr, walltime_accounting, "Test Particle Insertion");
 296
 297     /* The last charge group is the group to be inserted */
 298     const t_atoms &atomsToInsert = top_global->moltype[top_global->molblock.back().type].atoms;
 299     a_tp0 = top_global->natoms - atomsToInsert.nr;
 300     a_tp1 = top_global->natoms;
 301     if (debug)
 302     {
 303         fprintf(debug, "TPI atoms %d-%d\n", a_tp0, a_tp1);
 304     }
 305
 306     auto x = makeArrayRef(state_global->x);
 307
 308     if (EEL_PME(fr->ic->eeltype))
 309     {
 310         gmx_pme_reinit_atoms(fr->pmedata, a_tp0, nullptr);
 311     }
 312
 313     /* With reacion-field we have distance dependent potentials
 314      * between excluded atoms, we need to add these separately
 315      * for the inserted molecule.
 316      */
 317     real rfExclusionEnergy = 0;
 318     if (EEL_RF(fr->ic->eeltype))
 319     {
 320         rfExclusionEnergy = reactionFieldExclusionCorrection(x, *mdatoms, *fr->ic, a_tp0);
 321         if (debug)
 322         {
 323             fprintf(debug, "RF exclusion correction for inserted molecule: %f kJ/mol\n", rfExclusionEnergy);
 324         }
 325     }
 326
 327     snew(x_mol, a_tp1-a_tp0);
 328
 329     bDispCorr = (inputrec->eDispCorr != edispcNO);
 330     bCharge   = FALSE;
 331     for (i = a_tp0; i < a_tp1; i++)
 332     {
 333         /* Copy the coordinates of the molecule to be insterted */
 334         copy_rvec(x[i], x_mol[i-a_tp0]);
 335         /* Check if we need to print electrostatic energies */
 336         bCharge |= (mdatoms->chargeA[i] != 0 ||
 337                     ((mdatoms->chargeB != nullptr) && mdatoms->chargeB[i] != 0));
 338     }
 339     bRFExcl = (bCharge && EEL_RF(fr->ic->eeltype));
 340
 341     // Calculate the center of geometry of the molecule to insert
 342     rvec cog = { 0, 0, 0 };
 343     for (int a = a_tp0; a < a_tp1; a++)
 344     {
 345         rvec_inc(cog, x[a]);
 346     }
 347     svmul(1.0_real/(a_tp1 - a_tp0), cog, cog);
 348     real molRadius = 0;
 349     for (int a = a_tp0; a < a_tp1; a++)
 350     {
 351         molRadius = std::max(molRadius, distance2(x[a], cog));
 352     }
 353     molRadius = std::sqrt(molRadius);
 354
 355     const real maxCutoff =  std::max(inputrec->rvdw, inputrec->rcoulomb);
 356     if (bCavity)
 357     {
 358         if (norm(cog) > 0.5*maxCutoff && fplog)
 359         {
 360             fprintf(fplog, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
 361             fprintf(stderr, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
 362         }
 363     }
 364     else
 365     {
 366         /* Center the molecule to be inserted at zero */
 367         for (i = 0; i < a_tp1-a_tp0; i++)
 368         {
 369             rvec_dec(x_mol[i], cog);
 370         }
 371     }
 372
 373     if (fplog)
 374     {
 375         fprintf(fplog, "\nWill insert %d atoms %s partial charges\n",
 376                 a_tp1-a_tp0, bCharge ? "with" : "without");
 377
 378         fprintf(fplog, "\nWill insert %" PRId64 " times in each frame of %s\n",
 379                 nsteps, opt2fn("-rerun", nfile, fnm));
 380     }
 381
 382     if (!bCavity)
 383     {
 384         if (inputrec->nstlist > 1)
 385         {
 386
 387             /* With the same pair list we insert in a sphere of radius rtpi  in different orientations */
 388             if (drmax == 0 && a_tp1-a_tp0 == 1)
 389             {
 390                 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);
 391             }
 392             if (fplog)
 393             {
 394                 fprintf(fplog, "Will use the same neighborlist for %d insertions in a sphere of radius %f\n", inputrec->nstlist, drmax);
 395             }
 396         }
 397     }
 398     else
 399     {
 400         if (fplog)
 401         {
 402             fprintf(fplog, "Will insert randomly in a sphere of radius %f around the center of the cavity\n", drmax);
 403         }
 404     }
 405
 406     /* With the same pair list we insert in a sphere of radius rtpi
 407      * in different orientations. We generate the pairlist with all
 408      * inserted atoms located in the center of the sphere, so we need
 409      * a buffer of size of the sphere and molecule radius.
 410      */
 411     inputrec->rlist = maxCutoff + 2*inputrec->rtpi + 2*molRadius;
 412     fr->rlist       = inputrec->rlist;
 413     fr->nbv->changePairlistRadii(inputrec->rlist, inputrec->rlist);
 414
 415     ngid   = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
 416     gid_tp = GET_CGINFO_GID(fr->cginfo[a_tp0]);
 417     for (int a = a_tp0 + 1; a < a_tp1; a++)
 418     {
 419         if (GET_CGINFO_GID(fr->cginfo[a]) != gid_tp)
 420         {
 421             fprintf(fplog,
 422                     "NOTE: Atoms in the molecule to insert belong to different energy groups.\n"
 423                     "      Only contributions to the group of the first atom will be reported.\n");
 424             break;
 425         }
 426     }
 427     nener  = 1 + ngid;
 428     if (bDispCorr)
 429     {
 430         nener += 1;
 431     }
 432     if (bCharge)
 433     {
 434         nener += ngid;
 435         if (bRFExcl)
 436         {
 437             nener += 1;
 438         }
 439         if (EEL_FULL(fr->ic->eeltype))
 440         {
 441             nener += 1;
 442         }
 443     }
 444     snew(sum_UgembU, nener);
 445
 446     /* Copy the random seed set by the user */
 447     seed = inputrec->ld_seed;
 448
 449     gmx::ThreeFry2x64<16>                rng(seed, gmx::RandomDomain::TestParticleInsertion); // 16 bits internal counter => 2^16 * 2 = 131072 values per stream
 450     gmx::UniformRealDistribution<real>   dist;
 451
 452     if (MASTER(cr))
 453     {
 454         fp_tpi = xvgropen(opt2fn("-tpi", nfile, fnm),
 455                           "TPI energies", "Time (ps)",
 456                           "(kJ mol\\S-1\\N) / (nm\\S3\\N)", oenv);
 457         xvgr_subtitle(fp_tpi, "f. are averages over one frame", oenv);
 458         snew(leg, 4+nener);
 459         e = 0;
 460         sprintf(str, "-kT log(<Ve\\S-\\betaU\\N>/<V>)");
 461         leg[e++] = gmx_strdup(str);
 462         sprintf(str, "f. -kT log<e\\S-\\betaU\\N>");
 463         leg[e++] = gmx_strdup(str);
 464         sprintf(str, "f. <e\\S-\\betaU\\N>");
 465         leg[e++] = gmx_strdup(str);
 466         sprintf(str, "f. V");
 467         leg[e++] = gmx_strdup(str);
 468         sprintf(str, "f. <Ue\\S-\\betaU\\N>");
 469         leg[e++] = gmx_strdup(str);
 470         for (i = 0; i < ngid; i++)
 471         {
 472             sprintf(str, "f. <U\\sVdW %s\\Ne\\S-\\betaU\\N>",
 473                     *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][i]]));
 474             leg[e++] = gmx_strdup(str);
 475         }
 476         if (bDispCorr)
 477         {
 478             sprintf(str, "f. <U\\sdisp c\\Ne\\S-\\betaU\\N>");
 479             leg[e++] = gmx_strdup(str);
 480         }
 481         if (bCharge)
 482         {
 483             for (i = 0; i < ngid; i++)
 484             {
 485                 sprintf(str, "f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
 486                         *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][i]]));
 487                 leg[e++] = gmx_strdup(str);
 488             }
 489             if (bRFExcl)
 490             {
 491                 sprintf(str, "f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
 492                 leg[e++] = gmx_strdup(str);
 493             }
 494             if (EEL_FULL(fr->ic->eeltype))
 495             {
 496                 sprintf(str, "f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
 497                 leg[e++] = gmx_strdup(str);
 498             }
 499         }
 500         xvgr_legend(fp_tpi, 4+nener, leg, oenv);
 501         for (i = 0; i < 4+nener; i++)
 502         {
 503             sfree(leg[i]);
 504         }
 505         sfree(leg);
 506     }
 507     clear_rvec(x_init);
 508     V_all     = 0;
 509     VembU_all = 0;
 510
 511     invbinw = 10;
 512     nbin    = 10;
 513     snew(bin, nbin);
 514
 515     /* Avoid frame step numbers <= -1 */
 516     frame_step_prev = -1;
 517
 518     bNotLastFrame = read_first_frame(oenv, &status, opt2fn("-rerun", nfile, fnm),
 519                                      &rerun_fr, TRX_NEED_X);
 520     frame = 0;
 521
 522     if (rerun_fr.natoms - (bCavity ? nat_cavity : 0) !=
 523         mdatoms->nr - (a_tp1 - a_tp0))
 524     {
 525         gmx_fatal(FARGS, "Number of atoms in trajectory (%d)%s "
 526                   "is not equal the number in the run input file (%d) "
 527                   "minus the number of atoms to insert (%d)\n",
 528                   rerun_fr.natoms, bCavity ? " minus one" : "",
 529                   mdatoms->nr, a_tp1-a_tp0);
 530     }
 531
 532     refvolshift = log(det(rerun_fr.box));
 533
 534     switch (inputrec->eI)
 535     {
 536         case eiTPI:
 537             stepblocksize = inputrec->nstlist;
 538             break;
 539         case eiTPIC:
 540             stepblocksize = 1;
 541             break;
 542         default:
 543             gmx_fatal(FARGS, "Unknown integrator %s", ei_names[inputrec->eI]);
 544     }
 545
 546     while (bNotLastFrame)
 547     {
 548         frame_step      = rerun_fr.step;
 549         if (frame_step <= frame_step_prev)
 550         {
 551             /* We don't have step number in the trajectory file,
 552              * or we have constant or decreasing step numbers.
 553              * Ensure we have increasing step numbers, since we use
 554              * the step numbers as a counter for random numbers.
 555              */
 556             frame_step  = frame_step_prev + 1;
 557         }
 558         frame_step_prev = frame_step;
 559
 560         lambda = rerun_fr.lambda;
 561         t      = rerun_fr.time;
 562
 563         sum_embU = 0;
 564         for (e = 0; e < nener; e++)
 565         {
 566             sum_UgembU[e] = 0;
 567         }
 568
 569         /* Copy the coordinates from the input trajectory */
 570         auto x = makeArrayRef(state_global->x);
 571         for (i = 0; i < rerun_fr.natoms; i++)
 572         {
 573             copy_rvec(rerun_fr.x[i], x[i]);
 574         }
 575         copy_mat(rerun_fr.box, state_global->box);
 576         const matrix &box = state_global->box;
 577
 578         V    = det(box);
 579         logV = log(V);
 580
 581         bStateChanged = TRUE;
 582         bNS           = TRUE;
 583
 584         put_atoms_in_box(fr->ePBC, box, x);
 585
 586         /* Put all atoms except for the inserted ones on the grid */
 587         rvec vzero       = { 0, 0, 0 };
 588         rvec boxDiagonal = { box[XX][XX], box[YY][YY], box[ZZ][ZZ] };
 589         nbnxn_put_on_grid(fr->nbv.get(), box,
 590                           0, vzero, boxDiagonal,
 591                           nullptr, 0, a_tp0, -1,
 592                           fr->cginfo, x,
 593                           0, nullptr);
 594
 595         step = cr->nodeid*stepblocksize;
 596         while (step < nsteps)
 597         {
 598             /* Restart random engine using the frame and insertion step
 599              * as counters.
 600              * Note that we need to draw several random values per iteration,
 601              * but by using the internal subcounter functionality of ThreeFry2x64
 602              * we can draw 131072 unique 64-bit values before exhausting
 603              * the stream. This is a huge margin, and if something still goes
 604              * wrong you will get an exception when the stream is exhausted.
 605              */
 606             rng.restart(frame_step, step);
 607             dist.reset();  // erase any memory in the distribution
 608
 609             if (!bCavity)
 610             {
 611                 /* Random insertion in the whole volume */
 612                 bNS = (step % inputrec->nstlist == 0);
 613                 if (bNS)
 614                 {
 615                     /* Generate a random position in the box */
 616                     for (d = 0; d < DIM; d++)
 617                     {
 618                         x_init[d] = dist(rng)*state_global->box[d][d];
 619                     }
 620                 }
 621             }
 622             else
 623             {
 624                 /* Random insertion around a cavity location
 625                  * given by the last coordinate of the trajectory.
 626                  */
 627                 if (step == 0)
 628                 {
 629                     if (nat_cavity == 1)
 630                     {
 631                         /* Copy the location of the cavity */
 632                         copy_rvec(rerun_fr.x[rerun_fr.natoms-1], x_init);
 633                     }
 634                     else
 635                     {
 636                         /* Determine the center of mass of the last molecule */
 637                         clear_rvec(x_init);
 638                         mass_tot = 0;
 639                         for (i = 0; i < nat_cavity; i++)
 640                         {
 641                             for (d = 0; d < DIM; d++)
 642                             {
 643                                 x_init[d] +=
 644                                     mass_cavity[i]*rerun_fr.x[rerun_fr.natoms-nat_cavity+i][d];
 645                             }
 646                             mass_tot += mass_cavity[i];
 647                         }
 648                         for (d = 0; d < DIM; d++)
 649                         {
 650                             x_init[d] /= mass_tot;
 651                         }
 652                     }
 653                 }
 654             }
 655
 656             if (bNS)
 657             {
 658                 for (int a = a_tp0; a < a_tp1; a++)
 659                 {
 660                     x[a] = x_init;
 661                 }
 662
 663                 /* Put the inserted molecule on it's own search grid */
 664                 nbnxn_put_on_grid(fr->nbv.get(), box,
 665                                   1, x_init, x_init,
 666                                   nullptr, a_tp0, a_tp1, -1,
 667                                   fr->cginfo, x,
 668                                   0, nullptr);
 669
 670                 /* TODO: Avoid updating all atoms at every bNS step */
 671                 fr->nbv->setAtomProperties(*mdatoms, fr->cginfo);
 672
 673                 fr->nbv->constructPairlist(Nbnxm::InteractionLocality::Local,
 674                                            &top.excls, step, nrnb);
 675
 676                 bNS = FALSE;
 677             }
 678
 679             /* Add random displacement uniformly distributed in a sphere
 680              * of radius rtpi. We don't need to do this is we generate
 681              * a new center location every step.
 682              */
 683             rvec x_tp;
 684             if (bCavity || inputrec->nstlist > 1)
 685             {
 686                 /* Generate coordinates within |dx|=drmax of x_init */
 687                 do
 688                 {
 689                     for (d = 0; d < DIM; d++)
 690                     {
 691                         dx[d] = (2*dist(rng) - 1)*drmax;
 692                     }
 693                 }
 694                 while (norm2(dx) > drmax*drmax);
 695                 rvec_add(x_init, dx, x_tp);
 696             }
 697             else
 698             {
 699                 copy_rvec(x_init, x_tp);
 700             }
 701
 702             if (a_tp1 - a_tp0 == 1)
 703             {
 704                 /* Insert a single atom, just copy the insertion location */
 705                 copy_rvec(x_tp, x[a_tp0]);
 706             }
 707             else
 708             {
 709                 /* Copy the coordinates from the top file */
 710                 for (i = a_tp0; i < a_tp1; i++)
 711                 {
 712                     copy_rvec(x_mol[i-a_tp0], x[i]);
 713                 }
 714                 /* Rotate the molecule randomly */
 715                 real angleX = 2*M_PI*dist(rng);
 716                 real angleY = 2*M_PI*dist(rng);
 717                 real angleZ = 2*M_PI*dist(rng);
 718                 rotate_conf(a_tp1-a_tp0, state_global->x.rvec_array()+a_tp0, nullptr,
 719                             angleX, angleY, angleZ);
 720                 /* Shift to the insertion location */
 721                 for (i = a_tp0; i < a_tp1; i++)
 722                 {
 723                     rvec_inc(x[i], x_tp);
 724                 }
 725             }
 726
 727             /* Note: NonLocal refers to the inserted molecule */
 728             fr->nbv->setCoordinates(Nbnxm::AtomLocality::NonLocal, false,
 729                                     x, BufferOpsUseGpu::False, nullptr);
 730
 731             /* Clear some matrix variables  */
 732             clear_mat(force_vir);
 733             clear_mat(shake_vir);
 734             clear_mat(vir);
 735             clear_mat(pres);
 736
 737             /* Calc energy (no forces) on new positions.
 738              * Since we only need the intermolecular energy
 739              * and the RF exclusion terms of the inserted molecule occur
 740              * within a single charge group we can pass NULL for the graph.
 741              * This also avoids shifts that would move charge groups
 742              * out of the box. */
 743             /* Make do_force do a single node force calculation */
 744             cr->nnodes = 1;
 745
 746             // TPI might place a particle so close that the potential
 747             // is infinite. Since this is intended to happen, we
 748             // temporarily suppress any exceptions that the processor
 749             // might raise, then restore the old behaviour.
 750             std::fenv_t floatingPointEnvironment;
 751             std::feholdexcept(&floatingPointEnvironment);
 752             do_force(fplog, cr, ms, inputrec, nullptr, nullptr, imdSession,
 753                      pull_work,
 754                      step, nrnb, wcycle, &top,
 755                      state_global->box, state_global->x.arrayRefWithPadding(), &state_global->hist,
 756                      f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
 757                      state_global->lambda,
 758                      nullptr, fr, mdScheduleWork, nullptr, mu_tot, t, nullptr,
 759                      GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
 760                      (bStateChanged ? GMX_FORCE_STATECHANGED : 0),
 761                      DDBalanceRegionHandler(nullptr));
 762             std::feclearexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
 763             std::feupdateenv(&floatingPointEnvironment);
 764
 765             cr->nnodes    = nnodes;
 766             bStateChanged = FALSE;
 767
 768             if (fr->dispersionCorrection)
 769             {
 770                 /* Calculate long range corrections to pressure and energy */
 771                 const DispersionCorrection::Correction correction =
 772                     fr->dispersionCorrection->calculate(state_global->box, lambda);
 773                 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
 774                 enerd->term[F_DISPCORR] = correction.energy;
 775                 enerd->term[F_EPOT]    += correction.energy;
 776                 enerd->term[F_PRES]    += correction.pressure;
 777                 enerd->term[F_DVDL]    += correction.dvdl;
 778             }
 779             else
 780             {
 781                 enerd->term[F_DISPCORR]  = 0;
 782             }
 783             if (EEL_RF(fr->ic->eeltype))
 784             {
 785                 enerd->term[F_EPOT]    += rfExclusionEnergy;
 786             }
 787
 788             epot               = enerd->term[F_EPOT];
 789             bEnergyOutOfBounds = FALSE;
 790
 791             /* If the compiler doesn't optimize this check away
 792              * we catch the NAN energies.
 793              * The epot>GMX_REAL_MAX check catches inf values,
 794              * which should nicely result in embU=0 through the exp below,
 795              * but it does not hurt to check anyhow.
 796              */
 797             /* Non-bonded Interaction usually diverge at r=0.
 798              * With tabulated interaction functions the first few entries
 799              * should be capped in a consistent fashion between
 800              * repulsion, dispersion and Coulomb to avoid accidental
 801              * negative values in the total energy.
 802              * The table generation code in tables.c does this.
 803              * With user tbales the user should take care of this.
 804              */
 805             if (epot != epot || epot > GMX_REAL_MAX)
 806             {
 807                 bEnergyOutOfBounds = TRUE;
 808             }
 809             if (bEnergyOutOfBounds)
 810             {
 811                 if (debug)
 812                 {
 813                     fprintf(debug, "\n  time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n", t, static_cast<int>(step), epot);
 814                 }
 815                 embU = 0;
 816             }
 817             else
 818             {
 819                 // Exponent argument is fine in SP range, but output can be in DP range
 820                 embU      = exp(static_cast<double>(-beta*epot));
 821                 sum_embU += embU;
 822                 /* Determine the weighted energy contributions of each energy group */
 823                 e                = 0;
 824                 sum_UgembU[e++] += epot*embU;
 825                 if (fr->bBHAM)
 826                 {
 827                     for (i = 0; i < ngid; i++)
 828                     {
 829                         sum_UgembU[e++] +=
 830                             enerd->grpp.ener[egBHAMSR][GID(i, gid_tp, ngid)]*embU;
 831                     }
 832                 }
 833                 else
 834                 {
 835                     for (i = 0; i < ngid; i++)
 836                     {
 837                         sum_UgembU[e++] +=
 838                             enerd->grpp.ener[egLJSR][GID(i, gid_tp, ngid)]*embU;
 839                     }
 840                 }
 841                 if (bDispCorr)
 842                 {
 843                     sum_UgembU[e++] += enerd->term[F_DISPCORR]*embU;
 844                 }
 845                 if (bCharge)
 846                 {
 847                     for (i = 0; i < ngid; i++)
 848                     {
 849                         sum_UgembU[e++] += enerd->grpp.ener[egCOULSR][GID(i, gid_tp, ngid)] * embU;
 850                     }
 851                     if (bRFExcl)
 852                     {
 853                         sum_UgembU[e++] += rfExclusionEnergy*embU;
 854                     }
 855                     if (EEL_FULL(fr->ic->eeltype))
 856                     {
 857                         sum_UgembU[e++] += enerd->term[F_COUL_RECIP]*embU;
 858                     }
 859                 }
 860             }
 861
 862             if (embU == 0 || beta*epot > bU_bin_limit)
 863             {
 864                 bin[0]++;
 865             }
 866             else
 867             {
 868                 i = gmx::roundToInt((bU_logV_bin_limit
 869                                      - (beta*epot - logV + refvolshift))*invbinw);
 870                 if (i < 0)
 871                 {
 872                     i = 0;
 873                 }
 874                 if (i >= nbin)
 875                 {
 876                     realloc_bins(&bin, &nbin, i+10);
 877                 }
 878                 bin[i]++;
 879             }
 880
 881             if (debug)
 882             {
 883                 fprintf(debug, "TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
 884                         static_cast<int>(step), epot, x_tp[XX], x_tp[YY], x_tp[ZZ]);
 885             }
 886
 887             if (dump_pdb && epot <= dump_ener)
 888             {
 889                 sprintf(str, "t%g_step%d.pdb", t, static_cast<int>(step));
 890                 sprintf(str2, "t: %f step %d ener: %f", t, static_cast<int>(step), epot);
 891                 write_sto_conf_mtop(str, str2, top_global, state_global->x.rvec_array(), state_global->v.rvec_array(),
 892                                     inputrec->ePBC, state_global->box);
 893             }
 894
 895             step++;
 896             if ((step/stepblocksize) % cr->nnodes != cr->nodeid)
 897             {
 898                 /* Skip all steps assigned to the other MPI ranks */
 899                 step += (cr->nnodes - 1)*stepblocksize;
 900             }
 901         }
 902
 903         if (PAR(cr))
 904         {
 905             /* When running in parallel sum the energies over the processes */
 906             gmx_sumd(1,    &sum_embU, cr);
 907             gmx_sumd(nener, sum_UgembU, cr);
 908         }
 909
 910         frame++;
 911         V_all     += V;
 912         VembU_all += V*sum_embU/nsteps;
 913
 914         if (fp_tpi)
 915         {
 916             if (mdrunOptions.verbose || frame%10 == 0 || frame < 10)
 917             {
 918                 fprintf(stderr, "mu %10.3e <mu> %10.3e\n",
 919                         -log(sum_embU/nsteps)/beta, -log(VembU_all/V_all)/beta);
 920             }
 921
 922             fprintf(fp_tpi, "%10.3f %12.5e %12.5e %12.5e %12.5e",
 923                     t,
 924                     VembU_all == 0 ? 20/beta : -log(VembU_all/V_all)/beta,
 925                     sum_embU == 0  ? 20/beta : -log(sum_embU/nsteps)/beta,
 926                     sum_embU/nsteps, V);
 927             for (e = 0; e < nener; e++)
 928             {
 929                 fprintf(fp_tpi, " %12.5e", sum_UgembU[e]/nsteps);
 930             }
 931             fprintf(fp_tpi, "\n");
 932             fflush(fp_tpi);
 933         }
 934
 935         bNotLastFrame = read_next_frame(oenv, status, &rerun_fr);
 936     }   /* End of the loop  */
 937     walltime_accounting_end_time(walltime_accounting);
 938
 939     close_trx(status);
 940
 941     if (fp_tpi != nullptr)
 942     {
 943         xvgrclose(fp_tpi);
 944     }
 945
 946     if (fplog != nullptr)
 947     {
 948         fprintf(fplog, "\n");
 949         fprintf(fplog, "  <V>  = %12.5e nm^3\n", V_all/frame);
 950         const double mu = -log(VembU_all/V_all)/beta;
 951         fprintf(fplog, "  <mu> = %12.5e kJ/mol\n", mu);
 952
 953         if (!std::isfinite(mu))
 954         {
 955             fprintf(fplog, "\nThe computed chemical potential is not finite - consider increasing the number of steps and/or the number of frames to insert into.\n");
 956         }
 957     }
 958
 959     /* Write the Boltzmann factor histogram */
 960     if (PAR(cr))
 961     {
 962         /* When running in parallel sum the bins over the processes */
 963         i = nbin;
 964         global_max(cr, &i);
 965         realloc_bins(&bin, &nbin, i);
 966         gmx_sumd(nbin, bin, cr);
 967     }
 968     if (MASTER(cr))
 969     {
 970         fp_tpi = xvgropen(opt2fn("-tpid", nfile, fnm),
 971                           "TPI energy distribution",
 972                           "\\betaU - log(V/<V>)", "count", oenv);
 973         sprintf(str, "number \\betaU > %g: %9.3e", bU_bin_limit, bin[0]);
 974         xvgr_subtitle(fp_tpi, str, oenv);
 975         xvgr_legend(fp_tpi, 2, tpid_leg, oenv);
 976         for (i = nbin-1; i > 0; i--)
 977         {
 978             bUlogV = -i/invbinw + bU_logV_bin_limit - refvolshift + log(V_all/frame);
 979             fprintf(fp_tpi, "%6.2f %10d %12.5e\n",
 980                     bUlogV,
 981                     roundToInt(bin[i]),
 982                     bin[i]*exp(-bUlogV)*V_all/VembU_all);
 983         }
 984         xvgrclose(fp_tpi);
 985     }
 986     sfree(bin);
 987
 988     sfree(sum_UgembU);
 989
 990     walltime_accounting_set_nsteps_done(walltime_accounting, frame*inputrec->nsteps);
 991 }
 992
 993 } // namespace gmx