src/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) 2012,2013, by the GROMACS development team, led by
   6  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   7  * others, as listed in the AUTHORS file in the top-level source
   8  * 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 "maths.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             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 *               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){ n_cg++;}
 171             else if (wf[k0]==1){ n_ex++;}
 172             else {n_hyb++;}
 173         }
 174         else
 175         {
 176             mtot = 0.0;
 177             for(k=k0; (k<k1); k++)
 178             {
 179                 mtot += massT[k];
 180             }
 181             if (mtot > 0.0)
 182             {
 183                 inv_mtot = 1.0/mtot;
 184
 185                 clear_rvec(ix);
 186                 for(k=k0; (k<k1); k++)
 187                 {
 188                     for(d=0; (d<DIM); d++)
 189                     {
 190                         ix[d] += x[k][d]*massT[k];
 191                     }
 192                 }
 193                 for(d=0; (d<DIM); d++)
 194                 {
 195                     ix[d] *= inv_mtot;
 196                 }
 197             }
 198             /* Calculate the center of gravity if the charge group mtot=0 (only vsites) */
 199             else
 200             {
 201                 inv_ncg = 1.0/nrcg;
 202
 203                 clear_rvec(ix);
 204                 for(k=k0; (k<k1); k++)
 205                 {
 206                     for(d=0; (d<DIM); d++)
 207                     {
 208                         ix[d] += x[k][d];
 209                     }
 210                 }
 211                 for(d=0; (d<DIM); d++)
 212                 {
 213                     ix[d] *= inv_ncg;
 214                 }
 215             }
 216
 217             /* Set wf of all atoms in charge group equal to wf of com */
 218             wf[k0] = adress_weight(ix,adresstype,adressr,adressw,ref,pbc, fr);
 219
 220             if (wf[k0]==0){ n_cg++;}
 221             else if (wf[k0]==1){ n_ex++;}
 222             else {n_hyb++;}
 223
 224             for(k=(k0+1); (k<k1); k++)
 225             {
 226                 wf[k] = wf[k0];
 227             }
 228         }
 229     }
 230 }
 231
 232 void update_adress_weights_atom_per_atom(
 233                             int                  cg0,
 234                           int                  cg1,
 235                           t_block *            cgs,
 236                           rvec                 x[],
 237                           t_forcerec *         fr,
 238                           t_mdatoms *          mdatoms,
 239                           t_pbc *              pbc)
 240 {
 241     int            icg,k,k0,k1,d;
 242     real           nrcg,inv_ncg,mtot,inv_mtot;
 243     atom_id *      cgindex;
 244     rvec           ix;
 245     int            adresstype;
 246     real           adressr,adressw;
 247     rvec *         ref;
 248     real *         massT;
 249     real *         wf;
 250
 251
 252     int n_hyb, n_ex, n_cg;
 253
 254     n_hyb=0;
 255     n_cg=0;
 256     n_ex=0;
 257
 258     adresstype         = fr->adress_type;
 259     adressr            = fr->adress_ex_width;
 260     adressw            = fr->adress_hy_width;
 261     massT              = mdatoms->massT;
 262     wf                 = mdatoms->wf;
 263     ref                = &(fr->adress_refs);
 264
 265     cgindex = cgs->index;
 266
 267     /* Weighting function is determined for each atom individually.
 268      * This is an approximation
 269      * as in the theory requires an interpolation based on the center of masses.
 270      * Should be used with caution */
 271
 272     for (icg = cg0; (icg < cg1); icg++) {
 273         k0 = cgindex[icg];
 274         k1 = cgindex[icg + 1];
 275         nrcg = k1 - k0;
 276
 277         for (k = (k0); (k < k1); k++) {
 278             wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
 279             if (wf[k] == 0) {
 280                 n_cg++;
 281             } else if (wf[k] == 1) {
 282                 n_ex++;
 283             } else {
 284                 n_hyb++;
 285             }
 286         }
 287
 288     }
 289 }
 290
 291 void
 292 update_adress_weights_cog(t_iparams            ip[],
 293                           t_ilist              ilist[],
 294                           rvec                 x[],
 295                           t_forcerec *         fr,
 296                           t_mdatoms *          mdatoms,
 297                           t_pbc *              pbc)
 298 {
 299     int            i,j,k,nr,nra,inc;
 300     int            ftype,adresstype;
 301     t_iatom        avsite,ai,aj,ak,al;
 302     t_iatom *      ia;
 303     real           adressr,adressw;
 304     rvec *         ref;
 305     real *         wf;
 306     int            n_hyb, n_ex, n_cg;
 307
 308     adresstype         = fr->adress_type;
 309     adressr            = fr->adress_ex_width;
 310     adressw            = fr->adress_hy_width;
 311     wf                 = mdatoms->wf;
 312     ref                = &(fr->adress_refs);
 313
 314
 315     n_hyb=0;
 316     n_cg=0;
 317     n_ex=0;
 318
 319
 320     /* Since this is center of geometry AdResS, we know the vsite
 321      * is in the same charge group node as the constructing atoms.
 322      * Loop over vsite types, calculate the weight of the vsite,
 323      * then assign that weight to the constructing atoms. */
 324
 325     for(ftype=0; (ftype<F_NRE); ftype++)
 326     {
 327         if (interaction_function[ftype].flags & IF_VSITE)
 328         {
 329             nra    = interaction_function[ftype].nratoms;
 330             nr     = ilist[ftype].nr;
 331             ia     = ilist[ftype].iatoms;
 332
 333             for(i=0; (i<nr); )
 334             {
 335                 /* The vsite and first constructing atom */
 336                 avsite     = ia[1];
 337                 ai         = ia[2];
 338                 wf[avsite] = adress_weight(x[avsite],adresstype,adressr,adressw,ref,pbc,fr);
 339                 wf[ai]     = wf[avsite];
 340
 341                 if (wf[ai]  == 0) {
 342                     n_cg++;
 343                 } else if (wf[ai]  == 1) {
 344                     n_ex++;
 345                 } else {
 346                     n_hyb++;
 347                 }
 348
 349                 /* Assign the vsite wf to rest of constructing atoms depending on type */
 350                 inc = nra+1;
 351                 switch (ftype) {
 352                 case F_VSITE2:
 353                     aj     = ia[3];
 354                     wf[aj] = wf[avsite];
 355                     break;
 356                 case F_VSITE3:
 357                     aj     = ia[3];
 358                     wf[aj] = wf[avsite];
 359                     ak     = ia[4];
 360                     wf[ak] = wf[avsite];
 361                     break;
 362                 case F_VSITE3FD:
 363                     aj     = ia[3];
 364                     wf[aj] = wf[avsite];
 365                     ak     = ia[4];
 366                     wf[ak] = wf[avsite];
 367                     break;
 368                 case F_VSITE3FAD:
 369                     aj     = ia[3];
 370                     wf[aj] = wf[avsite];
 371                     ak     = ia[4];
 372                     wf[ak] = wf[avsite];
 373                     break;
 374                 case F_VSITE3OUT:
 375                     aj     = ia[3];
 376                     wf[aj] = wf[avsite];
 377                     ak     = ia[4];
 378                     wf[ak] = wf[avsite];
 379                     break;
 380                 case F_VSITE4FD:
 381                     aj     = ia[3];
 382                     wf[aj] = wf[avsite];
 383                     ak     = ia[4];
 384                     wf[ak] = wf[avsite];
 385                     al     = ia[5];
 386                     wf[al] = wf[avsite];
 387                     break;
 388                 case F_VSITE4FDN:
 389                     aj     = ia[3];
 390                     wf[aj] = wf[avsite];
 391                     ak     = ia[4];
 392                     wf[ak] = wf[avsite];
 393                     al     = ia[5];
 394                     wf[al] = wf[avsite];
 395                     break;
 396                 case F_VSITEN:
 397                     inc    = 3*ip[ia[0]].vsiten.n;
 398                     for(j=3; j<inc; j+=3)
 399                     {
 400                         ai = ia[j+2];
 401                         wf[ai] = wf[avsite];
 402                     }
 403                     break;
 404                 default:
 405                     gmx_fatal(FARGS,"No such vsite type %d in %s, line %d",
 406                               ftype,__FILE__,__LINE__);
 407                 }
 408
 409                 /* Increment loop variables */
 410                 i  += inc;
 411                 ia += inc;
 412             }
 413         }
 414     }
 415 }
 416
 417 void
 418 update_adress_weights_atom(int                  cg0,
 419                            int                  cg1,
 420                            t_block *            cgs,
 421                            rvec                 x[],
 422                            t_forcerec *         fr,
 423                            t_mdatoms *          mdatoms,
 424                            t_pbc *              pbc)
 425 {
 426     int            icg,k,k0,k1;
 427     atom_id *      cgindex;
 428     int            adresstype;
 429     real           adressr,adressw;
 430     rvec *         ref;
 431     real *         massT;
 432     real *         wf;
 433
 434     adresstype         = fr->adress_type;
 435     adressr            = fr->adress_ex_width;
 436     adressw            = fr->adress_hy_width;
 437     massT              = mdatoms->massT;
 438     wf                 = mdatoms->wf;
 439     ref                = &(fr->adress_refs);
 440     cgindex            = cgs->index;
 441
 442     /* Only use first atom in charge group.
 443      * We still can't be sure that the vsite and constructing
 444      * atoms are on the same processor, so we must calculate
 445      * in the same way as com adress. */
 446
 447     for(icg=cg0; (icg<cg1); icg++)
 448     {
 449         k0      = cgindex[icg];
 450         k1      = cgindex[icg+1];
 451         wf[k0] = adress_weight(x[k0],adresstype,adressr,adressw,ref,pbc,fr);
 452
 453         /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
 454         for(k=(k0+1); (k<k1); k++)
 455         {
 456             wf[k] = wf[k0];
 457         }
 458     }
 459 }
 460
 461 void
 462 adress_thermo_force(int                  start,
 463                     int                  homenr,
 464                     t_block *            cgs,
 465                     rvec                 x[],
 466                     rvec                 f[],
 467                     t_forcerec *         fr,
 468                     t_mdatoms *          mdatoms,
 469                     t_pbc *              pbc)
 470 {
 471     int              iatom,n0,nnn,nrcg, i;
 472     int              adresstype;
 473     real             adressw, adressr;
 474     atom_id *        cgindex;
 475     unsigned short * ptype;
 476     rvec *           ref;
 477     real *           wf;
 478     real             tabscale;
 479     real *           ATFtab;
 480     rvec             dr;
 481     real             w,wsq,wmin1,wmin1sq,wp,wt,rinv, sqr_dl, dl;
 482     real             eps,eps2,F,Geps,Heps2,Fp,dmu_dwp,dwp_dr,fscal;
 483
 484     adresstype       = fr->adress_type;
 485     adressw          = fr->adress_hy_width;
 486     adressr           = fr->adress_ex_width;
 487     cgindex          = cgs->index;
 488     ptype            = mdatoms->ptype;
 489     ref              = &(fr->adress_refs);
 490     wf               = mdatoms->wf;
 491
 492     for(iatom=start; (iatom<start+homenr); iatom++)
 493     {
 494         if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
 495         {
 496             if (ptype[iatom] == eptVSite)
 497             {
 498                 w    = wf[iatom];
 499                 /* is it hybrid or apply the thermodynamics force everywhere?*/
 500                 if ( mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
 501                 {
 502                     if (fr->n_adress_tf_grps > 0 ){
 503                         /* multi component tf is on, select the right table */
 504                         ATFtab = fr->atf_tabs[mdatoms->tf_table_index[iatom]].data;
 505                         tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
 506                     }
 507                     else {
 508                     /* just on component*/
 509                         ATFtab = fr->atf_tabs[DEFAULT_TF_TABLE].data;
 510                         tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
 511                     }
 512
 513                     fscal            = 0;
 514                     if (pbc)
 515                     {
 516                         pbc_dx(pbc,(*ref),x[iatom],dr);
 517                     }
 518                     else
 519                     {
 520                         rvec_sub((*ref),x[iatom],dr);
 521                     }
 522
 523
 524
 525
 526                     /* calculate distace to adress center again */
 527                     sqr_dl =0.0;
 528                     switch(adresstype)
 529                     {
 530                     case eAdressXSplit:
 531                         /* plane through center of ref, varies in x direction */
 532                         sqr_dl         = dr[0]*dr[0];
 533                         rinv             = gmx_invsqrt(dr[0]*dr[0]);
 534                         break;
 535                     case eAdressSphere:
 536                         /* point at center of ref, assuming cubic geometry */
 537                         for(i=0;i<3;i++){
 538                             sqr_dl    += dr[i]*dr[i];
 539                         }
 540                         rinv             = gmx_invsqrt(sqr_dl);
 541                         break;
 542                     default:
 543                         /* This case should not happen */
 544                         rinv = 0.0;
 545                     }
 546
 547                     dl=sqrt(sqr_dl);
 548                     /* table origin is adress center */
 549                     wt               = dl*tabscale;
 550                     n0               = wt;
 551                     eps              = wt-n0;
 552                     eps2             = eps*eps;
 553                     nnn              = 4*n0;
 554                     F                = ATFtab[nnn+1];
 555                     Geps             = eps*ATFtab[nnn+2];
 556                     Heps2            = eps2*ATFtab[nnn+3];
 557                     Fp               = F+Geps+Heps2;
 558                     F                = (Fp+Geps+2.0*Heps2)*tabscale;
 559
 560                     fscal            = F*rinv;
 561
 562                     f[iatom][0]        += fscal*dr[0];
 563                     if (adresstype != eAdressXSplit)
 564                     {
 565                         f[iatom][1]    += fscal*dr[1];
 566                         f[iatom][2]    += fscal*dr[2];
 567                     }
 568                 }
 569             }
 570         }
 571     }
 572 }
 573
 574 gmx_bool egp_explicit(t_forcerec *   fr, int egp_nr)
 575 {
 576     return fr->adress_group_explicit[egp_nr];
 577 }
 578
 579 gmx_bool egp_coarsegrained(t_forcerec *   fr, int egp_nr)
 580 {
 581    return !fr->adress_group_explicit[egp_nr];
 582 }