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, 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.
37 /* This file is completely threadsafe - keep it that way! */
50 #include "gromacs/gmxlib/network.h"
51 #include "gromacs/math/functions.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdtypes/commrec.h"
54 #include "gromacs/mdtypes/fcdata.h"
55 #include "gromacs/mdtypes/inputrec.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/mdtypes/state.h"
58 #include "gromacs/pbcutil/ishift.h"
59 #include "gromacs/pbcutil/mshift.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/topology/mtop_util.h"
62 #include "gromacs/topology/topology.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/pleasecite.h"
67 #include "gromacs/utility/smalloc.h"
69 void init_disres(FILE *fplog, const gmx_mtop_t *mtop,
70 t_inputrec *ir, const t_commrec *cr,
71 const gmx_multisim_t *ms,
72 t_fcdata *fcd, t_state *state, gmx_bool bIsREMD)
74 int fa, nmol, npair, np;
77 gmx_mtop_ilistloop_t iloop;
79 int type_min, type_max;
83 if (gmx_mtop_ftype_count(mtop, F_DISRES) == 0)
92 fprintf(fplog, "Initializing the distance restraints\n");
95 dd->dr_weighting = ir->eDisreWeighting;
96 dd->dr_fc = ir->dr_fc;
97 if (EI_DYNAMICS(ir->eI))
99 dd->dr_tau = ir->dr_tau;
105 if (dd->dr_tau == 0.0)
107 dd->dr_bMixed = FALSE;
112 /* We store the r^-6 time averages in an array that is indexed
113 * with the local disres iatom index, so this doesn't work with DD.
114 * Note that DD is not initialized yet here, so we check for PAR(cr),
115 * but there are probably also issues with e.g. NM MPI parallelization.
119 gmx_fatal(FARGS, "Time-averaged distance restraints are not supported with MPI parallelization. You can use OpenMP parallelization on a single node.");
122 dd->dr_bMixed = ir->bDisreMixed;
123 dd->ETerm = std::exp(-(ir->delta_t/ir->dr_tau));
125 dd->ETerm1 = 1.0 - dd->ETerm;
131 iloop = gmx_mtop_ilistloop_init(mtop);
132 while (const InteractionLists *il = gmx_mtop_ilistloop_next(iloop, &nmol))
134 if (nmol > 1 && (*il)[F_DISRES].size() > 0 && ir->eDisre != edrEnsemble)
136 gmx_fatal(FARGS, "NMR distance restraints with multiple copies of the same molecule are currently only supported with ensemble averaging. If you just want to restrain distances between atom pairs using a flat-bottomed potential, use a restraint potential (bonds type 10) instead.");
140 for (fa = 0; fa < (*il)[F_DISRES].size(); fa += 3)
144 type = (*il)[F_DISRES].iatoms[fa];
147 npair = mtop->ffparams.iparams[type].disres.npair;
150 dd->nres += (ir->eDisre == edrEnsemble ? 1 : nmol);
151 dd->npair += nmol*npair;
154 type_min = std::min(type_min, type);
155 type_max = std::max(type_max, type);
160 if (cr && PAR(cr) && ir->nstdisreout > 0)
162 /* With DD we currently only have local pair information available */
163 gmx_fatal(FARGS, "With MPI parallelization distance-restraint pair output is not supported. Use nstdisreout=0 or use OpenMP parallelization on a single node.");
166 /* For communicating and/or reducing (sums of) r^-6 for pairs over threads
167 * we use multiple arrays in t_disresdata. We need to have unique indices
168 * for each restraint that work over threads and MPI ranks. To this end
169 * we use the type index. These should all be in one block and there should
170 * be dd->nres types, but we check for this here.
171 * This setup currently does not allow for multiple copies of the same
172 * molecule without ensemble averaging, this is check for above.
174 GMX_RELEASE_ASSERT(type_max - type_min + 1 == dd->nres, "All distance restraint parameter entries in the topology should be consecutive");
176 dd->type_min = type_min;
178 snew(dd->rt, dd->npair);
180 if (dd->dr_tau != 0.0)
182 GMX_RELEASE_ASSERT(state != nullptr, "We need a valid state when using time-averaged distance restraints");
185 /* Set the "history lack" factor to 1 */
186 state->flags |= (1<<estDISRE_INITF);
187 hist->disre_initf = 1.0;
188 /* Allocate space for the r^-3 time averages */
189 state->flags |= (1<<estDISRE_RM3TAV);
190 hist->ndisrepairs = dd->npair;
191 snew(hist->disre_rm3tav, hist->ndisrepairs);
193 /* Allocate space for a copy of rm3tav,
194 * so we can call do_force without modifying the state.
196 snew(dd->rm3tav, dd->npair);
198 /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
199 * averaged over the processors in one call (in calc_disre_R_6)
201 snew(dd->Rt_6, 2*dd->nres);
202 dd->Rtav_6 = &(dd->Rt_6[dd->nres]);
204 ptr = getenv("GMX_DISRE_ENSEMBLE_SIZE");
205 if (cr && ms != nullptr && ptr != nullptr && !bIsREMD)
209 sscanf(ptr, "%d", &dd->nsystems);
212 fprintf(fplog, "Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n", dd->nsystems);
214 /* This check is only valid on MASTER(cr), so probably
215 * ensemble-averaged distance restraints are broken on more
216 * than one processor per simulation system. */
219 check_multi_int(fplog, ms, dd->nsystems,
220 "the number of systems per ensemble",
223 gmx_bcast_sim(sizeof(int), &dd->nsystems, cr);
225 /* We use to allow any value of nsystems which was a divisor
226 * of ms->nsim. But this required an extra communicator which
227 * was stored in t_fcdata. This pulled in mpi.h in nearly all C files.
229 if (!(ms->nsim == 1 || ms->nsim == dd->nsystems))
231 gmx_fatal(FARGS, "GMX_DISRE_ENSEMBLE_SIZE (%d) is not equal to 1 or the number of systems (option -multidir) %d", dd->nsystems, ms->nsim);
235 fprintf(fplog, "Our ensemble consists of systems:");
236 for (int i = 0; i < dd->nsystems; i++)
238 fprintf(fplog, " %d",
239 (ms->sim/dd->nsystems)*dd->nsystems+i);
241 fprintf(fplog, "\n");
250 if (dd->nsystems == 1)
252 dd->Rtl_6 = dd->Rt_6;
256 snew(dd->Rtl_6, dd->nres);
263 fprintf(fplog, "There are %d distance restraints involving %d atom pairs\n", dd->nres, dd->npair);
265 /* Have to avoid g_disre de-referencing cr blindly, mdrun not
266 * doing consistency checks for ensemble-averaged distance
267 * restraints when that's not happening, and only doing those
268 * checks from appropriate processes (since check_multi_int is
269 * too broken to check whether the communication will
271 if (cr && ms && dd->nsystems > 1 && MASTER(cr))
273 check_multi_int(fplog, ms, fcd->disres.nres,
274 "the number of distance restraints",
277 please_cite(fplog, "Tropp80a");
278 please_cite(fplog, "Torda89a");
282 void calc_disres_R_6(const t_commrec *cr,
283 const gmx_multisim_t *ms,
284 int nfa, const t_iatom forceatoms[],
285 const rvec x[], const t_pbc *pbc,
286 t_fcdata *fcd, history_t *hist)
289 real *rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6;
291 real ETerm, ETerm1, cf1 = 0, cf2 = 0;
295 bTav = (dd->dr_tau != 0);
306 /* scaling factor to smoothly turn on the restraint forces *
307 * when using time averaging */
308 dd->exp_min_t_tau = hist->disre_initf*ETerm;
310 cf1 = dd->exp_min_t_tau;
311 cf2 = 1.0/(1.0 - dd->exp_min_t_tau);
314 for (int res = 0; res < dd->nres; res++)
320 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
321 * the total number of atoms pairs is nfa/3 */
322 for (int fa = 0; fa < nfa; fa += 3)
324 int type = forceatoms[fa];
325 int res = type - dd->type_min;
327 int ai = forceatoms[fa+1];
328 int aj = forceatoms[fa+2];
332 pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
336 rvec_sub(x[ai], x[aj], dx);
338 real rt2 = iprod(dx, dx);
339 real rt_1 = gmx::invsqrt(rt2);
340 real rt_3 = rt_1*rt_1*rt_1;
345 /* Here we update rm3tav in t_fcdata using the data
347 * Thus the results stay correct when this routine
348 * is called multiple times.
350 rm3tav[pair] = cf2*((ETerm - cf1)*hist->disre_rm3tav[pair] +
358 /* Currently divide_bondeds_over_threads() ensures that pairs within
359 * the same restraint get assigned to the same thread, so we could
360 * run this loop thread-parallel.
362 Rt_6[res] += rt_3*rt_3;
363 Rtav_6[res] += rm3tav[pair]*rm3tav[pair];
366 /* NOTE: Rt_6 and Rtav_6 are stored consecutively in memory */
367 if (cr && DOMAINDECOMP(cr))
369 gmx_sum(2*dd->nres, dd->Rt_6, cr);
372 if (fcd->disres.nsystems > 1)
374 real invn = 1.0/dd->nsystems;
376 for (int res = 0; res < dd->nres; res++)
378 Rtl_6[res] = Rt_6[res];
383 GMX_ASSERT(cr != nullptr && ms != nullptr, "We need multisim with nsystems>1");
384 gmx_sum_sim(2*dd->nres, dd->Rt_6, ms);
386 if (DOMAINDECOMP(cr))
388 gmx_bcast(2*dd->nres, dd->Rt_6, cr);
392 /* Store the base forceatoms pointer, so we can re-calculate the pair
393 * index in ta_disres() for indexing pair data in t_disresdata when
394 * using thread parallelization.
396 dd->forceatomsStart = forceatoms;
401 real ta_disres(int nfa, const t_iatom forceatoms[], const t_iparams ip[],
402 const rvec x[], rvec4 f[], rvec fshift[],
403 const t_pbc *pbc, const t_graph *g,
404 real gmx_unused lambda, real gmx_unused *dvdlambda,
405 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
406 int gmx_unused *global_atom_index)
408 const real seven_three = 7.0/3.0;
412 real smooth_fc, Rt, Rtav, rt2, *Rtl_6, *Rt_6, *Rtav_6;
413 real k0, f_scal = 0, fmax_scal, fk_scal, fij;
414 real tav_viol, instant_viol, mixed_viol, violtot, vtot;
415 real tav_viol_Rtav7, instant_viol_Rtav7;
417 gmx_bool bConservative, bMixed, bViolation;
424 dr_weighting = dd->dr_weighting;
425 dr_bMixed = dd->dr_bMixed;
430 tav_viol = instant_viol = mixed_viol = tav_viol_Rtav7 = instant_viol_Rtav7 = 0;
432 smooth_fc = dd->dr_fc;
435 /* scaling factor to smoothly turn on the restraint forces *
436 * when using time averaging */
437 smooth_fc *= (1.0 - dd->exp_min_t_tau);
443 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
444 * the total number of atoms pairs is nfa/3 */
445 int faOffset = static_cast<int>(forceatoms - dd->forceatomsStart);
446 for (int fa = 0; fa < nfa; fa += 3)
448 int type = forceatoms[fa];
449 int npair = ip[type].disres.npair;
450 up1 = ip[type].disres.up1;
451 up2 = ip[type].disres.up2;
452 low = ip[type].disres.low;
453 k0 = smooth_fc*ip[type].disres.kfac;
455 int res = type - dd->type_min;
457 /* save some flops when there is only one pair */
458 if (ip[type].disres.type != 2)
460 bConservative = (dr_weighting == edrwConservative) && (npair > 1);
462 Rt = gmx::invsixthroot(Rt_6[res]);
463 Rtav = gmx::invsixthroot(Rtav_6[res]);
467 /* When rtype=2 use instantaneous not ensemble averaged distance */
468 bConservative = (npair > 1);
470 Rt = gmx::invsixthroot(Rtl_6[res]);
477 tav_viol = Rtav - up1;
482 tav_viol = Rtav - low;
491 /* Add 1/npair energy and violation for each of the npair pairs */
492 real pairFac = 1/static_cast<real>(npair);
495 * there is no real potential when time averaging is applied
497 vtot += 0.5*k0*gmx::square(tav_viol)*pairFac;
500 f_scal = -k0*tav_viol;
501 violtot += fabs(tav_viol)*pairFac;
509 instant_viol = Rt - up1;
520 instant_viol = Rt - low;
533 mixed_viol = std::sqrt(tav_viol*instant_viol);
534 f_scal = -k0*mixed_viol;
535 violtot += mixed_viol*pairFac;
542 fmax_scal = -k0*(up2-up1);
543 /* Correct the force for the number of restraints */
546 f_scal = std::max(f_scal, fmax_scal);
549 f_scal *= Rtav/Rtav_6[res];
553 f_scal /= 2*mixed_viol;
554 tav_viol_Rtav7 = tav_viol*Rtav/Rtav_6[res];
555 instant_viol_Rtav7 = instant_viol*Rt/Rt_6[res];
561 f_scal = std::max(f_scal, fmax_scal);
564 /* Exert the force ... */
566 int pair = (faOffset + fa)/3;
567 int ai = forceatoms[fa+1];
568 int aj = forceatoms[fa+2];
572 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
576 rvec_sub(x[ai], x[aj], dx);
580 weight_rt_1 = gmx::invsqrt(rt2);
586 weight_rt_1 *= std::pow(dd->rm3tav[pair], seven_three);
590 weight_rt_1 *= tav_viol_Rtav7*std::pow(dd->rm3tav[pair], seven_three)+
591 instant_viol_Rtav7/(dd->rt[pair]*gmx::power6(dd->rt[pair]));
595 fk_scal = f_scal*weight_rt_1;
599 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
603 for (int m = 0; m < DIM; m++)
609 fshift[ki][m] += fij;
610 fshift[CENTRAL][m] -= fij;
616 dd->sumviol += violtot;
622 void update_disres_history(const t_fcdata *fcd, history_t *hist)
624 const t_disresdata *dd;
630 /* Copy the new time averages that have been calculated
631 * in calc_disres_R_6.
633 hist->disre_initf = dd->exp_min_t_tau;
634 for (pair = 0; pair < dd->npair; pair++)
636 hist->disre_rm3tav[pair] = dd->rm3tav[pair];