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,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.
   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 "gromacs/utility/arrayref.h"
  49 #include "position_restraints.h"
  50
  51 #include <cassert>
  52 #include <cmath>
  53
  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"
  67
  68 struct gmx_wallcycle;
  69
  70 namespace
  71 {
  72
  73 /*! \brief returns dx, rdist, and dpdl for functions posres() and fbposres()
  74  */
  75 void posres_dx(const rvec      x,
  76                const rvec      pos0A,
  77                const rvec      pos0B,
  78                const rvec      comA_sc,
  79                const rvec      comB_sc,
  80                real            lambda,
  81                const t_pbc*    pbc,
  82                RefCoordScaling refcoord_scaling,
  83                int             npbcdim,
  84                rvec            dx,
  85                rvec            rdist,
  86                rvec            dpdl)
  87 {
  88     int  m, d;
  89     real posA, posB, L1, ref = 0.;
  90     rvec pos;
  91
  92     L1 = 1.0 - lambda;
  93
  94     for (m = 0; m < DIM; m++)
  95     {
  96         posA = pos0A[m];
  97         posB = pos0B[m];
  98         if (m < npbcdim)
  99         {
 100             switch (refcoord_scaling)
 101             {
 102                 case RefCoordScaling::No:
 103                     ref      = 0;
 104                     rdist[m] = L1 * posA + lambda * posB;
 105                     dpdl[m]  = posB - posA;
 106                     break;
 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++)
 113                     {
 114                         posA += pos0A[d] * pbc->box[d][m];
 115                         posB += pos0B[d] * pbc->box[d][m];
 116                     }
 117                     ref      = L1 * posA + lambda * posB;
 118                     rdist[m] = 0;
 119                     dpdl[m]  = posB - posA;
 120                     break;
 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;
 125                     break;
 126                 default: gmx_fatal(FARGS, "No such scaling method implemented");
 127             }
 128         }
 129         else
 130         {
 131             ref      = L1 * posA + lambda * posB;
 132             rdist[m] = 0;
 133             dpdl[m]  = posB - posA;
 134         }
 135
 136         /* We do pbc_dx with ref+rdist,
 137          * since with only ref we can be up to half a box vector wrong.
 138          */
 139         pos[m] = ref + rdist[m];
 140     }
 141
 142     if (pbc)
 143     {
 144         pbc_dx(pbc, x, pos, dx);
 145     }
 146     else
 147     {
 148         rvec_sub(x, pos, dx);
 149     }
 150 }
 151
 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)
 155 {
 156     int  d;
 157     real dr, dr2, invdr, v, rfb2;
 158
 159     dr2  = 0.0;
 160     rfb2 = gmx::square(rfb);
 161     v    = 0.0;
 162
 163     for (d = 0; d < DIM; d++)
 164     {
 165         if (d != fbdim)
 166         {
 167             dr2 += gmx::square(dx[d]);
 168         }
 169     }
 170
 171     if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
 172     {
 173         dr    = std::sqrt(dr2);
 174         invdr = 1. / dr;
 175         v     = 0.5 * kk * gmx::square(dr - rfb);
 176         for (d = 0; d < DIM; d++)
 177         {
 178             if (d != fbdim)
 179             {
 180                 fm[d] = -kk * (dr - rfb) * dx[d] * invdr; /* Force pointing to the center */
 181             }
 182         }
 183     }
 184
 185     return v;
 186 }
 187
 188 /*! \brief Compute energies and forces for flat-bottomed position restraints
 189  *
 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[],
 195               const rvec            x[],
 196               gmx::ForceWithVirial* forceWithVirial,
 197               const t_pbc*          pbc,
 198               RefCoordScaling       refcoord_scaling,
 199               PbcType               pbcType,
 200               const rvec            com)
 201 /* compute flat-bottomed positions restraints */
 202 {
 203     int              i, ai, m, d, type, npbcdim = 0, fbdim;
 204     const t_iparams* pr;
 205     real             kk, v;
 206     real             dr, dr2, rfb, rfb2, fact;
 207     rvec             com_sc, rdist, dx, dpdl, fm;
 208     gmx_bool         bInvert;
 209
 210     npbcdim = numPbcDimensions(pbcType);
 211     GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
 212     if (refcoord_scaling == RefCoordScaling::Com)
 213     {
 214         clear_rvec(com_sc);
 215         for (m = 0; m < npbcdim; m++)
 216         {
 217             assert(npbcdim <= DIM);
 218             for (d = m; d < npbcdim; d++)
 219             {
 220                 com_sc[m] += com[d] * pbc->box[d][m];
 221             }
 222         }
 223     }
 224
 225     rvec* f      = as_rvec_array(forceWithVirial->force_.data());
 226     real  vtot   = 0.0;
 227     rvec  virial = { 0 };
 228     for (i = 0; (i < nbonds);)
 229     {
 230         type = forceatoms[i++];
 231         ai   = forceatoms[i++];
 232         pr   = &forceparams[type];
 233
 234         /* same calculation as for normal posres, but with identical A and B states, and lambda==0 */
 235         posres_dx(x[ai],
 236                   forceparams[type].fbposres.pos0,
 237                   forceparams[type].fbposres.pos0,
 238                   com_sc,
 239                   com_sc,
 240                   0.0,
 241                   pbc,
 242                   refcoord_scaling,
 243                   npbcdim,
 244                   dx,
 245                   rdist,
 246                   dpdl);
 247
 248         clear_rvec(fm);
 249         v = 0.0;
 250
 251         kk   = pr->fbposres.k;
 252         rfb  = pr->fbposres.r;
 253         rfb2 = gmx::square(rfb);
 254
 255         /* with rfb<0, push particle out of the sphere/cylinder/layer */
 256         bInvert = FALSE;
 257         if (rfb < 0.)
 258         {
 259             bInvert = TRUE;
 260             rfb     = -rfb;
 261         }
 262
 263         switch (pr->fbposres.geom)
 264         {
 265             case efbposresSPHERE:
 266                 /* spherical flat-bottom posres */
 267                 dr2 = norm2(dx);
 268                 if (dr2 > 0.0 && ((dr2 > rfb2 && !bInvert) || (dr2 < rfb2 && bInvert)))
 269                 {
 270                     dr   = std::sqrt(dr2);
 271                     v    = 0.5 * kk * gmx::square(dr - rfb);
 272                     fact = -kk * (dr - rfb) / dr; /* Force pointing to the center pos0 */
 273                     svmul(fact, dx, fm);
 274                 }
 275                 break;
 276             case efbposresCYLINDERX:
 277                 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
 278                 fbdim = XX;
 279                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 280                 break;
 281             case efbposresCYLINDERY:
 282                 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
 283                 fbdim = YY;
 284                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 285                 break;
 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. */
 290                 fbdim = ZZ;
 291                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 292                 break;
 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;
 298                 dr    = dx[fbdim];
 299                 if ((dr > rfb && !bInvert) || (0 < dr && dr < rfb && bInvert))
 300                 {
 301                     v         = 0.5 * kk * gmx::square(dr - rfb);
 302                     fm[fbdim] = -kk * (dr - rfb);
 303                 }
 304                 else if ((dr < (-rfb) && !bInvert) || ((-rfb) < dr && dr < 0 && bInvert))
 305                 {
 306                     v         = 0.5 * kk * gmx::square(dr + rfb);
 307                     fm[fbdim] = -kk * (dr + rfb);
 308                 }
 309                 break;
 310         }
 311
 312         vtot += v;
 313
 314         for (m = 0; (m < DIM); m++)
 315         {
 316             f[ai][m] += fm[m];
 317             /* Here we correct for the pbc_dx which included rdist */
 318             virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm[m];
 319         }
 320     }
 321
 322     forceWithVirial->addVirialContribution(virial);
 323
 324     return vtot;
 325 }
 326
 327
 328 /*! \brief Compute energies and forces, when requested, for position restraints
 329  *
 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[],
 336             const rvec            x[],
 337             gmx::ForceWithVirial* forceWithVirial,
 338             const struct t_pbc*   pbc,
 339             real                  lambda,
 340             real*                 dvdlambda,
 341             RefCoordScaling       refcoord_scaling,
 342             PbcType               pbcType,
 343             const rvec            comA,
 344             const rvec            comB)
 345 {
 346     int              i, ai, m, d, type, npbcdim = 0;
 347     const t_iparams* pr;
 348     real             kk, fm;
 349     rvec             comA_sc, comB_sc, rdist, dpdl, dx;
 350
 351     npbcdim = numPbcDimensions(pbcType);
 352     GMX_ASSERT((pbcType == PbcType::No) == (npbcdim == 0), "");
 353     if (refcoord_scaling == RefCoordScaling::Com)
 354     {
 355         clear_rvec(comA_sc);
 356         clear_rvec(comB_sc);
 357         for (m = 0; m < npbcdim; m++)
 358         {
 359             assert(npbcdim <= DIM);
 360             for (d = m; d < npbcdim; d++)
 361             {
 362                 comA_sc[m] += comA[d] * pbc->box[d][m];
 363                 comB_sc[m] += comB[d] * pbc->box[d][m];
 364             }
 365         }
 366     }
 367
 368     const real L1 = 1.0 - lambda;
 369
 370     rvec* f;
 371     if (computeForce)
 372     {
 373         GMX_ASSERT(forceWithVirial != nullptr, "When forces are requested we need a force object");
 374         f = as_rvec_array(forceWithVirial->force_.data());
 375     }
 376     real vtot = 0.0;
 377     /* Use intermediate virial buffer to reduce reduction rounding errors */
 378     rvec virial = { 0 };
 379     for (i = 0; (i < nbonds);)
 380     {
 381         type = forceatoms[i++];
 382         ai   = forceatoms[i++];
 383         pr   = &forceparams[type];
 384
 385         /* return dx, rdist, and dpdl */
 386         posres_dx(x[ai],
 387                   forceparams[type].posres.pos0A,
 388                   forceparams[type].posres.pos0B,
 389                   comA_sc,
 390                   comB_sc,
 391                   lambda,
 392                   pbc,
 393                   refcoord_scaling,
 394                   npbcdim,
 395                   dx,
 396                   rdist,
 397                   dpdl);
 398
 399         for (m = 0; (m < DIM); m++)
 400         {
 401             kk = L1 * pr->posres.fcA[m] + lambda * pr->posres.fcB[m];
 402             fm = -kk * dx[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];
 405
 406             /* Here we correct for the pbc_dx which included rdist */
 407             if (computeForce)
 408             {
 409                 f[ai][m] += fm;
 410                 virial[m] -= 0.5 * (dx[m] + rdist[m]) * fm;
 411             }
 412         }
 413     }
 414
 415     if (computeForce)
 416     {
 417         forceWithVirial->addVirialContribution(virial);
 418     }
 419
 420     return vtot;
 421 }
 422
 423 } // namespace
 424
 425 void posres_wrapper(t_nrnb*                       nrnb,
 426                     const InteractionDefinitions& idef,
 427                     const struct t_pbc*           pbc,
 428                     const rvec*                   x,
 429                     gmx_enerdata_t*               enerd,
 430                     gmx::ArrayRef<const real>     lambda,
 431                     const t_forcerec*             fr,
 432                     gmx::ForceWithVirial*         forceWithVirial)
 433 {
 434     real v, dvdl;
 435
 436     dvdl = 0;
 437     v    = posres<true>(idef.il[F_POSRES].size(),
 438                      idef.il[F_POSRES].iatoms.data(),
 439                      idef.iparams_posres.data(),
 440                      x,
 441                      forceWithVirial,
 442                      fr->pbcType == PbcType::No ? nullptr : pbc,
 443                      lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Restraint)],
 444                      &dvdl,
 445                      fr->rc_scaling,
 446                      fr->pbcType,
 447                      fr->posres_com,
 448                      fr->posres_comB);
 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.
 452      */
 453     enerd->dvdl_nonlin[FreeEnergyPerturbationCouplingType::Restraint] += dvdl;
 454     inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef.il[F_POSRES].size(), 2));
 455 }
 456
 457 void posres_wrapper_lambda(struct gmx_wallcycle*         wcycle,
 458                            const t_lambda*               fepvals,
 459                            const InteractionDefinitions& idef,
 460                            const struct t_pbc*           pbc,
 461                            const rvec                    x[],
 462                            gmx_enerdata_t*               enerd,
 463                            gmx::ArrayRef<const real>     lambda,
 464                            const t_forcerec*             fr)
 465 {
 466     wallcycle_sub_start_nocount(wcycle, WallCycleSubCounter::Restraints);
 467
 468     auto& foreignTerms = enerd->foreignLambdaTerms;
 469     for (int i = 0; i < 1 + foreignTerms.numLambdas(); i++)
 470     {
 471         real dvdl = 0;
 472
 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(),
 479                                      x,
 480                                      nullptr,
 481                                      fr->pbcType == PbcType::No ? nullptr : pbc,
 482                                      lambda_dum,
 483                                      &dvdl,
 484                                      fr->rc_scaling,
 485                                      fr->pbcType,
 486                                      fr->posres_com,
 487                                      fr->posres_comB);
 488         foreignTerms.accumulate(i, v, dvdl);
 489     }
 490     wallcycle_sub_stop(wcycle, WallCycleSubCounter::Restraints);
 491 }
 492
 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,
 498                       const rvec*                   x,
 499                       gmx_enerdata_t*               enerd,
 500                       const t_forcerec*             fr,
 501                       gmx::ForceWithVirial*         forceWithVirial)
 502 {
 503     real v;
 504
 505     v = fbposres(idef.il[F_FBPOSRES].size(),
 506                  idef.il[F_FBPOSRES].iatoms.data(),
 507                  idef.iparams_fbposres.data(),
 508                  x,
 509                  forceWithVirial,
 510                  fr->pbcType == PbcType::No ? nullptr : pbc,
 511                  fr->rc_scaling,
 512                  fr->pbcType,
 513                  fr->posres_com);
 514     enerd->term[F_FBPOSRES] += v;
 515     inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef.il[F_FBPOSRES].size(), 2));
 516 }