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,2021, 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.
38 /* This file is completely threadsafe - keep it that way! */
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdrunutility/multisim.h"
55 #include "gromacs/mdtypes/commrec.h"
56 #include "gromacs/mdtypes/fcdata.h"
57 #include "gromacs/mdtypes/inputrec.h"
58 #include "gromacs/mdtypes/mdatom.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/mdtypes/state.h"
61 #include "gromacs/pbcutil/ishift.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/topology/mtop_util.h"
64 #include "gromacs/topology/topology.h"
65 #include "gromacs/utility/arrayref.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/futil.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/pleasecite.h"
70 #include "gromacs/utility/smalloc.h"
72 void init_disres(FILE* fplog,
73 const gmx_mtop_t& mtop,
75 DisResRunMode disResRunMode,
78 MPI_Comm communicator,
79 const gmx_multisim_t* ms,
86 int type_min, type_max;
88 if (gmx_mtop_ftype_count(mtop, F_DISRES) == 0)
97 fprintf(fplog, "Initializing the distance restraints\n");
100 dd->dr_weighting = ir->eDisreWeighting;
101 dd->dr_fc = ir->dr_fc;
102 if (EI_DYNAMICS(ir->eI))
104 dd->dr_tau = ir->dr_tau;
110 if (dd->dr_tau == 0.0)
112 dd->dr_bMixed = FALSE;
117 /* We store the r^-6 time averages in an array that is indexed
118 * with the local disres iatom index, so this doesn't work with DD.
119 * Note that DD is not initialized yet here, so we check that we are on multiple ranks,
120 * but there are probably also issues with e.g. NM MPI parallelization.
122 if ((disResRunMode == DisResRunMode::MDRun) && (numRanks == NumRanks::Multiple))
125 "Time-averaged distance restraints are not supported with MPI "
126 "parallelization. You can use OpenMP parallelization on a single node.");
129 dd->dr_bMixed = ir->bDisreMixed;
130 dd->ETerm = std::exp(-(ir->delta_t / ir->dr_tau));
132 dd->ETerm1 = 1.0 - dd->ETerm;
138 for (const auto il : IListRange(mtop))
140 if (il.nmol() > 1 && !il.list()[F_DISRES].empty()
141 && ir->eDisre != DistanceRestraintRefinement::Ensemble)
144 "NMR distance restraints with multiple copies of the same molecule are "
145 "currently only supported with ensemble averaging. If you just want to "
146 "restrain distances between atom pairs using a flat-bottomed potential, use "
147 "a restraint potential (bonds type 10) instead.");
151 for (int fa = 0; fa < il.list()[F_DISRES].size(); fa += 3)
153 int type = il.list()[F_DISRES].iatoms[fa];
156 int npair = mtop.ffparams.iparams[type].disres.npair;
159 dd->nres += (ir->eDisre == DistanceRestraintRefinement::Ensemble ? 1 : il.nmol());
160 dd->npair += il.nmol() * npair;
163 type_min = std::min(type_min, type);
164 type_max = std::max(type_max, type);
169 if ((disResRunMode == DisResRunMode::MDRun) && (numRanks == NumRanks::Multiple) && ir->nstdisreout > 0)
171 /* With DD we currently only have local pair information available */
173 "With MPI parallelization distance-restraint pair output is not supported. Use "
174 "nstdisreout=0 or use OpenMP parallelization on a single node.");
177 /* For communicating and/or reducing (sums of) r^-6 for pairs over threads
178 * we use multiple arrays in t_disresdata. We need to have unique indices
179 * for each restraint that work over threads and MPI ranks. To this end
180 * we use the type index. These should all be in one block and there should
181 * be dd->nres types, but we check for this here.
182 * This setup currently does not allow for multiple copies of the same
183 * molecule without ensemble averaging, this is check for above.
186 type_max - type_min + 1 == dd->nres,
187 "All distance restraint parameter entries in the topology should be consecutive");
189 dd->type_min = type_min;
191 snew(dd->rt, dd->npair);
193 if (dd->dr_tau != 0.0)
195 GMX_RELEASE_ASSERT(state != nullptr,
196 "We need a valid state when using time-averaged distance restraints");
199 /* Set the "history lack" factor to 1 */
200 state->flags |= enumValueToBitMask(StateEntry::DisreInitF);
201 hist->disre_initf = 1.0;
202 /* Allocate space for the r^-3 time averages */
203 state->flags |= enumValueToBitMask(StateEntry::DisreRm3Tav);
204 hist->disre_rm3tav.resize(dd->npair);
206 /* Allocate space for a copy of rm3tav,
207 * so we can call do_force without modifying the state.
209 snew(dd->rm3tav, dd->npair);
211 /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
212 * averaged over the processors in one call (in calc_disre_R_6)
214 snew(dd->Rt_6, 2 * dd->nres);
215 dd->Rtav_6 = &(dd->Rt_6[dd->nres]);
217 ptr = getenv("GMX_DISRE_ENSEMBLE_SIZE");
218 if ((disResRunMode == DisResRunMode::MDRun) && ms != nullptr && ptr != nullptr && !bIsREMD)
222 sscanf(ptr, "%d", &dd->nsystems);
225 fprintf(fplog, "Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n", dd->nsystems);
227 /* This check is only valid on MASTER(cr), so probably
228 * ensemble-averaged distance restraints are broken on more
229 * than one processor per simulation system. */
230 if (ddRole == DDRole::Master)
232 check_multi_int(fplog, ms, dd->nsystems, "the number of systems per ensemble", FALSE);
234 gmx_bcast(sizeof(int), &dd->nsystems, communicator);
236 /* We use to allow any value of nsystems which was a divisor
237 * of ms->numSimulations_. But this required an extra communicator which
238 * pulled in mpi.h in nearly all C files.
240 if (!(ms->numSimulations_ == 1 || ms->numSimulations_ == dd->nsystems))
243 "GMX_DISRE_ENSEMBLE_SIZE (%d) is not equal to 1 or the number of systems "
244 "(option -multidir) %d",
246 ms->numSimulations_);
250 fprintf(fplog, "Our ensemble consists of systems:");
251 for (int i = 0; i < dd->nsystems; i++)
253 fprintf(fplog, " %d", (ms->simulationIndex_ / dd->nsystems) * dd->nsystems + i);
255 fprintf(fplog, "\n");
258 GMX_UNUSED_VALUE(communicator);
266 if (dd->nsystems == 1)
268 dd->Rtl_6 = dd->Rt_6;
272 snew(dd->Rtl_6, dd->nres);
279 fprintf(fplog, "There are %d distance restraints involving %d atom pairs\n", dd->nres, dd->npair);
281 /* Have to avoid g_disre de-referencing cr blindly, mdrun not
282 * doing consistency checks for ensemble-averaged distance
283 * restraints when that's not happening, and only doing those
284 * checks from appropriate processes (since check_multi_int is
285 * too broken to check whether the communication will
287 if ((disResRunMode == DisResRunMode::MDRun) && ms && dd->nsystems > 1 && (ddRole == DDRole::Master))
289 check_multi_int(fplog, ms, dd->nres, "the number of distance restraints", FALSE);
291 please_cite(fplog, "Tropp80a");
292 please_cite(fplog, "Torda89a");
296 void calc_disres_R_6(const t_commrec* cr,
297 const gmx_multisim_t* ms,
299 const t_iatom forceatoms[],
303 const history_t* hist)
306 real * rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6;
307 real ETerm, ETerm1, cf1 = 0, cf2 = 0;
310 bTav = (dd->dr_tau != 0);
321 /* scaling factor to smoothly turn on the restraint forces *
322 * when using time averaging */
323 dd->exp_min_t_tau = hist->disre_initf * ETerm;
325 cf1 = dd->exp_min_t_tau;
326 cf2 = 1.0 / (1.0 - dd->exp_min_t_tau);
329 for (int res = 0; res < dd->nres; res++)
335 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
336 * the total number of atoms pairs is nfa/3 */
337 for (int fa = 0; fa < nfa; fa += 3)
339 int type = forceatoms[fa];
340 int res = type - dd->type_min;
342 int ai = forceatoms[fa + 1];
343 int aj = forceatoms[fa + 2];
347 pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
351 rvec_sub(x[ai], x[aj], dx);
353 real rt2 = iprod(dx, dx);
354 real rt_1 = gmx::invsqrt(rt2);
355 real rt_3 = rt_1 * rt_1 * rt_1;
357 rt[pair] = rt2 * rt_1;
360 /* Here we update rm3tav in t_disresdata using the data
362 * Thus the results stay correct when this routine
363 * is called multiple times.
365 rm3tav[pair] = cf2 * ((ETerm - cf1) * hist->disre_rm3tav[pair] + ETerm1 * rt_3);
372 /* Currently divide_bondeds_over_threads() ensures that pairs within
373 * the same restraint get assigned to the same thread, so we could
374 * run this loop thread-parallel.
376 Rt_6[res] += rt_3 * rt_3;
377 Rtav_6[res] += rm3tav[pair] * rm3tav[pair];
380 /* NOTE: Rt_6 and Rtav_6 are stored consecutively in memory */
381 if (cr && DOMAINDECOMP(cr))
383 gmx_sum(2 * dd->nres, dd->Rt_6, cr);
386 if (dd->nsystems > 1)
388 real invn = 1.0 / dd->nsystems;
390 for (int res = 0; res < dd->nres; res++)
392 Rtl_6[res] = Rt_6[res];
397 GMX_ASSERT(cr != nullptr && ms != nullptr, "We need multisim with nsystems>1");
398 gmx_sum_sim(2 * dd->nres, dd->Rt_6, ms);
400 if (DOMAINDECOMP(cr))
402 gmx_bcast(2 * dd->nres, dd->Rt_6, cr->mpi_comm_mygroup);
406 /* Store the base forceatoms pointer, so we can re-calculate the pair
407 * index in ta_disres() for indexing pair data in t_disresdata when
408 * using thread parallelization.
410 dd->forceatomsStart = forceatoms;
415 real ta_disres(int nfa,
416 const t_iatom* forceatoms,
422 real gmx_unused lambda,
423 real gmx_unused* dvdlambda,
424 gmx::ArrayRef<const real> /*charge*/,
425 t_fcdata gmx_unused* fcd,
426 t_disresdata* disresdata,
427 t_oriresdata gmx_unused* oriresdata,
428 int gmx_unused* global_atom_index)
430 const real seven_three = 7.0 / 3.0;
434 real smooth_fc, Rt, Rtav, rt2, *Rtl_6, *Rt_6, *Rtav_6;
435 real k0, f_scal = 0, fmax_scal, fk_scal, fij;
436 real tav_viol, instant_viol, mixed_viol, violtot, vtot;
437 real tav_viol_Rtav7, instant_viol_Rtav7;
439 gmx_bool bConservative, bMixed, bViolation;
442 DistanceRestraintWeighting dr_weighting = disresdata->dr_weighting;
443 dr_bMixed = disresdata->dr_bMixed;
444 Rtl_6 = disresdata->Rtl_6;
445 Rt_6 = disresdata->Rt_6;
446 Rtav_6 = disresdata->Rtav_6;
448 tav_viol = instant_viol = mixed_viol = tav_viol_Rtav7 = instant_viol_Rtav7 = 0;
450 smooth_fc = disresdata->dr_fc;
451 if (disresdata->dr_tau != 0)
453 /* scaling factor to smoothly turn on the restraint forces *
454 * when using time averaging */
455 smooth_fc *= (1.0 - disresdata->exp_min_t_tau);
461 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
462 * the total number of atoms pairs is nfa/3 */
463 int faOffset = static_cast<int>(forceatoms - disresdata->forceatomsStart);
464 for (int fa = 0; fa < nfa; fa += 3)
466 int type = forceatoms[fa];
467 int npair = ip[type].disres.npair;
468 up1 = ip[type].disres.up1;
469 up2 = ip[type].disres.up2;
470 low = ip[type].disres.low;
471 k0 = smooth_fc * ip[type].disres.kfac;
473 int res = type - disresdata->type_min;
475 /* save some flops when there is only one pair */
476 if (ip[type].disres.type != 2)
478 bConservative = (dr_weighting == DistanceRestraintWeighting::Conservative) && (npair > 1);
480 Rt = gmx::invsixthroot(Rt_6[res]);
481 Rtav = gmx::invsixthroot(Rtav_6[res]);
485 /* When rtype=2 use instantaneous not ensemble averaged distance */
486 bConservative = (npair > 1);
488 Rt = gmx::invsixthroot(Rtl_6[res]);
495 tav_viol = Rtav - up1;
500 tav_viol = Rtav - low;
509 /* Add 1/npair energy and violation for each of the npair pairs */
510 real pairFac = 1 / static_cast<real>(npair);
513 * there is no real potential when time averaging is applied
515 vtot += 0.5 * k0 * gmx::square(tav_viol) * pairFac;
518 f_scal = -k0 * tav_viol;
519 violtot += fabs(tav_viol) * pairFac;
527 instant_viol = Rt - up1;
538 instant_viol = Rt - low;
551 mixed_viol = std::sqrt(tav_viol * instant_viol);
552 f_scal = -k0 * mixed_viol;
553 violtot += mixed_viol * pairFac;
560 fmax_scal = -k0 * (up2 - up1);
561 /* Correct the force for the number of restraints */
564 f_scal = std::max(f_scal, fmax_scal);
567 f_scal *= Rtav / Rtav_6[res];
571 f_scal /= 2 * mixed_viol;
572 tav_viol_Rtav7 = tav_viol * Rtav / Rtav_6[res];
573 instant_viol_Rtav7 = instant_viol * Rt / Rt_6[res];
579 f_scal = std::max(f_scal, fmax_scal);
582 /* Exert the force ... */
584 int pair = (faOffset + fa) / 3;
585 int ai = forceatoms[fa + 1];
586 int aj = forceatoms[fa + 2];
587 int ki = gmx::c_centralShiftIndex;
590 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
594 rvec_sub(x[ai], x[aj], dx);
598 weight_rt_1 = gmx::invsqrt(rt2);
604 weight_rt_1 *= std::pow(disresdata->rm3tav[pair], seven_three);
609 tav_viol_Rtav7 * std::pow(disresdata->rm3tav[pair], seven_three)
611 / (disresdata->rt[pair] * gmx::power6(disresdata->rt[pair]));
615 fk_scal = f_scal * weight_rt_1;
617 for (int m = 0; m < DIM; m++)
619 fij = fk_scal * dx[m];
625 fshift[ki][m] += fij;
626 fshift[gmx::c_centralShiftIndex][m] -= fij;
633 disresdata->sumviol += violtot;
639 void update_disres_history(const t_disresdata& dd, history_t* hist)
643 /* Copy the new time averages that have been calculated
644 * in calc_disres_R_6.
646 hist->disre_initf = dd.exp_min_t_tau;
647 for (int pair = 0; pair < dd.npair; pair++)
649 hist->disre_rm3tav[pair] = dd.rm3tav[pair];