src/gromacs/mdlib/adress.cpp

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