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