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 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.
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.
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.
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.
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.
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.
40 * \brief This file defines the integrator for test particle insertion
42 * \author Berk Hess <hess@kth.se>
43 * \ingroup module_mdrun
56 #include "gromacs/commandline/filenm.h"
57 #include "gromacs/domdec/dlbtiming.h"
58 #include "gromacs/domdec/domdec.h"
59 #include "gromacs/ewald/pme.h"
60 #include "gromacs/fileio/confio.h"
61 #include "gromacs/fileio/trxio.h"
62 #include "gromacs/fileio/xvgr.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/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"
102 #include "legacysimulator.h"
104 //! Global max algorithm
105 static void global_max(t_commrec
* cr
, int* n
)
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
++)
114 *n
= std::max(*n
, sum
[i
]);
120 //! Reallocate arrays.
121 static void realloc_bins(double** bin
, int* nbin
, int nbin_new
)
125 if (nbin_new
!= *nbin
)
127 srenew(*bin
, nbin_new
);
128 for (i
= *nbin
; i
< nbin_new
; i
++)
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
,
144 for (int i
= beginAtom
; i
< mdatoms
.homenr
; i
++)
146 const real qi
= mdatoms
.chargeA
[i
];
147 energy
-= 0.5 * qi
* qi
* ic
.c_rf
;
149 for (int j
= i
+ 1; j
< mdatoms
.homenr
; j
++)
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
);
157 return ic
.epsfac
* energy
;
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()
166 GMX_RELEASE_ASSERT(gmx_omp_nthreads_get(emntDefault
) == 1, "TPI does not support OpenMP");
169 PaddedHostVector
<gmx::RVec
> f
{};
170 real lambda
, t
, temp
, beta
, drmax
, epot
;
171 double embU
, sum_embU
, *sum_UgembU
, V
, V_all
, VembU_all
;
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
;
178 rvec mu_tot
, x_init
, dx
;
180 int64_t frame_step_prev
, frame_step
;
181 int64_t nsteps
, stepblocksize
= 0, step
;
184 FILE* fp_tpi
= nullptr;
185 char * ptr
, *dump_pdb
, **leg
, str
[STRLEN
], str2
[STRLEN
];
186 double dbl
, dump_ener
;
188 int nat_cavity
= 0, d
;
189 real
* mass_cavity
= nullptr, mass_tot
;
191 double invbinw
, *bin
, refvolshift
, logV
, bUlogV
;
192 gmx_bool bEnergyOutOfBounds
;
193 const char* tpid_leg
[2] = { "direct", "reweighted" };
194 auto mdatoms
= mdAtoms
->mdatoms();
196 GMX_UNUSED_VALUE(outputProvider
);
198 if (EVDW_PME(inputrec
->vdwtype
))
200 gmx_fatal(FARGS
, "Test particle insertion not implemented with LJ-PME");
202 if (haveEwaldSurfaceContribution(*inputrec
))
205 "TPI with PME currently only works in a 3D geometry with tin-foil "
206 "boundary conditions");
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.");
216 /* Since there is no upper limit to the insertion energies,
217 * we need to set an upper limit for the distribution output.
219 real bU_bin_limit
= 50;
220 real bU_logV_bin_limit
= bU_bin_limit
+ 10;
224 gmx_mtop_generate_local_top(*top_global
, &top
, inputrec
->efep
!= efepNO
);
226 SimulationGroups
* groups
= &top_global
->groups
;
228 bCavity
= (inputrec
->eI
== eiTPIC
);
231 ptr
= getenv("GMX_TPIC_MASSES");
238 /* Read (multiple) masses from env var GMX_TPIC_MASSES,
239 * The center of mass of the last atoms is then used for TPIC.
242 while (sscanf(ptr
, "%20lf%n", &dbl
, &i
) > 0)
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
]);
252 gmx_fatal(FARGS
, "Found %d masses in GMX_TPIC_MASSES", nat_cavity
);
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];
266 for (i
= 1; (i
< inputrec
->opts
.ngtc
); i
++)
268 if (inputrec
->opts
.ref_t
[i
] != temp
)
271 "\nWARNING: The temperatures of the different temperature coupling groups "
272 "are not identical\n\n");
274 "\nWARNING: The temperatures of the different temperature coupling groups "
275 "are not identical\n\n");
278 fprintf(fplog
, "\n The temperature for test particle insertion is %.3f K\n\n", temp
);
280 beta
= 1.0 / (BOLTZ
* temp
);
282 /* Number of insertions per frame */
283 nsteps
= inputrec
->nsteps
;
285 /* Use the same neighborlist with more insertions points
286 * in a sphere of radius drmax around the initial point
288 /* This should be a proper mdp parameter */
289 drmax
= inputrec
->rtpi
;
291 /* An environment variable can be set to dump all configurations
292 * to pdb with an insertion energy <= this value.
294 dump_pdb
= getenv("GMX_TPI_DUMP");
298 sscanf(dump_pdb
, "%20lf", &dump_ener
);
301 atoms2md(top_global
, inputrec
, -1, nullptr, top_global
->natoms
, mdAtoms
);
302 update_mdatoms(mdatoms
, inputrec
->fepvals
->init_lambda
);
304 f
.resizeWithPadding(top_global
->natoms
);
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");
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
;
317 fprintf(debug
, "TPI atoms %d-%d\n", a_tp0
, a_tp1
);
320 auto x
= makeArrayRef(state_global
->x
);
322 if (EEL_PME(fr
->ic
->eeltype
))
324 gmx_pme_reinit_atoms(fr
->pmedata
, a_tp0
, nullptr);
327 /* With reacion-field we have distance dependent potentials
328 * between excluded atoms, we need to add these separately
329 * for the inserted molecule.
331 real rfExclusionEnergy
= 0;
332 if (EEL_RF(fr
->ic
->eeltype
))
334 rfExclusionEnergy
= reactionFieldExclusionCorrection(x
, *mdatoms
, *fr
->ic
, a_tp0
);
337 fprintf(debug
, "RF exclusion correction for inserted molecule: %f kJ/mol\n", rfExclusionEnergy
);
341 snew(x_mol
, a_tp1
- a_tp0
);
343 bDispCorr
= (inputrec
->eDispCorr
!= edispcNO
);
345 for (i
= a_tp0
; i
< a_tp1
; i
++)
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));
353 bRFExcl
= (bCharge
&& EEL_RF(fr
->ic
->eeltype
));
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
++)
361 svmul(1.0_real
/ (a_tp1
- a_tp0
), cog
, cog
);
363 for (int a
= a_tp0
; a
< a_tp1
; a
++)
365 molRadius
= std::max(molRadius
, distance2(x
[a
], cog
));
367 molRadius
= std::sqrt(molRadius
);
369 const real maxCutoff
= std::max(inputrec
->rvdw
, inputrec
->rcoulomb
);
372 if (norm(cog
) > 0.5 * maxCutoff
&& fplog
)
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");
380 /* Center the molecule to be inserted at zero */
381 for (i
= 0; i
< a_tp1
- a_tp0
; i
++)
383 rvec_dec(x_mol
[i
], cog
);
389 fprintf(fplog
, "\nWill insert %d atoms %s partial charges\n", a_tp1
- a_tp0
,
390 bCharge
? "with" : "without");
392 fprintf(fplog
, "\nWill insert %" PRId64
" times in each frame of %s\n", nsteps
,
393 opt2fn("-rerun", nfile
, fnm
));
398 if (inputrec
->nstlist
> 1)
401 /* With the same pair list we insert in a sphere of radius rtpi in different orientations */
402 if (drmax
== 0 && a_tp1
- a_tp0
== 1)
405 "Re-using the neighborlist %d times for insertions of a single atom in a "
406 "sphere of radius %f does not make sense",
407 inputrec
->nstlist
, drmax
);
412 "Will use the same neighborlist for %d insertions in a sphere of radius "
414 inputrec
->nstlist
, drmax
);
423 "Will insert randomly in a sphere of radius %f around the center of the "
429 /* With the same pair list we insert in a sphere of radius rtpi
430 * in different orientations. We generate the pairlist with all
431 * inserted atoms located in the center of the sphere, so we need
432 * a buffer of size of the sphere and molecule radius.
434 inputrec
->rlist
= maxCutoff
+ 2 * inputrec
->rtpi
+ 2 * molRadius
;
435 fr
->rlist
= inputrec
->rlist
;
436 fr
->nbv
->changePairlistRadii(inputrec
->rlist
, inputrec
->rlist
);
438 ngid
= groups
->groups
[SimulationAtomGroupType::EnergyOutput
].size();
439 gid_tp
= GET_CGINFO_GID(fr
->cginfo
[a_tp0
]);
440 for (int a
= a_tp0
+ 1; a
< a_tp1
; a
++)
442 if (GET_CGINFO_GID(fr
->cginfo
[a
]) != gid_tp
)
445 "NOTE: Atoms in the molecule to insert belong to different energy groups.\n"
446 " Only contributions to the group of the first atom will be reported.\n");
462 if (EEL_FULL(fr
->ic
->eeltype
))
467 snew(sum_UgembU
, nener
);
469 /* Copy the random seed set by the user */
470 seed
= inputrec
->ld_seed
;
472 gmx::ThreeFry2x64
<16> rng(
473 seed
, gmx::RandomDomain::TestParticleInsertion
); // 16 bits internal counter => 2^16 * 2 = 131072 values per stream
474 gmx::UniformRealDistribution
<real
> dist
;
478 fp_tpi
= xvgropen(opt2fn("-tpi", nfile
, fnm
), "TPI energies", "Time (ps)",
479 "(kJ mol\\S-1\\N) / (nm\\S3\\N)", oenv
);
480 xvgr_subtitle(fp_tpi
, "f. are averages over one frame", oenv
);
481 snew(leg
, 4 + nener
);
483 sprintf(str
, "-kT log(<Ve\\S-\\betaU\\N>/<V>)");
484 leg
[e
++] = gmx_strdup(str
);
485 sprintf(str
, "f. -kT log<e\\S-\\betaU\\N>");
486 leg
[e
++] = gmx_strdup(str
);
487 sprintf(str
, "f. <e\\S-\\betaU\\N>");
488 leg
[e
++] = gmx_strdup(str
);
489 sprintf(str
, "f. V");
490 leg
[e
++] = gmx_strdup(str
);
491 sprintf(str
, "f. <Ue\\S-\\betaU\\N>");
492 leg
[e
++] = gmx_strdup(str
);
493 for (i
= 0; i
< ngid
; i
++)
495 sprintf(str
, "f. <U\\sVdW %s\\Ne\\S-\\betaU\\N>",
496 *(groups
->groupNames
[groups
->groups
[SimulationAtomGroupType::EnergyOutput
][i
]]));
497 leg
[e
++] = gmx_strdup(str
);
501 sprintf(str
, "f. <U\\sdisp c\\Ne\\S-\\betaU\\N>");
502 leg
[e
++] = gmx_strdup(str
);
506 for (i
= 0; i
< ngid
; i
++)
508 sprintf(str
, "f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
509 *(groups
->groupNames
[groups
->groups
[SimulationAtomGroupType::EnergyOutput
][i
]]));
510 leg
[e
++] = gmx_strdup(str
);
514 sprintf(str
, "f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
515 leg
[e
++] = gmx_strdup(str
);
517 if (EEL_FULL(fr
->ic
->eeltype
))
519 sprintf(str
, "f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
520 leg
[e
++] = gmx_strdup(str
);
523 xvgr_legend(fp_tpi
, 4 + nener
, leg
, oenv
);
524 for (i
= 0; i
< 4 + nener
; i
++)
538 /* Avoid frame step numbers <= -1 */
539 frame_step_prev
= -1;
541 bNotLastFrame
= read_first_frame(oenv
, &status
, opt2fn("-rerun", nfile
, fnm
), &rerun_fr
, TRX_NEED_X
);
544 if (rerun_fr
.natoms
- (bCavity
? nat_cavity
: 0) != mdatoms
->nr
- (a_tp1
- a_tp0
))
547 "Number of atoms in trajectory (%d)%s "
548 "is not equal the number in the run input file (%d) "
549 "minus the number of atoms to insert (%d)\n",
550 rerun_fr
.natoms
, bCavity
? " minus one" : "", mdatoms
->nr
, a_tp1
- a_tp0
);
553 refvolshift
= log(det(rerun_fr
.box
));
555 switch (inputrec
->eI
)
557 case eiTPI
: stepblocksize
= inputrec
->nstlist
; break;
558 case eiTPIC
: stepblocksize
= 1; break;
559 default: gmx_fatal(FARGS
, "Unknown integrator %s", ei_names
[inputrec
->eI
]);
562 while (bNotLastFrame
)
564 frame_step
= rerun_fr
.step
;
565 if (frame_step
<= frame_step_prev
)
567 /* We don't have step number in the trajectory file,
568 * or we have constant or decreasing step numbers.
569 * Ensure we have increasing step numbers, since we use
570 * the step numbers as a counter for random numbers.
572 frame_step
= frame_step_prev
+ 1;
574 frame_step_prev
= frame_step
;
576 lambda
= rerun_fr
.lambda
;
580 for (e
= 0; e
< nener
; e
++)
585 /* Copy the coordinates from the input trajectory */
586 auto x
= makeArrayRef(state_global
->x
);
587 for (i
= 0; i
< rerun_fr
.natoms
; i
++)
589 copy_rvec(rerun_fr
.x
[i
], x
[i
]);
591 copy_mat(rerun_fr
.box
, state_global
->box
);
592 const matrix
& box
= state_global
->box
;
597 bStateChanged
= TRUE
;
600 put_atoms_in_box(fr
->pbcType
, box
, x
);
602 /* Put all atoms except for the inserted ones on the grid */
603 rvec vzero
= { 0, 0, 0 };
604 rvec boxDiagonal
= { box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
] };
605 nbnxn_put_on_grid(fr
->nbv
.get(), box
, 0, vzero
, boxDiagonal
, nullptr, { 0, a_tp0
}, -1,
606 fr
->cginfo
, x
, 0, nullptr);
608 step
= cr
->nodeid
* stepblocksize
;
609 while (step
< nsteps
)
611 /* Restart random engine using the frame and insertion step
613 * Note that we need to draw several random values per iteration,
614 * but by using the internal subcounter functionality of ThreeFry2x64
615 * we can draw 131072 unique 64-bit values before exhausting
616 * the stream. This is a huge margin, and if something still goes
617 * wrong you will get an exception when the stream is exhausted.
619 rng
.restart(frame_step
, step
);
620 dist
.reset(); // erase any memory in the distribution
624 /* Random insertion in the whole volume */
625 bNS
= (step
% inputrec
->nstlist
== 0);
628 /* Generate a random position in the box */
629 for (d
= 0; d
< DIM
; d
++)
631 x_init
[d
] = dist(rng
) * state_global
->box
[d
][d
];
637 /* Random insertion around a cavity location
638 * given by the last coordinate of the trajectory.
644 /* Copy the location of the cavity */
645 copy_rvec(rerun_fr
.x
[rerun_fr
.natoms
- 1], x_init
);
649 /* Determine the center of mass of the last molecule */
652 for (i
= 0; i
< nat_cavity
; i
++)
654 for (d
= 0; d
< DIM
; d
++)
656 x_init
[d
] += mass_cavity
[i
]
657 * rerun_fr
.x
[rerun_fr
.natoms
- nat_cavity
+ i
][d
];
659 mass_tot
+= mass_cavity
[i
];
661 for (d
= 0; d
< DIM
; d
++)
663 x_init
[d
] /= mass_tot
;
671 for (int a
= a_tp0
; a
< a_tp1
; a
++)
676 /* Put the inserted molecule on it's own search grid */
677 nbnxn_put_on_grid(fr
->nbv
.get(), box
, 1, x_init
, x_init
, nullptr, { a_tp0
, a_tp1
},
678 -1, fr
->cginfo
, x
, 0, nullptr);
680 /* TODO: Avoid updating all atoms at every bNS step */
681 fr
->nbv
->setAtomProperties(*mdatoms
, fr
->cginfo
);
683 fr
->nbv
->constructPairlist(InteractionLocality::Local
, top
.excls
, step
, nrnb
);
688 /* Add random displacement uniformly distributed in a sphere
689 * of radius rtpi. We don't need to do this is we generate
690 * a new center location every step.
693 if (bCavity
|| inputrec
->nstlist
> 1)
695 /* Generate coordinates within |dx|=drmax of x_init */
698 for (d
= 0; d
< DIM
; d
++)
700 dx
[d
] = (2 * dist(rng
) - 1) * drmax
;
702 } while (norm2(dx
) > drmax
* drmax
);
703 rvec_add(x_init
, dx
, x_tp
);
707 copy_rvec(x_init
, x_tp
);
710 if (a_tp1
- a_tp0
== 1)
712 /* Insert a single atom, just copy the insertion location */
713 copy_rvec(x_tp
, x
[a_tp0
]);
717 /* Copy the coordinates from the top file */
718 for (i
= a_tp0
; i
< a_tp1
; i
++)
720 copy_rvec(x_mol
[i
- a_tp0
], x
[i
]);
722 /* Rotate the molecule randomly */
723 real angleX
= 2 * M_PI
* dist(rng
);
724 real angleY
= 2 * M_PI
* dist(rng
);
725 real angleZ
= 2 * M_PI
* dist(rng
);
726 rotate_conf(a_tp1
- a_tp0
, state_global
->x
.rvec_array() + a_tp0
, nullptr, angleX
,
728 /* Shift to the insertion location */
729 for (i
= a_tp0
; i
< a_tp1
; i
++)
731 rvec_inc(x
[i
], x_tp
);
735 /* Note: NonLocal refers to the inserted molecule */
736 fr
->nbv
->convertCoordinates(AtomLocality::NonLocal
, false, x
);
738 /* Clear some matrix variables */
739 clear_mat(force_vir
);
740 clear_mat(shake_vir
);
744 /* Calc energy (no forces) on new positions.
745 * Since we only need the intermolecular energy
746 * and the RF exclusion terms of the inserted molecule occur
747 * within a single charge group we can pass NULL for the graph.
748 * This also avoids shifts that would move charge groups
750 /* Make do_force do a single node force calculation */
753 // TPI might place a particle so close that the potential
754 // is infinite. Since this is intended to happen, we
755 // temporarily suppress any exceptions that the processor
756 // might raise, then restore the old behaviour.
757 std::fenv_t floatingPointEnvironment
;
758 std::feholdexcept(&floatingPointEnvironment
);
759 do_force(fplog
, cr
, ms
, inputrec
, nullptr, nullptr, imdSession
, pull_work
, step
, nrnb
,
760 wcycle
, &top
, state_global
->box
, state_global
->x
.arrayRefWithPadding(),
761 &state_global
->hist
, f
.arrayRefWithPadding(), force_vir
, mdatoms
, enerd
, fcd
,
762 state_global
->lambda
, nullptr, fr
, runScheduleWork
, nullptr, mu_tot
, t
, nullptr,
763 GMX_FORCE_NONBONDED
| GMX_FORCE_ENERGY
| (bStateChanged
? GMX_FORCE_STATECHANGED
: 0),
764 DDBalanceRegionHandler(nullptr));
765 std::feclearexcept(FE_DIVBYZERO
| FE_INVALID
| FE_OVERFLOW
);
766 std::feupdateenv(&floatingPointEnvironment
);
769 bStateChanged
= FALSE
;
771 if (fr
->dispersionCorrection
)
773 /* Calculate long range corrections to pressure and energy */
774 const DispersionCorrection::Correction correction
=
775 fr
->dispersionCorrection
->calculate(state_global
->box
, lambda
);
776 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
777 enerd
->term
[F_DISPCORR
] = correction
.energy
;
778 enerd
->term
[F_EPOT
] += correction
.energy
;
779 enerd
->term
[F_PRES
] += correction
.pressure
;
780 enerd
->term
[F_DVDL
] += correction
.dvdl
;
784 enerd
->term
[F_DISPCORR
] = 0;
786 if (EEL_RF(fr
->ic
->eeltype
))
788 enerd
->term
[F_EPOT
] += rfExclusionEnergy
;
791 epot
= enerd
->term
[F_EPOT
];
792 bEnergyOutOfBounds
= FALSE
;
794 /* If the compiler doesn't optimize this check away
795 * we catch the NAN energies.
796 * The epot>GMX_REAL_MAX check catches inf values,
797 * which should nicely result in embU=0 through the exp below,
798 * but it does not hurt to check anyhow.
800 /* Non-bonded Interaction usually diverge at r=0.
801 * With tabulated interaction functions the first few entries
802 * should be capped in a consistent fashion between
803 * repulsion, dispersion and Coulomb to avoid accidental
804 * negative values in the total energy.
805 * The table generation code in tables.c does this.
806 * With user tbales the user should take care of this.
808 if (epot
!= epot
|| epot
> GMX_REAL_MAX
)
810 bEnergyOutOfBounds
= TRUE
;
812 if (bEnergyOutOfBounds
)
817 "\n time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n", t
,
818 static_cast<int>(step
), epot
);
824 // Exponent argument is fine in SP range, but output can be in DP range
825 embU
= exp(static_cast<double>(-beta
* epot
));
827 /* Determine the weighted energy contributions of each energy group */
829 sum_UgembU
[e
++] += epot
* embU
;
832 for (i
= 0; i
< ngid
; i
++)
834 sum_UgembU
[e
++] += enerd
->grpp
.ener
[egBHAMSR
][GID(i
, gid_tp
, ngid
)] * embU
;
839 for (i
= 0; i
< ngid
; i
++)
841 sum_UgembU
[e
++] += enerd
->grpp
.ener
[egLJSR
][GID(i
, gid_tp
, ngid
)] * embU
;
846 sum_UgembU
[e
++] += enerd
->term
[F_DISPCORR
] * embU
;
850 for (i
= 0; i
< ngid
; i
++)
852 sum_UgembU
[e
++] += enerd
->grpp
.ener
[egCOULSR
][GID(i
, gid_tp
, ngid
)] * embU
;
856 sum_UgembU
[e
++] += rfExclusionEnergy
* embU
;
858 if (EEL_FULL(fr
->ic
->eeltype
))
860 sum_UgembU
[e
++] += enerd
->term
[F_COUL_RECIP
] * embU
;
865 if (embU
== 0 || beta
* epot
> bU_bin_limit
)
871 i
= gmx::roundToInt((bU_logV_bin_limit
- (beta
* epot
- logV
+ refvolshift
)) * invbinw
);
878 realloc_bins(&bin
, &nbin
, i
+ 10);
885 fprintf(debug
, "TPI %7d %12.5e %12.5f %12.5f %12.5f\n", static_cast<int>(step
),
886 epot
, x_tp
[XX
], x_tp
[YY
], x_tp
[ZZ
]);
889 if (dump_pdb
&& epot
<= dump_ener
)
891 sprintf(str
, "t%g_step%d.pdb", t
, static_cast<int>(step
));
892 sprintf(str2
, "t: %f step %d ener: %f", t
, static_cast<int>(step
), epot
);
893 write_sto_conf_mtop(str
, str2
, top_global
, state_global
->x
.rvec_array(),
894 state_global
->v
.rvec_array(), inputrec
->pbcType
, state_global
->box
);
898 if ((step
/ stepblocksize
) % cr
->nnodes
!= cr
->nodeid
)
900 /* Skip all steps assigned to the other MPI ranks */
901 step
+= (cr
->nnodes
- 1) * stepblocksize
;
907 /* When running in parallel sum the energies over the processes */
908 gmx_sumd(1, &sum_embU
, cr
);
909 gmx_sumd(nener
, sum_UgembU
, cr
);
914 VembU_all
+= V
* sum_embU
/ nsteps
;
918 if (mdrunOptions
.verbose
|| frame
% 10 == 0 || frame
< 10)
920 fprintf(stderr
, "mu %10.3e <mu> %10.3e\n", -log(sum_embU
/ nsteps
) / beta
,
921 -log(VembU_all
/ V_all
) / beta
);
924 fprintf(fp_tpi
, "%10.3f %12.5e %12.5e %12.5e %12.5e", t
,
925 VembU_all
== 0 ? 20 / beta
: -log(VembU_all
/ V_all
) / beta
,
926 sum_embU
== 0 ? 20 / beta
: -log(sum_embU
/ nsteps
) / beta
, sum_embU
/ nsteps
, V
);
927 for (e
= 0; e
< nener
; e
++)
929 fprintf(fp_tpi
, " %12.5e", sum_UgembU
[e
] / nsteps
);
931 fprintf(fp_tpi
, "\n");
935 bNotLastFrame
= read_next_frame(oenv
, status
, &rerun_fr
);
936 } /* End of the loop */
937 walltime_accounting_end_time(walltime_accounting
);
941 if (fp_tpi
!= nullptr)
946 if (fplog
!= nullptr)
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
);
953 if (!std::isfinite(mu
))
956 "\nThe computed chemical potential is not finite - consider increasing the "
957 "number of steps and/or the number of frames to insert into.\n");
961 /* Write the Boltzmann factor histogram */
964 /* When running in parallel sum the bins over the processes */
967 realloc_bins(&bin
, &nbin
, i
);
968 gmx_sumd(nbin
, bin
, cr
);
972 fp_tpi
= xvgropen(opt2fn("-tpid", nfile
, fnm
), "TPI energy distribution",
973 "\\betaU - log(V/<V>)", "count", oenv
);
974 sprintf(str
, "number \\betaU > %g: %9.3e", bU_bin_limit
, bin
[0]);
975 xvgr_subtitle(fp_tpi
, str
, oenv
);
976 xvgr_legend(fp_tpi
, 2, tpid_leg
, oenv
);
977 for (i
= nbin
- 1; i
> 0; i
--)
979 bUlogV
= -i
/ invbinw
+ bU_logV_bin_limit
- refvolshift
+ log(V_all
/ frame
);
980 fprintf(fp_tpi
, "%6.2f %10d %12.5e\n", bUlogV
, roundToInt(bin
[i
]),
981 bin
[i
] * exp(-bUlogV
) * V_all
/ VembU_all
);
989 walltime_accounting_set_nsteps_done(walltime_accounting
, frame
* inputrec
->nsteps
);