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