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