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, 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! */
45 #include "gromacs/legacyheaders/typedefs.h"
46 #include "gromacs/legacyheaders/types/commrec.h"
47 #include "gromacs/legacyheaders/macros.h"
48 #include "gromacs/utility/futil.h"
49 #include "gromacs/legacyheaders/copyrite.h"
50 #include "gromacs/legacyheaders/disre.h"
51 #include "gromacs/legacyheaders/main.h"
52 #include "gromacs/topology/mtop_util.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/pbcutil/ishift.h"
56 #include "gromacs/pbcutil/mshift.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/smalloc.h"
61 void init_disres(FILE *fplog, const gmx_mtop_t *mtop,
62 t_inputrec *ir, const t_commrec *cr,
63 t_fcdata *fcd, t_state *state, gmx_bool bIsREMD)
65 int fa, nmol, i, npair, np;
69 gmx_mtop_ilistloop_t iloop;
75 if (gmx_mtop_ftype_count(mtop, F_DISRES) == 0)
84 fprintf(fplog, "Initializing the distance restraints\n");
88 if (ir->eDisre == edrEnsemble)
90 gmx_fatal(FARGS, "Sorry, distance restraints with ensemble averaging over multiple molecules in one system are not functional in this version of GROMACS");
93 dd->dr_weighting = ir->eDisreWeighting;
94 dd->dr_fc = ir->dr_fc;
95 if (EI_DYNAMICS(ir->eI))
97 dd->dr_tau = ir->dr_tau;
103 if (dd->dr_tau == 0.0)
105 dd->dr_bMixed = FALSE;
110 dd->dr_bMixed = ir->bDisreMixed;
111 dd->ETerm = exp(-(ir->delta_t/ir->dr_tau));
113 dd->ETerm1 = 1.0 - dd->ETerm;
115 ip = mtop->ffparams.iparams;
119 iloop = gmx_mtop_ilistloop_init(mtop);
120 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
123 for (fa = 0; fa < il[F_DISRES].nr; fa += 3)
126 npair = mtop->ffparams.iparams[il[F_DISRES].iatoms[fa]].disres.npair;
129 dd->nres += (ir->eDisre == edrEnsemble ? 1 : nmol)*npair;
130 dd->npair += nmol*npair;
138 /* Temporary check, will be removed when disre is implemented with DD */
139 const char *notestr = "NOTE: atoms involved in distance restraints should be within the same domain. If this is not the case mdrun generates a fatal error. If you encounter this, use a single MPI rank (Verlet+OpenMP+GPUs work fine).";
143 fprintf(stderr, "\n%s\n\n", notestr);
147 fprintf(fplog, "%s\n", notestr);
150 if (dd->dr_tau != 0 || ir->eDisre == edrEnsemble || cr->ms != NULL ||
151 dd->nres != dd->npair)
153 gmx_fatal(FARGS, "Time or ensemble averaged or multiple pair distance restraints do not work (yet) with domain decomposition, use a single MPI rank%s", cr->ms ? " per simulation" : "");
155 if (ir->nstdisreout != 0)
159 fprintf(fplog, "\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n\n");
163 fprintf(stderr, "\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n");
169 snew(dd->rt, dd->npair);
171 if (dd->dr_tau != 0.0)
174 /* Set the "history lack" factor to 1 */
175 state->flags |= (1<<estDISRE_INITF);
176 hist->disre_initf = 1.0;
177 /* Allocate space for the r^-3 time averages */
178 state->flags |= (1<<estDISRE_RM3TAV);
179 hist->ndisrepairs = dd->npair;
180 snew(hist->disre_rm3tav, hist->ndisrepairs);
182 /* Allocate space for a copy of rm3tav,
183 * so we can call do_force without modifying the state.
185 snew(dd->rm3tav, dd->npair);
187 /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
188 * averaged over the processors in one call (in calc_disre_R_6)
190 snew(dd->Rt_6, 2*dd->nres);
191 dd->Rtav_6 = &(dd->Rt_6[dd->nres]);
193 ptr = getenv("GMX_DISRE_ENSEMBLE_SIZE");
194 if (cr && cr->ms != NULL && ptr != NULL && !bIsREMD)
198 sscanf(ptr, "%d", &dd->nsystems);
201 fprintf(fplog, "Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n", dd->nsystems);
203 /* This check is only valid on MASTER(cr), so probably
204 * ensemble-averaged distance restraints are broken on more
205 * than one processor per simulation system. */
208 check_multi_int(fplog, cr->ms, dd->nsystems,
209 "the number of systems per ensemble",
212 gmx_bcast_sim(sizeof(int), &dd->nsystems, cr);
214 /* We use to allow any value of nsystems which was a divisor
215 * of ms->nsim. But this required an extra communicator which
216 * was stored in t_fcdata. This pulled in mpi.h in nearly all C files.
218 if (!(cr->ms->nsim == 1 || cr->ms->nsim == dd->nsystems))
220 gmx_fatal(FARGS, "GMX_DISRE_ENSEMBLE_SIZE (%d) is not equal to 1 or the number of systems (option -multi) %d", dd->nsystems, cr->ms->nsim);
224 fprintf(fplog, "Our ensemble consists of systems:");
225 for (i = 0; i < dd->nsystems; i++)
227 fprintf(fplog, " %d",
228 (cr->ms->sim/dd->nsystems)*dd->nsystems+i);
230 fprintf(fplog, "\n");
232 snew(dd->Rtl_6, dd->nres);
238 dd->Rtl_6 = dd->Rt_6;
245 fprintf(fplog, "There are %d distance restraints involving %d atom pairs\n", dd->nres, dd->npair);
247 /* Have to avoid g_disre de-referencing cr blindly, mdrun not
248 * doing consistency checks for ensemble-averaged distance
249 * restraints when that's not happening, and only doing those
250 * checks from appropriate processes (since check_multi_int is
251 * too broken to check whether the communication will
253 if (cr && cr->ms && dd->nsystems > 1 && MASTER(cr))
255 check_multi_int(fplog, cr->ms, fcd->disres.nres,
256 "the number of distance restraints",
259 please_cite(fplog, "Tropp80a");
260 please_cite(fplog, "Torda89a");
264 void calc_disres_R_6(int nfa, const t_iatom forceatoms[], const t_iparams ip[],
265 const rvec x[], const t_pbc *pbc,
266 t_fcdata *fcd, history_t *hist)
269 int fa, res, i, pair, ki, kj, m;
272 real *rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6;
273 real rt_1, rt_3, rt2;
276 real ETerm, ETerm1, cf1 = 0, cf2 = 0, invn = 0;
280 bTav = (dd->dr_tau != 0);
291 /* scaling factor to smoothly turn on the restraint forces *
292 * when using time averaging */
293 dd->exp_min_t_tau = hist->disre_initf*ETerm;
295 cf1 = dd->exp_min_t_tau;
296 cf2 = 1.0/(1.0 - dd->exp_min_t_tau);
299 if (dd->nsystems > 1)
301 invn = 1.0/dd->nsystems;
304 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
305 * the total number of atoms pairs is nfa/3 */
310 type = forceatoms[fa];
311 npair = ip[type].disres.npair;
316 /* Loop over the atom pairs of 'this' restraint */
318 while (fa < nfa && np < npair)
321 ai = forceatoms[fa+1];
322 aj = forceatoms[fa+2];
326 pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
330 rvec_sub(x[ai], x[aj], dx);
333 rt_1 = gmx_invsqrt(rt2);
334 rt_3 = rt_1*rt_1*rt_1;
336 rt[pair] = sqrt(rt2);
339 /* Here we update rm3tav in t_fcdata using the data
341 * Thus the results stay correct when this routine
342 * is called multiple times.
344 rm3tav[pair] = cf2*((ETerm - cf1)*hist->disre_rm3tav[pair] +
352 Rt_6[res] += rt_3*rt_3;
353 Rtav_6[res] += rm3tav[pair]*rm3tav[pair];
358 if (dd->nsystems > 1)
360 Rtl_6[res] = Rt_6[res];
369 real ta_disres(int nfa, const t_iatom forceatoms[], const t_iparams ip[],
370 const rvec x[], rvec f[], rvec fshift[],
371 const t_pbc *pbc, const t_graph *g,
372 real gmx_unused lambda, real gmx_unused *dvdlambda,
373 const t_mdatoms gmx_unused *md, t_fcdata *fcd,
374 int gmx_unused *global_atom_index)
376 const real sixth = 1.0/6.0;
377 const real seven_three = 7.0/3.0;
380 int fa, res, npair, p, pair, ki = CENTRAL, m;
384 real smooth_fc, Rt, Rtav, rt2, *Rtl_6, *Rt_6, *Rtav_6;
385 real k0, f_scal = 0, fmax_scal, fk_scal, fij;
386 real tav_viol, instant_viol, mixed_viol, violtot, vtot;
387 real tav_viol_Rtav7, instant_viol_Rtav7;
389 gmx_bool bConservative, bMixed, bViolation;
396 dr_weighting = dd->dr_weighting;
397 dr_bMixed = dd->dr_bMixed;
402 tav_viol = instant_viol = mixed_viol = tav_viol_Rtav7 = instant_viol_Rtav7 = 0;
404 smooth_fc = dd->dr_fc;
407 /* scaling factor to smoothly turn on the restraint forces *
408 * when using time averaging */
409 smooth_fc *= (1.0 - dd->exp_min_t_tau);
415 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
416 * the total number of atoms pairs is nfa/3 */
421 type = forceatoms[fa];
422 /* Take action depending on restraint, calculate scalar force */
423 npair = ip[type].disres.npair;
424 up1 = ip[type].disres.up1;
425 up2 = ip[type].disres.up2;
426 low = ip[type].disres.low;
427 k0 = smooth_fc*ip[type].disres.kfac;
429 /* save some flops when there is only one pair */
430 if (ip[type].disres.type != 2)
432 bConservative = (dr_weighting == edrwConservative) && (npair > 1);
434 Rt = pow(Rt_6[res], -sixth);
435 Rtav = pow(Rtav_6[res], -sixth);
439 /* When rtype=2 use instantaneous not ensemble avereged distance */
440 bConservative = (npair > 1);
442 Rt = pow(Rtl_6[res], -sixth);
449 tav_viol = Rtav - up1;
454 tav_viol = Rtav - low;
464 * there is no real potential when time averaging is applied
466 vtot += 0.5*k0*sqr(tav_viol);
469 printf("vtot is inf: %f\n", vtot);
473 f_scal = -k0*tav_viol;
474 violtot += fabs(tav_viol);
482 instant_viol = Rt - up1;
493 instant_viol = Rt - low;
506 mixed_viol = sqrt(tav_viol*instant_viol);
507 f_scal = -k0*mixed_viol;
508 violtot += mixed_viol;
515 fmax_scal = -k0*(up2-up1);
516 /* Correct the force for the number of restraints */
519 f_scal = max(f_scal, fmax_scal);
522 f_scal *= Rtav/Rtav_6[res];
526 f_scal /= 2*mixed_viol;
527 tav_viol_Rtav7 = tav_viol*Rtav/Rtav_6[res];
528 instant_viol_Rtav7 = instant_viol*Rt/Rt_6[res];
533 f_scal /= (real)npair;
534 f_scal = max(f_scal, fmax_scal);
537 /* Exert the force ... */
539 /* Loop over the atom pairs of 'this' restraint */
540 for (p = 0; p < npair; p++)
543 ai = forceatoms[fa+1];
544 aj = forceatoms[fa+2];
548 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
552 rvec_sub(x[ai], x[aj], dx);
556 weight_rt_1 = gmx_invsqrt(rt2);
562 weight_rt_1 *= pow(dd->rm3tav[pair], seven_three);
566 weight_rt_1 *= tav_viol_Rtav7*pow(dd->rm3tav[pair], seven_three)+
567 instant_viol_Rtav7*pow(dd->rt[pair], -7);
571 fk_scal = f_scal*weight_rt_1;
575 ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
579 for (m = 0; m < DIM; m++)
585 fshift[ki][m] += fij;
586 fshift[CENTRAL][m] -= fij;
593 /* No violation so force and potential contributions */
599 dd->sumviol = violtot;
605 void update_disres_history(t_fcdata *fcd, history_t *hist)
613 /* Copy the new time averages that have been calculated
614 * in calc_disres_R_6.
616 hist->disre_initf = dd->exp_min_t_tau;
617 for (pair = 0; pair < dd->npair; pair++)
619 hist->disre_rm3tav[pair] = dd->rm3tav[pair];