src/gromacs/listed-forces/listed-forces.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015, 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 high-level functions for mdrun to compute
  38  * energies and forces for listed interactions.
  39  *
  40  * \author Mark Abraham <mark.j.abraham@gmail.com>
  41  *
  42  * \ingroup module_listed-forces
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "listed-forces.h"
  47
  48 #include "config.h"
  49
  50 #include <assert.h>
  51
  52 #include <algorithm>
  53
  54 #include "gromacs/legacyheaders/disre.h"
  55 #include "gromacs/legacyheaders/force.h"
  56 #include "gromacs/legacyheaders/network.h"
  57 #include "gromacs/legacyheaders/nrnb.h"
  58 #include "gromacs/legacyheaders/orires.h"
  59 #include "gromacs/legacyheaders/types/force_flags.h"
  60 #include "gromacs/listed-forces/bonded.h"
  61 #include "gromacs/listed-forces/position-restraints.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/mdlib/forcerec-threading.h"
  64 #include "gromacs/pbcutil/ishift.h"
  65 #include "gromacs/pbcutil/pbc.h"
  66 #include "gromacs/simd/simd.h"
  67 #include "gromacs/timing/wallcycle.h"
  68 #include "gromacs/utility/smalloc.h"
  69
  70 #include "pairs.h"
  71
  72 namespace
  73 {
  74
  75 /*! \brief Return true if ftype is an explicit pair-listed LJ or
  76  * COULOMB interaction type: bonded LJ (usually 1-4), or special
  77  * listed non-bonded for FEP. */
  78 bool
  79 isPairInteraction(int ftype)
  80 {
  81     return ((ftype) >= F_LJ14 && (ftype) <= F_LJC_PAIRS_NB);
  82 }
  83
  84 /*! \brief Zero thread-local force-output buffers */
  85 void
  86 zero_thread_forces(f_thread_t *f_t, int n,
  87                    int nblock, int blocksize)
  88 {
  89     int b, a0, a1, a, i, j;
  90
  91     if (n > f_t->f_nalloc)
  92     {
  93         f_t->f_nalloc = over_alloc_large(n);
  94         srenew(f_t->f, f_t->f_nalloc);
  95     }
  96
  97     if (!bitmask_is_zero(f_t->red_mask))
  98     {
  99         for (b = 0; b < nblock; b++)
 100         {
 101             if (bitmask_is_set(f_t->red_mask, b))
 102             {
 103                 a0 = b*blocksize;
 104                 a1 = std::min((b+1)*blocksize, n);
 105                 for (a = a0; a < a1; a++)
 106                 {
 107                     clear_rvec(f_t->f[a]);
 108                 }
 109             }
 110         }
 111     }
 112     for (i = 0; i < SHIFTS; i++)
 113     {
 114         clear_rvec(f_t->fshift[i]);
 115     }
 116     for (i = 0; i < F_NRE; i++)
 117     {
 118         f_t->ener[i] = 0;
 119     }
 120     for (i = 0; i < egNR; i++)
 121     {
 122         for (j = 0; j < f_t->grpp.nener; j++)
 123         {
 124             f_t->grpp.ener[i][j] = 0;
 125         }
 126     }
 127     for (i = 0; i < efptNR; i++)
 128     {
 129         f_t->dvdl[i] = 0;
 130     }
 131 }
 132
 133 /*! \brief The max thread number is arbitrary, we used a fixed number
 134  * to avoid memory management.  Using more than 16 threads is probably
 135  * never useful performance wise. */
 136 #define MAX_BONDED_THREADS 256
 137
 138 /*! \brief Reduce thread-local force buffers */
 139 void
 140 reduce_thread_force_buffer(int n, rvec *f,
 141                            int nthreads, f_thread_t *f_t,
 142                            int nblock, int block_size)
 143 {
 144     int b;
 145
 146     if (nthreads > MAX_BONDED_THREADS)
 147     {
 148         gmx_fatal(FARGS, "Can not reduce bonded forces on more than %d threads",
 149                   MAX_BONDED_THREADS);
 150     }
 151
 152     /* This reduction can run on any number of threads,
 153      * independently of nthreads.
 154      */
 155 #pragma omp parallel for num_threads(nthreads) schedule(static)
 156     for (b = 0; b < nblock; b++)
 157     {
 158         rvec *fp[MAX_BONDED_THREADS];
 159         int   nfb, ft, fb;
 160         int   a0, a1, a;
 161
 162         /* Determine which threads contribute to this block */
 163         nfb = 0;
 164         for (ft = 1; ft < nthreads; ft++)
 165         {
 166             if (bitmask_is_set(f_t[ft].red_mask, b))
 167             {
 168                 fp[nfb++] = f_t[ft].f;
 169             }
 170         }
 171         if (nfb > 0)
 172         {
 173             /* Reduce force buffers for threads that contribute */
 174             a0 =  b   *block_size;
 175             a1 = (b+1)*block_size;
 176             a1 = std::min(a1, n);
 177             for (a = a0; a < a1; a++)
 178             {
 179                 for (fb = 0; fb < nfb; fb++)
 180                 {
 181                     rvec_inc(f[a], fp[fb][a]);
 182                 }
 183             }
 184         }
 185     }
 186 }
 187
 188 /*! \brief Reduce thread-local forces */
 189 void
 190 reduce_thread_forces(int n, rvec *f, rvec *fshift,
 191                      real *ener, gmx_grppairener_t *grpp, real *dvdl,
 192                      int nthreads, f_thread_t *f_t,
 193                      int nblock, int block_size,
 194                      gmx_bool bCalcEnerVir,
 195                      gmx_bool bDHDL)
 196 {
 197     if (nblock > 0)
 198     {
 199         /* Reduce the bonded force buffer */
 200         reduce_thread_force_buffer(n, f, nthreads, f_t, nblock, block_size);
 201     }
 202
 203     /* When necessary, reduce energy and virial using one thread only */
 204     if (bCalcEnerVir)
 205     {
 206         int t, i, j;
 207
 208         for (i = 0; i < SHIFTS; i++)
 209         {
 210             for (t = 1; t < nthreads; t++)
 211             {
 212                 rvec_inc(fshift[i], f_t[t].fshift[i]);
 213             }
 214         }
 215         for (i = 0; i < F_NRE; i++)
 216         {
 217             for (t = 1; t < nthreads; t++)
 218             {
 219                 ener[i] += f_t[t].ener[i];
 220             }
 221         }
 222         for (i = 0; i < egNR; i++)
 223         {
 224             for (j = 0; j < f_t[1].grpp.nener; j++)
 225             {
 226                 for (t = 1; t < nthreads; t++)
 227                 {
 228
 229                     grpp->ener[i][j] += f_t[t].grpp.ener[i][j];
 230                 }
 231             }
 232         }
 233         if (bDHDL)
 234         {
 235             for (i = 0; i < efptNR; i++)
 236             {
 237
 238                 for (t = 1; t < nthreads; t++)
 239                 {
 240                     dvdl[i] += f_t[t].dvdl[i];
 241                 }
 242             }
 243         }
 244     }
 245 }
 246
 247 /*! \brief Calculate one element of the list of bonded interactions
 248     for this thread */
 249 real
 250 calc_one_bond(int thread,
 251               int ftype, const t_idef *idef,
 252               const rvec x[], rvec f[], rvec fshift[],
 253               t_forcerec *fr,
 254               const t_pbc *pbc, const t_graph *g,
 255               gmx_grppairener_t *grpp,
 256               t_nrnb *nrnb,
 257               real *lambda, real *dvdl,
 258               const t_mdatoms *md, t_fcdata *fcd,
 259               gmx_bool bCalcEnerVir,
 260               int *global_atom_index)
 261 {
 262     int      nat1, nbonds, efptFTYPE;
 263     real     v = 0;
 264     t_iatom *iatoms;
 265     int      nb0, nbn;
 266
 267     if (IS_RESTRAINT_TYPE(ftype))
 268     {
 269         efptFTYPE = efptRESTRAINT;
 270     }
 271     else
 272     {
 273         efptFTYPE = efptBONDED;
 274     }
 275
 276     nat1      = interaction_function[ftype].nratoms + 1;
 277     nbonds    = idef->il[ftype].nr/nat1;
 278     iatoms    = idef->il[ftype].iatoms;
 279
 280     nb0 = idef->il_thread_division[ftype*(idef->nthreads+1)+thread];
 281     nbn = idef->il_thread_division[ftype*(idef->nthreads+1)+thread+1] - nb0;
 282
 283     if (!isPairInteraction(ftype))
 284     {
 285         if (ftype == F_CMAP)
 286         {
 287             v = cmap_dihs(nbn, iatoms+nb0,
 288                           idef->iparams, &idef->cmap_grid,
 289                           x, f, fshift,
 290                           pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 291                           md, fcd, global_atom_index);
 292         }
 293 #ifdef GMX_SIMD_HAVE_REAL
 294         else if (ftype == F_ANGLES &&
 295                  !bCalcEnerVir && fr->efep == efepNO)
 296         {
 297             /* No energies, shift forces, dvdl */
 298             angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 299                                idef->iparams,
 300                                x, f,
 301                                pbc, g, lambda[efptFTYPE], md, fcd,
 302                                global_atom_index);
 303             v = 0;
 304         }
 305 #endif
 306         else if (ftype == F_PDIHS &&
 307                  !bCalcEnerVir && fr->efep == efepNO)
 308         {
 309             /* No energies, shift forces, dvdl */
 310 #ifdef GMX_SIMD_HAVE_REAL
 311             pdihs_noener_simd
 312 #else
 313             pdihs_noener
 314 #endif
 315                 (nbn, idef->il[ftype].iatoms+nb0,
 316                 idef->iparams,
 317                 x, f,
 318                 pbc, g, lambda[efptFTYPE], md, fcd,
 319                 global_atom_index);
 320             v = 0;
 321         }
 322 #ifdef GMX_SIMD_HAVE_REAL
 323         else if (ftype == F_RBDIHS &&
 324                  !bCalcEnerVir && fr->efep == efepNO)
 325         {
 326             /* No energies, shift forces, dvdl */
 327             rbdihs_noener_simd(nbn, idef->il[ftype].iatoms+nb0,
 328                                idef->iparams,
 329                                (const rvec*)x, f,
 330                                pbc, g, lambda[efptFTYPE], md, fcd,
 331                                global_atom_index);
 332             v = 0;
 333         }
 334 #endif
 335         else
 336         {
 337             v = interaction_function[ftype].ifunc(nbn, iatoms+nb0,
 338                                                   idef->iparams,
 339                                                   x, f, fshift,
 340                                                   pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]),
 341                                                   md, fcd, global_atom_index);
 342         }
 343     }
 344     else
 345     {
 346         v = do_pairs(ftype, nbn, iatoms+nb0, idef->iparams, x, f, fshift,
 347                      pbc, g, lambda, dvdl, md, fr, grpp, global_atom_index);
 348     }
 349
 350     if (thread == 0)
 351     {
 352         inc_nrnb(nrnb, interaction_function[ftype].nrnb_ind, nbonds);
 353     }
 354
 355     return v;
 356 }
 357
 358 } // namespace
 359
 360 gmx_bool
 361 ftype_is_bonded_potential(int ftype)
 362 {
 363     return
 364         (interaction_function[ftype].flags & IF_BOND) &&
 365         !(ftype == F_CONNBONDS || ftype == F_POSRES || ftype == F_FBPOSRES) &&
 366         (ftype < F_GB12 || ftype > F_GB14);
 367 }
 368
 369 void calc_listed(const gmx_multisim_t *ms,
 370                  const t_idef *idef,
 371                  const rvec x[], history_t *hist,
 372                  rvec f[], t_forcerec *fr,
 373                  const struct t_pbc *pbc,
 374                  const struct t_pbc *pbc_full,
 375                  const struct t_graph *g,
 376                  gmx_enerdata_t *enerd, t_nrnb *nrnb,
 377                  real *lambda,
 378                  const t_mdatoms *md,
 379                  t_fcdata *fcd, int *global_atom_index,
 380                  int force_flags)
 381 {
 382     gmx_bool      bCalcEnerVir;
 383     int           i;
 384     real          dvdl[efptNR]; /* The dummy array is to have a place to store the dhdl at other values
 385                                                         of lambda, which will be thrown away in the end*/
 386     const  t_pbc *pbc_null;
 387     int           thread;
 388
 389     assert(fr->nthreads == idef->nthreads);
 390
 391     bCalcEnerVir = (force_flags & (GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY));
 392
 393     for (i = 0; i < efptNR; i++)
 394     {
 395         dvdl[i] = 0.0;
 396     }
 397     if (fr->bMolPBC)
 398     {
 399         pbc_null = pbc;
 400     }
 401     else
 402     {
 403         pbc_null = NULL;
 404     }
 405
 406 #ifdef DEBUG
 407     if (g && debug)
 408     {
 409         p_graph(debug, "Bondage is fun", g);
 410     }
 411 #endif
 412
 413     if (idef->il[F_POSRES].nr > 0)
 414     {
 415         posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr);
 416     }
 417
 418     if (idef->il[F_FBPOSRES].nr > 0)
 419     {
 420         fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr);
 421     }
 422
 423     /* Do pre force calculation stuff which might require communication */
 424     if (idef->il[F_ORIRES].nr)
 425     {
 426         enerd->term[F_ORIRESDEV] =
 427             calc_orires_dev(ms, idef->il[F_ORIRES].nr,
 428                             idef->il[F_ORIRES].iatoms,
 429                             idef->iparams, md, x,
 430                             pbc_null, fcd, hist);
 431     }
 432     if (idef->il[F_DISRES].nr)
 433     {
 434         calc_disres_R_6(idef->il[F_DISRES].nr,
 435                         idef->il[F_DISRES].iatoms,
 436                         idef->iparams, x, pbc_null,
 437                         fcd, hist);
 438 #ifdef GMX_MPI
 439         if (fcd->disres.nsystems > 1)
 440         {
 441             gmx_sum_sim(2*fcd->disres.nres, fcd->disres.Rt_6, ms);
 442         }
 443 #endif
 444     }
 445
 446 #pragma omp parallel for num_threads(fr->nthreads) schedule(static)
 447     for (thread = 0; thread < fr->nthreads; thread++)
 448     {
 449         int                ftype;
 450         real              *epot, v;
 451         /* thread stuff */
 452         rvec              *ft, *fshift;
 453         real              *dvdlt;
 454         gmx_grppairener_t *grpp;
 455
 456         if (thread == 0)
 457         {
 458             ft     = f;
 459             fshift = fr->fshift;
 460             epot   = enerd->term;
 461             grpp   = &enerd->grpp;
 462             dvdlt  = dvdl;
 463         }
 464         else
 465         {
 466             zero_thread_forces(&fr->f_t[thread], fr->natoms_force,
 467                                fr->red_nblock, 1<<fr->red_ashift);
 468
 469             ft     = fr->f_t[thread].f;
 470             fshift = fr->f_t[thread].fshift;
 471             epot   = fr->f_t[thread].ener;
 472             grpp   = &fr->f_t[thread].grpp;
 473             dvdlt  = fr->f_t[thread].dvdl;
 474         }
 475         /* Loop over all bonded force types to calculate the bonded forces */
 476         for (ftype = 0; (ftype < F_NRE); ftype++)
 477         {
 478             if (idef->il[ftype].nr > 0 && ftype_is_bonded_potential(ftype))
 479             {
 480                 v = calc_one_bond(thread, ftype, idef, x,
 481                                   ft, fshift, fr, pbc_null, g, grpp,
 482                                   nrnb, lambda, dvdlt,
 483                                   md, fcd, bCalcEnerVir,
 484                                   global_atom_index);
 485                 epot[ftype] += v;
 486             }
 487         }
 488     }
 489     if (fr->nthreads > 1)
 490     {
 491         reduce_thread_forces(fr->natoms_force, f, fr->fshift,
 492                              enerd->term, &enerd->grpp, dvdl,
 493                              fr->nthreads, fr->f_t,
 494                              fr->red_nblock, 1<<fr->red_ashift,
 495                              bCalcEnerVir,
 496                              force_flags & GMX_FORCE_DHDL);
 497     }
 498     if (force_flags & GMX_FORCE_DHDL)
 499     {
 500         for (i = 0; i < efptNR; i++)
 501         {
 502             enerd->dvdl_nonlin[i] += dvdl[i];
 503         }
 504     }
 505
 506     /* Copy the sum of violations for the distance restraints from fcd */
 507     if (fcd)
 508     {
 509         enerd->term[F_DISRESVIOL] = fcd->disres.sumviol;
 510
 511     }
 512 }
 513
 514 void calc_listed_lambda(const t_idef *idef,
 515                         const rvec x[],
 516                         t_forcerec *fr,
 517                         const struct t_pbc *pbc, const struct t_graph *g,
 518                         gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb,
 519                         real *lambda,
 520                         const t_mdatoms *md,
 521                         t_fcdata *fcd,
 522                         int *global_atom_index)
 523 {
 524     int           ftype, nr_nonperturbed, nr;
 525     real          v;
 526     real          dvdl_dum[efptNR] = {0};
 527     rvec         *f, *fshift;
 528     const  t_pbc *pbc_null;
 529     t_idef        idef_fe;
 530
 531     if (fr->bMolPBC)
 532     {
 533         pbc_null = pbc;
 534     }
 535     else
 536     {
 537         pbc_null = NULL;
 538     }
 539
 540     /* Copy the whole idef, so we can modify the contents locally */
 541     idef_fe          = *idef;
 542     idef_fe.nthreads = 1;
 543     snew(idef_fe.il_thread_division, F_NRE*(idef_fe.nthreads+1));
 544
 545     /* We already have the forces, so we use temp buffers here */
 546     snew(f, fr->natoms_force);
 547     snew(fshift, SHIFTS);
 548
 549     /* Loop over all bonded force types to calculate the bonded energies */
 550     for (ftype = 0; (ftype < F_NRE); ftype++)
 551     {
 552         if (ftype_is_bonded_potential(ftype))
 553         {
 554             /* Set the work range of thread 0 to the perturbed bondeds only */
 555             nr_nonperturbed                       = idef->il[ftype].nr_nonperturbed;
 556             nr                                    = idef->il[ftype].nr;
 557             idef_fe.il_thread_division[ftype*2+0] = nr_nonperturbed;
 558             idef_fe.il_thread_division[ftype*2+1] = nr;
 559
 560             /* This is only to get the flop count correct */
 561             idef_fe.il[ftype].nr = nr - nr_nonperturbed;
 562
 563             if (nr - nr_nonperturbed > 0)
 564             {
 565                 v = calc_one_bond(0, ftype, &idef_fe,
 566                                   x, f, fshift, fr, pbc_null, g,
 567                                   grpp, nrnb, lambda, dvdl_dum,
 568                                   md, fcd, TRUE,
 569                                   global_atom_index);
 570                 epot[ftype] += v;
 571             }
 572         }
 573     }
 574
 575     sfree(fshift);
 576     sfree(f);
 577
 578     sfree(idef_fe.il_thread_division);
 579 }
 580
 581 void
 582 do_force_listed(gmx_wallcycle        *wcycle,
 583                 matrix                box,
 584                 const t_lambda       *fepvals,
 585                 const gmx_multisim_t *ms,
 586                 const t_idef         *idef,
 587                 const rvec            x[],
 588                 history_t            *hist,
 589                 rvec                  f[],
 590                 t_forcerec           *fr,
 591                 const struct t_pbc   *pbc,
 592                 const struct t_graph *graph,
 593                 gmx_enerdata_t       *enerd,
 594                 t_nrnb               *nrnb,
 595                 real                 *lambda,
 596                 const t_mdatoms      *md,
 597                 t_fcdata             *fcd,
 598                 int                  *global_atom_index,
 599                 int                   flags)
 600 {
 601     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 602
 603     if (!(flags & GMX_FORCE_LISTED))
 604     {
 605         return;
 606     }
 607
 608     wallcycle_sub_start(wcycle, ewcsLISTED);
 609
 610     if ((idef->il[F_POSRES].nr > 0) ||
 611         (idef->il[F_FBPOSRES].nr > 0))
 612     {
 613         set_pbc(&pbc_full, fr->ePBC, box);
 614     }
 615     calc_listed(ms, idef, x, hist, f, fr, pbc, &pbc_full,
 616                 graph, enerd, nrnb, lambda, md, fcd,
 617                 global_atom_index, flags);
 618
 619     /* Check if we have to determine energy differences
 620      * at foreign lambda's.
 621      */
 622     if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL))
 623     {
 624         posres_wrapper_lambda(fepvals, idef, &pbc_full, x, enerd, lambda, fr);
 625
 626         if (idef->ilsort != ilsortNO_FE)
 627         {
 628             if (idef->ilsort != ilsortFE_SORTED)
 629             {
 630                 gmx_incons("The bonded interactions are not sorted for free energy");
 631             }
 632             for (int i = 0; i < enerd->n_lambda; i++)
 633             {
 634                 real lam_i[efptNR];
 635
 636                 reset_foreign_enerdata(enerd);
 637                 for (int j = 0; j < efptNR; j++)
 638                 {
 639                     lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);
 640                 }
 641                 calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,
 642                                    fcd, global_atom_index);
 643                 sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
 644                 enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
 645             }
 646         }
 647     }
 648     debug_gmx();
 649
 650     wallcycle_sub_stop(wcycle, ewcsLISTED);
 651 }