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