src/gromacs/listed_forces/position_restraints.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   9  *
  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.
  14  *
  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.
  19  *
  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.
  24  *
  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.
  32  *
  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.
  35  */
  36 /*! \internal \file
  37  *
  38  * \brief This file defines low-level functions necessary for
  39  * computing energies and forces for position restraints.
  40  *
  41  * \author Mark Abraham <mark.j.abraham@gmail.com>
  42  *
  43  * \ingroup module_listed_forces
  44  */
  45
  46 #include "gmxpre.h"
  47
  48 #include "position_restraints.h"
  49
  50 #include <cassert>
  51 #include <cmath>
  52
  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"
  66
  67 struct gmx_wallcycle;
  68
  69 namespace
  70 {
  71
  72 /*! \brief returns dx, rdist, and dpdl for functions posres() and fbposres()
  73  */
  74 void posres_dx(const rvec   x,
  75                const rvec   pos0A,
  76                const rvec   pos0B,
  77                const rvec   comA_sc,
  78                const rvec   comB_sc,
  79                real         lambda,
  80                const t_pbc* pbc,
  81                int          refcoord_scaling,
  82                int          npbcdim,
  83                rvec         dx,
  84                rvec         rdist,
  85                rvec         dpdl)
  86 {
  87     int  m, d;
  88     real posA, posB, L1, ref = 0.;
  89     rvec pos;
  90
  91     L1 = 1.0 - lambda;
  92
  93     for (m = 0; m < DIM; m++)
  94     {
  95         posA = pos0A[m];
  96         posB = pos0B[m];
  97         if (m < npbcdim)
  98         {
  99             switch (refcoord_scaling)
 100             {
 101                 case erscNO:
 102                     ref      = 0;
 103                     rdist[m] = L1 * posA + lambda * posB;
 104                     dpdl[m]  = posB - posA;
 105                     break;
 106                 case erscALL:
 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++)
 112                     {
 113                         posA += pos0A[d] * pbc->box[d][m];
 114                         posB += pos0B[d] * pbc->box[d][m];
 115                     }
 116                     ref      = L1 * posA + lambda * posB;
 117                     rdist[m] = 0;
 118                     dpdl[m]  = posB - posA;
 119                     break;
 120                 case erscCOM:
 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;
 124                     break;
 125                 default: gmx_fatal(FARGS, "No such scaling method implemented");
 126             }
 127         }
 128         else
 129         {
 130             ref      = L1 * posA + lambda * posB;
 131             rdist[m] = 0;
 132             dpdl[m]  = posB - posA;
 133         }
 134
 135         /* We do pbc_dx with ref+rdist,
 136          * since with only ref we can be up to half a box vector wrong.
 137          */
 138         pos[m] = ref + rdist[m];
 139     }
 140
 141     if (pbc)
 142     {
 143         pbc_dx(pbc, x, pos, dx);
 144     }
 145     else
 146     {
 147         rvec_sub(x, pos, dx);
 148     }
 149 }
 150
 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)
 154 {
 155     int  d;
 156     real dr, dr2, invdr, v, rfb2;
 157
 158     dr2  = 0.0;
 159     rfb2 = gmx::square(rfb);
 160     v    = 0.0;
 161
 162     for (d = 0; d < DIM; d++)
 163     {
 164         if (d != fbdim)
 165         {
 166             dr2 += gmx::square(dx[d]);
 167         }
 168     }
 169
 170     if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
 171     {
 172         dr    = std::sqrt(dr2);
 173         invdr = 1. / dr;
 174         v     = 0.5 * kk * gmx::square(dr - rfb);
 175         for (d = 0; d < DIM; d++)
 176         {
 177             if (d != fbdim)
 178             {
 179                 fm[d] = -kk * (dr - rfb) * dx[d] * invdr; /* Force pointing to the center */
 180             }
 181         }
 182     }
 183
 184     return v;
 185 }
 186
 187 /*! \brief Compute energies and forces for flat-bottomed position restraints
 188  *
 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[],
 194               const rvec            x[],
 195               gmx::ForceWithVirial* forceWithVirial,
 196               const t_pbc*          pbc,
 197               int                   refcoord_scaling,
 198               PbcType               pbcType,
 199               const rvec            com)
 200 /* compute flat-bottomed positions restraints */
 201 {
 202     int              i, ai, m, d, type, npbcdim = 0, fbdim;
 203     const t_iparams* pr;
 204     real             kk, v;
 205     real             dr, dr2, rfb, rfb2, fact;
 206     rvec             com_sc, rdist, dx, dpdl, fm;
 207     gmx_bool         bInvert;
 208
 209     npbcdim = numPbcDimensions(pbcType);
 210     GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
 211     if (refcoord_scaling == erscCOM)
 212     {
 213         clear_rvec(com_sc);
 214         for (m = 0; m < npbcdim; m++)
 215         {
 216             assert(npbcdim <= DIM);
 217             for (d = m; d < npbcdim; d++)
 218             {
 219                 com_sc[m] += com[d] * pbc->box[d][m];
 220             }
 221         }
 222     }
 223
 224     rvec* f      = as_rvec_array(forceWithVirial->force_.data());
 225     real  vtot   = 0.0;
 226     rvec  virial = { 0 };
 227     for (i = 0; (i < nbonds);)
 228     {
 229         type = forceatoms[i++];
 230         ai   = forceatoms[i++];
 231         pr   = &forceparams[type];
 232
 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);
 236
 237         clear_rvec(fm);
 238         v = 0.0;
 239
 240         kk   = pr->fbposres.k;
 241         rfb  = pr->fbposres.r;
 242         rfb2 = gmx::square(rfb);
 243
 244         /* with rfb<0, push particle out of the sphere/cylinder/layer */
 245         bInvert = FALSE;
 246         if (rfb < 0.)
 247         {
 248             bInvert = TRUE;
 249             rfb     = -rfb;
 250         }
 251
 252         switch (pr->fbposres.geom)
 253         {
 254             case efbposresSPHERE:
 255                 /* spherical flat-bottom posres */
 256                 dr2 = norm2(dx);
 257                 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
 258                 {
 259                     dr   = std::sqrt(dr2);
 260                     v    = 0.5 * kk * gmx::square(dr - rfb);
 261                     fact = -kk * (dr - rfb) / dr; /* Force pointing to the center pos0 */
 262                     svmul(fact, dx, fm);
 263                 }
 264                 break;
 265             case efbposresCYLINDERX:
 266                 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
 267                 fbdim = XX;
 268                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 269                 break;
 270             case efbposresCYLINDERY:
 271                 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
 272                 fbdim = YY;
 273                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 274                 break;
 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. */
 279                 fbdim = ZZ;
 280                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 281                 break;
 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;
 287                 dr    = dx[fbdim];
 288                 if ((dr > rfb && !bInvert) || (0 < dr && dr < rfb && bInvert))
 289                 {
 290                     v         = 0.5 * kk * gmx::square(dr - rfb);
 291                     fm[fbdim] = -kk * (dr - rfb);
 292                 }
 293                 else if ((dr < (-rfb) && !bInvert) || ((-rfb) < dr && dr < 0 && bInvert))
 294                 {
 295                     v         = 0.5 * kk * gmx::square(dr + rfb);
 296                     fm[fbdim] = -kk * (dr + rfb);
 297                 }
 298                 break;
 299         }
 300
 301         vtot += v;
 302
 303         for (m = 0; (m < DIM); m++)
 304         {
 305             f[ai][m] += fm[m];
 306             /* Here we correct for the pbc_dx which included rdist */
 307             virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm[m];
 308         }
 309     }
 310
 311     forceWithVirial->addVirialContribution(virial);
 312
 313     return vtot;
 314 }
 315
 316
 317 /*! \brief Compute energies and forces, when requested, for position restraints
 318  *
 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[],
 325             const rvec            x[],
 326             gmx::ForceWithVirial* forceWithVirial,
 327             const struct t_pbc*   pbc,
 328             real                  lambda,
 329             real*                 dvdlambda,
 330             int                   refcoord_scaling,
 331             PbcType               pbcType,
 332             const rvec            comA,
 333             const rvec            comB)
 334 {
 335     int              i, ai, m, d, type, npbcdim = 0;
 336     const t_iparams* pr;
 337     real             kk, fm;
 338     rvec             comA_sc, comB_sc, rdist, dpdl, dx;
 339
 340     npbcdim = numPbcDimensions(pbcType);
 341     GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
 342     if (refcoord_scaling == erscCOM)
 343     {
 344         clear_rvec(comA_sc);
 345         clear_rvec(comB_sc);
 346         for (m = 0; m < npbcdim; m++)
 347         {
 348             assert(npbcdim <= DIM);
 349             for (d = m; d < npbcdim; d++)
 350             {
 351                 comA_sc[m] += comA[d] * pbc->box[d][m];
 352                 comB_sc[m] += comB[d] * pbc->box[d][m];
 353             }
 354         }
 355     }
 356
 357     const real L1 = 1.0 - lambda;
 358
 359     rvec* f;
 360     if (computeForce)
 361     {
 362         GMX_ASSERT(forceWithVirial != nullptr, "When forces are requested we need a force object");
 363         f = as_rvec_array(forceWithVirial->force_.data());
 364     }
 365     real vtot = 0.0;
 366     /* Use intermediate virial buffer to reduce reduction rounding errors */
 367     rvec virial = { 0 };
 368     for (i = 0; (i < nbonds);)
 369     {
 370         type = forceatoms[i++];
 371         ai   = forceatoms[i++];
 372         pr   = &forceparams[type];
 373
 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);
 377
 378         for (m = 0; (m < DIM); m++)
 379         {
 380             kk = L1 * pr->posres.fcA[m] + lambda * pr->posres.fcB[m];
 381             fm = -kk * dx[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];
 384
 385             /* Here we correct for the pbc_dx which included rdist */
 386             if (computeForce)
 387             {
 388                 f[ai][m] += fm;
 389                 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm;
 390             }
 391         }
 392     }
 393
 394     if (computeForce)
 395     {
 396         forceWithVirial->addVirialContribution(virial);
 397     }
 398
 399     return vtot;
 400 }
 401
 402 } // namespace
 403
 404 void posres_wrapper(t_nrnb*                       nrnb,
 405                     const InteractionDefinitions& idef,
 406                     const struct t_pbc*           pbc,
 407                     const rvec*                   x,
 408                     gmx_enerdata_t*               enerd,
 409                     const real*                   lambda,
 410                     const t_forcerec*             fr,
 411                     gmx::ForceWithVirial*         forceWithVirial)
 412 {
 413     real v, dvdl;
 414
 415     dvdl = 0;
 416     v    = posres<true>(idef.il[F_POSRES].size(), idef.il[F_POSRES].iatoms.data(),
 417                      idef.iparams_posres.data(), x, forceWithVirial,
 418                      fr->pbcType == PbcType::No ? nullptr : pbc, lambda[efptRESTRAINT], &dvdl,
 419                      fr->rc_scaling, fr->pbcType, 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.
 423      */
 424     enerd->dvdl_nonlin[efptRESTRAINT] += dvdl;
 425     inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef.il[F_POSRES].size(), 2));
 426 }
 427
 428 void posres_wrapper_lambda(struct gmx_wallcycle*         wcycle,
 429                            const t_lambda*               fepvals,
 430                            const InteractionDefinitions& idef,
 431                            const struct t_pbc*           pbc,
 432                            const rvec                    x[],
 433                            gmx_enerdata_t*               enerd,
 434                            const real*                   lambda,
 435                            const t_forcerec*             fr)
 436 {
 437     real v;
 438
 439     if (idef.il[F_POSRES].empty())
 440     {
 441         return;
 442     }
 443
 444     wallcycle_sub_start_nocount(wcycle, ewcsRESTRAINTS);
 445     for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++)
 446     {
 447         real dvdl = 0;
 448
 449         const real lambda_dum =
 450                 (i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i - 1]);
 451         v = posres<false>(idef.il[F_POSRES].size(), idef.il[F_POSRES].iatoms.data(),
 452                           idef.iparams_posres.data(), x, nullptr,
 453                           fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, &dvdl,
 454                           fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB);
 455         enerd->enerpart_lambda[i] += v;
 456         enerd->dhdlLambda[i] += dvdl;
 457     }
 458     wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 459 }
 460
 461 /*! \brief Helper function that wraps calls to fbposres for
 462     free-energy perturbation */
 463 void fbposres_wrapper(t_nrnb*                       nrnb,
 464                       const InteractionDefinitions& idef,
 465                       const struct t_pbc*           pbc,
 466                       const rvec*                   x,
 467                       gmx_enerdata_t*               enerd,
 468                       const t_forcerec*             fr,
 469                       gmx::ForceWithVirial*         forceWithVirial)
 470 {
 471     real v;
 472
 473     v = fbposres(idef.il[F_FBPOSRES].size(), idef.il[F_FBPOSRES].iatoms.data(),
 474                  idef.iparams_fbposres.data(), x, forceWithVirial,
 475                  fr->pbcType == PbcType::No ? nullptr : pbc, fr->rc_scaling, fr->pbcType, fr->posres_com);
 476     enerd->term[F_FBPOSRES] += v;
 477     inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef.il[F_FBPOSRES].size(), 2));
 478 }