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, 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/mdlib/main.h"
54 #include "gromacs/mdtypes/commrec.h"
55 #include "gromacs/mdtypes/fcdata.h"
56 #include "gromacs/mdtypes/inputrec.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/mdtypes/state.h"
59 #include "gromacs/pbcutil/ishift.h"
60 #include "gromacs/pbcutil/mshift.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/topology/mtop_util.h"
63 #include "gromacs/topology/topology.h"
64 #include "gromacs/utility/fatalerror.h"
65 #include "gromacs/utility/futil.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/pleasecite.h"
68 #include "gromacs/utility/smalloc.h"
70 void init_disres(FILE *fplog
, const gmx_mtop_t
*mtop
,
71 t_inputrec
*ir
, const t_commrec
*cr
,
72 t_fcdata
*fcd
, t_state
*state
, gmx_bool bIsREMD
)
74 int fa
, nmol
, npair
, np
;
77 gmx_mtop_ilistloop_t iloop
;
80 int type_min
, type_max
;
84 if (gmx_mtop_ftype_count(mtop
, F_DISRES
) == 0)
93 fprintf(fplog
, "Initializing the distance restraints\n");
96 dd
->dr_weighting
= ir
->eDisreWeighting
;
97 dd
->dr_fc
= ir
->dr_fc
;
98 if (EI_DYNAMICS(ir
->eI
))
100 dd
->dr_tau
= ir
->dr_tau
;
106 if (dd
->dr_tau
== 0.0)
108 dd
->dr_bMixed
= FALSE
;
113 /* We store the r^-6 time averages in an array that is indexed
114 * with the local disres iatom index, so this doesn't work with DD.
115 * Note that DD is not initialized yet here, so we check for PAR(cr),
116 * but there are probably also issues with e.g. NM MPI parallelization.
120 gmx_fatal(FARGS
, "Time-averaged distance restraints are not supported with MPI parallelization. You can use OpenMP parallelization on a single node.");
123 dd
->dr_bMixed
= ir
->bDisreMixed
;
124 dd
->ETerm
= std::exp(-(ir
->delta_t
/ir
->dr_tau
));
126 dd
->ETerm1
= 1.0 - dd
->ETerm
;
132 iloop
= gmx_mtop_ilistloop_init(mtop
);
133 while (gmx_mtop_ilistloop_next(iloop
, &il
, &nmol
))
135 if (nmol
> 1 && il
[F_DISRES
].nr
> 0 && ir
->eDisre
!= edrEnsemble
)
137 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.");
141 for (fa
= 0; fa
< il
[F_DISRES
].nr
; fa
+= 3)
145 type
= il
[F_DISRES
].iatoms
[fa
];
148 npair
= mtop
->ffparams
.iparams
[type
].disres
.npair
;
151 dd
->nres
+= (ir
->eDisre
== edrEnsemble
? 1 : nmol
);
152 dd
->npair
+= nmol
*npair
;
155 type_min
= std::min(type_min
, type
);
156 type_max
= std::max(type_max
, type
);
161 if (cr
&& PAR(cr
) && ir
->nstdisreout
> 0)
163 /* With DD we currently only have local pair information available */
164 gmx_fatal(FARGS
, "With MPI parallelization distance-restraint pair output is not supported. Use nstdisreout=0 or use OpenMP parallelization on a single node.");
167 /* For communicating and/or reducing (sums of) r^-6 for pairs over threads
168 * we use multiple arrays in t_disresdata. We need to have unique indices
169 * for each restraint that work over threads and MPI ranks. To this end
170 * we use the type index. These should all be in one block and there should
171 * be dd->nres types, but we check for this here.
172 * This setup currently does not allow for multiple copies of the same
173 * molecule without ensemble averaging, this is check for above.
175 GMX_RELEASE_ASSERT(type_max
- type_min
+ 1 == dd
->nres
, "All distance restraint parameter entries in the topology should be consecutive");
177 dd
->type_min
= type_min
;
179 snew(dd
->rt
, dd
->npair
);
181 if (dd
->dr_tau
!= 0.0)
184 /* Set the "history lack" factor to 1 */
185 state
->flags
|= (1<<estDISRE_INITF
);
186 hist
->disre_initf
= 1.0;
187 /* Allocate space for the r^-3 time averages */
188 state
->flags
|= (1<<estDISRE_RM3TAV
);
189 hist
->ndisrepairs
= dd
->npair
;
190 snew(hist
->disre_rm3tav
, hist
->ndisrepairs
);
192 /* Allocate space for a copy of rm3tav,
193 * so we can call do_force without modifying the state.
195 snew(dd
->rm3tav
, dd
->npair
);
197 /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
198 * averaged over the processors in one call (in calc_disre_R_6)
200 snew(dd
->Rt_6
, 2*dd
->nres
);
201 dd
->Rtav_6
= &(dd
->Rt_6
[dd
->nres
]);
203 ptr
= getenv("GMX_DISRE_ENSEMBLE_SIZE");
204 if (cr
&& cr
->ms
!= nullptr && ptr
!= nullptr && !bIsREMD
)
208 sscanf(ptr
, "%d", &dd
->nsystems
);
211 fprintf(fplog
, "Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n", dd
->nsystems
);
213 /* This check is only valid on MASTER(cr), so probably
214 * ensemble-averaged distance restraints are broken on more
215 * than one processor per simulation system. */
218 check_multi_int(fplog
, cr
->ms
, dd
->nsystems
,
219 "the number of systems per ensemble",
222 gmx_bcast_sim(sizeof(int), &dd
->nsystems
, cr
);
224 /* We use to allow any value of nsystems which was a divisor
225 * of ms->nsim. But this required an extra communicator which
226 * was stored in t_fcdata. This pulled in mpi.h in nearly all C files.
228 if (!(cr
->ms
->nsim
== 1 || cr
->ms
->nsim
== dd
->nsystems
))
230 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
);
234 fprintf(fplog
, "Our ensemble consists of systems:");
235 for (int i
= 0; i
< dd
->nsystems
; i
++)
237 fprintf(fplog
, " %d",
238 (cr
->ms
->sim
/dd
->nsystems
)*dd
->nsystems
+i
);
240 fprintf(fplog
, "\n");
249 if (dd
->nsystems
== 1)
251 dd
->Rtl_6
= dd
->Rt_6
;
255 snew(dd
->Rtl_6
, dd
->nres
);
262 fprintf(fplog
, "There are %d distance restraints involving %d atom pairs\n", dd
->nres
, dd
->npair
);
264 /* Have to avoid g_disre de-referencing cr blindly, mdrun not
265 * doing consistency checks for ensemble-averaged distance
266 * restraints when that's not happening, and only doing those
267 * checks from appropriate processes (since check_multi_int is
268 * too broken to check whether the communication will
270 if (cr
&& cr
->ms
&& dd
->nsystems
> 1 && MASTER(cr
))
272 check_multi_int(fplog
, cr
->ms
, fcd
->disres
.nres
,
273 "the number of distance restraints",
276 please_cite(fplog
, "Tropp80a");
277 please_cite(fplog
, "Torda89a");
281 void calc_disres_R_6(const t_commrec
*cr
,
282 int nfa
, const t_iatom forceatoms
[],
283 const rvec x
[], const t_pbc
*pbc
,
284 t_fcdata
*fcd
, history_t
*hist
)
287 real
*rt
, *rm3tav
, *Rtl_6
, *Rt_6
, *Rtav_6
;
289 real ETerm
, ETerm1
, cf1
= 0, cf2
= 0;
293 bTav
= (dd
->dr_tau
!= 0);
304 /* scaling factor to smoothly turn on the restraint forces *
305 * when using time averaging */
306 dd
->exp_min_t_tau
= hist
->disre_initf
*ETerm
;
308 cf1
= dd
->exp_min_t_tau
;
309 cf2
= 1.0/(1.0 - dd
->exp_min_t_tau
);
312 for (int res
= 0; res
< dd
->nres
; res
++)
318 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
319 * the total number of atoms pairs is nfa/3 */
320 for (int fa
= 0; fa
< nfa
; fa
+= 3)
322 int type
= forceatoms
[fa
];
323 int res
= type
- dd
->type_min
;
325 int ai
= forceatoms
[fa
+1];
326 int aj
= forceatoms
[fa
+2];
330 pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], dx
);
334 rvec_sub(x
[ai
], x
[aj
], dx
);
336 real rt2
= iprod(dx
, dx
);
337 real rt_1
= gmx::invsqrt(rt2
);
338 real rt_3
= rt_1
*rt_1
*rt_1
;
343 /* Here we update rm3tav in t_fcdata using the data
345 * Thus the results stay correct when this routine
346 * is called multiple times.
348 rm3tav
[pair
] = cf2
*((ETerm
- cf1
)*hist
->disre_rm3tav
[pair
] +
356 /* Currently divide_bondeds_over_threads() ensures that pairs within
357 * the same restraint get assigned to the same thread, so we could
358 * run this loop thread-parallel.
360 Rt_6
[res
] += rt_3
*rt_3
;
361 Rtav_6
[res
] += rm3tav
[pair
]*rm3tav
[pair
];
364 /* NOTE: Rt_6 and Rtav_6 are stored consecutively in memory */
365 if (cr
&& DOMAINDECOMP(cr
))
367 gmx_sum(2*dd
->nres
, dd
->Rt_6
, cr
);
370 if (fcd
->disres
.nsystems
> 1)
372 real invn
= 1.0/dd
->nsystems
;
374 for (int res
= 0; res
< dd
->nres
; res
++)
376 Rtl_6
[res
] = Rt_6
[res
];
381 GMX_ASSERT(cr
!= NULL
&& cr
->ms
!= NULL
, "We need multisim with nsystems>1");
382 gmx_sum_sim(2*dd
->nres
, dd
->Rt_6
, cr
->ms
);
384 if (DOMAINDECOMP(cr
))
386 gmx_bcast(2*dd
->nres
, dd
->Rt_6
, cr
);
390 /* Store the base forceatoms pointer, so we can re-calculate the pair
391 * index in ta_disres() for indexing pair data in t_disresdata when
392 * using thread parallelization.
394 dd
->forceatomsStart
= forceatoms
;
399 real
ta_disres(int nfa
, const t_iatom forceatoms
[], const t_iparams ip
[],
400 const rvec x
[], rvec4 f
[], rvec fshift
[],
401 const t_pbc
*pbc
, const t_graph
*g
,
402 real gmx_unused lambda
, real gmx_unused
*dvdlambda
,
403 const t_mdatoms gmx_unused
*md
, t_fcdata
*fcd
,
404 int gmx_unused
*global_atom_index
)
406 const real seven_three
= 7.0/3.0;
410 real smooth_fc
, Rt
, Rtav
, rt2
, *Rtl_6
, *Rt_6
, *Rtav_6
;
411 real k0
, f_scal
= 0, fmax_scal
, fk_scal
, fij
;
412 real tav_viol
, instant_viol
, mixed_viol
, violtot
, vtot
;
413 real tav_viol_Rtav7
, instant_viol_Rtav7
;
415 gmx_bool bConservative
, bMixed
, bViolation
;
422 dr_weighting
= dd
->dr_weighting
;
423 dr_bMixed
= dd
->dr_bMixed
;
428 tav_viol
= instant_viol
= mixed_viol
= tav_viol_Rtav7
= instant_viol_Rtav7
= 0;
430 smooth_fc
= dd
->dr_fc
;
433 /* scaling factor to smoothly turn on the restraint forces *
434 * when using time averaging */
435 smooth_fc
*= (1.0 - dd
->exp_min_t_tau
);
441 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
442 * the total number of atoms pairs is nfa/3 */
443 int faOffset
= static_cast<int>(forceatoms
- dd
->forceatomsStart
);
444 for (int fa
= 0; fa
< nfa
; fa
+= 3)
446 int type
= forceatoms
[fa
];
447 int npair
= ip
[type
].disres
.npair
;
448 up1
= ip
[type
].disres
.up1
;
449 up2
= ip
[type
].disres
.up2
;
450 low
= ip
[type
].disres
.low
;
451 k0
= smooth_fc
*ip
[type
].disres
.kfac
;
453 int res
= type
- dd
->type_min
;
455 /* save some flops when there is only one pair */
456 if (ip
[type
].disres
.type
!= 2)
458 bConservative
= (dr_weighting
== edrwConservative
) && (npair
> 1);
460 Rt
= gmx::invsixthroot(Rt_6
[res
]);
461 Rtav
= gmx::invsixthroot(Rtav_6
[res
]);
465 /* When rtype=2 use instantaneous not ensemble averaged distance */
466 bConservative
= (npair
> 1);
468 Rt
= gmx::invsixthroot(Rtl_6
[res
]);
475 tav_viol
= Rtav
- up1
;
480 tav_viol
= Rtav
- low
;
489 /* Add 1/npair energy and violation for each of the npair pairs */
490 real pairFac
= 1/static_cast<real
>(npair
);
493 * there is no real potential when time averaging is applied
495 vtot
+= 0.5*k0
*gmx::square(tav_viol
)*pairFac
;
498 f_scal
= -k0
*tav_viol
;
499 violtot
+= fabs(tav_viol
)*pairFac
;
507 instant_viol
= Rt
- up1
;
518 instant_viol
= Rt
- low
;
531 mixed_viol
= std::sqrt(tav_viol
*instant_viol
);
532 f_scal
= -k0
*mixed_viol
;
533 violtot
+= mixed_viol
*pairFac
;
540 fmax_scal
= -k0
*(up2
-up1
);
541 /* Correct the force for the number of restraints */
544 f_scal
= std::max(f_scal
, fmax_scal
);
547 f_scal
*= Rtav
/Rtav_6
[res
];
551 f_scal
/= 2*mixed_viol
;
552 tav_viol_Rtav7
= tav_viol
*Rtav
/Rtav_6
[res
];
553 instant_viol_Rtav7
= instant_viol
*Rt
/Rt_6
[res
];
559 f_scal
= std::max(f_scal
, fmax_scal
);
562 /* Exert the force ... */
564 int pair
= (faOffset
+ fa
)/3;
565 int ai
= forceatoms
[fa
+1];
566 int aj
= forceatoms
[fa
+2];
570 ki
= pbc_dx_aiuc(pbc
, x
[ai
], x
[aj
], dx
);
574 rvec_sub(x
[ai
], x
[aj
], dx
);
578 weight_rt_1
= gmx::invsqrt(rt2
);
584 weight_rt_1
*= std::pow(dd
->rm3tav
[pair
], seven_three
);
588 weight_rt_1
*= tav_viol_Rtav7
*std::pow(dd
->rm3tav
[pair
], seven_three
)+
589 instant_viol_Rtav7
/(dd
->rt
[pair
]*gmx::power6(dd
->rt
[pair
]));
593 fk_scal
= f_scal
*weight_rt_1
;
597 ivec_sub(SHIFT_IVEC(g
, ai
), SHIFT_IVEC(g
, aj
), dt
);
601 for (int m
= 0; m
< DIM
; m
++)
607 fshift
[ki
][m
] += fij
;
608 fshift
[CENTRAL
][m
] -= fij
;
614 dd
->sumviol
+= violtot
;
620 void update_disres_history(t_fcdata
*fcd
, history_t
*hist
)
628 /* Copy the new time averages that have been calculated
629 * in calc_disres_R_6.
631 hist
->disre_initf
= dd
->exp_min_t_tau
;
632 for (pair
= 0; pair
< dd
->npair
; pair
++)
634 hist
->disre_rm3tav
[pair
] = dd
->rm3tav
[pair
];