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->ndisrepairs = dd->npair;
205 snew(hist->disre_rm3tav, hist->ndisrepairs);
207 /* Allocate space for a copy of rm3tav,
208 * so we can call do_force without modifying the state.
210 snew(dd->rm3tav, dd->npair);
212 /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
213 * averaged over the processors in one call (in calc_disre_R_6)
215 snew(dd->Rt_6, 2 * dd->nres);
216 dd->Rtav_6 = &(dd->Rt_6[dd->nres]);
218 ptr = getenv("GMX_DISRE_ENSEMBLE_SIZE");
219 if ((disResRunMode == DisResRunMode::MDRun) && ms != nullptr && ptr != nullptr && !bIsREMD)
223 sscanf(ptr, "%d", &dd->nsystems);
226 fprintf(fplog, "Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n", dd->nsystems);
228 /* This check is only valid on MASTER(cr), so probably
229 * ensemble-averaged distance restraints are broken on more
230 * than one processor per simulation system. */
231 if (ddRole == DDRole::Master)
233 check_multi_int(fplog, ms, dd->nsystems, "the number of systems per ensemble", FALSE);
235 gmx_bcast(sizeof(int), &dd->nsystems, communicator);
237 /* We use to allow any value of nsystems which was a divisor
238 * of ms->numSimulations_. But this required an extra communicator which
239 * pulled in mpi.h in nearly all C files.
241 if (!(ms->numSimulations_ == 1 || ms->numSimulations_ == dd->nsystems))
244 "GMX_DISRE_ENSEMBLE_SIZE (%d) is not equal to 1 or the number of systems "
245 "(option -multidir) %d",
247 ms->numSimulations_);
251 fprintf(fplog, "Our ensemble consists of systems:");
252 for (int i = 0; i < dd->nsystems; i++)
254 fprintf(fplog, " %d", (ms->simulationIndex_ / dd->nsystems) * dd->nsystems + i);
256 fprintf(fplog, "\n");
259 GMX_UNUSED_VALUE(communicator);
267 if (dd->nsystems == 1)
269 dd->Rtl_6 = dd->Rt_6;
273 snew(dd->Rtl_6, dd->nres);
280 fprintf(fplog, "There are %d distance restraints involving %d atom pairs\n", dd->nres, dd->npair);
282 /* Have to avoid g_disre de-referencing cr blindly, mdrun not
283 * doing consistency checks for ensemble-averaged distance
284 * restraints when that's not happening, and only doing those
285 * checks from appropriate processes (since check_multi_int is
286 * too broken to check whether the communication will
288 if ((disResRunMode == DisResRunMode::MDRun) && ms && dd->nsystems > 1 && (ddRole == DDRole::Master))
290 check_multi_int(fplog, ms, dd->nres, "the number of distance restraints", FALSE);
292 please_cite(fplog, "Tropp80a");
293 please_cite(fplog, "Torda89a");
297 void calc_disres_R_6(const t_commrec* cr,
298 const gmx_multisim_t* ms,
300 const t_iatom forceatoms[],
304 const history_t* hist)
307 real * rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6;
308 real ETerm, ETerm1, cf1 = 0, cf2 = 0;
311 bTav = (dd->dr_tau != 0);
322 /* scaling factor to smoothly turn on the restraint forces *
323 * when using time averaging */
324 dd->exp_min_t_tau = hist->disre_initf * ETerm;
326 cf1 = dd->exp_min_t_tau;
327 cf2 = 1.0 / (1.0 - dd->exp_min_t_tau);
330 for (int res = 0; res < dd->nres; res++)
336 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
337 * the total number of atoms pairs is nfa/3 */
338 for (int fa = 0; fa < nfa; fa += 3)
340 int type = forceatoms[fa];
341 int res = type - dd->type_min;
343 int ai = forceatoms[fa + 1];
344 int aj = forceatoms[fa + 2];
348 pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
352 rvec_sub(x[ai], x[aj], dx);
354 real rt2 = iprod(dx, dx);
355 real rt_1 = gmx::invsqrt(rt2);
356 real rt_3 = rt_1 * rt_1 * rt_1;
358 rt[pair] = rt2 * rt_1;
361 /* Here we update rm3tav in t_disresdata using the data
363 * Thus the results stay correct when this routine
364 * is called multiple times.
366 rm3tav[pair] = cf2 * ((ETerm - cf1) * hist->disre_rm3tav[pair] + ETerm1 * rt_3);
373 /* Currently divide_bondeds_over_threads() ensures that pairs within
374 * the same restraint get assigned to the same thread, so we could
375 * run this loop thread-parallel.
377 Rt_6[res] += rt_3 * rt_3;
378 Rtav_6[res] += rm3tav[pair] * rm3tav[pair];
381 /* NOTE: Rt_6 and Rtav_6 are stored consecutively in memory */
382 if (cr && DOMAINDECOMP(cr))
384 gmx_sum(2 * dd->nres, dd->Rt_6, cr);
387 if (dd->nsystems > 1)
389 real invn = 1.0 / dd->nsystems;
391 for (int res = 0; res < dd->nres; res++)
393 Rtl_6[res] = Rt_6[res];
398 GMX_ASSERT(cr != nullptr && ms != nullptr, "We need multisim with nsystems>1");
399 gmx_sum_sim(2 * dd->nres, dd->Rt_6, ms);
401 if (DOMAINDECOMP(cr))
403 gmx_bcast(2 * dd->nres, dd->Rt_6, cr->mpi_comm_mygroup);
407 /* Store the base forceatoms pointer, so we can re-calculate the pair
408 * index in ta_disres() for indexing pair data in t_disresdata when
409 * using thread parallelization.
411 dd->forceatomsStart = forceatoms;
416 real ta_disres(int nfa,
417 const t_iatom* forceatoms,
423 real gmx_unused lambda,
424 real gmx_unused* dvdlambda,
425 gmx::ArrayRef<const real> /*charge*/,
426 t_fcdata gmx_unused* fcd,
427 t_disresdata* disresdata,
428 t_oriresdata gmx_unused* oriresdata,
429 int gmx_unused* global_atom_index)
431 const real seven_three = 7.0 / 3.0;
435 real smooth_fc, Rt, Rtav, rt2, *Rtl_6, *Rt_6, *Rtav_6;
436 real k0, f_scal = 0, fmax_scal, fk_scal, fij;
437 real tav_viol, instant_viol, mixed_viol, violtot, vtot;
438 real tav_viol_Rtav7, instant_viol_Rtav7;
440 gmx_bool bConservative, bMixed, bViolation;
443 DistanceRestraintWeighting dr_weighting = disresdata->dr_weighting;
444 dr_bMixed = disresdata->dr_bMixed;
445 Rtl_6 = disresdata->Rtl_6;
446 Rt_6 = disresdata->Rt_6;
447 Rtav_6 = disresdata->Rtav_6;
449 tav_viol = instant_viol = mixed_viol = tav_viol_Rtav7 = instant_viol_Rtav7 = 0;
451 smooth_fc = disresdata->dr_fc;
452 if (disresdata->dr_tau != 0)
454 /* scaling factor to smoothly turn on the restraint forces *
455 * when using time averaging */
456 smooth_fc *= (1.0 - disresdata->exp_min_t_tau);
462 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
463 * the total number of atoms pairs is nfa/3 */
464 int faOffset = static_cast<int>(forceatoms - disresdata->forceatomsStart);
465 for (int fa = 0; fa < nfa; fa += 3)
467 int type = forceatoms[fa];
468 int npair = ip[type].disres.npair;
469 up1 = ip[type].disres.up1;
470 up2 = ip[type].disres.up2;
471 low = ip[type].disres.low;
472 k0 = smooth_fc * ip[type].disres.kfac;
474 int res = type - disresdata->type_min;
476 /* save some flops when there is only one pair */
477 if (ip[type].disres.type != 2)
479 bConservative = (dr_weighting == DistanceRestraintWeighting::Conservative) && (npair > 1);
481 Rt = gmx::invsixthroot(Rt_6[res]);
482 Rtav = gmx::invsixthroot(Rtav_6[res]);
486 /* When rtype=2 use instantaneous not ensemble averaged distance */
487 bConservative = (npair > 1);
489 Rt = gmx::invsixthroot(Rtl_6[res]);
496 tav_viol = Rtav - up1;
501 tav_viol = Rtav - low;
510 /* Add 1/npair energy and violation for each of the npair pairs */
511 real pairFac = 1 / static_cast<real>(npair);
514 * there is no real potential when time averaging is applied
516 vtot += 0.5 * k0 * gmx::square(tav_viol) * pairFac;
519 f_scal = -k0 * tav_viol;
520 violtot += fabs(tav_viol) * pairFac;
528 instant_viol = Rt - up1;
539 instant_viol = Rt - low;
552 mixed_viol = std::sqrt(tav_viol * instant_viol);
553 f_scal = -k0 * mixed_viol;
554 violtot += mixed_viol * pairFac;
561 fmax_scal = -k0 * (up2 - up1);
562 /* Correct the force for the number of restraints */
565 f_scal = std::max(f_scal, fmax_scal);
568 f_scal *= Rtav / Rtav_6[res];
572 f_scal /= 2 * mixed_viol;
573 tav_viol_Rtav7 = tav_viol * Rtav / Rtav_6[res];
574 instant_viol_Rtav7 = instant_viol * Rt / Rt_6[res];
580 f_scal = std::max(f_scal, fmax_scal);
583 /* Exert the force ... */
585 int pair = (faOffset + fa) / 3;
586 int ai = forceatoms[fa + 1];
587 int aj = forceatoms[fa + 2];
591 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
595 rvec_sub(x[ai], x[aj], dx);
599 weight_rt_1 = gmx::invsqrt(rt2);
605 weight_rt_1 *= std::pow(disresdata->rm3tav[pair], seven_three);
610 tav_viol_Rtav7 * std::pow(disresdata->rm3tav[pair], seven_three)
612 / (disresdata->rt[pair] * gmx::power6(disresdata->rt[pair]));
616 fk_scal = f_scal * weight_rt_1;
618 for (int m = 0; m < DIM; m++)
620 fij = fk_scal * dx[m];
626 fshift[ki][m] += fij;
627 fshift[CENTRAL][m] -= fij;
634 disresdata->sumviol += violtot;
640 void update_disres_history(const t_disresdata& dd, history_t* hist)
644 /* Copy the new time averages that have been calculated
645 * in calc_disres_R_6.
647 hist->disre_initf = dd.exp_min_t_tau;
648 for (int pair = 0; pair < dd.npair; pair++)
650 hist->disre_rm3tav[pair] = dd.rm3tav[pair];