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