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