2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
55 void make_wall_tables(FILE *fplog,const output_env_t oenv,
56 const t_inputrec *ir,const char *tabfn,
57 const gmx_groups_t *groups,
60 int w,negp_pp,egp,i,j;
65 negp_pp = ir->opts.ngener - ir->nwall;
66 nm_ind = groups->grps[egcENER].nm_ind;
69 fprintf(fplog,"Reading user tables for %d energy groups with %d walls\n",
73 snew(fr->wall_tab,ir->nwall);
74 for(w=0; w<ir->nwall; w++) {
75 snew(fr->wall_tab[w],negp_pp);
76 for(egp=0; egp<negp_pp; egp++) {
77 /* If the energy group pair is excluded, we don't need a table */
78 if (!(fr->egp_flags[egp*ir->opts.ngener+negp_pp+w] & EGP_EXCL)) {
79 tab = &fr->wall_tab[w][egp];
80 sprintf(buf,"%s",tabfn);
81 sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1,"_%s_%s.%s",
82 *groups->grpname[nm_ind[egp]],
83 *groups->grpname[nm_ind[negp_pp+w]],
85 *tab = make_tables(fplog,oenv,fr,FALSE,buf,0,GMX_MAKETABLES_FORCEUSER);
86 /* Since wall have no charge, we can compress the table */
87 for(i=0; i<=tab->n; i++)
89 tab->data[8*i+j] = tab->data[12*i+4+j];
95 static void wall_error(int a,rvec *x,real r)
98 "An atom is beyond the wall: coordinates %f %f %f, distance %f\n"
99 "You might want to use the mdp option wall_r_linpot",
100 x[a][XX],x[a][YY],x[a][ZZ],r);
103 real do_walls(t_inputrec *ir,t_forcerec *fr,matrix box,t_mdatoms *md,
104 rvec x[],rvec f[],real lambda,real Vlj[],t_nrnb *nrnb)
107 int ntw[2],at,ntype,ngid,ggid,*egp_flags,*type;
108 real *nbfp,lamfac,fac_d[2],fac_r[2],Cd,Cr,Vtot,Fwall[2];
109 real wall_z[2],r,mr,r1,r2,r4,Vd,Vr,V=0,Fd,Fr,F=0,dvdlambda;
112 real tabscale,*VFtab,rt,eps,eps2,Yt,Ft,Geps,Heps,Heps2,Fp,VV,FF;
113 unsigned short *gid=md->cENER;
117 ngid = ir->opts.ngener;
120 egp_flags = fr->egp_flags;
122 for(w=0; w<nwall; w++)
124 ntw[w] = 2*ntype*ir->wall_atomtype[w];
125 switch (ir->wall_type)
128 fac_d[w] = ir->wall_density[w]*M_PI/6;
129 fac_r[w] = ir->wall_density[w]*M_PI/45;
132 fac_d[w] = ir->wall_density[w]*M_PI/2;
133 fac_r[w] = ir->wall_density[w]*M_PI/5;
141 wall_z[1] = box[ZZ][ZZ];
146 for(lam=0; lam<(md->nPerturbed ? 2 : 1); lam++)
166 for(i=md->start; i<md->start+md->homenr; i++)
168 for(w=0; w<nwall; w++)
170 /* The wall energy groups are always at the end of the list */
171 ggid = gid[i]*ngid + ngid - nwall + w;
173 /* nbfp now includes the 6.0/12.0 derivative prefactors */
174 Cd = nbfp[ntw[w]+2*at]/6.0;
175 Cr = nbfp[ntw[w]+2*at+1]/12.0;
176 if (!((Cd==0 && Cr==0) || (egp_flags[ggid] & EGP_EXCL)))
184 r = wall_z[1] - x[i][ZZ];
186 if (r < ir->wall_r_linpot)
188 mr = ir->wall_r_linpot - r;
189 r = ir->wall_r_linpot;
195 switch (ir->wall_type)
202 tab = &(fr->wall_tab[w][gid[i]]);
203 tabscale = tab->scale;
210 /* Beyond the table range, set V and F to zero */
222 Geps = VFtab[nnn+2]*eps;
223 Heps2 = VFtab[nnn+3]*eps2;
224 Fp = Ft + Geps + Heps2;
226 FF = Fp + Geps + 2.0*Heps2;
233 Geps = VFtab[nnn+2]*eps;
234 Heps2 = VFtab[nnn+3]*eps2;
235 Fp = Ft + Geps + Heps2;
237 FF = Fp + Geps + 2.0*Heps2;
241 F = -lamfac*(Fd + Fr)*tabscale;
252 Vd = fac_d[w]*Cd*r2*r1;
253 Vr = fac_r[w]*Cr*r4*r4*r1;
255 F = lamfac*(9*Vr - 3*Vd)*r1;
266 Vr = fac_r[w]*Cr*r4*r4*r2;
268 F = lamfac*(10*Vr - 4*Vd)*r1;
281 F = lamfac*(12*Vr - 6*Vd)*r1;
294 Vlj[ggid] += lamfac*V;
297 /* Because of the single sum virial calculation we need
298 * to add the full virial contribution of the walls.
299 * Since the force only has a z-component, there is only
300 * a contribution to the z component of the virial tensor.
301 * We could also determine the virial contribution directly,
302 * which would be cheaper here, but that would require extra
303 * communication for f_novirsum for with virtual sites
306 xf_z[XX] -= x[i][XX]*F;
307 xf_z[YY] -= x[i][YY]*F;
308 xf_z[ZZ] -= wall_z[w]*F;
314 dvdlambda += (lam==0 ? -1 : 1)*Vtot;
317 inc_nrnb(nrnb,eNR_WALLS,md->homenr);
322 fr->vir_wall_z[i] = -0.5*xf_z[i];