1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
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,
57 int w,negp_pp,egp,i,j;
62 negp_pp = ir->opts.ngener - ir->nwall;
63 nm_ind = groups->grps[egcENER].nm_ind;
66 fprintf(fplog,"Reading user tables for %d energy groups with %d walls\n",
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]],
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++)
86 tab->tab[8*i+j] = tab->tab[12*i+4+j];
92 static void wall_error(int a,rvec *x,real r)
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);
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)
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;
109 real tabscale,*VFtab,rt,eps,eps2,Yt,Ft,Geps,Heps,Heps2,Fp,VV,FF;
110 unsigned short *gid=md->cENER;
114 ngid = ir->opts.ngener;
117 egp_flags = fr->egp_flags;
119 for(w=0; w<nwall; w++)
121 ntw[w] = 2*ntype*ir->wall_atomtype[w];
122 switch (ir->wall_type)
125 fac_d[w] = ir->wall_density[w]*M_PI/6;
126 fac_r[w] = ir->wall_density[w]*M_PI/45;
129 fac_d[w] = ir->wall_density[w]*M_PI/2;
130 fac_r[w] = ir->wall_density[w]*M_PI/5;
138 wall_z[1] = box[ZZ][ZZ];
143 for(lam=0; lam<(md->nPerturbed ? 2 : 1); lam++)
163 for(i=md->start; i<md->start+md->homenr; i++)
165 for(w=0; w<nwall; w++)
167 /* The wall energy groups are always at the end of the list */
168 ggid = gid[i]*ngid + ngid - nwall + w;
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)))
180 r = wall_z[1] - x[i][ZZ];
182 if (r < ir->wall_r_linpot)
184 mr = ir->wall_r_linpot - r;
185 r = ir->wall_r_linpot;
191 switch (ir->wall_type)
198 tab = &(fr->wall_tab[w][gid[i]]);
199 tabscale = tab->scale;
206 /* Beyond the table range, set V and F to zero */
218 Geps = VFtab[nnn+2]*eps;
219 Heps2 = VFtab[nnn+3]*eps2;
220 Fp = Ft + Geps + Heps2;
222 FF = Fp + Geps + 2.0*Heps2;
229 Geps = VFtab[nnn+2]*eps;
230 Heps2 = VFtab[nnn+3]*eps2;
231 Fp = Ft + Geps + Heps2;
233 FF = Fp + Geps + 2.0*Heps2;
237 F = -lamfac*(Fd + Fr)*tabscale;
248 Vd = fac_d[w]*Cd*r2*r1;
249 Vr = fac_r[w]*Cr*r4*r4*r1;
251 F = lamfac*(9*Vr - 3*Vd)*r1;
262 Vr = fac_r[w]*Cr*r4*r4*r2;
264 F = lamfac*(10*Vr - 4*Vd)*r1;
277 F = lamfac*(12*Vr - 6*Vd)*r1;
290 Vlj[ggid] += lamfac*V;
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
302 xf_z[XX] -= x[i][XX]*F;
303 xf_z[YY] -= x[i][YY]*F;
304 xf_z[ZZ] -= wall_z[w]*F;
310 dvdlambda += (lam==0 ? -1 : 1)*Vtot;
313 inc_nrnb(nrnb,eNR_WALLS,md->homenr);
318 fr->vir_wall_z[i] = -0.5*xf_z[i];