src/gromacs/mdlib/adress.c

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