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->data[8*i+j] = tab->data[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 /* nbfp now includes the 6.0/12.0 derivative prefactors */
171 Cd = nbfp[ntw[w]+2*at]/6.0;
172 Cr = nbfp[ntw[w]+2*at+1]/12.0;
173 if (!((Cd==0 && Cr==0) || (egp_flags[ggid] & EGP_EXCL)))
181 r = wall_z[1] - x[i][ZZ];
183 if (r < ir->wall_r_linpot)
185 mr = ir->wall_r_linpot - r;
186 r = ir->wall_r_linpot;
192 switch (ir->wall_type)
199 tab = &(fr->wall_tab[w][gid[i]]);
200 tabscale = tab->scale;
207 /* Beyond the table range, set V and F to zero */
219 Geps = VFtab[nnn+2]*eps;
220 Heps2 = VFtab[nnn+3]*eps2;
221 Fp = Ft + Geps + Heps2;
223 FF = Fp + Geps + 2.0*Heps2;
230 Geps = VFtab[nnn+2]*eps;
231 Heps2 = VFtab[nnn+3]*eps2;
232 Fp = Ft + Geps + Heps2;
234 FF = Fp + Geps + 2.0*Heps2;
238 F = -lamfac*(Fd + Fr)*tabscale;
249 Vd = fac_d[w]*Cd*r2*r1;
250 Vr = fac_r[w]*Cr*r4*r4*r1;
252 F = lamfac*(9*Vr - 3*Vd)*r1;
263 Vr = fac_r[w]*Cr*r4*r4*r2;
265 F = lamfac*(10*Vr - 4*Vd)*r1;
278 F = lamfac*(12*Vr - 6*Vd)*r1;
291 Vlj[ggid] += lamfac*V;
294 /* Because of the single sum virial calculation we need
295 * to add the full virial contribution of the walls.
296 * Since the force only has a z-component, there is only
297 * a contribution to the z component of the virial tensor.
298 * We could also determine the virial contribution directly,
299 * which would be cheaper here, but that would require extra
300 * communication for f_novirsum for with virtual sites
303 xf_z[XX] -= x[i][XX]*F;
304 xf_z[YY] -= x[i][YY]*F;
305 xf_z[ZZ] -= wall_z[w]*F;
311 dvdlambda += (lam==0 ? -1 : 1)*Vtot;
314 inc_nrnb(nrnb,eNR_WALLS,md->homenr);
319 fr->vir_wall_z[i] = -0.5*xf_z[i];