2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017,2018, 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"
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, const rvec pos0A, const rvec pos0B,
75 const rvec comA_sc, const rvec comB_sc,
77 const t_pbc *pbc, int refcoord_scaling, int npbcdim,
78 rvec dx, rvec rdist, rvec dpdl)
81 real posA, posB, L1, ref = 0.;
86 for (m = 0; m < DIM; m++)
92 switch (refcoord_scaling)
96 rdist[m] = L1*posA + lambda*posB;
97 dpdl[m] = posB - posA;
100 /* Box relative coordinates are stored for dimensions with pbc */
101 posA *= pbc->box[m][m];
102 posB *= pbc->box[m][m];
103 assert(npbcdim <= DIM);
104 for (d = m+1; d < npbcdim; d++)
106 posA += pos0A[d]*pbc->box[d][m];
107 posB += pos0B[d]*pbc->box[d][m];
109 ref = L1*posA + lambda*posB;
111 dpdl[m] = posB - posA;
114 ref = L1*comA_sc[m] + lambda*comB_sc[m];
115 rdist[m] = L1*posA + lambda*posB;
116 dpdl[m] = comB_sc[m] - comA_sc[m] + posB - posA;
119 gmx_fatal(FARGS, "No such scaling method implemented");
124 ref = L1*posA + lambda*posB;
126 dpdl[m] = posB - posA;
129 /* We do pbc_dx with ref+rdist,
130 * since with only ref we can be up to half a box vector wrong.
132 pos[m] = ref + rdist[m];
137 pbc_dx(pbc, x, pos, dx);
141 rvec_sub(x, pos, dx);
145 /*! \brief Computes forces and potential for flat-bottom cylindrical restraints.
146 * Returns the flat-bottom potential. */
147 real do_fbposres_cylinder(int fbdim, rvec fm, rvec dx, real rfb, real kk, gmx_bool bInvert)
150 real dr, dr2, invdr, v, rfb2;
153 rfb2 = gmx::square(rfb);
156 for (d = 0; d < DIM; d++)
160 dr2 += gmx::square(dx[d]);
165 ( (dr2 > rfb2 && bInvert == FALSE ) || (dr2 < rfb2 && bInvert == TRUE ) )
170 v = 0.5*kk*gmx::square(dr - rfb);
171 for (d = 0; d < DIM; d++)
175 fm[d] = -kk*(dr-rfb)*dx[d]*invdr; /* Force pointing to the center */
183 /*! \brief Compute energies and forces for flat-bottomed position restraints
185 * Returns the flat-bottomed potential. Same PBC treatment as in
186 * normal position restraints */
187 real fbposres(int nbonds,
188 const t_iatom forceatoms[], const t_iparams forceparams[],
190 gmx::ForceWithVirial *forceWithVirial,
192 int refcoord_scaling, int ePBC, const rvec com)
193 /* compute flat-bottomed positions restraints */
195 int i, ai, m, d, type, npbcdim = 0, fbdim;
198 real dr, dr2, rfb, rfb2, fact;
199 rvec com_sc, rdist, dx, dpdl, fm;
202 npbcdim = ePBC2npbcdim(ePBC);
204 if (refcoord_scaling == erscCOM)
207 for (m = 0; m < npbcdim; m++)
209 assert(npbcdim <= DIM);
210 for (d = m; d < npbcdim; d++)
212 com_sc[m] += com[d]*pbc->box[d][m];
217 rvec *f = as_rvec_array(forceWithVirial->force_.data());
220 for (i = 0; (i < nbonds); )
222 type = forceatoms[i++];
223 ai = forceatoms[i++];
224 pr = &forceparams[type];
226 /* same calculation as for normal posres, but with identical A and B states, and lambda==0 */
227 posres_dx(x[ai], forceparams[type].fbposres.pos0, forceparams[type].fbposres.pos0,
229 pbc, refcoord_scaling, npbcdim,
236 rfb = pr->fbposres.r;
237 rfb2 = gmx::square(rfb);
239 /* with rfb<0, push particle out of the sphere/cylinder/layer */
247 switch (pr->fbposres.geom)
249 case efbposresSPHERE:
250 /* spherical flat-bottom posres */
253 ( (dr2 > rfb2 && bInvert == FALSE ) || (dr2 < rfb2 && bInvert == TRUE ) )
257 v = 0.5*kk*gmx::square(dr - rfb);
258 fact = -kk*(dr-rfb)/dr; /* Force pointing to the center pos0 */
262 case efbposresCYLINDERX:
263 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
265 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
267 case efbposresCYLINDERY:
268 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
270 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
272 case efbposresCYLINDER:
273 /* equivalent to efbposresCYLINDERZ for backwards compatibility */
274 case efbposresCYLINDERZ:
275 /* cylindrical flat-bottom posres in x-y plane. fm[ZZ] = 0. */
277 v = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
279 case efbposresX: /* fbdim=XX */
280 case efbposresY: /* fbdim=YY */
281 case efbposresZ: /* fbdim=ZZ */
282 /* 1D flat-bottom potential */
283 fbdim = pr->fbposres.geom - efbposresX;
285 if ( ( dr > rfb && bInvert == FALSE ) || ( 0 < dr && dr < rfb && bInvert == TRUE ) )
287 v = 0.5*kk*gmx::square(dr - rfb);
288 fm[fbdim] = -kk*(dr - rfb);
290 else if ( (dr < (-rfb) && bInvert == FALSE ) || ( (-rfb) < dr && dr < 0 && bInvert == TRUE ))
292 v = 0.5*kk*gmx::square(dr + rfb);
293 fm[fbdim] = -kk*(dr + rfb);
300 for (m = 0; (m < DIM); m++)
303 /* Here we correct for the pbc_dx which included rdist */
304 virial[m] -= 0.5*(dx[m] + rdist[m])*fm[m];
308 forceWithVirial->addVirialContribution(virial);
314 /*! \brief Compute energies and forces, when requested, for position restraints
316 * Note that position restraints require a different pbc treatment
317 * from other bondeds */
318 template<bool computeForce>
319 real posres(int nbonds,
320 const t_iatom forceatoms[], const t_iparams forceparams[],
322 gmx::ForceWithVirial *forceWithVirial,
323 const struct t_pbc *pbc,
324 real lambda, real *dvdlambda,
325 int refcoord_scaling, int ePBC, const rvec comA, const rvec comB)
327 int i, ai, m, d, type, npbcdim = 0;
330 rvec comA_sc, comB_sc, rdist, dpdl, dx;
332 npbcdim = ePBC2npbcdim(ePBC);
334 if (refcoord_scaling == erscCOM)
338 for (m = 0; m < npbcdim; m++)
340 assert(npbcdim <= DIM);
341 for (d = m; d < npbcdim; d++)
343 comA_sc[m] += comA[d]*pbc->box[d][m];
344 comB_sc[m] += comB[d]*pbc->box[d][m];
349 const real L1 = 1.0 - lambda;
354 GMX_ASSERT(forceWithVirial != nullptr, "When forces are requested we need a force object");
355 f = as_rvec_array(forceWithVirial->force_.data());
358 /* Use intermediate virial buffer to reduce reduction rounding errors */
360 for (i = 0; (i < nbonds); )
362 type = forceatoms[i++];
363 ai = forceatoms[i++];
364 pr = &forceparams[type];
366 /* return dx, rdist, and dpdl */
367 posres_dx(x[ai], forceparams[type].posres.pos0A, forceparams[type].posres.pos0B,
368 comA_sc, comB_sc, lambda,
369 pbc, refcoord_scaling, npbcdim,
372 for (m = 0; (m < DIM); m++)
374 kk = L1*pr->posres.fcA[m] + lambda*pr->posres.fcB[m];
376 vtot += 0.5*kk*dx[m]*dx[m];
378 0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]
381 /* Here we correct for the pbc_dx which included rdist */
385 virial[m] -= 0.5*(dx[m] + rdist[m])*fm;
392 forceWithVirial->addVirialContribution(virial);
401 posres_wrapper(t_nrnb *nrnb,
403 const struct t_pbc *pbc,
405 gmx_enerdata_t *enerd,
407 const t_forcerec *fr,
408 gmx::ForceWithVirial *forceWithVirial)
413 v = posres<true>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms,
414 idef->iparams_posres,
417 fr->ePBC == epbcNONE ? nullptr : pbc,
418 lambda[efptRESTRAINT], &dvdl,
419 fr->rc_scaling, fr->ePBC, fr->posres_com, fr->posres_comB);
420 enerd->term[F_POSRES] += v;
421 /* If just the force constant changes, the FEP term is linear,
422 * but if k changes, it is not.
424 enerd->dvdl_nonlin[efptRESTRAINT] += dvdl;
425 inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef->il[F_POSRES].nr, 2));
429 posres_wrapper_lambda(struct gmx_wallcycle *wcycle,
430 const t_lambda *fepvals,
432 const struct t_pbc *pbc,
434 gmx_enerdata_t *enerd,
436 const t_forcerec *fr)
441 if (0 == idef->il[F_POSRES].nr)
446 wallcycle_sub_start_nocount(wcycle, ewcsRESTRAINTS);
447 for (i = 0; i < enerd->n_lambda; i++)
449 real dvdl_dum = 0, lambda_dum;
451 lambda_dum = (i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i-1]);
452 v = posres<false>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms,
453 idef->iparams_posres,
455 fr->ePBC == epbcNONE ? nullptr : pbc, lambda_dum, &dvdl_dum,
456 fr->rc_scaling, fr->ePBC, fr->posres_com, fr->posres_comB);
457 enerd->enerpart_lambda[i] += v;
459 wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
462 /*! \brief Helper function that wraps calls to fbposres for
463 free-energy perturbation */
464 void fbposres_wrapper(t_nrnb *nrnb,
466 const struct t_pbc *pbc,
468 gmx_enerdata_t *enerd,
469 const t_forcerec *fr,
470 gmx::ForceWithVirial *forceWithVirial)
474 v = fbposres(idef->il[F_FBPOSRES].nr, idef->il[F_FBPOSRES].iatoms,
475 idef->iparams_fbposres,
478 fr->ePBC == epbcNONE ? nullptr : pbc,
479 fr->rc_scaling, fr->ePBC, fr->posres_com);
480 enerd->term[F_FBPOSRES] += v;
481 inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef->il[F_FBPOSRES].nr, 2));