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