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, 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     adress_set_kernel_flags(n_ex, n_hyb, n_cg, mdatoms);
 233
 234
 235 }
 236 void update_adress_weights_atom_per_atom(
 237                             int                  cg0,
 238                           int                  cg1,
 239                           t_block *            cgs,
 240                           rvec                 x[],
 241                           t_forcerec *         fr,
 242                           t_mdatoms *          mdatoms,
 243                           t_pbc *              pbc)
 244 {
 245     int            icg,k,k0,k1,d;
 246     real           nrcg,inv_ncg,mtot,inv_mtot;
 247     atom_id *      cgindex;
 248     rvec           ix;
 249     int            adresstype;
 250     real           adressr,adressw;
 251     rvec *         ref;
 252     real *         massT;
 253     real *         wf;
 254
 255
 256     int n_hyb, n_ex, n_cg;
 257
 258     n_hyb=0;
 259     n_cg=0;
 260     n_ex=0;
 261
 262     adresstype         = fr->adress_type;
 263     adressr            = fr->adress_ex_width;
 264     adressw            = fr->adress_hy_width;
 265     massT              = mdatoms->massT;
 266     wf                 = mdatoms->wf;
 267     ref                = &(fr->adress_refs);
 268
 269     cgindex = cgs->index;
 270
 271     /* Weighting function is determined for each atom individually.
 272      * This is an approximation
 273      * as in the theory requires an interpolation based on the center of masses.
 274      * Should be used with caution */
 275
 276     for (icg = cg0; (icg < cg1); icg++) {
 277         k0 = cgindex[icg];
 278         k1 = cgindex[icg + 1];
 279         nrcg = k1 - k0;
 280
 281         for (k = (k0); (k < k1); k++) {
 282             wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
 283             if (wf[k] == 0) {
 284                 n_cg++;
 285             } else if (wf[k] == 1) {
 286                 n_ex++;
 287             } else {
 288                 n_hyb++;
 289             }
 290         }
 291
 292     }
 293     adress_set_kernel_flags(n_ex, n_hyb, n_cg, mdatoms);
 294 }
 295
 296 void
 297 update_adress_weights_cog(t_iparams            ip[],
 298                           t_ilist              ilist[],
 299                           rvec                 x[],
 300                           t_forcerec *         fr,
 301                           t_mdatoms *          mdatoms,
 302                           t_pbc *              pbc)
 303 {
 304     int            i,j,k,nr,nra,inc;
 305     int            ftype,adresstype;
 306     t_iatom        avsite,ai,aj,ak,al;
 307     t_iatom *      ia;
 308     real           adressr,adressw;
 309     rvec *         ref;
 310     real *         wf;
 311     int            n_hyb, n_ex, n_cg;
 312
 313     adresstype         = fr->adress_type;
 314     adressr            = fr->adress_ex_width;
 315     adressw            = fr->adress_hy_width;
 316     wf                 = mdatoms->wf;
 317     ref                = &(fr->adress_refs);
 318
 319
 320     n_hyb=0;
 321     n_cg=0;
 322     n_ex=0;
 323
 324
 325     /* Since this is center of geometry AdResS, we know the vsite
 326      * is in the same charge group node as the constructing atoms.
 327      * Loop over vsite types, calculate the weight of the vsite,
 328      * then assign that weight to the constructing atoms. */
 329
 330     for(ftype=0; (ftype<F_NRE); ftype++)
 331     {
 332         if (interaction_function[ftype].flags & IF_VSITE)
 333         {
 334             nra    = interaction_function[ftype].nratoms;
 335             nr     = ilist[ftype].nr;
 336             ia     = ilist[ftype].iatoms;
 337
 338             for(i=0; (i<nr); )
 339             {
 340                 /* The vsite and first constructing atom */
 341                 avsite     = ia[1];
 342                 ai         = ia[2];
 343                 wf[avsite] = adress_weight(x[avsite],adresstype,adressr,adressw,ref,pbc,fr);
 344                 wf[ai]     = wf[avsite];
 345
 346                 if (wf[ai]  == 0) {
 347                     n_cg++;
 348                 } else if (wf[ai]  == 1) {
 349                     n_ex++;
 350                 } else {
 351                     n_hyb++;
 352                 }
 353
 354                 /* Assign the vsite wf to rest of constructing atoms depending on type */
 355                 inc = nra+1;
 356                 switch (ftype) {
 357                 case F_VSITE2:
 358                     aj     = ia[3];
 359                     wf[aj] = wf[avsite];
 360                     break;
 361                 case F_VSITE3:
 362                     aj     = ia[3];
 363                     wf[aj] = wf[avsite];
 364                     ak     = ia[4];
 365                     wf[ak] = wf[avsite];
 366                     break;
 367                 case F_VSITE3FD:
 368                     aj     = ia[3];
 369                     wf[aj] = wf[avsite];
 370                     ak     = ia[4];
 371                     wf[ak] = wf[avsite];
 372                     break;
 373                 case F_VSITE3FAD:
 374                     aj     = ia[3];
 375                     wf[aj] = wf[avsite];
 376                     ak     = ia[4];
 377                     wf[ak] = wf[avsite];
 378                     break;
 379                 case F_VSITE3OUT:
 380                     aj     = ia[3];
 381                     wf[aj] = wf[avsite];
 382                     ak     = ia[4];
 383                     wf[ak] = wf[avsite];
 384                     break;
 385                 case F_VSITE4FD:
 386                     aj     = ia[3];
 387                     wf[aj] = wf[avsite];
 388                     ak     = ia[4];
 389                     wf[ak] = wf[avsite];
 390                     al     = ia[5];
 391                     wf[al] = wf[avsite];
 392                     break;
 393                 case F_VSITE4FDN:
 394                     aj     = ia[3];
 395                     wf[aj] = wf[avsite];
 396                     ak     = ia[4];
 397                     wf[ak] = wf[avsite];
 398                     al     = ia[5];
 399                     wf[al] = wf[avsite];
 400                     break;
 401                 case F_VSITEN:
 402                     inc    = 3*ip[ia[0]].vsiten.n;
 403                     for(j=3; j<inc; j+=3)
 404                     {
 405                         ai = ia[j+2];
 406                         wf[ai] = wf[avsite];
 407                     }
 408                     break;
 409                 default:
 410                     gmx_fatal(FARGS,"No such vsite type %d in %s, line %d",
 411                               ftype,__FILE__,__LINE__);
 412                 }
 413
 414                 /* Increment loop variables */
 415                 i  += inc;
 416                 ia += inc;
 417             }
 418         }
 419     }
 420
 421     adress_set_kernel_flags(n_ex, n_hyb, n_cg, mdatoms);
 422 }
 423
 424 void
 425 update_adress_weights_atom(int                  cg0,
 426                            int                  cg1,
 427                            t_block *            cgs,
 428                            rvec                 x[],
 429                            t_forcerec *         fr,
 430                            t_mdatoms *          mdatoms,
 431                            t_pbc *              pbc)
 432 {
 433     int            icg,k,k0,k1;
 434     atom_id *      cgindex;
 435     int            adresstype;
 436     real           adressr,adressw;
 437     rvec *         ref;
 438     real *         massT;
 439     real *         wf;
 440
 441     adresstype         = fr->adress_type;
 442     adressr            = fr->adress_ex_width;
 443     adressw            = fr->adress_hy_width;
 444     massT              = mdatoms->massT;
 445     wf                 = mdatoms->wf;
 446     ref                = &(fr->adress_refs);
 447     cgindex            = cgs->index;
 448
 449     /* Only use first atom in charge group.
 450      * We still can't be sure that the vsite and constructing
 451      * atoms are on the same processor, so we must calculate
 452      * in the same way as com adress. */
 453
 454     for(icg=cg0; (icg<cg1); icg++)
 455     {
 456         k0      = cgindex[icg];
 457         k1      = cgindex[icg+1];
 458         wf[k0] = adress_weight(x[k0],adresstype,adressr,adressw,ref,pbc,fr);
 459
 460         /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
 461         for(k=(k0+1); (k<k1); k++)
 462         {
 463             wf[k] = wf[k0];
 464         }
 465     }
 466 }
 467
 468 void adress_set_kernel_flags(int n_ex, int n_hyb, int n_cg, t_mdatoms * mdatoms){
 469
 470     /* With domain decomposition we can check weather a cpu calculates only
 471      * coarse-grained or explicit interactions. If so we use standard gromacs kernels
 472      * on this proc. See also nonbonded.c */
 473
 474     if (n_hyb ==0 && n_ex == 0){
 475      /* all particles on this proc are coarse-grained, use standard gromacs kernels */
 476         if (!mdatoms->purecg){
 477             mdatoms->purecg = TRUE;
 478            if (debug) fprintf (debug, "adress.c: pure cg kernels on this proc\n");
 479         }
 480     }
 481     else
 482     {
 483         if (mdatoms->purecg){
 484          /* now this processor has hybrid particles again, call the hybrid kernels */
 485             mdatoms->purecg = FALSE;
 486         }
 487     }
 488
 489     if (n_hyb ==0 && n_cg == 0){
 490     /* all particles on this proc are atomistic, use standard gromacs kernels */
 491         if (!mdatoms->pureex){
 492              mdatoms->pureex = TRUE;
 493              if (debug) fprintf (debug, "adress.c: pure ex kernels on this proc\n");
 494         }
 495     }
 496     else
 497     {
 498         if (mdatoms->pureex){
 499             mdatoms->pureex = FALSE;
 500         }
 501     }
 502 }
 503
 504 void
 505 adress_thermo_force(int                  start,
 506                     int                  homenr,
 507                     t_block *            cgs,
 508                     rvec                 x[],
 509                     rvec                 f[],
 510                     t_forcerec *         fr,
 511                     t_mdatoms *          mdatoms,
 512                     t_pbc *              pbc)
 513 {
 514     int              iatom,n0,nnn,nrcg, i;
 515     int              adresstype;
 516     real             adressw, adressr;
 517     atom_id *        cgindex;
 518     unsigned short * ptype;
 519     rvec *           ref;
 520     real *           wf;
 521     real             tabscale;
 522     real *           ATFtab;
 523     rvec             dr;
 524     real             w,wsq,wmin1,wmin1sq,wp,wt,rinv, sqr_dl, dl;
 525     real             eps,eps2,F,Geps,Heps2,Fp,dmu_dwp,dwp_dr,fscal;
 526
 527     adresstype       = fr->adress_type;
 528     adressw          = fr->adress_hy_width;
 529     adressr           = fr->adress_ex_width;
 530     cgindex          = cgs->index;
 531     ptype            = mdatoms->ptype;
 532     ref              = &(fr->adress_refs);
 533     wf               = mdatoms->wf;
 534
 535     for(iatom=start; (iatom<start+homenr); iatom++)
 536     {
 537         if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
 538         {
 539             if (ptype[iatom] == eptVSite)
 540             {
 541                 w    = wf[iatom];
 542                 /* is it hybrid or apply the thermodynamics force everywhere?*/
 543                 if ( mdatoms->tf_table_index[iatom] != NO_TF_TABLE)
 544                 {
 545                     if (fr->n_adress_tf_grps > 0 ){
 546                         /* multi component tf is on, select the right table */
 547                         ATFtab = fr->atf_tabs[mdatoms->tf_table_index[iatom]].data;
 548                         tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
 549                     }
 550                     else {
 551                     /* just on component*/
 552                         ATFtab = fr->atf_tabs[DEFAULT_TF_TABLE].data;
 553                         tabscale = fr->atf_tabs[DEFAULT_TF_TABLE].scale;
 554                     }
 555
 556                     fscal            = 0;
 557                     if (pbc)
 558                     {
 559                         pbc_dx(pbc,(*ref),x[iatom],dr);
 560                     }
 561                     else
 562                     {
 563                         rvec_sub((*ref),x[iatom],dr);
 564                     }
 565
 566
 567
 568
 569                     /* calculate distace to adress center again */
 570                     sqr_dl =0.0;
 571                     switch(adresstype)
 572                     {
 573                     case eAdressXSplit:
 574                         /* plane through center of ref, varies in x direction */
 575                         sqr_dl         = dr[0]*dr[0];
 576                         rinv             = gmx_invsqrt(dr[0]*dr[0]);
 577                         break;
 578                     case eAdressSphere:
 579                         /* point at center of ref, assuming cubic geometry */
 580                         for(i=0;i<3;i++){
 581                             sqr_dl    += dr[i]*dr[i];
 582                         }
 583                         rinv             = gmx_invsqrt(sqr_dl);
 584                         break;
 585                     default:
 586                         /* This case should not happen */
 587                         rinv = 0.0;
 588                     }
 589
 590                     dl=sqrt(sqr_dl);
 591                     /* table origin is adress center */
 592                     wt               = dl*tabscale;
 593                     n0               = wt;
 594                     eps              = wt-n0;
 595                     eps2             = eps*eps;
 596                     nnn              = 4*n0;
 597                     F                = ATFtab[nnn+1];
 598                     Geps             = eps*ATFtab[nnn+2];
 599                     Heps2            = eps2*ATFtab[nnn+3];
 600                     Fp               = F+Geps+Heps2;
 601                     F                = (Fp+Geps+2.0*Heps2)*tabscale;
 602
 603                     fscal            = F*rinv;
 604
 605                     f[iatom][0]        += fscal*dr[0];
 606                     if (adresstype != eAdressXSplit)
 607                     {
 608                         f[iatom][1]    += fscal*dr[1];
 609                         f[iatom][2]    += fscal*dr[2];
 610                     }
 611                 }
 612             }
 613         }
 614     }
 615 }
 616
 617 gmx_bool egp_explicit(t_forcerec *   fr, int egp_nr)
 618 {
 619     return fr->adress_group_explicit[egp_nr];
 620 }
 621
 622 gmx_bool egp_coarsegrained(t_forcerec *   fr, int egp_nr)
 623 {
 624    return !fr->adress_group_explicit[egp_nr];
 625 }