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