src/gromacs/mdlib/adress.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
   5  * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36
  37 #include "adress.h"
  38 #include "gromacs/math/utilities.h"
  39 #include "pbc.h"
  40 #include "types/simple.h"
  41 #include "typedefs.h"
  42 #include "vec.h"
  43
  44 real
  45 adress_weight(rvec                 x,
  46               int                  adresstype,
  47               real                 adressr,
  48               real                 adressw,
  49               rvec      *          ref,
  50               t_pbc      *         pbc,
  51               t_forcerec *         fr )
  52 {
  53     int  i;
  54     real l2 = adressr+adressw;
  55     real sqr_dl, dl;
  56     real tmp;
  57     rvec dx;
  58
  59     sqr_dl = 0.0;
  60
  61     if (pbc)
  62     {
  63         pbc_dx(pbc, (*ref), x, dx);
  64     }
  65     else
  66     {
  67         rvec_sub((*ref), x, dx);
  68     }
  69
  70     switch (adresstype)
  71     {
  72         case eAdressOff:
  73             /* default to explicit simulation */
  74             return 1;
  75         case eAdressConst:
  76             /* constant value for weighting function = adressw */
  77             return fr->adress_const_wf;
  78         case eAdressXSplit:
  79             /* plane through center of ref, varies in x direction */
  80             sqr_dl         = dx[0]*dx[0];
  81             break;
  82         case eAdressSphere:
  83             /* point at center of ref, assuming cubic geometry */
  84             for (i = 0; i < 3; i++)
  85             {
  86                 sqr_dl    += dx[i]*dx[i];
  87             }
  88             break;
  89         default:
  90             /* default to explicit simulation */
  91             return 1;
  92     }
  93
  94     dl = sqrt(sqr_dl);
  95
  96     /* molecule is coarse grained */
  97     if (dl > l2)
  98     {
  99         return 0;
 100     }
 101     /* molecule is explicit */
 102     else if (dl < adressr)
 103     {
 104         return 1;
 105     }
 106     /* hybrid region */
 107     else
 108     {
 109         tmp = cos((dl-adressr)*M_PI/2/adressw);
 110         return tmp*tmp;
 111     }
 112 }
 113
 114 void
 115 update_adress_weights_com(FILE gmx_unused    * fplog,
 116                           int                  cg0,
 117                           int                  cg1,
 118                           t_block *            cgs,
 119                           rvec                 x[],
 120                           t_forcerec *         fr,
 121                           t_mdatoms *          mdatoms,
 122                           t_pbc *              pbc)
 123 {
 124     int            icg, k, k0, k1, d;
 125     real           nrcg, inv_ncg, mtot, inv_mtot;
 126     atom_id *      cgindex;
 127     rvec           ix;
 128     int            adresstype;
 129     real           adressr, adressw;
 130     rvec *         ref;
 131     real *         massT;
 132     real *         wf;
 133
 134
 135     int n_hyb, n_ex, n_cg;
 136
 137     n_hyb = 0;
 138     n_cg  = 0;
 139     n_ex  = 0;
 140
 141     adresstype         = fr->adress_type;
 142     adressr            = fr->adress_ex_width;
 143     adressw            = fr->adress_hy_width;
 144     massT              = mdatoms->massT;
 145     wf                 = mdatoms->wf;
 146     ref                = &(fr->adress_refs);
 147
 148
 149     /* Since this is center of mass AdResS, the vsite is not guaranteed
 150      * to be on the same node as the constructing atoms.  Therefore we
 151      * loop over the charge groups, calculate their center of mass,
 152      * then use this to calculate wf for each atom.  This wastes vsite
 153      * construction, but it's the only way to assure that the explicit
 154      * atoms have the same wf as their vsite. */
 155
 156 #ifdef DEBUG
 157     fprintf(fplog, "Calculating center of mass for charge groups %d to %d\n",
 158             cg0, cg1);
 159 #endif
 160     cgindex = cgs->index;
 161
 162     /* Compute the center of mass for all charge groups */
 163     for (icg = cg0; (icg < cg1); icg++)
 164     {
 165         k0      = cgindex[icg];
 166         k1      = cgindex[icg+1];
 167         nrcg    = k1-k0;
 168         if (nrcg == 1)
 169         {
 170             wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
 171             if (wf[k0] == 0)
 172             {
 173                 n_cg++;
 174             }
 175             else if (wf[k0] == 1)
 176             {
 177                 n_ex++;
 178             }
 179             else
 180             {
 181                 n_hyb++;
 182             }
 183         }
 184         else
 185         {
 186             mtot = 0.0;
 187             for (k = k0; (k < k1); k++)
 188             {
 189                 mtot += massT[k];
 190             }
 191             if (mtot > 0.0)
 192             {
 193                 inv_mtot = 1.0/mtot;
 194
 195                 clear_rvec(ix);
 196                 for (k = k0; (k < k1); k++)
 197                 {
 198                     for (d = 0; (d < DIM); d++)
 199                     {
 200                         ix[d] += x[k][d]*massT[k];
 201                     }
 202                 }
 203                 for (d = 0; (d < DIM); d++)
 204                 {
 205                     ix[d] *= inv_mtot;
 206                 }
 207             }
 208             /* Calculate the center of gravity if the charge group mtot=0 (only vsites) */
 209             else
 210             {
 211                 inv_ncg = 1.0/nrcg;
 212
 213                 clear_rvec(ix);
 214                 for (k = k0; (k < k1); k++)
 215                 {
 216                     for (d = 0; (d < DIM); d++)
 217                     {
 218                         ix[d] += x[k][d];
 219                     }
 220                 }
 221                 for (d = 0; (d < DIM); d++)
 222                 {
 223                     ix[d] *= inv_ncg;
 224                 }
 225             }
 226
 227             /* Set wf of all atoms in charge group equal to wf of com */
 228             wf[k0] = adress_weight(ix, adresstype, adressr, adressw, ref, pbc, fr);
 229
 230             if (wf[k0] == 0)
 231             {
 232                 n_cg++;
 233             }
 234             else if (wf[k0] == 1)
 235             {
 236                 n_ex++;
 237             }
 238             else
 239             {
 240                 n_hyb++;
 241             }
 242
 243             for (k = (k0+1); (k < k1); k++)
 244             {
 245                 wf[k] = wf[k0];
 246             }
 247         }
 248     }
 249 }
 250
 251 void update_adress_weights_atom_per_atom(
 252         int                  cg0,
 253         int                  cg1,
 254         t_block *            cgs,
 255         rvec                 x[],
 256         t_forcerec *         fr,
 257         t_mdatoms *          mdatoms,
 258         t_pbc *              pbc)
 259 {
 260     int            icg, k, k0, k1, d;
 261     real           nrcg, inv_ncg, mtot, inv_mtot;
 262     atom_id *      cgindex;
 263     rvec           ix;
 264     int            adresstype;
 265     real           adressr, adressw;
 266     rvec *         ref;
 267     real *         massT;
 268     real *         wf;
 269
 270
 271     int n_hyb, n_ex, n_cg;
 272
 273     n_hyb = 0;
 274     n_cg  = 0;
 275     n_ex  = 0;
 276
 277     adresstype         = fr->adress_type;
 278     adressr            = fr->adress_ex_width;
 279     adressw            = fr->adress_hy_width;
 280     massT              = mdatoms->massT;
 281     wf                 = mdatoms->wf;
 282     ref                = &(fr->adress_refs);
 283
 284     cgindex = cgs->index;
 285
 286     /* Weighting function is determined for each atom individually.
 287      * This is an approximation
 288      * as in the theory requires an interpolation based on the center of masses.
 289      * Should be used with caution */
 290
 291     for (icg = cg0; (icg < cg1); icg++)
 292     {
 293         k0   = cgindex[icg];
 294         k1   = cgindex[icg + 1];
 295         nrcg = k1 - k0;
 296
 297         for (k = (k0); (k < k1); k++)
 298         {
 299             wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
 300             if (wf[k] == 0)
 301             {
 302                 n_cg++;
 303             }
 304             else if (wf[k] == 1)
 305             {
 306                 n_ex++;
 307             }
 308             else
 309             {
 310                 n_hyb++;
 311             }
 312         }
 313
 314     }
 315 }
 316
 317 void
 318 update_adress_weights_cog(t_iparams            ip[],
 319                           t_ilist              ilist[],
 320                           rvec                 x[],
 321                           t_forcerec *         fr,
 322                           t_mdatoms *          mdatoms,
 323                           t_pbc *              pbc)
 324 {
 325     int            i, j, k, nr, nra, inc;
 326     int            ftype, adresstype;
 327     t_iatom        avsite, ai, aj, ak, al;
 328     t_iatom *      ia;
 329     real           adressr, adressw;
 330     rvec *         ref;
 331     real *         wf;
 332     int            n_hyb, n_ex, n_cg;
 333
 334     adresstype         = fr->adress_type;
 335     adressr            = fr->adress_ex_width;
 336     adressw            = fr->adress_hy_width;
 337     wf                 = mdatoms->wf;
 338     ref                = &(fr->adress_refs);
 339
 340
 341     n_hyb = 0;
 342     n_cg  = 0;
 343     n_ex  = 0;
 344
 345
 346     /* Since this is center of geometry AdResS, we know the vsite
 347      * is in the same charge group node as the constructing atoms.
 348      * Loop over vsite types, calculate the weight of the vsite,
 349      * then assign that weight to the constructing atoms. */
 350
 351     for (ftype = 0; (ftype < F_NRE); ftype++)
 352     {
 353         if (interaction_function[ftype].flags & IF_VSITE)
 354         {
 355             nra    = interaction_function[ftype].nratoms;
 356             nr     = ilist[ftype].nr;
 357             ia     = ilist[ftype].iatoms;
 358
 359             for (i = 0; (i < nr); )
 360             {
 361                 /* The vsite and first constructing atom */
 362                 avsite     = ia[1];
 363                 ai         = ia[2];
 364                 wf[avsite] = adress_weight(x[avsite], adresstype, adressr, adressw, ref, pbc, fr);
 365                 wf[ai]     = wf[avsite];
 366
 367                 if (wf[ai]  == 0)
 368                 {
 369                     n_cg++;
 370                 }
 371                 else if (wf[ai]  == 1)
 372                 {
 373                     n_ex++;
 374                 }
 375                 else
 376                 {
 377                     n_hyb++;
 378                 }
 379
 380                 /* Assign the vsite wf to rest of constructing atoms depending on type */
 381                 inc = nra+1;
 382                 switch (ftype)
 383                 {
 384                     case F_VSITE2:
 385                         aj     = ia[3];
 386                         wf[aj] = wf[avsite];
 387                         break;
 388                     case F_VSITE3:
 389                         aj     = ia[3];
 390                         wf[aj] = wf[avsite];
 391                         ak     = ia[4];
 392                         wf[ak] = wf[avsite];
 393                         break;
 394                     case F_VSITE3FD:
 395                         aj     = ia[3];
 396                         wf[aj] = wf[avsite];
 397                         ak     = ia[4];
 398                         wf[ak] = wf[avsite];
 399                         break;
 400                     case F_VSITE3FAD:
 401                         aj     = ia[3];
 402                         wf[aj] = wf[avsite];
 403                         ak     = ia[4];
 404                         wf[ak] = wf[avsite];
 405                         break;
 406                     case F_VSITE3OUT:
 407                         aj     = ia[3];
 408                         wf[aj] = wf[avsite];
 409                         ak     = ia[4];
 410                         wf[ak] = wf[avsite];
 411                         break;
 412                     case F_VSITE4FD:
 413                         aj     = ia[3];
 414                         wf[aj] = wf[avsite];
 415                         ak     = ia[4];
 416                         wf[ak] = wf[avsite];
 417                         al     = ia[5];
 418                         wf[al] = wf[avsite];
 419                         break;
 420                     case F_VSITE4FDN:
 421                         aj     = ia[3];
 422                         wf[aj] = wf[avsite];
 423                         ak     = ia[4];
 424                         wf[ak] = wf[avsite];
 425                         al     = ia[5];
 426                         wf[al] = wf[avsite];
 427                         break;
 428                     case F_VSITEN:
 429                         inc    = 3*ip[ia[0]].vsiten.n;
 430                         for (j = 3; j < inc; j += 3)
 431                         {
 432                             ai     = ia[j+2];
 433                             wf[ai] = wf[avsite];
 434                         }
 435                         break;
 436                     default:
 437                         gmx_fatal(FARGS, "No such vsite type %d in %s, line %d",
 438                                   ftype, __FILE__, __LINE__);
 439                 }
 440
 441                 /* Increment loop variables */
 442                 i  += inc;
 443                 ia += inc;
 444             }
 445         }
 446     }
 447 }
 448
 449 void
 450 update_adress_weights_atom(int                  cg0,
 451                            int                  cg1,
 452                            t_block *            cgs,
 453                            rvec                 x[],
 454                            t_forcerec *         fr,
 455                            t_mdatoms *          mdatoms,
 456                            t_pbc *              pbc)
 457 {
 458     int            icg, k, k0, k1;
 459     atom_id *      cgindex;
 460     int            adresstype;
 461     real           adressr, adressw;
 462     rvec *         ref;
 463     real *         massT;
 464     real *         wf;
 465
 466     adresstype         = fr->adress_type;
 467     adressr            = fr->adress_ex_width;
 468     adressw            = fr->adress_hy_width;
 469     massT              = mdatoms->massT;
 470     wf                 = mdatoms->wf;
 471     ref                = &(fr->adress_refs);
 472     cgindex            = cgs->index;
 473
 474     /* Only use first atom in charge group.
 475      * We still can't be sure that the vsite and constructing
 476      * atoms are on the same processor, so we must calculate
 477      * in the same way as com adress. */
 478
 479     for (icg = cg0; (icg < cg1); icg++)
 480     {
 481         k0      = cgindex[icg];
 482         k1      = cgindex[icg+1];
 483         wf[k0]  = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
 484
 485         /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
 486         for (k = (k0+1); (k < k1); k++)
 487         {
 488             wf[k] = wf[k0];
 489         }
 490     }
 491 }
 492
 493 void
 494 adress_thermo_force(int                  start,
 495                     int                  homenr,
 496                     t_block *            cgs,
 497                     rvec                 x[],
 498                     rvec                 f[],
 499                     t_forcerec *         fr,
 500                     t_mdatoms *          mdatoms,
 501                     t_pbc *              pbc)
 502 {
 503     int              iatom, n0, nnn, nrcg, i;
 504     int              adresstype;
 505     real             adressw, adressr;
 506     atom_id *        cgindex;
 507     unsigned short * ptype;
 508     rvec *           ref;
 509     real *           wf;
 510     real             tabscale;
 511     real *           ATFtab;
 512     rvec             dr;
 513     real             w, wsq, wmin1, wmin1sq, wp, wt, rinv, sqr_dl, dl;
 514     real             eps, eps2, F, Geps, Heps2, Fp, dmu_dwp, dwp_dr, fscal;
 515
 516     adresstype        = fr->adress_type;
 517     adressw           = fr->adress_hy_width;
 518     adressr           = fr->adress_ex_width;
 519     cgindex           = cgs->index;
 520     ptype             = mdatoms->ptype;
 521     ref               = &(fr->adress_refs);
 522     wf                = mdatoms->wf;
 523
 524     for (iatom = start; (iatom < start+homenr); iatom++)
 525     {
 526         if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
 527         {
 528             if (ptype[iatom] == eptVSite)
 529             {
 530                 w    = wf[iatom];
 531                 /* is it hybrid or apply the thermodynamics force everywhere?*/
 532                 if (mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
 533                 {
 534                     if (fr->n_adress_tf_grps > 0)
 535                     {
 536                         /* multi component tf is on, select the right table */
 537                         ATFtab   = fr->atf_tabs[mdatoms->tf_table_index[iatom]].data;
 538                         tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
 539                     }
 540                     else
 541                     {
 542                         /* just on component*/
 543                         ATFtab   = fr->atf_tabs[DEFAULT_TF_TABLE].data;
 544                         tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
 545                     }
 546
 547                     fscal            = 0;
 548                     if (pbc)
 549                     {
 550                         pbc_dx(pbc, (*ref), x[iatom], dr);
 551                     }
 552                     else
 553                     {
 554                         rvec_sub((*ref), x[iatom], dr);
 555                     }
 556
 557
 558
 559
 560                     /* calculate distace to adress center again */
 561                     sqr_dl = 0.0;
 562                     switch (adresstype)
 563                     {
 564                         case eAdressXSplit:
 565                             /* plane through center of ref, varies in x direction */
 566                             sqr_dl           = dr[0]*dr[0];
 567                             rinv             = gmx_invsqrt(dr[0]*dr[0]);
 568                             break;
 569                         case eAdressSphere:
 570                             /* point at center of ref, assuming cubic geometry */
 571                             for (i = 0; i < 3; i++)
 572                             {
 573                                 sqr_dl    += dr[i]*dr[i];
 574                             }
 575                             rinv             = gmx_invsqrt(sqr_dl);
 576                             break;
 577                         default:
 578                             /* This case should not happen */
 579                             rinv = 0.0;
 580                     }
 581
 582                     dl = sqrt(sqr_dl);
 583                     /* table origin is adress center */
 584                     wt               = dl*tabscale;
 585                     n0               = wt;
 586                     eps              = wt-n0;
 587                     eps2             = eps*eps;
 588                     nnn              = 4*n0;
 589                     F                = ATFtab[nnn+1];
 590                     Geps             = eps*ATFtab[nnn+2];
 591                     Heps2            = eps2*ATFtab[nnn+3];
 592                     Fp               = F+Geps+Heps2;
 593                     F                = (Fp+Geps+2.0*Heps2)*tabscale;
 594
 595                     fscal            = F*rinv;
 596
 597                     f[iatom][0]        += fscal*dr[0];
 598                     if (adresstype != eAdressXSplit)
 599                     {
 600                         f[iatom][1]    += fscal*dr[1];
 601                         f[iatom][2]    += fscal*dr[2];
 602                     }
 603                 }
 604             }
 605         }
 606     }
 607 }
 608
 609 gmx_bool egp_explicit(t_forcerec *   fr, int egp_nr)
 610 {
 611     return fr->adress_group_explicit[egp_nr];
 612 }
 613
 614 gmx_bool egp_coarsegrained(t_forcerec *   fr, int egp_nr)
 615 {
 616     return !fr->adress_group_explicit[egp_nr];
 617 }