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