2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * \brief This file defines low-level functions necessary for
38 * computing energies and forces for position restraints.
40 * \author Mark Abraham <mark.j.abraham@gmail.com>
42 * \ingroup module_listed_forces
47 #include "position_restraints.h"
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/mdtypes/enerdata.h"
56 #include "gromacs/mdtypes/forceoutput.h"
57 #include "gromacs/mdtypes/forcerec.h"
58 #include "gromacs/mdtypes/inputrec.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/timing/wallcycle.h"
62 #include "gromacs/topology/idef.h"
63 #include "gromacs/utility/basedefinitions.h"
64 #include "gromacs/utility/fatalerror.h"
71 /*! \brief returns dx, rdist, and dpdl for functions posres() and fbposres()
73 void posres_dx(const rvec x,
87 real posA, posB, L1, ref = 0.;
92 for (m = 0; m < DIM; m++)
98 switch (refcoord_scaling)
102 rdist[m] = L1 * posA + lambda * posB;
103 dpdl[m] = posB - posA;
106 /* Box relative coordinates are stored for dimensions with pbc */
107 posA *= pbc->box[m][m];
108 posB *= pbc->box[m][m];
109 assert(npbcdim <= DIM);
110 for (d = m + 1; d < npbcdim && d < DIM; d++)
112 posA += pos0A[d] * pbc->box[d][m];
113 posB += pos0B[d] * pbc->box[d][m];
115 ref = L1 * posA + lambda * posB;
117 dpdl[m] = posB - posA;
120 ref = L1 * comA_sc[m] + lambda * comB_sc[m];
121 rdist[m] = L1 * posA + lambda * posB;
122 dpdl[m] = comB_sc[m] - comA_sc[m] + posB - posA;
124 default: gmx_fatal(FARGS, "No such scaling method implemented");
129 ref = L1 * posA + lambda * posB;
131 dpdl[m] = posB - posA;
134 /* We do pbc_dx with ref+rdist,
135 * since with only ref we can be up to half a box vector wrong.
137 pos[m] = ref + rdist[m];
142 pbc_dx(pbc, x, pos, dx);
146 rvec_sub(x, pos, dx);
150 /*! \brief Computes forces and potential for flat-bottom cylindrical restraints.
151 * Returns the flat-bottom potential. */
152 real do_fbposres_cylinder(int fbdim, rvec fm, rvec dx, real rfb, real kk, gmx_bool bInvert)
155 real dr, dr2, invdr, v, rfb2;
158 rfb2 = gmx::square(rfb);
161 for (d = 0; d < DIM; d++)
165 dr2 += gmx::square(dx[d]);
169 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
173 v = 0.5 * kk * gmx::square(dr - rfb);
174 for (d = 0; d < DIM; d++)
178 fm[d] = -kk * (dr - rfb) * dx[d] * invdr; /* Force pointing to the center */
186 /*! \brief Compute energies and forces for flat-bottomed position restraints
188 * Returns the flat-bottomed potential. Same PBC treatment as in
189 * normal position restraints */
190 real fbposres(int nbonds,
191 const t_iatom forceatoms[],
192 const t_iparams forceparams[],
194 gmx::ForceWithVirial* forceWithVirial,
196 int refcoord_scaling,
199 /* compute flat-bottomed positions restraints */
201 int i, ai, m, d, type, npbcdim = 0, fbdim;
204 real dr, dr2, rfb, rfb2, fact;
205 rvec com_sc, rdist, dx, dpdl, fm;
208 npbcdim = ePBC2npbcdim(ePBC);
209 GMX_ASSERT((ePBC == epbcNONE) == (npbcdim == 0), "");
210 if (refcoord_scaling == erscCOM)
213 for (m = 0; m < npbcdim; m++)
215 assert(npbcdim <= DIM);
216 for (d = m; d < npbcdim; d++)
218 com_sc[m] += com[d] * pbc->box[d][m];
223 rvec* f = as_rvec_array(forceWithVirial->force_.data());
226 for (i = 0; (i < nbonds);)
228 type = forceatoms[i++];
229 ai = forceatoms[i++];
230 pr = &forceparams[type];
232 /* same calculation as for normal posres, but with identical A and B states, and lambda==0 */
233 posres_dx(x[ai], forceparams[type].fbposres.pos0, forceparams[type].fbposres.pos0, com_sc,
234 com_sc, 0.0, pbc, refcoord_scaling, npbcdim, dx, rdist, dpdl);
240 rfb = pr->fbposres.r;
241 rfb2 = gmx::square(rfb);
243 /* with rfb<0, push particle out of the sphere/cylinder/layer */
251 switch (pr->fbposres.geom)
253 case efbposresSPHERE:
254 /* spherical flat-bottom posres */
256 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
259 v = 0.5 * kk * gmx::square(dr - rfb);
260 fact = -kk * (dr - rfb) / dr; /* Force pointing to the center pos0 */
264 case efbposresCYLINDERX:
265 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
267 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
269 case efbposresCYLINDERY:
270 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
272 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
274 case efbposresCYLINDER:
275 /* equivalent to efbposresCYLINDERZ for backwards compatibility */
276 case efbposresCYLINDERZ:
277 /* cylindrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */
279 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
281 case efbposresX: /* fbdim=XX */
282 case efbposresY: /* fbdim=YY */
283 case efbposresZ: /* fbdim=ZZ */
284 /* 1D flat-bottom potential */
285 fbdim = pr->fbposres.geom - efbposresX;
287 if ((dr > rfb && !bInvert) || (0 < dr && dr < rfb && bInvert))
289 v = 0.5 * kk * gmx::square(dr - rfb);
290 fm[fbdim] = -kk * (dr - rfb);
292 else if ((dr < (-rfb) && !bInvert) || ((-rfb) < dr && dr < 0 && bInvert))
294 v = 0.5 * kk * gmx::square(dr + rfb);
295 fm[fbdim] = -kk * (dr + rfb);
302 for (m = 0; (m < DIM); m++)
305 /* Here we correct for the pbc_dx which included rdist */
306 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm[m];
310 forceWithVirial->addVirialContribution(virial);
316 /*! \brief Compute energies and forces, when requested, for position restraints
318 * Note that position restraints require a different pbc treatment
319 * from other bondeds */
320 template<bool computeForce>
321 real posres(int nbonds,
322 const t_iatom forceatoms[],
323 const t_iparams forceparams[],
325 gmx::ForceWithVirial* forceWithVirial,
326 const struct t_pbc* pbc,
329 int refcoord_scaling,
334 int i, ai, m, d, type, npbcdim = 0;
337 rvec comA_sc, comB_sc, rdist, dpdl, dx;
339 npbcdim = ePBC2npbcdim(ePBC);
340 GMX_ASSERT((ePBC == epbcNONE) == (npbcdim == 0), "");
341 if (refcoord_scaling == erscCOM)
345 for (m = 0; m < npbcdim; m++)
347 assert(npbcdim <= DIM);
348 for (d = m; d < npbcdim; d++)
350 comA_sc[m] += comA[d] * pbc->box[d][m];
351 comB_sc[m] += comB[d] * pbc->box[d][m];
356 const real L1 = 1.0 - lambda;
361 GMX_ASSERT(forceWithVirial != nullptr, "When forces are requested we need a force object");
362 f = as_rvec_array(forceWithVirial->force_.data());
365 /* Use intermediate virial buffer to reduce reduction rounding errors */
367 for (i = 0; (i < nbonds);)
369 type = forceatoms[i++];
370 ai = forceatoms[i++];
371 pr = &forceparams[type];
373 /* return dx, rdist, and dpdl */
374 posres_dx(x[ai], forceparams[type].posres.pos0A, forceparams[type].posres.pos0B, comA_sc,
375 comB_sc, lambda, pbc, refcoord_scaling, npbcdim, dx, rdist, dpdl);
377 for (m = 0; (m < DIM); m++)
379 kk = L1 * pr->posres.fcA[m] + lambda * pr->posres.fcB[m];
381 vtot += 0.5 * kk * dx[m] * dx[m];
382 *dvdlambda += 0.5 * (pr->posres.fcB[m] - pr->posres.fcA[m]) * dx[m] * dx[m] + fm * dpdl[m];
384 /* Here we correct for the pbc_dx which included rdist */
388 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm;
395 forceWithVirial->addVirialContribution(virial);
403 void posres_wrapper(t_nrnb* nrnb,
405 const struct t_pbc* pbc,
407 gmx_enerdata_t* enerd,
409 const t_forcerec* fr,
410 gmx::ForceWithVirial* forceWithVirial)
415 v = posres<true>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x,
416 forceWithVirial, fr->ePBC == epbcNONE ? nullptr : pbc, lambda[efptRESTRAINT],
417 &dvdl, fr->rc_scaling, fr->ePBC, fr->posres_com, fr->posres_comB);
418 enerd->term[F_POSRES] += v;
419 /* If just the force constant changes, the FEP term is linear,
420 * but if k changes, it is not.
422 enerd->dvdl_nonlin[efptRESTRAINT] += dvdl;
423 inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef->il[F_POSRES].nr, 2));
426 void posres_wrapper_lambda(struct gmx_wallcycle* wcycle,
427 const t_lambda* fepvals,
429 const struct t_pbc* pbc,
431 gmx_enerdata_t* enerd,
433 const t_forcerec* fr)
437 if (0 == idef->il[F_POSRES].nr)
442 wallcycle_sub_start_nocount(wcycle, ewcsRESTRAINTS);
443 for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++)
445 real dvdl_dum = 0, lambda_dum;
447 lambda_dum = (i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i - 1]);
448 v = posres<false>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x,
449 nullptr, fr->ePBC == epbcNONE ? nullptr : pbc, lambda_dum, &dvdl_dum,
450 fr->rc_scaling, fr->ePBC, fr->posres_com, fr->posres_comB);
451 enerd->enerpart_lambda[i] += v;
453 wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
456 /*! \brief Helper function that wraps calls to fbposres for
457 free-energy perturbation */
458 void fbposres_wrapper(t_nrnb* nrnb,
460 const struct t_pbc* pbc,
462 gmx_enerdata_t* enerd,
463 const t_forcerec* fr,
464 gmx::ForceWithVirial* forceWithVirial)
468 v = fbposres(idef->il[F_FBPOSRES].nr, idef->il[F_FBPOSRES].iatoms, idef->iparams_fbposres, x,
469 forceWithVirial, fr->ePBC == epbcNONE ? nullptr : pbc, fr->rc_scaling, fr->ePBC,
471 enerd->term[F_FBPOSRES] += v;
472 inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef->il[F_FBPOSRES].nr, 2));