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, 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 <assert.h>
  50
  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, const rvec pos0A, const rvec pos0B,
  75                const rvec comA_sc, const rvec comB_sc,
  76                real lambda,
  77                const t_pbc *pbc, int refcoord_scaling, int npbcdim,
  78                rvec dx, rvec rdist, rvec dpdl)
  79 {
  80     int  m, d;
  81     real posA, posB, L1, ref = 0.;
  82     rvec pos;
  83
  84     L1 = 1.0-lambda;
  85
  86     for (m = 0; m < DIM; m++)
  87     {
  88         posA = pos0A[m];
  89         posB = pos0B[m];
  90         if (m < npbcdim)
  91         {
  92             switch (refcoord_scaling)
  93             {
  94                 case erscNO:
  95                     ref      = 0;
  96                     rdist[m] = L1*posA + lambda*posB;
  97                     dpdl[m]  = posB - posA;
  98                     break;
  99                 case erscALL:
 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++)
 105                     {
 106                         posA += pos0A[d]*pbc->box[d][m];
 107                         posB += pos0B[d]*pbc->box[d][m];
 108                     }
 109                     ref      = L1*posA + lambda*posB;
 110                     rdist[m] = 0;
 111                     dpdl[m]  = posB - posA;
 112                     break;
 113                 case erscCOM:
 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;
 117                     break;
 118                 default:
 119                     gmx_fatal(FARGS, "No such scaling method implemented");
 120             }
 121         }
 122         else
 123         {
 124             ref      = L1*posA + lambda*posB;
 125             rdist[m] = 0;
 126             dpdl[m]  = posB - posA;
 127         }
 128
 129         /* We do pbc_dx with ref+rdist,
 130          * since with only ref we can be up to half a box vector wrong.
 131          */
 132         pos[m] = ref + rdist[m];
 133     }
 134
 135     if (pbc)
 136     {
 137         pbc_dx(pbc, x, pos, dx);
 138     }
 139     else
 140     {
 141         rvec_sub(x, pos, dx);
 142     }
 143 }
 144
 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)
 148 {
 149     int     d;
 150     real    dr, dr2, invdr, v, rfb2;
 151
 152     dr2  = 0.0;
 153     rfb2 = gmx::square(rfb);
 154     v    = 0.0;
 155
 156     for (d = 0; d < DIM; d++)
 157     {
 158         if (d != fbdim)
 159         {
 160             dr2 += gmx::square(dx[d]);
 161         }
 162     }
 163
 164     if  (dr2 > 0.0 &&
 165          ( (dr2 > rfb2 && bInvert == FALSE ) || (dr2 < rfb2 && bInvert == TRUE ) )
 166          )
 167     {
 168         dr     = std::sqrt(dr2);
 169         invdr  = 1./dr;
 170         v      = 0.5*kk*gmx::square(dr - rfb);
 171         for (d = 0; d < DIM; d++)
 172         {
 173             if (d != fbdim)
 174             {
 175                 fm[d] = -kk*(dr-rfb)*dx[d]*invdr; /* Force pointing to the center */
 176             }
 177         }
 178     }
 179
 180     return v;
 181 }
 182
 183 /*! \brief Compute energies and forces for flat-bottomed position restraints
 184  *
 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[],
 189               const rvec x[],
 190               gmx::ForceWithVirial *forceWithVirial,
 191               const t_pbc *pbc,
 192               int refcoord_scaling, int ePBC, const rvec com)
 193 /* compute flat-bottomed positions restraints */
 194 {
 195     int              i, ai, m, d, type, npbcdim = 0, fbdim;
 196     const t_iparams *pr;
 197     real             kk, v;
 198     real             dr, dr2, rfb, rfb2, fact;
 199     rvec             com_sc, rdist, dx, dpdl, fm;
 200     gmx_bool         bInvert;
 201
 202     npbcdim = ePBC2npbcdim(ePBC);
 203
 204     if (refcoord_scaling == erscCOM)
 205     {
 206         clear_rvec(com_sc);
 207         for (m = 0; m < npbcdim; m++)
 208         {
 209             assert(npbcdim <= DIM);
 210             for (d = m; d < npbcdim; d++)
 211             {
 212                 com_sc[m] += com[d]*pbc->box[d][m];
 213             }
 214         }
 215     }
 216
 217     rvec *f      = as_rvec_array(forceWithVirial->force_.data());
 218     real  vtot   = 0.0;
 219     rvec  virial = { 0 };
 220     for (i = 0; (i < nbonds); )
 221     {
 222         type = forceatoms[i++];
 223         ai   = forceatoms[i++];
 224         pr   = &forceparams[type];
 225
 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,
 228                   com_sc, com_sc, 0.0,
 229                   pbc, refcoord_scaling, npbcdim,
 230                   dx, rdist, dpdl);
 231
 232         clear_rvec(fm);
 233         v = 0.0;
 234
 235         kk   = pr->fbposres.k;
 236         rfb  = pr->fbposres.r;
 237         rfb2 = gmx::square(rfb);
 238
 239         /* with rfb<0, push particle out of the sphere/cylinder/layer */
 240         bInvert = FALSE;
 241         if (rfb < 0.)
 242         {
 243             bInvert = TRUE;
 244             rfb     = -rfb;
 245         }
 246
 247         switch (pr->fbposres.geom)
 248         {
 249             case efbposresSPHERE:
 250                 /* spherical flat-bottom posres */
 251                 dr2 = norm2(dx);
 252                 if (dr2 > 0.0 &&
 253                     ( (dr2 > rfb2 && bInvert == FALSE ) || (dr2 < rfb2 && bInvert == TRUE ) )
 254                     )
 255                 {
 256                     dr   = std::sqrt(dr2);
 257                     v    = 0.5*kk*gmx::square(dr - rfb);
 258                     fact = -kk*(dr-rfb)/dr; /* Force pointing to the center pos0 */
 259                     svmul(fact, dx, fm);
 260                 }
 261                 break;
 262             case efbposresCYLINDERX:
 263                 /* cylindrical flat-bottom posres in y-z plane. fm[XX] = 0. */
 264                 fbdim = XX;
 265                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 266                 break;
 267             case efbposresCYLINDERY:
 268                 /* cylindrical flat-bottom posres in x-z plane. fm[YY] = 0. */
 269                 fbdim = YY;
 270                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 271                 break;
 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. */
 276                 fbdim = ZZ;
 277                 v     = do_fbposres_cylinder(fbdim, fm, dx, rfb, kk, bInvert);
 278                 break;
 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;
 284                 dr    = dx[fbdim];
 285                 if ( ( dr > rfb && bInvert == FALSE ) || ( 0 < dr && dr < rfb && bInvert == TRUE )  )
 286                 {
 287                     v         = 0.5*kk*gmx::square(dr - rfb);
 288                     fm[fbdim] = -kk*(dr - rfb);
 289                 }
 290                 else if ( (dr < (-rfb) && bInvert == FALSE ) || ( (-rfb) < dr && dr < 0 && bInvert == TRUE ))
 291                 {
 292                     v         = 0.5*kk*gmx::square(dr + rfb);
 293                     fm[fbdim] = -kk*(dr + rfb);
 294                 }
 295                 break;
 296         }
 297
 298         vtot += v;
 299
 300         for (m = 0; (m < DIM); m++)
 301         {
 302             f[ai][m]  += fm[m];
 303             /* Here we correct for the pbc_dx which included rdist */
 304             virial[m] -= 0.5*(dx[m] + rdist[m])*fm[m];
 305         }
 306     }
 307
 308     forceWithVirial->addVirialContribution(virial);
 309
 310     return vtot;
 311 }
 312
 313
 314 /*! \brief Compute energies and forces, when requested, for position restraints
 315  *
 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[],
 321             const rvec x[],
 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)
 326 {
 327     int              i, ai, m, d, type, npbcdim = 0;
 328     const t_iparams *pr;
 329     real             kk, fm;
 330     rvec             comA_sc, comB_sc, rdist, dpdl, dx;
 331
 332     npbcdim = ePBC2npbcdim(ePBC);
 333
 334     if (refcoord_scaling == erscCOM)
 335     {
 336         clear_rvec(comA_sc);
 337         clear_rvec(comB_sc);
 338         for (m = 0; m < npbcdim; m++)
 339         {
 340             assert(npbcdim <= DIM);
 341             for (d = m; d < npbcdim; d++)
 342             {
 343                 comA_sc[m] += comA[d]*pbc->box[d][m];
 344                 comB_sc[m] += comB[d]*pbc->box[d][m];
 345             }
 346         }
 347     }
 348
 349     const real  L1     = 1.0 - lambda;
 350
 351     rvec       *f;
 352     if (computeForce)
 353     {
 354         GMX_ASSERT(forceWithVirial != nullptr, "When forces are requested we need a force object");
 355         f              = as_rvec_array(forceWithVirial->force_.data());
 356     }
 357     real        vtot   = 0.0;
 358     /* Use intermediate virial buffer to reduce reduction rounding errors */
 359     rvec        virial = { 0 };
 360     for (i = 0; (i < nbonds); )
 361     {
 362         type = forceatoms[i++];
 363         ai   = forceatoms[i++];
 364         pr   = &forceparams[type];
 365
 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,
 370                   dx, rdist, dpdl);
 371
 372         for (m = 0; (m < DIM); m++)
 373         {
 374             kk          = L1*pr->posres.fcA[m] + lambda*pr->posres.fcB[m];
 375             fm          = -kk*dx[m];
 376             vtot       += 0.5*kk*dx[m]*dx[m];
 377             *dvdlambda +=
 378                 0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]
 379                 + fm*dpdl[m];
 380
 381             /* Here we correct for the pbc_dx which included rdist */
 382             if (computeForce)
 383             {
 384                 f[ai][m]  += fm;
 385                 virial[m] -= 0.5*(dx[m] + rdist[m])*fm;
 386             }
 387         }
 388     }
 389
 390     if (computeForce)
 391     {
 392         forceWithVirial->addVirialContribution(virial);
 393     }
 394
 395     return vtot;
 396 }
 397
 398 } // namespace
 399
 400 void
 401 posres_wrapper(t_nrnb               *nrnb,
 402                const t_idef         *idef,
 403                const struct t_pbc   *pbc,
 404                const rvec           *x,
 405                gmx_enerdata_t       *enerd,
 406                const real           *lambda,
 407                const t_forcerec     *fr,
 408                gmx::ForceWithVirial *forceWithVirial)
 409 {
 410     real  v, dvdl;
 411
 412     dvdl = 0;
 413     v    = posres<true>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms,
 414                         idef->iparams_posres,
 415                         x,
 416                         forceWithVirial,
 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.
 423      */
 424     enerd->dvdl_nonlin[efptRESTRAINT] += dvdl;
 425     inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef->il[F_POSRES].nr, 2));
 426 }
 427
 428 void
 429 posres_wrapper_lambda(struct gmx_wallcycle *wcycle,
 430                       const t_lambda       *fepvals,
 431                       const t_idef         *idef,
 432                       const struct t_pbc   *pbc,
 433                       const rvec            x[],
 434                       gmx_enerdata_t       *enerd,
 435                       const real           *lambda,
 436                       const t_forcerec     *fr)
 437 {
 438     real  v;
 439     int   i;
 440
 441     if (0 == idef->il[F_POSRES].nr)
 442     {
 443         return;
 444     }
 445
 446     wallcycle_sub_start_nocount(wcycle, ewcsRESTRAINTS);
 447     for (i = 0; i < enerd->n_lambda; i++)
 448     {
 449         real dvdl_dum = 0, lambda_dum;
 450
 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,
 454                                    x, nullptr,
 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;
 458     }
 459     wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
 460 }
 461
 462 /*! \brief Helper function that wraps calls to fbposres for
 463     free-energy perturbation */
 464 void fbposres_wrapper(t_nrnb               *nrnb,
 465                       const t_idef         *idef,
 466                       const struct t_pbc   *pbc,
 467                       const rvec           *x,
 468                       gmx_enerdata_t       *enerd,
 469                       const t_forcerec     *fr,
 470                       gmx::ForceWithVirial *forceWithVirial)
 471 {
 472     real  v;
 473
 474     v = fbposres(idef->il[F_FBPOSRES].nr, idef->il[F_FBPOSRES].iatoms,
 475                  idef->iparams_fbposres,
 476                  x,
 477                  forceWithVirial,
 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));
 482 }