2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 * \brief This file defines the integrator for test particle insertion
41 * \author Berk Hess <hess@kth.se>
42 * \ingroup module_mdrun
55 #include "gromacs/commandline/filenm.h"
56 #include "gromacs/domdec/domdec.h"
57 #include "gromacs/ewald/pme.h"
58 #include "gromacs/fileio/confio.h"
59 #include "gromacs/fileio/trxio.h"
60 #include "gromacs/fileio/xvgr.h"
61 #include "gromacs/gmxlib/chargegroup.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/energyoutput.h"
69 #include "gromacs/mdlib/force.h"
70 #include "gromacs/mdlib/force_flags.h"
71 #include "gromacs/mdlib/mdatoms.h"
72 #include "gromacs/mdlib/mdrun.h"
73 #include "gromacs/mdlib/ns.h"
74 #include "gromacs/mdlib/sim_util.h"
75 #include "gromacs/mdlib/tgroup.h"
76 #include "gromacs/mdlib/update.h"
77 #include "gromacs/mdlib/vsite.h"
78 #include "gromacs/mdtypes/commrec.h"
79 #include "gromacs/mdtypes/forcerec.h"
80 #include "gromacs/mdtypes/group.h"
81 #include "gromacs/mdtypes/inputrec.h"
82 #include "gromacs/mdtypes/md_enums.h"
83 #include "gromacs/mdtypes/state.h"
84 #include "gromacs/pbcutil/pbc.h"
85 #include "gromacs/random/threefry.h"
86 #include "gromacs/random/uniformrealdistribution.h"
87 #include "gromacs/timing/wallcycle.h"
88 #include "gromacs/timing/walltime_accounting.h"
89 #include "gromacs/topology/mtop_util.h"
90 #include "gromacs/trajectory/trajectoryframe.h"
91 #include "gromacs/utility/cstringutil.h"
92 #include "gromacs/utility/fatalerror.h"
93 #include "gromacs/utility/gmxassert.h"
94 #include "gromacs/utility/logger.h"
95 #include "gromacs/utility/smalloc.h"
97 #include "integrator.h"
99 //! Global max algorithm
100 static void global_max(t_commrec
*cr
, int *n
)
104 snew(sum
, cr
->nnodes
);
105 sum
[cr
->nodeid
] = *n
;
106 gmx_sumi(cr
->nnodes
, sum
, cr
);
107 for (i
= 0; i
< cr
->nnodes
; i
++)
109 *n
= std::max(*n
, sum
[i
]);
115 //! Reallocate arrays.
116 static void realloc_bins(double **bin
, int *nbin
, int nbin_new
)
120 if (nbin_new
!= *nbin
)
122 srenew(*bin
, nbin_new
);
123 for (i
= *nbin
; i
< nbin_new
; i
++)
138 gmx_groups_t
*groups
;
139 gmx_enerdata_t
*enerd
;
140 PaddedVector
<gmx::RVec
> f
{};
141 real lambda
, t
, temp
, beta
, drmax
, epot
;
142 double embU
, sum_embU
, *sum_UgembU
, V
, V_all
, VembU_all
;
145 gmx_bool bDispCorr
, bCharge
, bRFExcl
, bNotLastFrame
, bStateChanged
, bNS
;
146 tensor force_vir
, shake_vir
, vir
, pres
;
147 int cg_tp
, a_tp0
, a_tp1
, ngid
, gid_tp
, nener
, e
;
149 rvec mu_tot
, x_init
, dx
, x_tp
;
151 int64_t frame_step_prev
, frame_step
;
152 int64_t nsteps
, stepblocksize
= 0, step
;
155 FILE *fp_tpi
= nullptr;
156 char *ptr
, *dump_pdb
, **leg
, str
[STRLEN
], str2
[STRLEN
];
157 double dbl
, dump_ener
;
159 int nat_cavity
= 0, d
;
160 real
*mass_cavity
= nullptr, mass_tot
;
162 double invbinw
, *bin
, refvolshift
, logV
, bUlogV
;
163 real prescorr
, enercorr
, dvdlcorr
;
164 gmx_bool bEnergyOutOfBounds
;
165 const char *tpid_leg
[2] = {"direct", "reweighted"};
166 auto mdatoms
= mdAtoms
->mdatoms();
168 GMX_UNUSED_VALUE(outputProvider
);
170 GMX_LOG(mdlog
.info
).asParagraph().
171 appendText("Note that it is planned to change the command gmx mdrun -tpi "
172 "(and -tpic) to make the functionality available in a different "
173 "form in a future version of GROMACS, e.g. gmx test-particle-insertion.");
175 /* Since there is no upper limit to the insertion energies,
176 * we need to set an upper limit for the distribution output.
178 real bU_bin_limit
= 50;
179 real bU_logV_bin_limit
= bU_bin_limit
+ 10;
181 if (inputrec
->cutoff_scheme
== ecutsVERLET
)
183 gmx_fatal(FARGS
, "TPI does not work (yet) with the Verlet cut-off scheme");
188 gmx_mtop_generate_local_top(*top_global
, &top
, inputrec
->efep
!= efepNO
);
190 groups
= &top_global
->groups
;
192 bCavity
= (inputrec
->eI
== eiTPIC
);
195 ptr
= getenv("GMX_TPIC_MASSES");
202 /* Read (multiple) masses from env var GMX_TPIC_MASSES,
203 * The center of mass of the last atoms is then used for TPIC.
206 while (sscanf(ptr
, "%20lf%n", &dbl
, &i
) > 0)
208 srenew(mass_cavity
, nat_cavity
+1);
209 mass_cavity
[nat_cavity
] = dbl
;
210 fprintf(fplog
, "mass[%d] = %f\n",
211 nat_cavity
+1, mass_cavity
[nat_cavity
]);
217 gmx_fatal(FARGS
, "Found %d masses in GMX_TPIC_MASSES", nat_cavity
);
223 init_em(fplog,TPI,inputrec,&lambda,nrnb,mu_tot,
224 state_global->box,fr,mdatoms,top,cr,nfile,fnm,NULL,NULL);*/
225 /* We never need full pbc for TPI */
227 /* Determine the temperature for the Boltzmann weighting */
228 temp
= inputrec
->opts
.ref_t
[0];
231 for (i
= 1; (i
< inputrec
->opts
.ngtc
); i
++)
233 if (inputrec
->opts
.ref_t
[i
] != temp
)
235 fprintf(fplog
, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
236 fprintf(stderr
, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
240 "\n The temperature for test particle insertion is %.3f K\n\n",
243 beta
= 1.0/(BOLTZ
*temp
);
245 /* Number of insertions per frame */
246 nsteps
= inputrec
->nsteps
;
248 /* Use the same neighborlist with more insertions points
249 * in a sphere of radius drmax around the initial point
251 /* This should be a proper mdp parameter */
252 drmax
= inputrec
->rtpi
;
254 /* An environment variable can be set to dump all configurations
255 * to pdb with an insertion energy <= this value.
257 dump_pdb
= getenv("GMX_TPI_DUMP");
261 sscanf(dump_pdb
, "%20lf", &dump_ener
);
264 atoms2md(top_global
, inputrec
, -1, nullptr, top_global
->natoms
, mdAtoms
);
265 update_mdatoms(mdatoms
, inputrec
->fepvals
->init_lambda
);
268 init_enerdata(groups
->grps
[egcENER
].nr
, inputrec
->fepvals
->n_lambda
, enerd
);
269 f
.resizeWithPadding(top_global
->natoms
);
271 /* Print to log file */
272 walltime_accounting_start_time(walltime_accounting
);
273 wallcycle_start(wcycle
, ewcRUN
);
274 print_start(fplog
, cr
, walltime_accounting
, "Test Particle Insertion");
276 /* The last charge group is the group to be inserted */
277 cg_tp
= top
.cgs
.nr
- 1;
278 a_tp0
= top
.cgs
.index
[cg_tp
];
279 a_tp1
= top
.cgs
.index
[cg_tp
+1];
282 fprintf(debug
, "TPI cg %d, atoms %d-%d\n", cg_tp
, a_tp0
, a_tp1
);
285 GMX_RELEASE_ASSERT(inputrec
->rcoulomb
<= inputrec
->rlist
&& inputrec
->rvdw
<= inputrec
->rlist
, "Twin-range interactions are not supported with TPI");
287 snew(x_mol
, a_tp1
-a_tp0
);
289 bDispCorr
= (inputrec
->eDispCorr
!= edispcNO
);
291 auto x
= makeArrayRef(state_global
->x
);
292 for (i
= a_tp0
; i
< a_tp1
; i
++)
294 /* Copy the coordinates of the molecule to be insterted */
295 copy_rvec(x
[i
], x_mol
[i
-a_tp0
]);
296 /* Check if we need to print electrostatic energies */
297 bCharge
|= (mdatoms
->chargeA
[i
] != 0 ||
298 ((mdatoms
->chargeB
!= nullptr) && mdatoms
->chargeB
[i
] != 0));
300 bRFExcl
= (bCharge
&& EEL_RF(fr
->ic
->eeltype
));
302 calc_cgcm(fplog
, cg_tp
, cg_tp
+1, &(top
.cgs
), state_global
->x
.rvec_array(), fr
->cg_cm
);
305 if (norm(fr
->cg_cm
[cg_tp
]) > 0.5*inputrec
->rlist
&& fplog
)
307 fprintf(fplog
, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
308 fprintf(stderr
, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
313 /* Center the molecule to be inserted at zero */
314 for (i
= 0; i
< a_tp1
-a_tp0
; i
++)
316 rvec_dec(x_mol
[i
], fr
->cg_cm
[cg_tp
]);
322 fprintf(fplog
, "\nWill insert %d atoms %s partial charges\n",
323 a_tp1
-a_tp0
, bCharge
? "with" : "without");
325 fprintf(fplog
, "\nWill insert %" PRId64
" times in each frame of %s\n",
326 nsteps
, opt2fn("-rerun", nfile
, fnm
));
331 if (inputrec
->nstlist
> 1)
333 if (drmax
== 0 && a_tp1
-a_tp0
== 1)
335 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
);
339 fprintf(fplog
, "Will use the same neighborlist for %d insertions in a sphere of radius %f\n", inputrec
->nstlist
, drmax
);
347 fprintf(fplog
, "Will insert randomly in a sphere of radius %f around the center of the cavity\n", drmax
);
351 ngid
= groups
->grps
[egcENER
].nr
;
352 gid_tp
= GET_CGINFO_GID(fr
->cginfo
[cg_tp
]);
365 if (EEL_FULL(fr
->ic
->eeltype
))
370 snew(sum_UgembU
, nener
);
372 /* Copy the random seed set by the user */
373 seed
= inputrec
->ld_seed
;
375 gmx::ThreeFry2x64
<16> rng(seed
, gmx::RandomDomain::TestParticleInsertion
); // 16 bits internal counter => 2^16 * 2 = 131072 values per stream
376 gmx::UniformRealDistribution
<real
> dist
;
380 fp_tpi
= xvgropen(opt2fn("-tpi", nfile
, fnm
),
381 "TPI energies", "Time (ps)",
382 "(kJ mol\\S-1\\N) / (nm\\S3\\N)", oenv
);
383 xvgr_subtitle(fp_tpi
, "f. are averages over one frame", oenv
);
386 sprintf(str
, "-kT log(<Ve\\S-\\betaU\\N>/<V>)");
387 leg
[e
++] = gmx_strdup(str
);
388 sprintf(str
, "f. -kT log<e\\S-\\betaU\\N>");
389 leg
[e
++] = gmx_strdup(str
);
390 sprintf(str
, "f. <e\\S-\\betaU\\N>");
391 leg
[e
++] = gmx_strdup(str
);
392 sprintf(str
, "f. V");
393 leg
[e
++] = gmx_strdup(str
);
394 sprintf(str
, "f. <Ue\\S-\\betaU\\N>");
395 leg
[e
++] = gmx_strdup(str
);
396 for (i
= 0; i
< ngid
; i
++)
398 sprintf(str
, "f. <U\\sVdW %s\\Ne\\S-\\betaU\\N>",
399 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
400 leg
[e
++] = gmx_strdup(str
);
404 sprintf(str
, "f. <U\\sdisp c\\Ne\\S-\\betaU\\N>");
405 leg
[e
++] = gmx_strdup(str
);
409 for (i
= 0; i
< ngid
; i
++)
411 sprintf(str
, "f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
412 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
413 leg
[e
++] = gmx_strdup(str
);
417 sprintf(str
, "f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
418 leg
[e
++] = gmx_strdup(str
);
420 if (EEL_FULL(fr
->ic
->eeltype
))
422 sprintf(str
, "f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
423 leg
[e
++] = gmx_strdup(str
);
426 xvgr_legend(fp_tpi
, 4+nener
, leg
, oenv
);
427 for (i
= 0; i
< 4+nener
; i
++)
441 /* Avoid frame step numbers <= -1 */
442 frame_step_prev
= -1;
444 bNotLastFrame
= read_first_frame(oenv
, &status
, opt2fn("-rerun", nfile
, fnm
),
445 &rerun_fr
, TRX_NEED_X
);
448 if (rerun_fr
.natoms
- (bCavity
? nat_cavity
: 0) !=
449 mdatoms
->nr
- (a_tp1
- a_tp0
))
451 gmx_fatal(FARGS
, "Number of atoms in trajectory (%d)%s "
452 "is not equal the number in the run input file (%d) "
453 "minus the number of atoms to insert (%d)\n",
454 rerun_fr
.natoms
, bCavity
? " minus one" : "",
455 mdatoms
->nr
, a_tp1
-a_tp0
);
458 refvolshift
= log(det(rerun_fr
.box
));
460 switch (inputrec
->eI
)
463 stepblocksize
= inputrec
->nstlist
;
469 gmx_fatal(FARGS
, "Unknown integrator %s", ei_names
[inputrec
->eI
]);
472 while (bNotLastFrame
)
474 frame_step
= rerun_fr
.step
;
475 if (frame_step
<= frame_step_prev
)
477 /* We don't have step number in the trajectory file,
478 * or we have constant or decreasing step numbers.
479 * Ensure we have increasing step numbers, since we use
480 * the step numbers as a counter for random numbers.
482 frame_step
= frame_step_prev
+ 1;
484 frame_step_prev
= frame_step
;
486 lambda
= rerun_fr
.lambda
;
490 for (e
= 0; e
< nener
; e
++)
495 /* Copy the coordinates from the input trajectory */
496 auto x
= makeArrayRef(state_global
->x
);
497 for (i
= 0; i
< rerun_fr
.natoms
; i
++)
499 copy_rvec(rerun_fr
.x
[i
], x
[i
]);
501 copy_mat(rerun_fr
.box
, state_global
->box
);
503 V
= det(state_global
->box
);
506 bStateChanged
= TRUE
;
509 step
= cr
->nodeid
*stepblocksize
;
510 while (step
< nsteps
)
512 /* Restart random engine using the frame and insertion step
514 * Note that we need to draw several random values per iteration,
515 * but by using the internal subcounter functionality of ThreeFry2x64
516 * we can draw 131072 unique 64-bit values before exhausting
517 * the stream. This is a huge margin, and if something still goes
518 * wrong you will get an exception when the stream is exhausted.
520 rng
.restart(frame_step
, step
);
521 dist
.reset(); // erase any memory in the distribution
525 /* Random insertion in the whole volume */
526 bNS
= (step
% inputrec
->nstlist
== 0);
529 /* Generate a random position in the box */
530 for (d
= 0; d
< DIM
; d
++)
532 x_init
[d
] = dist(rng
)*state_global
->box
[d
][d
];
536 if (inputrec
->nstlist
== 1)
538 copy_rvec(x_init
, x_tp
);
542 /* Generate coordinates within |dx|=drmax of x_init */
545 for (d
= 0; d
< DIM
; d
++)
547 dx
[d
] = (2*dist(rng
) - 1)*drmax
;
550 while (norm2(dx
) > drmax
*drmax
);
551 rvec_add(x_init
, dx
, x_tp
);
556 /* Random insertion around a cavity location
557 * given by the last coordinate of the trajectory.
563 /* Copy the location of the cavity */
564 copy_rvec(rerun_fr
.x
[rerun_fr
.natoms
-1], x_init
);
568 /* Determine the center of mass of the last molecule */
571 for (i
= 0; i
< nat_cavity
; i
++)
573 for (d
= 0; d
< DIM
; d
++)
576 mass_cavity
[i
]*rerun_fr
.x
[rerun_fr
.natoms
-nat_cavity
+i
][d
];
578 mass_tot
+= mass_cavity
[i
];
580 for (d
= 0; d
< DIM
; d
++)
582 x_init
[d
] /= mass_tot
;
586 /* Generate coordinates within |dx|=drmax of x_init */
589 for (d
= 0; d
< DIM
; d
++)
591 dx
[d
] = (2*dist(rng
) - 1)*drmax
;
594 while (norm2(dx
) > drmax
*drmax
);
595 rvec_add(x_init
, dx
, x_tp
);
598 if (a_tp1
- a_tp0
== 1)
600 /* Insert a single atom, just copy the insertion location */
601 copy_rvec(x_tp
, x
[a_tp0
]);
605 /* Copy the coordinates from the top file */
606 for (i
= a_tp0
; i
< a_tp1
; i
++)
608 copy_rvec(x_mol
[i
-a_tp0
], x
[i
]);
610 /* Rotate the molecule randomly */
611 real angleX
= 2*M_PI
*dist(rng
);
612 real angleY
= 2*M_PI
*dist(rng
);
613 real angleZ
= 2*M_PI
*dist(rng
);
614 rotate_conf(a_tp1
-a_tp0
, state_global
->x
.rvec_array()+a_tp0
, nullptr,
615 angleX
, angleY
, angleZ
);
616 /* Shift to the insertion location */
617 for (i
= a_tp0
; i
< a_tp1
; i
++)
619 rvec_inc(x
[i
], x_tp
);
623 /* Clear some matrix variables */
624 clear_mat(force_vir
);
625 clear_mat(shake_vir
);
629 /* Set the charge group center of mass of the test particle */
630 copy_rvec(x_init
, fr
->cg_cm
[top
.cgs
.nr
-1]);
632 /* Calc energy (no forces) on new positions.
633 * Since we only need the intermolecular energy
634 * and the RF exclusion terms of the inserted molecule occur
635 * within a single charge group we can pass NULL for the graph.
636 * This also avoids shifts that would move charge groups
638 /* Make do_force do a single node force calculation */
641 // TPI might place a particle so close that the potential
642 // is infinite. Since this is intended to happen, we
643 // temporarily suppress any exceptions that the processor
644 // might raise, then restore the old behaviour.
645 std::fenv_t floatingPointEnvironment
;
646 std::feholdexcept(&floatingPointEnvironment
);
647 do_force(fplog
, cr
, ms
, inputrec
, nullptr, nullptr,
648 step
, nrnb
, wcycle
, &top
, &top_global
->groups
,
649 state_global
->box
, state_global
->x
.arrayRefWithPadding(), &state_global
->hist
,
650 f
.arrayRefWithPadding(), force_vir
, mdatoms
, enerd
, fcd
,
651 state_global
->lambda
,
652 nullptr, fr
, ppForceWorkload
, nullptr, mu_tot
, t
, nullptr,
653 GMX_FORCE_NONBONDED
| GMX_FORCE_ENERGY
|
654 (bNS
? GMX_FORCE_DYNAMICBOX
| GMX_FORCE_NS
: 0) |
655 (bStateChanged
? GMX_FORCE_STATECHANGED
: 0),
656 DdOpenBalanceRegionBeforeForceComputation::no
,
657 DdCloseBalanceRegionAfterForceComputation::no
);
658 std::feclearexcept(FE_DIVBYZERO
| FE_INVALID
| FE_OVERFLOW
);
659 std::feupdateenv(&floatingPointEnvironment
);
662 bStateChanged
= FALSE
;
665 /* Calculate long range corrections to pressure and energy */
666 calc_dispcorr(inputrec
, fr
, state_global
->box
,
667 lambda
, pres
, vir
, &prescorr
, &enercorr
, &dvdlcorr
);
668 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
669 enerd
->term
[F_DISPCORR
] = enercorr
;
670 enerd
->term
[F_EPOT
] += enercorr
;
671 enerd
->term
[F_PRES
] += prescorr
;
672 enerd
->term
[F_DVDL_VDW
] += dvdlcorr
;
674 epot
= enerd
->term
[F_EPOT
];
675 bEnergyOutOfBounds
= FALSE
;
677 /* If the compiler doesn't optimize this check away
678 * we catch the NAN energies.
679 * The epot>GMX_REAL_MAX check catches inf values,
680 * which should nicely result in embU=0 through the exp below,
681 * but it does not hurt to check anyhow.
683 /* Non-bonded Interaction usually diverge at r=0.
684 * With tabulated interaction functions the first few entries
685 * should be capped in a consistent fashion between
686 * repulsion, dispersion and Coulomb to avoid accidental
687 * negative values in the total energy.
688 * The table generation code in tables.c does this.
689 * With user tbales the user should take care of this.
691 if (epot
!= epot
|| epot
> GMX_REAL_MAX
)
693 bEnergyOutOfBounds
= TRUE
;
695 if (bEnergyOutOfBounds
)
699 fprintf(debug
, "\n time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n", t
, static_cast<int>(step
), epot
);
705 // Exponent argument is fine in SP range, but output can be in DP range
706 embU
= exp(static_cast<double>(-beta
*epot
));
708 /* Determine the weighted energy contributions of each energy group */
710 sum_UgembU
[e
++] += epot
*embU
;
713 for (i
= 0; i
< ngid
; i
++)
716 enerd
->grpp
.ener
[egBHAMSR
][GID(i
, gid_tp
, ngid
)]*embU
;
721 for (i
= 0; i
< ngid
; i
++)
724 enerd
->grpp
.ener
[egLJSR
][GID(i
, gid_tp
, ngid
)]*embU
;
729 sum_UgembU
[e
++] += enerd
->term
[F_DISPCORR
]*embU
;
733 for (i
= 0; i
< ngid
; i
++)
735 sum_UgembU
[e
++] += enerd
->grpp
.ener
[egCOULSR
][GID(i
, gid_tp
, ngid
)] * embU
;
739 sum_UgembU
[e
++] += enerd
->term
[F_RF_EXCL
]*embU
;
741 if (EEL_FULL(fr
->ic
->eeltype
))
743 sum_UgembU
[e
++] += enerd
->term
[F_COUL_RECIP
]*embU
;
748 if (embU
== 0 || beta
*epot
> bU_bin_limit
)
754 i
= gmx::roundToInt((bU_logV_bin_limit
755 - (beta
*epot
- logV
+ refvolshift
))*invbinw
);
762 realloc_bins(&bin
, &nbin
, i
+10);
769 fprintf(debug
, "TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
770 static_cast<int>(step
), epot
, x_tp
[XX
], x_tp
[YY
], x_tp
[ZZ
]);
773 if (dump_pdb
&& epot
<= dump_ener
)
775 sprintf(str
, "t%g_step%d.pdb", t
, static_cast<int>(step
));
776 sprintf(str2
, "t: %f step %d ener: %f", t
, static_cast<int>(step
), epot
);
777 write_sto_conf_mtop(str
, str2
, top_global
, state_global
->x
.rvec_array(), state_global
->v
.rvec_array(),
778 inputrec
->ePBC
, state_global
->box
);
782 if ((step
/stepblocksize
) % cr
->nnodes
!= cr
->nodeid
)
784 /* Skip all steps assigned to the other MPI ranks */
785 step
+= (cr
->nnodes
- 1)*stepblocksize
;
791 /* When running in parallel sum the energies over the processes */
792 gmx_sumd(1, &sum_embU
, cr
);
793 gmx_sumd(nener
, sum_UgembU
, cr
);
798 VembU_all
+= V
*sum_embU
/nsteps
;
802 if (mdrunOptions
.verbose
|| frame
%10 == 0 || frame
< 10)
804 fprintf(stderr
, "mu %10.3e <mu> %10.3e\n",
805 -log(sum_embU
/nsteps
)/beta
, -log(VembU_all
/V_all
)/beta
);
808 fprintf(fp_tpi
, "%10.3f %12.5e %12.5e %12.5e %12.5e",
810 VembU_all
== 0 ? 20/beta
: -log(VembU_all
/V_all
)/beta
,
811 sum_embU
== 0 ? 20/beta
: -log(sum_embU
/nsteps
)/beta
,
813 for (e
= 0; e
< nener
; e
++)
815 fprintf(fp_tpi
, " %12.5e", sum_UgembU
[e
]/nsteps
);
817 fprintf(fp_tpi
, "\n");
821 bNotLastFrame
= read_next_frame(oenv
, status
, &rerun_fr
);
822 } /* End of the loop */
823 walltime_accounting_end_time(walltime_accounting
);
827 if (fp_tpi
!= nullptr)
832 if (fplog
!= nullptr)
834 fprintf(fplog
, "\n");
835 fprintf(fplog
, " <V> = %12.5e nm^3\n", V_all
/frame
);
836 const double mu
= -log(VembU_all
/V_all
)/beta
;
837 fprintf(fplog
, " <mu> = %12.5e kJ/mol\n", mu
);
839 if (!std::isfinite(mu
))
841 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");
845 /* Write the Boltzmann factor histogram */
848 /* When running in parallel sum the bins over the processes */
851 realloc_bins(&bin
, &nbin
, i
);
852 gmx_sumd(nbin
, bin
, cr
);
856 fp_tpi
= xvgropen(opt2fn("-tpid", nfile
, fnm
),
857 "TPI energy distribution",
858 "\\betaU - log(V/<V>)", "count", oenv
);
859 sprintf(str
, "number \\betaU > %g: %9.3e", bU_bin_limit
, bin
[0]);
860 xvgr_subtitle(fp_tpi
, str
, oenv
);
861 xvgr_legend(fp_tpi
, 2, tpid_leg
, oenv
);
862 for (i
= nbin
-1; i
> 0; i
--)
864 bUlogV
= -i
/invbinw
+ bU_logV_bin_limit
- refvolshift
+ log(V_all
/frame
);
865 fprintf(fp_tpi
, "%6.2f %10d %12.5e\n",
868 bin
[i
]*exp(-bUlogV
)*V_all
/VembU_all
);
876 walltime_accounting_set_nsteps_done(walltime_accounting
, frame
*inputrec
->nsteps
);