src/mdlib/wall.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  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2008, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  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  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <string.h>
  42 #include "maths.h"
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "macros.h"
  46 #include "smalloc.h"
  47 #include "force.h"
  48 #include "filenm.h"
  49 #include "nrnb.h"
  50 #include "vec.h"
  51
  52 void make_wall_tables(FILE *fplog,const output_env_t oenv,
  53                       const t_inputrec *ir,const char *tabfn,
  54                       const gmx_groups_t *groups,
  55                       t_forcerec *fr)
  56 {
  57   int  w,negp_pp,egp,i,j;
  58   int  *nm_ind;
  59   char buf[STRLEN];
  60   t_forcetable *tab;
  61
  62   negp_pp = ir->opts.ngener - ir->nwall;
  63   nm_ind  = groups->grps[egcENER].nm_ind;
  64
  65   if (fplog) {
  66     fprintf(fplog,"Reading user tables for %d energy groups with %d walls\n",
  67             negp_pp,ir->nwall);
  68   }
  69
  70   snew(fr->wall_tab,ir->nwall);
  71   for(w=0; w<ir->nwall; w++) {
  72     snew(fr->wall_tab[w],negp_pp);
  73     for(egp=0; egp<negp_pp; egp++) {
  74       /* If the energy group pair is excluded, we don't need a table */
  75       if (!(fr->egp_flags[egp*ir->opts.ngener+negp_pp+w] & EGP_EXCL)) {
  76         tab = &fr->wall_tab[w][egp];
  77         sprintf(buf,"%s",tabfn);
  78         sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1,"_%s_%s.%s",
  79                 *groups->grpname[nm_ind[egp]],
  80                 *groups->grpname[nm_ind[negp_pp+w]],
  81                 ftp2ext(efXVG));
  82         *tab = make_tables(fplog,oenv,fr,FALSE,buf,0,GMX_MAKETABLES_FORCEUSER);
  83         /* Since wall have no charge, we can compress the table */
  84         for(i=0; i<=tab->n; i++)
  85           for(j=0; j<8; j++)
  86             tab->tab[8*i+j] = tab->tab[12*i+4+j];
  87       }
  88     }
  89   }
  90 }
  91
  92 static void wall_error(int a,rvec *x,real r)
  93 {
  94     gmx_fatal(FARGS,
  95               "An atom is beyond the wall: coordinates %f %f %f, distance %f\n"
  96               "You might want to use the mdp option wall_r_linpot",
  97               x[a][XX],x[a][YY],x[a][ZZ],r);
  98 }
  99
 100 real do_walls(t_inputrec *ir,t_forcerec *fr,matrix box,t_mdatoms *md,
 101               rvec x[],rvec f[],real lambda,real Vlj[],t_nrnb *nrnb)
 102 {
 103     int  nwall,w,lam,i;
 104     int  ntw[2],at,ntype,ngid,ggid,*egp_flags,*type;
 105     real *nbfp,lamfac,fac_d[2],fac_r[2],Cd,Cr,Vtot,Fwall[2];
 106     real wall_z[2],r,mr,r1,r2,r4,Vd,Vr,V=0,Fd,Fr,F=0,dvdlambda;
 107     dvec xf_z;
 108     int  n0,nnn;
 109     real tabscale,*VFtab,rt,eps,eps2,Yt,Ft,Geps,Heps,Heps2,Fp,VV,FF;
 110     unsigned short *gid=md->cENER;
 111     t_forcetable *tab;
 112
 113     nwall     = ir->nwall;
 114     ngid      = ir->opts.ngener;
 115     ntype     = fr->ntype;
 116     nbfp      = fr->nbfp;
 117     egp_flags = fr->egp_flags;
 118
 119     for(w=0; w<nwall; w++)
 120     {
 121         ntw[w] = 2*ntype*ir->wall_atomtype[w];
 122         switch (ir->wall_type)
 123         {
 124         case ewt93:
 125             fac_d[w] = ir->wall_density[w]*M_PI/6;
 126             fac_r[w] = ir->wall_density[w]*M_PI/45;
 127             break;
 128         case ewt104:
 129             fac_d[w] = ir->wall_density[w]*M_PI/2;
 130             fac_r[w] = ir->wall_density[w]*M_PI/5;
 131             break;
 132         default:
 133             break;
 134         }
 135         Fwall[w] = 0;
 136     }
 137     wall_z[0] = 0;
 138     wall_z[1] = box[ZZ][ZZ];
 139
 140     Vtot = 0;
 141     dvdlambda = 0;
 142     clear_dvec(xf_z);
 143     for(lam=0; lam<(md->nPerturbed ? 2 : 1); lam++)
 144     {
 145         if (md->nPerturbed)
 146         {
 147             if (lam == 0)
 148             {
 149                 lamfac = 1 - lambda;
 150                 type = md->typeA;
 151             }
 152             else
 153             {
 154                 lamfac = 0;
 155                 type = md->typeB;
 156             }
 157         }
 158         else
 159         {
 160             lamfac = 1;
 161             type = md->typeA;
 162         }
 163         for(i=md->start; i<md->start+md->homenr; i++)
 164         {
 165             for(w=0; w<nwall; w++)
 166             {
 167                 /* The wall energy groups are always at the end of the list */
 168                 ggid = gid[i]*ngid + ngid - nwall + w;
 169                 at = type[i];
 170                 Cd = nbfp[ntw[w]+2*at];
 171                 Cr = nbfp[ntw[w]+2*at+1];
 172                 if (!((Cd==0 && Cr==0) || egp_flags[ggid] & EGP_EXCL))
 173                 {
 174                     if (w == 0)
 175                     {
 176                         r = x[i][ZZ];
 177                     }
 178                     else
 179                     {
 180                         r = wall_z[1] - x[i][ZZ];
 181                     }
 182                     if (r < ir->wall_r_linpot)
 183                     {
 184                         mr = ir->wall_r_linpot - r;
 185                         r  = ir->wall_r_linpot;
 186                     }
 187                     else
 188                     {
 189                         mr = 0;
 190                     }
 191                     switch (ir->wall_type)
 192                     {
 193                     case ewtTABLE:
 194                         if (r < 0)
 195                         {
 196                             wall_error(i,x,r);
 197                         }
 198                         tab = &(fr->wall_tab[w][gid[i]]);
 199                         tabscale = tab->scale;
 200                         VFtab    = tab->tab;
 201
 202                         rt    = r*tabscale;
 203                         n0    = rt;
 204                         if (n0 >= tab->n)
 205                         {
 206                             /* Beyond the table range, set V and F to zero */
 207                             V     = 0;
 208                             F     = 0;
 209                         }
 210                         else
 211                         {
 212                             eps   = rt - n0;
 213                             eps2  = eps*eps;
 214                             /* Dispersion */
 215                             nnn   = 8*n0;
 216                             Yt    = VFtab[nnn];
 217                             Ft    = VFtab[nnn+1];
 218                             Geps  = VFtab[nnn+2]*eps;
 219                             Heps2 = VFtab[nnn+3]*eps2;
 220                             Fp    = Ft + Geps + Heps2;
 221                             VV    = Yt + Fp*eps;
 222                             FF    = Fp + Geps + 2.0*Heps2;
 223                             Vd    = Cd*VV;
 224                             Fd    = Cd*FF;
 225                             /* Repulsion */
 226                             nnn   = nnn + 4;
 227                             Yt    = VFtab[nnn];
 228                             Ft    = VFtab[nnn+1];
 229                             Geps  = VFtab[nnn+2]*eps;
 230                             Heps2 = VFtab[nnn+3]*eps2;
 231                             Fp    = Ft + Geps + Heps2;
 232                             VV    = Yt + Fp*eps;
 233                             FF    = Fp + Geps + 2.0*Heps2;
 234                             Vr    = Cr*VV;
 235                             Fr    = Cr*FF;
 236                             V     = Vd + Vr;
 237                             F     = -lamfac*(Fd + Fr)*tabscale;
 238                         }
 239                         break;
 240                     case ewt93:
 241                         if (r <= 0)
 242                         {
 243                             wall_error(i,x,r);
 244                         }
 245                         r1 = 1/r;
 246                         r2 = r1*r1;
 247                         r4 = r2*r2;
 248                         Vd = fac_d[w]*Cd*r2*r1;
 249                         Vr = fac_r[w]*Cr*r4*r4*r1;
 250                         V  = Vr - Vd;
 251                         F  = lamfac*(9*Vr - 3*Vd)*r1;
 252                         break;
 253                     case ewt104:
 254                         if (r <= 0)
 255                         {
 256                             wall_error(i,x,r);
 257                         }
 258                         r1 = 1/r;
 259                         r2 = r1*r1;
 260                         r4 = r2*r2;
 261                         Vd = fac_d[w]*Cd*r4;
 262                         Vr = fac_r[w]*Cr*r4*r4*r2;
 263                         V  = Vr - Vd;
 264                         F  = lamfac*(10*Vr - 4*Vd)*r1;
 265                         break;
 266                     case ewt126:
 267                         if (r <= 0)
 268                         {
 269                             wall_error(i,x,r);
 270                         }
 271                         r1 = 1/r;
 272                         r2 = r1*r1;
 273                         r4 = r2*r2;
 274                         Vd = Cd*r4*r2;
 275                         Vr = Cr*r4*r4*r4;
 276                         V  = Vr - Vd;
 277                         F  = lamfac*(12*Vr - 6*Vd)*r1;
 278                         break;
 279                     default:
 280                         break;
 281                     }
 282                     if (mr > 0)
 283                     {
 284                         V += mr*F;
 285                     }
 286                     if (w == 1)
 287                     {
 288                         F = -F;
 289                     }
 290                     Vlj[ggid] += lamfac*V;
 291                     Vtot      += V;
 292                     f[i][ZZ]  += F;
 293                     /* Because of the single sum virial calculation we need
 294                      * to add  the full virial contribution of the walls.
 295                      * Since the force only has a z-component, there is only
 296                      * a contribution to the z component of the virial tensor.
 297                      * We could also determine the virial contribution directly,
 298                      * which would be cheaper here, but that would require extra
 299                      * communication for f_novirsum for with virtual sites
 300                      * in parallel.
 301                      */
 302                     xf_z[XX]  -= x[i][XX]*F;
 303                     xf_z[YY]  -= x[i][YY]*F;
 304                     xf_z[ZZ]  -= wall_z[w]*F;
 305                 }
 306             }
 307         }
 308         if (md->nPerturbed)
 309         {
 310             dvdlambda += (lam==0 ? -1 : 1)*Vtot;
 311         }
 312
 313         inc_nrnb(nrnb,eNR_WALLS,md->homenr);
 314     }
 315
 316     for(i=0; i<DIM; i++)
 317     {
 318         fr->vir_wall_z[i] = -0.5*xf_z[i];
 319     }
 320
 321     return dvdlambda;
 322 }