2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
5 * Copyright (c) 2019,2020, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * \brief This file defines low-level functions necessary for
39 * computing energies and forces for position restraints.
41 * \author Mark Abraham <mark.j.abraham@gmail.com>
43 * \ingroup module_listed_forces
48 #include "position_restraints.h"
53 #include "gromacs/gmxlib/nrnb.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/mdtypes/enerdata.h"
57 #include "gromacs/mdtypes/forceoutput.h"
58 #include "gromacs/mdtypes/forcerec.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/timing/wallcycle.h"
63 #include "gromacs/topology/idef.h"
64 #include "gromacs/utility/basedefinitions.h"
65 #include "gromacs/utility/fatalerror.h"
72 /*! \brief returns dx, rdist, and dpdl for functions posres() and fbposres()
74 void posres_dx(const rvec x,
88 real posA, posB, L1, ref = 0.;
93 for (m = 0; m < DIM; m++)
99 switch (refcoord_scaling)
103 rdist[m] = L1 * posA + lambda * posB;
104 dpdl[m] = posB - posA;
107 /* Box relative coordinates are stored for dimensions with pbc */
108 posA *= pbc->box[m][m];
109 posB *= pbc->box[m][m];
110 assert(npbcdim <= DIM);
111 for (d = m + 1; d < npbcdim && d < DIM; d++)
113 posA += pos0A[d] * pbc->box[d][m];
114 posB += pos0B[d] * pbc->box[d][m];
116 ref = L1 * posA + lambda * posB;
118 dpdl[m] = posB - posA;
121 ref = L1 * comA_sc[m] + lambda * comB_sc[m];
122 rdist[m] = L1 * posA + lambda * posB;
123 dpdl[m] = comB_sc[m] - comA_sc[m] + posB - posA;
125 default: gmx_fatal(FARGS, "No such scaling method implemented");
130 ref = L1 * posA + lambda * posB;
132 dpdl[m] = posB - posA;
135 /* We do pbc_dx with ref+rdist,
136 * since with only ref we can be up to half a box vector wrong.
138 pos[m] = ref + rdist[m];
143 pbc_dx(pbc, x, pos, dx);
147 rvec_sub(x, pos, dx);
151 /*! \brief Computes forces and potential for flat-bottom cylindrical restraints.
152 * Returns the flat-bottom potential. */
153 real do_fbposres_cylinder(int fbdim, rvec fm, rvec dx, real rfb, real kk, gmx_bool bInvert)
156 real dr, dr2, invdr, v, rfb2;
159 rfb2 = gmx::square(rfb);
162 for (d = 0; d < DIM; d++)
166 dr2 += gmx::square(dx[d]);
170 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
174 v = 0.5 * kk * gmx::square(dr - rfb);
175 for (d = 0; d < DIM; d++)
179 fm[d] = -kk * (dr - rfb) * dx[d] * invdr; /* Force pointing to the center */
187 /*! \brief Compute energies and forces for flat-bottomed position restraints
189 * Returns the flat-bottomed potential. Same PBC treatment as in
190 * normal position restraints */
191 real fbposres(int nbonds,
192 const t_iatom forceatoms[],
193 const t_iparams forceparams[],
195 gmx::ForceWithVirial* forceWithVirial,
197 int refcoord_scaling,
200 /* compute flat-bottomed positions restraints */
202 int i, ai, m, d, type, npbcdim = 0, fbdim;
205 real dr, dr2, rfb, rfb2, fact;
206 rvec com_sc, rdist, dx, dpdl, fm;
209 npbcdim = numPbcDimensions(pbcType);
210 GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
211 if (refcoord_scaling == erscCOM)
214 for (m = 0; m < npbcdim; m++)
216 assert(npbcdim <= DIM);
217 for (d = m; d < npbcdim; d++)
219 com_sc[m] += com[d] * pbc->box[d][m];
224 rvec* f = as_rvec_array(forceWithVirial->force_.data());
227 for (i = 0; (i < nbonds);)
229 type = forceatoms[i++];
230 ai = forceatoms[i++];
231 pr = &forceparams[type];
233 /* same calculation as for normal posres, but with identical A and B states, and lambda==0 */
234 posres_dx(x[ai], forceparams[type].fbposres.pos0, forceparams[type].fbposres.pos0, com_sc,
235 com_sc, 0.0, pbc, refcoord_scaling, npbcdim, dx, rdist, dpdl);
241 rfb = pr->fbposres.r;
242 rfb2 = gmx::square(rfb);
244 /* with rfb<0, push particle out of the sphere/cylinder/layer */
252 switch (pr->fbposres.geom)
254 case efbposresSPHERE:
255 /* spherical flat-bottom posres */
257 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
260 v = 0.5 * kk * gmx::square(dr - rfb);
261 fact = -kk * (dr - rfb) / dr; /* Force pointing to the center pos0 */
265 case efbposresCYLINDERX:
266 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
268 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
270 case efbposresCYLINDERY:
271 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
273 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
275 case efbposresCYLINDER:
276 /* equivalent to efbposresCYLINDERZ for backwards compatibility */
277 case efbposresCYLINDERZ:
278 /* cylindrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */
280 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
282 case efbposresX: /* fbdim=XX */
283 case efbposresY: /* fbdim=YY */
284 case efbposresZ: /* fbdim=ZZ */
285 /* 1D flat-bottom potential */
286 fbdim = pr->fbposres.geom - efbposresX;
288 if ((dr > rfb && !bInvert) || (0 < dr && dr < rfb && bInvert))
290 v = 0.5 * kk * gmx::square(dr - rfb);
291 fm[fbdim] = -kk * (dr - rfb);
293 else if ((dr < (-rfb) && !bInvert) || ((-rfb) < dr && dr < 0 && bInvert))
295 v = 0.5 * kk * gmx::square(dr + rfb);
296 fm[fbdim] = -kk * (dr + rfb);
303 for (m = 0; (m < DIM); m++)
306 /* Here we correct for the pbc_dx which included rdist */
307 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm[m];
311 forceWithVirial->addVirialContribution(virial);
317 /*! \brief Compute energies and forces, when requested, for position restraints
319 * Note that position restraints require a different pbc treatment
320 * from other bondeds */
321 template<bool computeForce>
322 real posres(int nbonds,
323 const t_iatom forceatoms[],
324 const t_iparams forceparams[],
326 gmx::ForceWithVirial* forceWithVirial,
327 const struct t_pbc* pbc,
330 int refcoord_scaling,
335 int i, ai, m, d, type, npbcdim = 0;
338 rvec comA_sc, comB_sc, rdist, dpdl, dx;
340 npbcdim = numPbcDimensions(pbcType);
341 GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
342 if (refcoord_scaling == erscCOM)
346 for (m = 0; m < npbcdim; m++)
348 assert(npbcdim <= DIM);
349 for (d = m; d < npbcdim; d++)
351 comA_sc[m] += comA[d] * pbc->box[d][m];
352 comB_sc[m] += comB[d] * pbc->box[d][m];
357 const real L1 = 1.0 - lambda;
362 GMX_ASSERT(forceWithVirial != nullptr, "When forces are requested we need a force object");
363 f = as_rvec_array(forceWithVirial->force_.data());
366 /* Use intermediate virial buffer to reduce reduction rounding errors */
368 for (i = 0; (i < nbonds);)
370 type = forceatoms[i++];
371 ai = forceatoms[i++];
372 pr = &forceparams[type];
374 /* return dx, rdist, and dpdl */
375 posres_dx(x[ai], forceparams[type].posres.pos0A, forceparams[type].posres.pos0B, comA_sc,
376 comB_sc, lambda, pbc, refcoord_scaling, npbcdim, dx, rdist, dpdl);
378 for (m = 0; (m < DIM); m++)
380 kk = L1 * pr->posres.fcA[m] + lambda * pr->posres.fcB[m];
382 vtot += 0.5 * kk * dx[m] * dx[m];
383 *dvdlambda += 0.5 * (pr->posres.fcB[m] - pr->posres.fcA[m]) * dx[m] * dx[m] + fm * dpdl[m];
385 /* Here we correct for the pbc_dx which included rdist */
389 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm;
396 forceWithVirial->addVirialContribution(virial);
404 void posres_wrapper(t_nrnb* nrnb,
406 const struct t_pbc* pbc,
408 gmx_enerdata_t* enerd,
410 const t_forcerec* fr,
411 gmx::ForceWithVirial* forceWithVirial)
416 v = posres<true>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x,
417 forceWithVirial, fr->pbcType == PbcType::No ? nullptr : pbc, lambda[efptRESTRAINT],
418 &dvdl, fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB);
419 enerd->term[F_POSRES] += v;
420 /* If just the force constant changes, the FEP term is linear,
421 * but if k changes, it is not.
423 enerd->dvdl_nonlin[efptRESTRAINT] += dvdl;
424 inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef->il[F_POSRES].nr, 2));
427 void posres_wrapper_lambda(struct gmx_wallcycle* wcycle,
428 const t_lambda* fepvals,
430 const struct t_pbc* pbc,
432 gmx_enerdata_t* enerd,
434 const t_forcerec* fr)
438 if (0 == idef->il[F_POSRES].nr)
443 wallcycle_sub_start_nocount(wcycle, ewcsRESTRAINTS);
444 for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++)
448 const real lambda_dum =
449 (i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i - 1]);
450 v = posres<false>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x,
451 nullptr, fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, &dvdl,
452 fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB);
453 enerd->enerpart_lambda[i] += v;
454 enerd->dhdlLambda[i] += dvdl;
456 wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
459 /*! \brief Helper function that wraps calls to fbposres for
460 free-energy perturbation */
461 void fbposres_wrapper(t_nrnb* nrnb,
463 const struct t_pbc* pbc,
465 gmx_enerdata_t* enerd,
466 const t_forcerec* fr,
467 gmx::ForceWithVirial* forceWithVirial)
471 v = fbposres(idef->il[F_FBPOSRES].nr, idef->il[F_FBPOSRES].iatoms, idef->iparams_fbposres, x,
472 forceWithVirial, fr->pbcType == PbcType::No ? nullptr : pbc, fr->rc_scaling,
473 fr->pbcType, fr->posres_com);
474 enerd->term[F_FBPOSRES] += v;
475 inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef->il[F_FBPOSRES].nr, 2));