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,2021, 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 "gromacs/utility/arrayref.h"
49 #include "position_restraints.h"
54 #include "gromacs/gmxlib/nrnb.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdtypes/enerdata.h"
58 #include "gromacs/mdtypes/forceoutput.h"
59 #include "gromacs/mdtypes/forcerec.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/timing/wallcycle.h"
64 #include "gromacs/topology/idef.h"
65 #include "gromacs/utility/basedefinitions.h"
66 #include "gromacs/utility/fatalerror.h"
73 /*! \brief returns dx, rdist, and dpdl for functions posres() and fbposres()
75 void posres_dx(const rvec x,
82 RefCoordScaling refcoord_scaling,
89 real posA, posB, L1, ref = 0.;
94 for (m = 0; m < DIM; m++)
100 switch (refcoord_scaling)
102 case RefCoordScaling::No:
104 rdist[m] = L1 * posA + lambda * posB;
105 dpdl[m] = posB - posA;
107 case RefCoordScaling::All:
108 /* Box relative coordinates are stored for dimensions with pbc */
109 posA *= pbc->box[m][m];
110 posB *= pbc->box[m][m];
111 assert(npbcdim <= DIM);
112 for (d = m + 1; d < npbcdim && d < DIM; d++)
114 posA += pos0A[d] * pbc->box[d][m];
115 posB += pos0B[d] * pbc->box[d][m];
117 ref = L1 * posA + lambda * posB;
119 dpdl[m] = posB - posA;
121 case RefCoordScaling::Com:
122 ref = L1 * comA_sc[m] + lambda * comB_sc[m];
123 rdist[m] = L1 * posA + lambda * posB;
124 dpdl[m] = comB_sc[m] - comA_sc[m] + posB - posA;
126 default: gmx_fatal(FARGS, "No such scaling method implemented");
131 ref = L1 * posA + lambda * posB;
133 dpdl[m] = posB - posA;
136 /* We do pbc_dx with ref+rdist,
137 * since with only ref we can be up to half a box vector wrong.
139 pos[m] = ref + rdist[m];
144 pbc_dx(pbc, x, pos, dx);
148 rvec_sub(x, pos, dx);
152 /*! \brief Computes forces and potential for flat-bottom cylindrical restraints.
153 * Returns the flat-bottom potential. */
154 real do_fbposres_cylinder(int fbdim, rvec fm, rvec dx, real rfb, real kk, gmx_bool bInvert)
157 real dr, dr2, invdr, v, rfb2;
160 rfb2 = gmx::square(rfb);
163 for (d = 0; d < DIM; d++)
167 dr2 += gmx::square(dx[d]);
171 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
175 v = 0.5 * kk * gmx::square(dr - rfb);
176 for (d = 0; d < DIM; d++)
180 fm[d] = -kk * (dr - rfb) * dx[d] * invdr; /* Force pointing to the center */
188 /*! \brief Compute energies and forces for flat-bottomed position restraints
190 * Returns the flat-bottomed potential. Same PBC treatment as in
191 * normal position restraints */
192 real fbposres(int nbonds,
193 const t_iatom forceatoms[],
194 const t_iparams forceparams[],
196 gmx::ForceWithVirial* forceWithVirial,
198 RefCoordScaling refcoord_scaling,
201 /* compute flat-bottomed positions restraints */
203 int i, ai, m, d, type, npbcdim = 0, fbdim;
206 real dr, dr2, rfb, rfb2, fact;
207 rvec com_sc, rdist, dx, dpdl, fm;
210 npbcdim = numPbcDimensions(pbcType);
211 GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
212 if (refcoord_scaling == RefCoordScaling::Com)
215 for (m = 0; m < npbcdim; m++)
217 assert(npbcdim <= DIM);
218 for (d = m; d < npbcdim; d++)
220 com_sc[m] += com[d] * pbc->box[d][m];
225 rvec* f = as_rvec_array(forceWithVirial->force_.data());
228 for (i = 0; (i < nbonds);)
230 type = forceatoms[i++];
231 ai = forceatoms[i++];
232 pr = &forceparams[type];
234 /* same calculation as for normal posres, but with identical A and B states, and lambda==0 */
236 forceparams[type].fbposres.pos0,
237 forceparams[type].fbposres.pos0,
252 rfb = pr->fbposres.r;
253 rfb2 = gmx::square(rfb);
255 /* with rfb<0, push particle out of the sphere/cylinder/layer */
263 switch (pr->fbposres.geom)
265 case efbposresSPHERE:
266 /* spherical flat-bottom posres */
268 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
271 v = 0.5 * kk * gmx::square(dr - rfb);
272 fact = -kk * (dr - rfb) / dr; /* Force pointing to the center pos0 */
276 case efbposresCYLINDERX:
277 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
279 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
281 case efbposresCYLINDERY:
282 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
284 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
286 case efbposresCYLINDER:
287 /* equivalent to efbposresCYLINDERZ for backwards compatibility */
288 case efbposresCYLINDERZ:
289 /* cylindrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */
291 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
293 case efbposresX: /* fbdim=XX */
294 case efbposresY: /* fbdim=YY */
295 case efbposresZ: /* fbdim=ZZ */
296 /* 1D flat-bottom potential */
297 fbdim = pr->fbposres.geom - efbposresX;
299 if ((dr > rfb && !bInvert) || (0 < dr && dr < rfb && bInvert))
301 v = 0.5 * kk * gmx::square(dr - rfb);
302 fm[fbdim] = -kk * (dr - rfb);
304 else if ((dr < (-rfb) && !bInvert) || ((-rfb) < dr && dr < 0 && bInvert))
306 v = 0.5 * kk * gmx::square(dr + rfb);
307 fm[fbdim] = -kk * (dr + rfb);
314 for (m = 0; (m < DIM); m++)
317 /* Here we correct for the pbc_dx which included rdist */
318 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm[m];
322 forceWithVirial->addVirialContribution(virial);
328 /*! \brief Compute energies and forces, when requested, for position restraints
330 * Note that position restraints require a different pbc treatment
331 * from other bondeds */
332 template<bool computeForce>
333 real posres(int nbonds,
334 const t_iatom forceatoms[],
335 const t_iparams forceparams[],
337 gmx::ForceWithVirial* forceWithVirial,
338 const struct t_pbc* pbc,
341 RefCoordScaling refcoord_scaling,
346 int i, ai, m, d, type, npbcdim = 0;
349 rvec comA_sc, comB_sc, rdist, dpdl, dx;
351 npbcdim = numPbcDimensions(pbcType);
352 GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
353 if (refcoord_scaling == RefCoordScaling::Com)
357 for (m = 0; m < npbcdim; m++)
359 assert(npbcdim <= DIM);
360 for (d = m; d < npbcdim; d++)
362 comA_sc[m] += comA[d] * pbc->box[d][m];
363 comB_sc[m] += comB[d] * pbc->box[d][m];
368 const real L1 = 1.0 - lambda;
373 GMX_ASSERT(forceWithVirial != nullptr, "When forces are requested we need a force object");
374 f = as_rvec_array(forceWithVirial->force_.data());
377 /* Use intermediate virial buffer to reduce reduction rounding errors */
379 for (i = 0; (i < nbonds);)
381 type = forceatoms[i++];
382 ai = forceatoms[i++];
383 pr = &forceparams[type];
385 /* return dx, rdist, and dpdl */
387 forceparams[type].posres.pos0A,
388 forceparams[type].posres.pos0B,
399 for (m = 0; (m < DIM); m++)
401 kk = L1 * pr->posres.fcA[m] + lambda * pr->posres.fcB[m];
403 vtot += 0.5 * kk * dx[m] * dx[m];
404 *dvdlambda += 0.5 * (pr->posres.fcB[m] - pr->posres.fcA[m]) * dx[m] * dx[m] + fm * dpdl[m];
406 /* Here we correct for the pbc_dx which included rdist */
410 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm;
417 forceWithVirial->addVirialContribution(virial);
425 void posres_wrapper(t_nrnb* nrnb,
426 const InteractionDefinitions& idef,
427 const struct t_pbc* pbc,
429 gmx_enerdata_t* enerd,
430 gmx::ArrayRef<const real> lambda,
431 const t_forcerec* fr,
432 gmx::ForceWithVirial* forceWithVirial)
437 v = posres<true>(idef.il[F_POSRES].size(),
438 idef.il[F_POSRES].iatoms.data(),
439 idef.iparams_posres.data(),
442 fr->pbcType == PbcType::No ? nullptr : pbc,
443 lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Restraint)],
449 enerd->term[F_POSRES] += v;
450 /* If just the force constant changes, the FEP term is linear,
451 * but if k changes, it is not.
453 enerd->dvdl_nonlin[FreeEnergyPerturbationCouplingType::Restraint] += dvdl;
454 inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef.il[F_POSRES].size(), 2));
457 void posres_wrapper_lambda(struct gmx_wallcycle* wcycle,
458 const t_lambda* fepvals,
459 const InteractionDefinitions& idef,
460 const struct t_pbc* pbc,
462 gmx_enerdata_t* enerd,
463 gmx::ArrayRef<const real> lambda,
464 const t_forcerec* fr)
466 wallcycle_sub_start_nocount(wcycle, ewcsRESTRAINTS);
468 auto& foreignTerms = enerd->foreignLambdaTerms;
469 for (int i = 0; i < 1 + foreignTerms.numLambdas(); i++)
473 const real lambda_dum =
474 (i == 0 ? lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Restraint)]
475 : fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Restraint][i - 1]);
476 const real v = posres<false>(idef.il[F_POSRES].size(),
477 idef.il[F_POSRES].iatoms.data(),
478 idef.iparams_posres.data(),
481 fr->pbcType == PbcType::No ? nullptr : pbc,
488 foreignTerms.accumulate(i, v, dvdl);
490 wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
493 /*! \brief Helper function that wraps calls to fbposres for
494 free-energy perturbation */
495 void fbposres_wrapper(t_nrnb* nrnb,
496 const InteractionDefinitions& idef,
497 const struct t_pbc* pbc,
499 gmx_enerdata_t* enerd,
500 const t_forcerec* fr,
501 gmx::ForceWithVirial* forceWithVirial)
505 v = fbposres(idef.il[F_FBPOSRES].size(),
506 idef.il[F_FBPOSRES].iatoms.data(),
507 idef.iparams_fbposres.data(),
510 fr->pbcType == PbcType::No ? nullptr : pbc,
514 enerd->term[F_FBPOSRES] += v;
515 inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef.il[F_FBPOSRES].size(), 2));